用FLUENT中的FW-H噪声板块计算建筑的气动噪声，现在利用大涡模拟得到了监测点的瞬时风压时程，接下来要计算这个监测点的总声压级和A计权声压级，在网上找了个MATLAB程序，自己修改整合了之后的程序如下：
%输入需要分析的时程数据
x=importdata('xx.txt');
fs=1000;%采样频率
p0=2*10^(-5);%基准声压
y=x(:,2);%声压
n=length(y);%输入数据的长度
t=0:1/fs:(n-1)/fs;%时间
%%%%%%绘制声压时程图
figure(1);
plot(t,y);
xlabel('时间(s)');
ylabel('声压(pa)');
title('声压时程图');
grid on;
%定义三分之一倍频程的中心频率
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,12500,16000,20000];
%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,-9.3];
%加汉宁窗
w=hanning(n);
yy=1.633*y.*w;
%大于并接近n的2的幂次方长度
nfft=2^nextpow2(n);
%FFT变换
a=fft(yy,nfft);
%0~fs/2之间的频率分布
f=fs/2*linspace(0,1,nfft/2);
for i=1:nfft/2
p(i)=2*abs(a(i))/n;
end
%%%%%绘制声压频谱图
figure(2);
plot(f,p);
xlabel('频率(Hz)');
ylabel('声压级(dB)');
title('声压频谱图');
grid on;
%将声压转换成声压级
for i=1:nfft/2
Lp1(i)=20*log10((p(i)/p0));
end
%%%%%绘制声压级频谱图
figure(3);
plot(f,Lp1);
xlabel('频率(Hz)');
ylabel('声压级(dB)');
title('声压级频谱图');
grid on;
Lp_sum1=0;
for j=1:(n+1)/2
Lp_sum1=Lp_sum1+10^(0.1*Lp1(j));
end
%输出总声压级
fprintf('总声压级为：')
OSPL1=10*log10(Lp_sum1)
%%%%%%%%%%%%%三分之一倍频%%%%%%%%%%%
%中心频率与下限频率的比值
oc6=2^(1/6);
%中心频率长度
nc=length(fc);
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
%三分之一倍频对应的声压级
for i=1:m
Lp2(i)=20*log10(yc(i)/p0);
end
Lp_sum2=0;
for j=1:m
Lp_sum2=Lp_sum2+10^(0.1*Lp2(j));
end
%三分之一倍频对应的未计权总声压级
fprintf('三分之一倍频对应的未计权总声压级：')
OSPL2=10*log10(Lp_sum2)
%三分之一倍频对应的A计权声压级
Lp_A=Lp2+cf(1:m);
Lp_sumA=0;
for j=1:m
Lp_sumA=Lp_sumA+10^(0.1*Lp_A(j));
end
%A计权总声压级
fprintf('三分之一倍频对应的A计权总声压级：')
OSPL_A=10*log10(Lp_sumA)
%%%%%%绘制A计权声压级频谱图
figure(4);
bar(Lp_A(1:m));
gg=zeros(1,m);
for i=1:m
gg(1:m)=fc(1:m);
end
ggg=1:m;
set(gca,'xtick',ggg);%划分x轴的刻度
set(gca,'xticklabel',gg);%设置x轴刻度对应的数值
xlabel('A计权中心频率(Hz)');
ylabel('A计权声压级(dB)');
title('A计权声压级频谱图');
grid on;
现在的问题是，我的来流风速是7m/s，用这个程序计算得到的总声压级达到了110dB以上，感觉这个数太大了(貌似分贝数达到145人就会晕过去。。)，不知道程序有没有什么问题？另，总声压级和A计权声压级相差多少比较正常呢？(看文献上两者相差20+，但是我算出来的仅相差10+)麻烦高手帮忙看下我的程序！急！急！急！ 小女子先行谢过！！
