基于Matlab的1/3倍频程计算

倍频程定义

倍频程表示的是相应倍频程带内的功率（能量）之和，其详细理论课参考其他文章，在计算倍频程之前需要将时域信号进行傅里叶变化，本文基于功率谱密度进行倍频程计算。

计算自功率谱密度

计算一维数组x的自功率谱密度：

n=length(x);
if mod(n,2)==0
	%偶数
else 
	n=n-1;
end
windows=hann(n);
x1=detrend(x); %去除线性趋势项
x1=x1(1:n).*windows; %时域数据加汉宁窗
y=fft(x1,n);
f=(0:(n-1)/2)/n*fs; %根据对称性取频率一半，fs为采样周期
P1=y(1:n/2);
P1(2:end-1)=2*P1(2:end-1); %得到单边复数频谱
P1=P1/n; %单边复数频谱的真实值
AP=P1.*conj(P1); %计算自功率谱
APL=AP/(fs/n); %自功率谱密度
APL=APL*1.6333; %加窗造成的衰减，需要乘以汉宁窗恢复系数。

计算倍频程

ajrms=zeros(1,23);
for i=1:23
	fdnum=find(f>=param.fid(i));
	%param.fid为各倍频程带的下限频率
	funum=find(f<param.fiu(i));
	%param.fid为各倍频程带的上限频率
	api=APL(fdnum(1):funum(end));
	aptrapz=trapz(fs/n,api);
	ajrms(1,i)=aptrapz^0.5;
end
%ajrms即为各倍频程带的倍频程

倍频程的计算和噪声倍频程或许存在出入，本案例参照GB 4970汽车平顺性的倍频程计算方法计算。