方法一：%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=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:n))));

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(Lap(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);

方法二：

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=2^nextpow2(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 fu>fs/2

m=j-1;

break;

end

p=w(nl:nu);

lp=length(p);

k=0;

for ii=1:lp

if ii+2>lp

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);

end

Lp=10*log10(p2/p0^2);

for jj=1:length(Lp)

Lp1(jj)=10^(Lp(jj)/10);

end

Lpt=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);

end

ggg=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=2^nextpow2(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 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(b(1:n))));

end

%绘制输入时程曲线图形

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('有效值');

grid on;

%保存倍频程数据

fid=fopen(fno,’w ’);

for k=1:m;

fprintf(fid,’%f %f\n’,fc(k),yc(k));

end

status=fclose(fid);