%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];

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

總結

以上是生活随笔為你收集整理的matlab1 3倍频程,瞬时声压时域数据怎么用matlab进行1/3倍频程声压级分析的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。