%产生含有噪声的序列
xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n)); nfft=1024;
cxn=xcorr(xn,'unbiased'); %计算序列的自相关函数 CXk=fft(cxn,nfft); Pxx=abs(CXk); index=0:round(nfft/2-1); k=index*Fs/nfft;
plot_Pxx=10*log10(Pxx(index+1)); plot(k,plot_Pxx); xlabel('频率') ylabel('功率/DB')
35302520功率/DB151050-5050100150200250频率300350400450500
图7 直接法功率谱图
当M=N-1时,BT法与周期图法估计出的功率谱是一样的;当 M < N-1 时,BT 法的偏差大于周期图法,在窗函数满足一定条件时是渐进无偏估计;方差小于周期图的
11
方差;分辨率比周期图法低,与窗函数的选择有关。
BT法的缺点在于当 M →N 时, R(m)的方差很大,使谱估计质量下降;由R(m) 得到的S(w)不一定为正值,从而可能失去功率谱的物理意义。 3.3 改进的周期图法仿真
对于直接法的功率谱估计，当数据长度N太大时，谱曲线起伏加剧，若N太小，谱的分辨率又不好，因此需要改进。
对周期图法改进的思想是将信号分段进行估计，然后再将这些估计结果进行平均，从而减小估计的协方差，是功率谱图变得比直接法更平滑。增加分段数可以进一步降低估计的协方差，然而若每段中的数据点数太少，就会使估计的频率分辨率下降很多[8]。从样本信号序列总点数一定的条件下，可以采用使分段相互重叠的方法来增加分段数，从而保持每段信号点数不变，这样就在保证频率分辨率的前提下进一步降低估计协方差。主要的改进方法有Bartlett法和Welch法 3.3.1 Bartlett法
Bartlett平均周期图的方法是将N点的有限长序列x(n)分段求周期图再平均。功率谱图如图8所示。 Matlab代码示例： clear； Fs=1000; n=0:1/Fs:1;
xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n)); nfft=1024;
window=boxcar(length(n)); %矩形窗 noverlap=0; %数据无重叠 p=0.9; %置信概率
[Pxx,Pxxc]=psd(xn,nfft,Fs,window,noverlap,p); index=0:round(nfft/2-1); k=index*Fs/nfft;
plot_Pxx=10*log10(Pxx(index+1));
12
???
plot_Pxxc=10*log10(Pxxc(index+1)); figure(1) plot(k,plot_Pxx); pause; figure(2)
plot(k,[plot_Pxx plot_Pxx-plot_Pxxc plot_Pxx+plot_Pxxc]); xlabel('频率') ylabel('功率/DB')
35302520功率/DB151050-5050100150200250频率300350400450500
图8 bartlett法估计功率谱图
3.3.2 Welch法
现在比较常用的改进方法是Welch法，又叫加权交叠平均法[9]，简记为WOSA法，这种方法以加窗（加权）求取平滑，以分段重叠求得平均，因此集平均与平滑的优点于一体，同时也不可避免带有两者的缺点，因此归根到底是一种折中。其主要步骤是：
（1）将N长的数据段分成L个小段，每小段M点，相邻小段间交叠M/2点（即2:1分段）。因为L（M/2）+M/2=N,所以段数
L?N?M/2 (3-1)
M/213
（2）对各小段加同样的平滑窗后起傅氏变换
Xi(e)??xi(n)ax(n)e?jwn,i?1,...L, （3-2）
jwn?0M?1（3）用下式求各小段功率谱的平均
1l1jw2Si(e)??Xi(e) （3-3）
Li?1MU?jw这里，代表窗函数平均功率，MU是M长窗函数的能量。仍以上
述平稳随机信号为例,采样频率、采样点数、FFT点数、窗长度及重叠数据不变,窗函数采用矩形窗、Blackman 窗、 Hamming窗，仿真结果如图所示。图9是加矩形窗得到的功率谱图；图10是加海明窗所得的功率谱图；图11是加blackman窗所得到的功率谱图。
Matlab代码示例： clear; Fs=1000; n=0:1/Fs:1;
xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n)); nfft=1024;
window=boxcar(100); %矩形窗 window1=hamming(100); %海明窗 window2=blackman(100); %blackman窗 noverlap=20; %数据无重叠
range='half'; %频率间隔为[0 Fs/2]，只计算一半的频率 [Pxx,f]=pwelch(xn,window,noverlap,nfft,Fs,range); [Pxx1,f]=pwelch(xn,window1,noverlap,nfft,Fs,range); [Pxx2,f]=pwelch(xn,window2,noverlap,nfft,Fs,range); plot_Pxx=10*log10(Pxx);
14
plot_Pxx1=10*log10(Pxx1); plot_Pxx2=10*log10(Pxx2); figure(1) plot(f,plot_Pxx); xlabel('频率') ylabel('功率/DB') pause; figure(2) plot(f,plot_Pxx1); xlabel('频率') ylabel('功率/DB') pause; figure(3) plot(f,plot_Pxx2); xlabel('频率') ylabel('功率/DB')
0-5-10功率/DB-15-20-25-30-35050100150200250频率300350400450500
图9 矩形窗
15