A声级及1/3倍频程计算matlab程序

simplebinbin

小弟最近入门声学分析，用声级计采集了瞬时声压数据，自己按照振动分析的代码写了个程序，
贴出来供大家参考，一起讨论学习交流！
通过验证，与声级计分析出来的A计权总声压级在有的时间段差别很小0.5dB（A),有的时间段差别比较大，好几个dB(A)，
希望大家一起完善！

1. %A计权声压级频谱分析
   
2. clc;
   
3. clear;
   
4. close all;
   
5. %时域分析
   
6. y=wavread('abc.wav');
   
7. %频域分析
   
8. 
   
9. fs=51200;%采样频率
   
10. p0=2e-5;%参考声压
    
11. 
    
12. f=[1.00 1.25 1.600 2.00 2.50 3.15 4.00 5.00 6.30 8.0]; %基准中心频率
    
13. f1=[20.00 25.0 31.5 40.0 50.0 63.0 80];
    
14. fc=[f1,100*f,1000*f,10000*f]; %%%%%%%%%中心频率%%%%%%%%
    
15. %20-16000Hz A声级计权值
    
16. 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];
    
17. 
    
18. x=y(t1*fs:t2*fs);%截取需要处理的数据段
    
19. n=length(x);
    
20. t=(0:1/fs:(n-1)/fs);
    
21. 
    
22. subplot(221);
    
23. plot(t,x);%瞬时声压时程图
    
24. 
    
25. w=hanning(n); %汉宁窗
    
26. xx=1.633*x.*w; %加汉宁窗(恢复系数为1.633)
    
27. 
    
28. 
    
29. nfft=2^nextpow2(n);
    
30. %nextpow2(n)-取最接近的较大2次幂
    
31. a = fft(xx,nfft);
    
32. f = fs/2*linspace(0,1,nfft/2);
    
33. 
    
34. w=2*abs(a(1:nfft/2)/n);
    
35. subplot(222）；
    
36. plot(f,w);%绘制频谱图
    
37. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
38. %1/3倍频程计算
    
39. oc6=2^(1/6);
    
40. nc=length(cf);
    
41. %下面这个求1/3倍频程的程序是按照振动振级计算那个来的
    
42. for j=1:nc
    
43. fl=fc(j)/oc6;
    
44. fu=fc(j)*oc6;
    
45. nl=round(fl*nfft/fs+1);
    
46. nu=round(fu*nfft/fs+1);
    
47. if fu>fs/2
    
48. m=j-1;
    
49. break;
    
50. end
    
51. 
    
52. b=zeros(1,nfft);
    
53. b(nl:nu)=a(nl:nu);
    
54. b(nfft-nu+1:nfft-nl+1)=a(nfft-nu+1:nfft-nl+1);
    
55. c=ifft(b,nfft);
    
56. yc(j)=sqrt(var(real(c(1:nnn))));
    
57. 
    
58. end
    
59. 
    
60. aj_sumn=0;
    
61. for i=1:nc
    
62. Lp1(i)=20*log10(yc(i)/p0);%未计权1/3倍频程声压级
    
63. end
    
64. 
    
65. %%%%%
    
66. for jj=1:nc
    
67. aj_sumn=aj_sumn+10^(0.1*Lp1(j));
    
68. end
    
69. Lp=10*log10(aj_sumn);%未计权总声压级
    
70. subplot(223);%绘制未计权1/3倍频程声压级图谱
    
71. 
    
72. bar(Lp1(1:nc));
    
73. gg=zeros(1,nc);
    
74. for i=1:nc
    
75. gg(1:nc)=fc(1:nc);
    
76. end
    
77. ggg=1:nc;
    
78. set(gca,'xtick',ggg);
    
79. set(gca,'xticklabel',gg);
    
80. 
    
81. %%%%%A计权1/3倍频程声压级
    
82. Lap=Lp1+cf;
    
83. aj_sum=0;
    
84. for j=1:nc
    
85. aj_sum=aj_sum+10^(0.1*Lap(j));
    
86. end
    
87. LA=10*log10(aj_sum);%Aa计权总声压级
    
88. 
    
89. subplot(224);%绘制A计权1/3倍频程声压级图谱
    
90. bar(Lp1(1:nc));
    
91. gg=zeros(1,nc);
    
92. for i=1:nc
    
93. gg(1:nc)=fc(1:nc);
    
94. end
    
95. ggg=1:nc;
    
96. set(gca,'xtick',ggg);
    
97. set(gca,'xticklabel',gg);
    
98. 
    
99. 
    

复制代码

simplebinbin

代码第90行应该是：90.bar(Lap(1:nc));

mywl

谢谢分享哦

mxlzhenzhu

总体思路我看了，多谢分享；

关于其中一个问题：


w=hanning(n); %汉宁窗

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





nfft=2^nextpow2(n);

%nextpow2(n)-取最接近的较大2次幂

a = fft(xx,nfft);

f = fs/2*linspace(0,1,nfft/2);


顺序好像有问题；你得先计算nfft，然后再计算窗函数的权重；
即25行应该为：

w=hanning(nfft); %汉宁窗

mxlzhenzhu

这样才是整周期采样

simplebinbin

mxlzhenzhu 发表于 2014-7-30 00:13
这样才是整周期采样

与整周期采样关系不大
你这样只能处理特定的数据，即采样频率为f=2^n,数据长度为T*f（T=2^m），m、n均为正整数
缺失了程序的通用性

mxlzhenzhu

simplebinbin 发表于 2014-8-1 11:52
与整周期采样关系不大
你这样只能处理特定的数据，即采样频率为f=2^n,数据长度为T*f（T=2^m），m、n均为 ...

matlab 内部采样点个数就是2^n，这样也和FFT的要求一致；

如果你的数据长度不是2^n，matlab会自动截取，或者自动补充零以满足FFT对数据的要求；

加窗会衰减为零，如果长度不同，就会不合理乘以一个权系数，你觉得数据不会受到很大影响么？

simplebinbin

mxlzhenzhu 发表于 2014-8-1 20:21
matlab 内部采样频率就是2^n，这样也和FFT的要求一致；

如果你的数据长度不是2^n，matlab会自动截取 ...

我分析了一下，发现这个取值长度在对频谱分析幅值有一定影响（似乎是数据越短影响越大），但目前没有看到相关文献说对这个长度有严格的规定，为了使数据频谱分析准确，尽量使截取数据长度为2^n。

junethtwo

不错，加油哟

lq12801914

你好，请问代码中cf具体是什么意思？

mxlzhenzhu

lq12801914 发表于 2014-8-20 17:42
你好，请问代码中cf具体是什么意思？

81&82代码显示，是用来计权的，得到A计权结果。

小海豚zc

感谢分享，请教个问题：上面程序中 求1/3倍频程的程序 跟《MATLAB在振动信号处理中的应用》的一样，我用简单的正弦信号验证，为什么幅值是不准的？（不加窗且整周期采样）
程序如下

1. %三分之一倍频程处理
   
2. clear
   
3. %clc
   
4. %close all hidden
   
5. format long
   
6. 
   
7. sf=2000;
   
8. N=100000;
   
9. nn=1:N;
   
10. t=nn./sf;
    
11. x=2*sin(2*pi*2*t)+2*sqrt(2)*sin(2*pi*10*t)+2*sin(2*pi*30*t);
    
12. %定义三分之一倍频程的中心频率
    
13. f=[1.00 1.25 1.60 2.00 2.50 3.15 4.00 5.00 6.30 8.00];
    
14. fc=[f,10*f,100*f,1000*f,10000*f];
    
15. %中心频率与下限频率的比值
    
16. oc6=2^(1/6);
    
17. nc=length(fc);
    
18. n=length(x);
    
19. nfft=2^nextpow2(n);
    
20. a=fft(x,nfft);
    
21. for j=1:nc
    
22.     fl=fc(j)/oc6;
    
23.     fu=fc(j)*oc6;
    
24.     nl=round(fl*nfft/sf+1);
    
25.     nu=round(fu*nfft/sf+1);
    
26.     if fu>sf/2
    
27.         m=j-1;break
    
28.     end
    
29.     b=zeros(1,nfft);
    
30.     b(nl:nu)=a(nl:nu);
    
31.     b(nfft-nu+1:nfft-nl+1)=a(nfft-nu+1:nfft-nl+1);
    
32.     c=ifft(b,nfft);
    
33.     yc(j)=sqrt(var(real(c)));
    
34. end
    
35. figure
    
36. subplot(2,1,1);
    
37. t=0:1/sf:(n-1)/sf;
    
38. plot(t,x);
    
39. xlabel('时间(s)');
    
40. ylabel('加速度(g)');
    
41. grid on;
    
42. subplot(2,1,2);
    
43. plot(fc(1:m),yc(1:m));
    
44. xlabel('频率(Hz)');
    
45. ylabel('有效值');
    
46. xlim([0 100])
    
47. grid on;
    

复制代码

simplebinbin

小海豚zc 发表于 2015-8-3 14:14
感谢分享，请教个问题：上面程序中 求1/3倍频程的程序 跟《MATLAB在振动信号处理中的应用》的一样，我用简 ...

07.sf=2000;
08.N=100000;
你的采样频率不符合FFT分析法则，
采样频率sf=2^n
N=2^nn
这样就可以了
例如sf=128;N=1024;
采样频率一般设置为：2^n

小海豚zc

simplebinbin 发表于 2015-8-11 16:35
07.sf=2000;
08.N=100000;
你的采样频率不符合FFT分析法则，

确实，如您所说，改了之后幅值没有误差了。谢谢！！
那在实际应用时，有必要将样本长度取成2^n吗？还是说只是在fft时控制点数，例如 fft（xn,2^n）？当然前者精度是最好的吧

simplebinbin

小海豚zc 发表于 2015-8-17 09:20
确实，如您所说，改了之后幅值没有误差了。谢谢！！
那在实际应用时，有必要将样本长度取成2^n吗？还是 ...

对，如果条件允许情况下，采样长度取2^n比较好，“fft时控制点数”是截取数据的，会造成谱泄露。
但我目前常用后面一种方法，因为我实际处理信号时，没有遇到这种连续的周期振动/噪声信号，让我进行整周期采样。
就拿我程序中的声压分析来讲，如果分析时间长度为1s或者2s,则都是符合2^n,如果分析长度为3s（需要分析的噪声，刚好出现了3s),则就不符合2^n，在进行fft分析时，只能让其自动加0了，数据截取时，一般都需要加窗修正，减小谱泄露的影响
所以实际应用和理论还是有差别的