A声级及1/3倍频程计算matlab程序
小弟最近入门声学分析，用声级计采集了瞬时声压数据，自己按照振动分析的代码写了个程序，贴出来供大家参考，一起讨论学习交流！
通过验证，与声级计分析出来的A计权总声压级在有的时间段差别很小0.5dB（A),有的时间段差别比较大，好几个dB(A)，
希望大家一起完善！
%A计权声压级频谱分析
clc;
clear;
close all;
%时域分析
y=wavread('abc.wav');
%频域分析
fs=51200;%采样频率
p0=2e-5;%参考声压
f=; %基准中心频率
f1=;
fc=; %%%%%%%%%中心频率%%%%%%%%
%20-16000Hz A声级计权值
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];
x=y(t1*fs:t2*fs);%截取需要处理的数据段
n=length(x);
t=(0:1/fs:(n-1)/fs);
subplot(221);
plot(t,x);%瞬时声压时程图
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);
w=2*abs(a(1:nfft/2)/n);
subplot(222）；
plot(f,w);%绘制频谱图
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%1/3倍频程计算
oc6=2^(1/6);
nc=length(cf);
%下面这个求1/3倍频程的程序是按照振动振级计算那个来的
for j=1:nc
fl=fc(j)/oc6;
fu=fc(j)*oc6;
nl=round(fl*nfft/fs+1);
nu=round(fu*nfft/fs+1);
if fu>fs/2
m=j-1;
break;
end
b=zeros(1,nfft);
b(nl:nu)=a(nl:nu);
b(nfft-nu+1:nfft-nl+1)=a(nfft-nu+1:nfft-nl+1);
c=ifft(b,nfft);
yc(j)=sqrt(var(real(c(1:nnn))));
end
aj_sumn=0;
for i=1:nc
Lp1(i)=20*log10(yc(i)/p0);%未计权1/3倍频程声压级
end
%%%%%
for jj=1:nc
aj_sumn=aj_sumn+10^(0.1*Lp1(j));
end
Lp=10*log10(aj_sumn);%未计权总声压级
subplot(223);%绘制未计权1/3倍频程声压级图谱
bar(Lp1(1:nc));
gg=zeros(1,nc);
for i=1:nc
gg(1:nc)=fc(1:nc);
end
ggg=1:nc;
set(gca,'xtick',ggg);
set(gca,'xticklabel',gg);
%%%%%A计权1/3倍频程声压级
Lap=Lp1+cf;
aj_sum=0;
for j=1:nc
aj_sum=aj_sum+10^(0.1*Lap(j));
end
LA=10*log10(aj_sum);%Aa计权总声压级
subplot(224);%绘制A计权1/3倍频程声压级图谱
bar(Lp1(1:nc));
gg=zeros(1,nc);
for i=1:nc
gg(1:nc)=fc(1:nc);
end
ggg=1:nc;
set(gca,'xtick',ggg);
set(gca,'xticklabel',gg);
代码第90行应该是：90.bar(Lap(1:nc)); 谢谢分享哦 {:{03}:} 本帖最后由 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); %汉宁窗
这样才是整周期采样 mxlzhenzhu 发表于 2014-7-30 00:13
这样才是整周期采样
与整周期采样关系不大
你这样只能处理特定的数据，即采样频率为f=2^n,数据长度为T*f（T=2^m），m、n均为正整数
缺失了程序的通用性 本帖最后由 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对数据的要求；
加窗会衰减为零，如果长度不同，就会不合理乘以一个权系数，你觉得数据不会受到很大影响么？
mxlzhenzhu 发表于 2014-8-1 20:21
matlab 内部采样频率就是2^n，这样也和FFT的要求一致；
如果你的数据长度不是2^n，matlab会自动截取 ...
我分析了一下，发现这个取值长度在对频谱分析幅值有一定影响（似乎是数据越短影响越大），但目前没有看到相关文献说对这个长度有严格的规定，为了使数据频谱分析准确，尽量使截取数据长度为2^n。
不错，加油哟 你好，请问代码中cf具体是什么意思？ lq12801914 发表于 2014-8-20 17:42
你好，请问代码中cf具体是什么意思？
81&82代码显示，是用来计权的，得到A计权结果。 感谢分享，请教个问题：上面程序中 求1/3倍频程的程序 跟《MATLAB在振动信号处理中的应用》的一样，我用简单的正弦信号验证，为什么幅值是不准的？（不加窗且整周期采样）
程序如下
%三分之一倍频程处理
clear
%clc
%close all hidden
format long
sf=2000;
N=100000;
nn=1:N;
t=nn./sf;
x=2*sin(2*pi*2*t)+2*sqrt(2)*sin(2*pi*10*t)+2*sin(2*pi*30*t);
%定义三分之一倍频程的中心频率
f=;
fc=;
%中心频率与下限频率的比值
oc6=2^(1/6);
nc=length(fc);
n=length(x);
nfft=2^nextpow2(n);
a=fft(x,nfft);
for j=1:nc
fl=fc(j)/oc6;
fu=fc(j)*oc6;
nl=round(fl*nfft/sf+1);
nu=round(fu*nfft/sf+1);
if fu>sf/2
m=j-1;break
end
b=zeros(1,nfft);
b(nl:nu)=a(nl:nu);
b(nfft-nu+1:nfft-nl+1)=a(nfft-nu+1:nfft-nl+1);
c=ifft(b,nfft);
yc(j)=sqrt(var(real(c)));
end
figure
subplot(2,1,1);
t=0:1/sf:(n-1)/sf;
plot(t,x);
xlabel('时间(s)');
ylabel('加速度(g)');
grid on;
subplot(2,1,2);
plot(fc(1:m),yc(1:m));
xlabel('频率(Hz)');
ylabel('有效值');
xlim()
grid on;
小海豚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 simplebinbin 发表于 2015-8-11 16:35
07.sf=2000;
08.N=100000;
你的采样频率不符合FFT分析法则，
确实，如您所说，改了之后幅值没有误差了。谢谢！！
那在实际应用时，有必要将样本长度取成2^n吗？还是说只是在fft时控制点数，例如 fft（xn,2^n）？当然前者精度是最好的吧 小海豚zc 发表于 2015-8-17 09:20
确实，如您所说，改了之后幅值没有误差了。谢谢！！
那在实际应用时，有必要将样本长度取成2^n吗？还是 ...
对，如果条件允许情况下，采样长度取2^n比较好，“fft时控制点数”是截取数据的，会造成谱泄露。
但我目前常用后面一种方法，因为我实际处理信号时，没有遇到这种连续的周期振动/噪声信号，让我进行整周期采样。
就拿我程序中的声压分析来讲，如果分析时间长度为1s或者2s,则都是符合2^n,如果分析长度为3s（需要分析的噪声，刚好出现了3s),则就不符合2^n，在进行fft分析时，只能让其自动加0了，数据截取时，一般都需要加窗修正，减小谱泄露的影响
所以实际应用和理论还是有差别的
页:
[1]
2