谱估计的分类及应用.doc

一、谱估计的分类 1 经典功率谱估计 1.1 相关函数法(BT法) 该方法先由序列 x(n)估计出自相关函数 R(n), 然后对 R(n)进行傅立叶变换, 便得到 x(n)的功率谱估计。 当延迟与数据长度相比很小时,可以有良好的估计精度。 Matlab 代码示例 1: Fs=500;%采样频率 n=0:1/Fs:1; %产生含有噪声的序列 xn=cos(2*pi*40*n)+3*cos(2*pi*90*n)+randn(size(n)); nfft=512; 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)); figure(1)% plot(k,plot_Pxx); 1.2 周期图法( periodogram) 周期图法是把随机序列 x(n)的 N个观测数据视为一能量有限的序列, 直接计算 x(n)的离散傅立叶变换,得 X(k), 然后再取其幅值的平方, 并除以 N, 作为序列x(n)真实功率谱的估计。 Matlab 代码示例 2: Fs=600; %采样频率 n=0:1/Fs:1;%产生含有噪声的序列 xn=cos(2*pi*40*n)+3*cos(2*pi*90*n)+0.1*randn(size(n)); window=boxcar(length(xn));%矩形窗 nfft=512; [Pxx,f]=periodogram(xn,window,nfft,Fs);%直接法 plot(f,10*log10(Pxx)); window=boxcar(length(xn));%矩形窗 nfft=1024; [Pxx,f]=periodogram(xn,window,nfft,Fs); %直接法 figure(1) plot(f,10*log10(Pxx)); 对于周期图的功率谱估计, 当数据长度 N 太大时, 谱曲线起伏加剧, 若 N 太小, 谱的分辨率又不好,因此需要改进。两种改进的估计法是平均周期图法和平滑平均周期图法。 1.3 平均周期图法（Bartlett）: Bartlett 平均周期图的方法是将 N 点的有限长序列 x(n)分段求周期图再平均。 Matlab 代码示例 3: fs=600; n=0:1/fs:1; xn=cos(2*pi*20*n)+3*cos(2*pi*90*n)+randn(size(n)); nfft=512; window=hamming(nfft); %矩形窗 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)); plot_Pxxc=10*log10(Pxxc(index+1)); figure(1) plot(k,plot_Pxx); figure(2) plot(k,[plot_Pxx plot_Pxx- plot_Pxxc plot_Pxx+plot_Pxxc]); 1.4平滑平均周期图法（Welch 法）:对 Bartlett 法进行了两方面的修正, 一是选择适当的窗函数 w(n), 并在周期图计算前直接加进去, 加窗的优点是无论什么样的窗函数均可使谱估计非负。二是在分段时, 可使各段之间有重叠,这样会使方差减小。 Matlab 代码示例 4: Fs=600; n=0:1/Fs:1; xn=cos(2*pi*40*n)+3*cos(2*pi*90*n)+randn(size(n)); nfft=512; 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); plot_Pxx1=10*log10(Pxx1); plot_Pxx