三分之一倍频程程序.doc

方法一：%A计权声压级频谱分析clc;clear;close all;y=wavread(abc.wav);fs=51200;%采样频率p0=2e-5;%参考声压f=1.00 1.25 1.600 2.00 2.50 3.15 4.00 5.00 6.30 8.0; %基准中心频率f1=20.00 25.0 31.5 40.0 50.0 63.0 80;fc=f1,100*f,1000*f,10000*f; %中心频率%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; t1=1;t2=2;x=y(t1*fs+1: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=2nextpow2(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 fufs/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:n); endaj_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);endLp=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);endggg=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);endLA=10*log10(aj_sum);%Aa计权总声压级subplot(224);%绘制A计权1/3倍频程声压级图谱bar(Lap(1:nc);gg=zeros(1,nc);for i=1:nc gg(1:nc)=fc(1:nc);endggg=1:nc; set(gca,xtick,ggg);set(gca,xticklabel,gg); 方法二：clc;clear;close all;%时域分析y=wavread(abc.wav);%频域分析fs=51200;%采样频率p0=2e-5;%参考声压f=1.00 1.25 1.600 2.00 2.50 3.15 4.00 5.00 6.30 8.0; %基准中心频率f1=20.00 25.0 31.5 40.0 50.0 63.0 80;fc=f1,100*f,1000*f,10000*f; %中心频率%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; n=length(y);t=(0:1/fs:(n-1)/fs);h1=figure;plot(t,y);title(瞬时声压时程);xlabel(Time(s);ylabel(Sound Presure Value(Pa);%t1=0;t2=4;x=y(t1*fs+1:t2*fs);%截取需要处理的数据段n=length(x);%t=(0:1/fs:(n-1)/fs);%plot(t,x);%瞬时声压时程图w=1.633*hanning(n); %汉宁窗(恢复系数为1.633)%w=1.812*blackmanharris(n); %布拉克曼窗(功率相等恢复系数1.812)xx=x.*w; %加汉宁窗nfft=2nextpow2(n); %nextpow2(n)-取最接近的较大2次幂a = fft(xx,nfft)/n;%f = fs/2*linspace(0,1,nfft/2);w=2*abs(a(1:nfft/2);oc6=2(1/6);nc=length(cf);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 fufs/2 m=j-1; break; end p=w(nl:nu); lp=length(p); k=0; for ii=1:lp if ii+2lp break end if p(ii+1)p(ii)&p(ii+1)p(ii+2) k=k+1; pp(k)=p(ii+1)/sqrt(2); end end p2(j)=sum(pp.*pp); endLp=10*log10(p2/p02);for jj=1:length(Lp) Lp1(jj)=10(Lp(jj)/10);endLpt=10*log10(sum(Lp1)h2=figure;mm=nc;bar(Lp(1:mm);gg=zeros(1,mm);for i=1:mm gg(1:mm)=fc(1:mm);endggg=1:mm; set(gca,xtick,ggg);set(gca,xticklabel,gg);set(gcf,PaperPosition,1,1,40,20)set(gca,fontsize,10)xlabel(Frequency(Hz);ylabel(SPL(dB);title(未计权声压级);grid on;方法三：% 三分之一倍频程处理clear;clc;close all;%s = xlsread(ay.xls);%输入时程数据sf=256; %采样频率x=s(:,2); %定义三分之一倍频程的中心频率f=1.00 1.25 1.60 2.00 2.50 3.15 4.00 5.00 6.30 8.00;fc=f,10*f,100*f,1000*f,10000*f;%中心频率与下限频率的比值oc6=2(1/6);%取中心频率总的长度nc=length(fc);%输入数据的长度n=length(x);%大于并接近n的2的幂次方长度nfft=2nextpow2(n);%FFT变换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 fusf/2m=j-1;breakend%以每个中心频率段为通带进行带通频率滤波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(b(1:n);end%绘制输入时程曲线图形subplot(2,1,1);t=0:1/sf:(n-1)/sf;plot(t,x);xlabel(时间(s);ylabe