杭州锐达数字技术有限公司
查看: 7343|回复: 19

[声学测量] A声级及1/3倍频程计算matlab程序

[复制链接]
发表于 2014-5-14 10:05 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有帐号?我要加入

x
小弟最近入门声学分析,用声级计采集了瞬时声压数据,自己按照振动分析的代码写了个程序,
贴出来供大家参考,一起讨论学习交流!
通过验证,与声级计分析出来的A计权总声压级在有的时间段差别很小0.5dB(A),有的时间段差别比较大,好几个dB(A),
希望大家一起完善!
  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. b=zeros(1,nfft);
  44. b(nl:nu)=a(nl:nu);
  45. b(nfft-nu+1:nfft-nl+1)=a(nfft-nu+1:nfft-nl+1);
  46. c=ifft(b,nfft);
  47. yc(j)=sqrt(var(real(c(1:nnn))));

  48. end

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

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

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

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

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


复制代码

评分

1

查看全部评分

本帖被以下淘专辑推荐:

回复
分享到:

使用道具 举报

 楼主| 发表于 2014-5-14 10:13 | 显示全部楼层
代码第90行应该是:90.bar(Lap(1:nc));
发表于 2014-7-28 15:55 | 显示全部楼层
谢谢分享哦
发表于 2014-7-30 00:10 | 显示全部楼层
本帖最后由 mxlzhenzhu 于 2014-7-30 00:12 编辑

总体思路我看了,多谢分享;

关于其中一个问题:


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); %汉宁窗


发表于 2014-7-30 00:13 | 显示全部楼层
这样才是整周期采样
 楼主| 发表于 2014-8-1 11:52 | 显示全部楼层
mxlzhenzhu 发表于 2014-7-30 00:13
这样才是整周期采样

与整周期采样关系不大
你这样只能处理特定的数据,即采样频率为f=2^n,数据长度为T*f(T=2^m),m、n均为正整数
缺失了程序的通用性
发表于 2014-8-1 20:21 | 显示全部楼层
本帖最后由 mxlzhenzhu 于 2014-8-5 00:22 编辑
simplebinbin 发表于 2014-8-1 11:52
与整周期采样关系不大
你这样只能处理特定的数据,即采样频率为f=2^n,数据长度为T*f(T=2^m),m、n均为 ...


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

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

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

 楼主| 发表于 2014-8-4 14:18 | 显示全部楼层
mxlzhenzhu 发表于 2014-8-1 20:21
matlab 内部采样频率就是2^n,这样也和FFT的要求一致;

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

我分析了一下,发现这个取值长度在对频谱分析幅值有一定影响(似乎是数据越短影响越大),但目前没有看到相关文献说对这个长度有严格的规定,为了使数据频谱分析准确,尽量使截取数据长度为2^n。
发表于 2014-8-12 19:32 | 显示全部楼层
不错,加油哟
发表于 2014-8-20 17:42 | 显示全部楼层
你好,请问代码中cf具体是什么意思?
发表于 2014-8-21 15:59 | 显示全部楼层
lq12801914 发表于 2014-8-20 17:42
你好,请问代码中cf具体是什么意思?

81&82代码显示,是用来计权的,得到A计权结果。
发表于 2015-8-3 14:14 | 显示全部楼层
感谢分享,请教个问题:上面程序中 求1/3倍频程的程序 跟《MATLAB在振动信号处理中的应用》的一样,我用简单的正弦信号验证,为什么幅值是不准的?(不加窗且整周期采样)
程序如下

  1. %三分之一倍频程处理
  2. clear
  3. %clc
  4. %close all hidden
  5. format long

  6. sf=2000;
  7. N=100000;
  8. nn=1:N;
  9. t=nn./sf;
  10. x=2*sin(2*pi*2*t)+2*sqrt(2)*sin(2*pi*10*t)+2*sin(2*pi*30*t);
  11. %定义三分之一倍频程的中心频率
  12. f=[1.00 1.25 1.60 2.00 2.50 3.15 4.00 5.00 6.30 8.00];
  13. fc=[f,10*f,100*f,1000*f,10000*f];
  14. %中心频率与下限频率的比值
  15. oc6=2^(1/6);
  16. nc=length(fc);
  17. n=length(x);
  18. nfft=2^nextpow2(n);
  19. a=fft(x,nfft);
  20. for j=1:nc
  21.     fl=fc(j)/oc6;
  22.     fu=fc(j)*oc6;
  23.     nl=round(fl*nfft/sf+1);
  24.     nu=round(fu*nfft/sf+1);
  25.     if fu>sf/2
  26.         m=j-1;break
  27.     end
  28.     b=zeros(1,nfft);
  29.     b(nl:nu)=a(nl:nu);
  30.     b(nfft-nu+1:nfft-nl+1)=a(nfft-nu+1:nfft-nl+1);
  31.     c=ifft(b,nfft);
  32.     yc(j)=sqrt(var(real(c)));
  33. end
  34. figure
  35. subplot(2,1,1);
  36. t=0:1/sf:(n-1)/sf;
  37. plot(t,x);
  38. xlabel('时间(s)');
  39. ylabel('加速度(g)');
  40. grid on;
  41. subplot(2,1,2);
  42. plot(fc(1:m),yc(1:m));
  43. xlabel('频率(Hz)');
  44. ylabel('有效值');
  45. xlim([0 100])
  46. grid on;
复制代码
 楼主| 发表于 2015-8-11 16:35 | 显示全部楼层
小海豚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
发表于 2015-8-17 09:20 | 显示全部楼层
simplebinbin 发表于 2015-8-11 16:35
07.sf=2000;
08.N=100000;
你的采样频率不符合FFT分析法则,

确实,如您所说,改了之后幅值没有误差了。谢谢!!
那在实际应用时,有必要将样本长度取成2^n吗?还是说只是在fft时控制点数,例如 fft(xn,2^n)?当然前者精度是最好的吧
 楼主| 发表于 2015-8-18 14:09 | 显示全部楼层
小海豚zc 发表于 2015-8-17 09:20
确实,如您所说,改了之后幅值没有误差了。谢谢!!
那在实际应用时,有必要将样本长度取成2^n吗?还是 ...

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

本版积分规则

快速回复 返回顶部 返回列表