[已解决] 瞬时声压时域数据怎么用matlab进行1/3倍频程声压级分析

  [复制链接]
simplebinbin 发表于 2014-5-8 14:02:49
求大神指导,小弟采集了一列瞬时声压数据,怎么样用matlab做1/3倍频声压级呢?
哪位大神有程序,给小弟一个,感激不尽啊,我自己也正在编程序,编号了发上来请各位指导。

最佳答案


songzy41 发表于 2014-5-16 17:42:48
那可以按照LZ在2层的程序运行。但程序中fs可以从wav文件中读出,是48000,同LZ设置的不对,又t1和t2没有设定,又nnn没有设定。

103 条回复


simplebinbin 发表于 2014-5-14 11:12:33
  1. %A计权声压级频谱分析
  2. clc;
  3. clear;
  4. close all;
  5. %时域分析
  6. y=wavread('abc.wav');
  7. %频域分析

  8. fs=51200;%采样频率
  9. p0=2e-5;%参考声压

  10. f=[1.00 1.25 1.600 2.00 2.50 3.15 4.00 5.00 6.30 8.0]; %基准中心频率
  11. f1=[20.00 25.0 31.5 40.0 50.0 63.0 80];
  12. fc=[f1,100*f,1000*f,10000*f]; %%%%%%%%%中心频率%%%%%%%%
  13. %20-16000Hz A声级计权值
  14. cf=[-50.5,-44.7,-39.4,-34.6,-30.2,-26.2,-22.5,-19.1,-16.1,-13.4,-10.9,-8.6,-6.6,-4.8,-3.2,-1.9,-0.8,0,0.6,1.0,1.2,1.3,1.2,1.0,0.5,-0.1,-1.1,-2.5,-4.3,-6.6];

  15. x=y(t1*fs:t2*fs);%截取需要处理的数据段
  16. n=length(x);
  17. t=(0:1/fs:(n-1)/fs);

  18. subplot(221);
  19. plot(t,x);%瞬时声压时程图

  20. w=hanning(n);    %汉宁窗
  21. xx=1.633*x.*w;         %加汉宁窗(恢复系数为1.633)


  22. nfft=2^nextpow2(n);
  23. %nextpow2(n)-取最接近的较大2次幂
  24. a = fft(xx,nfft);
  25. f = fs/2*linspace(0,1,nfft/2);

  26. w=2*abs(a(1:nfft/2)/n);
  27. subplot(222);
  28. plot(f,w);%绘制频谱图
  29. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  30. %1/3倍频程计算
  31. oc6=2^(1/6);
  32. nc=length(cf);
  33. %下面这个求1/3倍频程的程序是按照振动振级计算那个来的
  34. for j=1:nc
  35.     fl=fc(j)/oc6;
  36.     fu=fc(j)*oc6;
  37.     nl=round(fl*nfft/fs+1);
  38.     nu=round(fu*nfft/fs+1);
  39.     if fu>fs/2
  40.         m=j-1;
  41.         break;
  42.     end
  43.    
  44.     b=zeros(1,nfft);
  45.     b(nl:nu)=a(nl:nu);
  46.     b(nfft-nu+1:nfft-nl+1)=a(nfft-nu+1:nfft-nl+1);
  47.     c=ifft(b,nfft);
  48.     yc(j)=sqrt(var(real(c(1:nnn))));
  49.    
  50. end

  51. aj_sumn=0;
  52. for i=1:nc
  53.     Lp1(i)=20*log10(yc(i)/p0);%未计权1/3倍频程声压级
  54. end

  55. %%%%%
  56. for jj=1:nc
  57.      aj_sumn=aj_sumn+10^(0.1*Lp1(j));
  58. end
  59. Lp=10*log10(aj_sumn);%未计权总声压级
  60. subplot(223);%绘制未计权1/3倍频程声压级图谱

  61. bar(Lp1(1:nc));
  62. gg=zeros(1,nc);
  63. for i=1:nc
  64.     gg(1:nc)=fc(1:nc);
  65. end
  66. ggg=1:nc;   
  67. set(gca,'xtick',ggg);
  68. set(gca,'xticklabel',gg);

  69. %%%%%A计权1/3倍频程声压级
  70. Lap=Lp1+cf;
  71. aj_sum=0;
  72. for j=1:nc
  73.     aj_sum=aj_sum+10^(0.1*Lap(j));
  74. end
  75. LA=10*log10(aj_sum);%Aa计权总声压级

  76. subplot(224);%绘制A计权1/3倍频程声压级图谱
  77. bar(Lap(1:nc));
  78. gg=zeros(1,nc);
  79. for i=1:nc
  80.     gg(1:nc)=fc(1:nc);
  81. end
  82. ggg=1:nc;   
  83. set(gca,'xtick',ggg);
  84. set(gca,'xticklabel',gg);

  85.   
复制代码

simplebinbin 发表于 2014-5-14 11:14:54
按照振动的代码写了个,与声级计计算的有的时段差别较小(0.5dB(A)左右),有的时段差别较大5dB(A)左右,大家有好的建议可以提,一起讨论

songzy41 发表于 2014-5-14 14:40:33
LZ在声压计中对于p0肯定不是一个设置的值:p0=2e-5;%参考声压,而是通过标准信号发生器对了传声器发出一个标准的1000HZ信号得到的。同样在计算中也不能用设置的办法来给出p0。p0=2e-5帕是没有错,它是一个声压单位,但LZ的声信号已经过传声器、放大器、A/D变换和量化等,输入到计算机是电压信号。LZ对p0的设置应设置为与2e-5帕对应的输入计算机的值。

simplebinbin 发表于 2014-5-14 15:08:07
songzy41 发表于 2014-5-14 14:40
LZ在声压计中对于p0肯定不是一个设置的值:p0=2e-5;%参考声压,而是通过标准信号发生器对了传声器发出一个 ...

版主,您好,我采集的信号是瞬时声压信号,y=wavread('abc.wav'),abc.wav就是保存的原始瞬时声压值p(1秒内有Fs个声压值),单位是Pa(您说这个p是电压信号,我觉得这就跟振动信号差不多,比如加速度信号,得到的是加速度幅值)
程序的核心思想是根据书上的公式:
将瞬时声压值进行1/3倍频程滤波,在逆变换,求中心频率的有效值,用声压级公式:Lpi=20log(p/p0);
计算出中心频带的声压级,然后用HJ 453-2008中的公式,计算总声级:Ltp=10log∑10^(Lpi/10);
A计权就在对应频带上减去计权值就可以了,总A声级计算公式一致。

songzy41 发表于 2014-5-14 17:00:42
不论声压、振动或其它物理都是模拟的连续信号,但经过换能器、又通过放大器(不管是调理线路中的前置放大,或把信号放大的放大器),此时相关的物 理量已变成统一的电压(或电流)信号,再经A/D变换,连续量转换成离散量,计算机中就是这种离散值。除了对离散值进行校正,例如对声信号来说,离散值在频率1000HZ时对应于声压是多少帕,或者p0对应于离散值是多少,这样才有可能计算绝对声压级,否则计算的都是相对声压级。
LZ读入的信号y=wavread('abc.wav'),本身并不带有单位帕,如果没有进行上述的校正,加上单位帕只是人为的行为,而不是实际的帕,没有太大的意义。

simplebinbin 发表于 2014-5-15 15:02:54
songzy41 发表于 2014-5-14 17:00
不论声压、振动或其它物理都是模拟的连续信号,但经过换能器、又通过放大器(不管是调理线路中的前置放大, ...

版主,您好!我在声学方面的分析刚开始,目前能得到的数据就只有声级计导出的原始数据'abc.wav',我知道这些数据经过放大器放大及A/D变换得到,按照您的意思,需要对'abc.wav'的数据进行修正,请问怎样进行呢,“离散值在频率1000HZ时对应于声压是多少帕”,这个怎么弄呢?p0对应一个离散值,多个离散值就需要多个参考声压吗?
谢谢!
目前网上对这种计算的程序资料很少,无法参看,我只能按照自己的理解来弄这个程序,很难验证。

simplebinbin 发表于 2014-5-15 15:06:50
abc.wav文件,声级计导出的原始数据。

abc.rar

384.99 KB, 下载次数: 36055

abc.wav文件


songzy41 发表于 2014-5-15 20:21:48
声级机导出的数据是wav格式?它的幅值是在-1~1之间,这幅值的单位是帕?
现在测量问题属于声学测量范畴,建议查一下声压计的有关手册,或问一下你的指导老师。

simplebinbin 发表于 2014-5-16 10:06:06
songzy41 发表于 2014-5-15 20:21
声级机导出的数据是wav格式?它的幅值是在-1~1之间,这幅值的单位是帕?
现在测量问题属于声学测量范畴,建 ...

宋老师,不好意思,刚接触声学这块,不太了解。我刚买了本您的书《MATLAB在语音信号分析与合成中的应用》先慢慢学习着。
声级计导出的是wav格式文件,声级计厂家说这个数据是原始声压数据,可以导出来直接做频谱分析及A声级计算,abc.wav的幅值如下图

原始声压数据

原始声压数据


songzy41 发表于 2014-5-16 14:16:16
我也在网上查了一些资料,知道声级计可以用wav格式输出。但我昨日问LZ,读入的wav格式的数据单位是什么?是伏还是帕?从声级计校正时能得到声级计的灵敏度,如下图所示。如果输出单位是伏就能按这灵敏度计算出对应的声压和声压级了。
calibration.JPG

simplebinbin 发表于 2014-5-16 15:27:55
songzy41 发表于 2014-5-16 14:16
我也在网上查了一些资料,知道声级计可以用wav格式输出。但我昨日问LZ,读入的wav格式的数据单位是什么?是 ...

宋老师,您好,传感器厂家说这个wav格式的数据是原始的瞬时声压值,单位是“pa"。

songzy41 发表于 2014-5-16 17:42:48
那可以按照LZ在2层的程序运行。但程序中fs可以从wav文件中读出,是48000,同LZ设置的不对,又t1和t2没有设定,又nnn没有设定。
回复此楼

simplebinbin 发表于 2014-5-19 14:00:58
songzy41 发表于 2014-5-16 17:42
那可以按照LZ在2层的程序运行。但程序中fs可以从wav文件中读出,是48000,同LZ设置的不对,又t1和t2没有设 ...

宋老师,您好,我问了传感器厂家,原始采样频率是48000,程序中的fs是我自己计算假定的,t1,t2根据想处理的数据任意设定,nnn应该是n,我对程序做了修改,放在楼下,可以直接运行了。
但是还有一个问题,在高频段,未计权声压级出现了负值,正常情况是不应该含有负值,是吗?

simplebinbin 发表于 2014-5-19 14:01:43
  1. %A计权声压级频谱分析
  2. clc;
  3. clear;
  4. close all;

  5. y=wavread('abc.wav');

  6. fs=51200;%采样频率
  7. p0=2e-5;%参考声压

  8. f=[1.00 1.25 1.600 2.00 2.50 3.15 4.00 5.00 6.30 8.0]; %基准中心频率
  9. f1=[20.00 25.0 31.5 40.0 50.0 63.0 80];
  10. fc=[f1,100*f,1000*f,10000*f]; %%%%%%%%%中心频率%%%%%%%%
  11. %20-16000Hz A声级计权值
  12. cf=[-50.5,-44.7,-39.4,-34.6,-30.2,-26.2,-22.5,-19.1,-16.1,-13.4,-10.9,-8.6,-6.6,-4.8,-3.2,-1.9,-0.8,0,0.6,1.0,1.2,1.3,1.2,1.0,0.5,-0.1,-1.1,-2.5,-4.3,-6.6];
  13. t1=1;
  14. t2=2;

  15. x=y(t1*fs+1:t2*fs);%截取需要处理的数据段
  16. n=length(x);
  17. t=(0:1/fs:(n-1)/fs);

  18. subplot(221);
  19. plot(t,x);%瞬时声压时程图

  20. w=hanning(n);    %汉宁窗
  21. xx=1.633*x.*w;         %加汉宁窗(恢复系数为1.633)


  22. nfft=2^nextpow2(n);
  23. %nextpow2(n)-取最接近的较大2次幂
  24. a = fft(xx,nfft);
  25. f = fs/2*linspace(0,1,nfft/2);

  26. w=2*abs(a(1:nfft/2)/n);
  27. subplot(222);
  28. plot(f,w);%绘制频谱图
  29. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  30. %1/3倍频程计算
  31. oc6=2^(1/6);
  32. nc=length(cf);
  33. %下面这个求1/3倍频程的程序是按照振动振级计算那个来的
  34. for j=1:nc
  35.     fl=fc(j)/oc6;
  36.     fu=fc(j)*oc6;
  37.     nl=round(fl*nfft/fs+1);
  38.     nu=round(fu*nfft/fs+1);
  39.     if fu>fs/2
  40.         m=j-1;
  41.         break;
  42.     end
  43.    
  44.     b=zeros(1,nfft);
  45.     b(nl:nu)=a(nl:nu);
  46.     b(nfft-nu+1:nfft-nl+1)=a(nfft-nu+1:nfft-nl+1);
  47.     c=ifft(b,nfft);
  48.     yc(j)=sqrt(var(real(c(1:n))));
  49.    
  50. end

  51. aj_sumn=0;
  52. for i=1:nc
  53.     Lp1(i)=20*log10(yc(i)/p0);%未计权1/3倍频程声压级
  54. end

  55. %%%%%
  56. for jj=1:nc
  57.      aj_sumn=aj_sumn+10^(0.1*Lp1(j));
  58. end
  59. Lp=10*log10(aj_sumn);%未计权总声压级
  60. subplot(223);%绘制未计权1/3倍频程声压级图谱

  61. bar(Lp1(1:nc));
  62. gg=zeros(1,nc);
  63. for i=1:nc
  64.     gg(1:nc)=fc(1:nc);
  65. end
  66. ggg=1:nc;   
  67. set(gca,'xtick',ggg);
  68. set(gca,'xticklabel',gg);

  69. %%%%%A计权1/3倍频程声压级
  70. Lap=Lp1+cf;
  71. aj_sum=0;
  72. for j=1:nc
  73.     aj_sum=aj_sum+10^(0.1*Lap(j));
  74. end
  75. LA=10*log10(aj_sum);%Aa计权总声压级

  76. subplot(224);%绘制A计权1/3倍频程声压级图谱
  77. bar(Lap(1:nc));
  78. gg=zeros(1,nc);
  79. for i=1:nc
  80.     gg(1:nc)=fc(1:nc);
  81. end
  82. ggg=1:nc;   
  83. set(gca,'xtick',ggg);
  84. set(gca,'xticklabel',gg);

  85.   
复制代码

simplebinbin 发表于 2014-5-21 10:00:06

67行中67.for jj=1:nc

jj改为j

simplebinbin 发表于 2014-5-21 10:39:02
参考别的文献,发现有新的方法计算噪声1/3倍频程计权声级的方法,参考文献:张登攀《噪声1/3倍频程计权声级算法》[J].
编写了部分程序:
  1. clc;
  2. clear;
  3. close all;
  4. %时域分析
  5. y=wavread('abc.wav');
  6. %频域分析
  7. fs=51200;%采样频率
  8. p0=2e-5;%参考声压
  9. f=[1.00 1.25 1.600 2.00 2.50 3.15 4.00 5.00 6.30 8.0]; %基准中心频率
  10. f1=[20.00 25.0 31.5 40.0 50.0 63.0 80];
  11. fc=[f1,100*f,1000*f,10000*f]; %%%%%%%%%中心频率%%%%%%%%
  12. %20-16000Hz A声级计权值
  13. cf=[-50.5,-44.7,-39.4,-34.6,-30.2,-26.2,-22.5,-19.1,-16.1,-13.4,-10.9,-8.6,-6.6,-4.8,-3.2,-1.9,-0.8,0,0.6,1.0,1.2,1.3,1.2,1.0,0.5,-0.1,-1.1,-2.5,-4.3,-6.6];

  14. n=length(y);
  15. t=(0:1/fs:(n-1)/fs);
  16. h1=figure;
  17. plot(t,y);
  18. title('瞬时声压时程');
  19. xlabel('Time(s)');
  20. ylabel('Sound Presure Value(Pa)');

  21. %
  22. t1=0;
  23. t2=4;

  24. x=y(t1*fs+1:t2*fs);%截取需要处理的数据段
  25. n=length(x);
  26. %t=(0:1/fs:(n-1)/fs);
  27. %plot(t,x);%瞬时声压时程图
  28. w=1.633*hanning(n);    %汉宁窗(恢复系数为1.633)
  29. %w=1.812*blackmanharris(n);    %布拉克曼窗(功率相等恢复系数1.812)
  30. xx=x.*w;         %加汉宁窗
  31. nfft=2^nextpow2(n);    %nextpow2(n)-取最接近的较大2次幂
  32. a = fft(xx,nfft)/n;
  33. %f = fs/2*linspace(0,1,nfft/2);
  34. w=2*abs(a(1:nfft/2));

  35. oc6=2^(1/6);
  36. nc=length(cf);

  37. for j=1:nc
  38.     fl=fc(j)/oc6;
  39.     fu=fc(j)*oc6;
  40.     nl=round(fl*nfft/fs+1);
  41.     nu=round(fu*nfft/fs+1);
  42.     if fu>fs/2
  43.         m=j-1;
  44.         break;
  45.     end
  46.    
  47.     p=w(nl:nu);
  48.     lp=length(p);
  49.     k=0;
  50.     for ii=1:lp
  51.         if ii+2>lp
  52.         break
  53.         end
  54.         if p(ii+1)>p(ii)&&p(ii+1)>p(ii+2)
  55.             k=k+1;
  56.             pp(k)=p(ii+1)/sqrt(2);
  57.         end
  58.     end
  59.      p2(j)=sum(pp.*pp);
  60.      
  61. end
  62. Lp=10*log10(p2/p0^2);
  63. for jj=1:length(Lp)
  64.     Lp1(jj)=10^(Lp(jj)/10);
  65. end
  66. Lpt=10*log10(sum(Lp1))

  67. h2=figure;
  68. mm=nc;
  69. bar(Lp(1:mm));
  70. gg=zeros(1,mm);
  71. for i=1:mm
  72.     gg(1:mm)=fc(1:mm);
  73. end
  74. ggg=1:mm;   
  75. set(gca,'xtick',ggg);
  76. set(gca,'xticklabel',gg);
  77. set(gcf,'PaperPosition',[1,1,40,20])
  78. set(gca,'fontsize',10)
  79. xlabel('Frequency(Hz)');
  80. ylabel('SPL(dB)');
  81. title('未计权声压级');
  82. grid on;
复制代码

simplebinbin 发表于 2014-5-21 10:54:19
根据这两种方法计算结果还是有一些差别的,前一种方法称为方法1,后一种为方法2。计算abc.wav数据中前4s时程。总未计权声级分别为39.7dB、36.3dB。1/3倍频程图谱见附件

方法1

方法1

方法2

方法2

Themoment 发表于 2014-6-26 09:49:37
学习了好多知识,感谢大家的分享。
倍频程后面的负值是不是有问题?

simplebinbin 发表于 2014-6-26 16:43:07
Themoment 发表于 2014-6-26 09:49
学习了好多知识,感谢大家的分享。
倍频程后面的负值是不是有问题?

没有问题,负值说明这个频段有效声压值小于参考声压值

其行浩浩 发表于 2015-1-12 11:07:00
很好的帖子,学习啦!

simplebinbin 发表于 2015-1-13 10:47:52
最近在分析wav数据时发现,wav文件一般都做了归一化处理,MATLAB读取时需要进行修正(利用校准进行幅值修正)
y=wavread('abc.wav');读取数据时,matlab进行归一化
y=wavread('abc.wav','native');读取原始的wav文件数据,与瞬时声压值有差别

大家在进行声级计算时需要修正,频谱分析不修正也可以,对频谱峰值分布没有影响

laoxuzi 发表于 2015-2-20 23:50:26
学习中。。。。。

mushroom不止 发表于 2015-3-1 21:28:31
simplebinbin 发表于 2014-5-21 10:54
根据这两种方法计算结果还是有一些差别的,前一种方法称为方法1,后一种为方法2。计算abc.wav数据中前4s时 ...

楼主您好,仔细阅读了您的程序,收货颇丰。但是有一点疑问,希望指教。根据声级计分析的结果,是方法1更接近真实的A声压级分布,还是方法2?

simplebinbin 发表于 2015-3-23 09:12:35
mushroom不止 发表于 2015-3-1 21:28
楼主您好,仔细阅读了您的程序,收货颇丰。但是有一点疑问,希望指教。根据声级计分析的结果,是方法1更 ...

方法1与声级计计算结果差不多
您需要登录后才可以回帖 登录 | 注册

本版积分规则

相关帖子
相关文章
热门教程
站长推荐
快速回复 返回顶部 返回列表