一、谱估计的分类 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_Pxx2=10*log10(Pxx2); figure(1)
plot(f,plot_Pxx); figure(2)
plot(f,plot_Pxx1); figure(3)
plot(f,plot_Pxx2);
2 现代功率谱估计
2.1 基于AR模型的功率谱
AR 模型又称为自回归模型, 它是一个全极点模型, 要利用AR 模型进行功率谱估计须通过 levinson _dubin 递推算法由正则方程求得 AR 的参数: a1 , a2 , ..., ap。在Mat lab 仿真中可调用 Pburg 函数直接画出基于 burg算法的功率谱估计的曲线图如图1所示。用周期图法求出的功率谱曲线和 burg 算法求出的AR 功率谱曲线( p= 50) , 其程序如下: fs= 200;
n= 0: 1/ fs: 1;
xn= cos( 2* pi* 40* n) + cos( 2* pi* 41* n) +34* co s( 2* pi+ 90* n) + 0. 1* randn( size( n) ) ; window = bo xcar ( 1eng th( x n) ) ; nfft= 512;
[ pxx , f] = perio do gr am( x n, w indow , nfft, fs) ; subplot( 121) ;
plot( f, 10+ log 10( pxx ) ) ; xlabel(‘ ffequency( hz)’) ;
ylabel(‘pow er spectr al density( Db/ H z)’) ; title(‘perio do gram PSE estimate ‘) ; or der l= 50; range= ‘ half ‘; magunits= ‘db ‘; subplot( 122) ;
pburg ( x n, order l, nfft, fs, range) ;
2.2基于MA模型的功率谱
基于MA 模型算法原理是由N 点数据x ( n) 建立一个 p 阶的AR模型, 从而求出p 阶AR系数a, 再利用AR系数建立线性预测,等效q阶AR模型,再利用求解AR 模型的出b,从而实现MA 功率谱估计. 程序如下: A1
A14= [ 1, a4 ??]
B14= fliplr ( conv ( fliplr ( B1) , fliplr ( A14 ) ) ) ; y 24 = filter
( B14, A1, randn( 1, N) ) ; %. * [ zer os( 1, 200) , ones( 1, 256) ] ; [ Ama4, Ema4] = arburg ( y24, 32) , B1
b4= arburg ( Ama4, 4) %- - - 求功率谱- - - % w= linspace( 0, pi, 512) ; %H1= fr eqz( B1, A1, w) H14= f reqz( b4, A14, w ) ; %Ps1= abs( H 1) . ^2; Py 14= abs( H14) . ^2; %if Py14> 200 % PPy14= 200; %elseif Py14< 200 % PPy14= Py14; %end
SPy 14= SPy14+ Py14;
VSPy14= VSPy14+ abs( Py 14) . ^2; fig ure( 4)
plot( w . / ( 2* pi) , Ps1, w. / ( 2* pi) , Py 14) ;
leg end(‘真实功率谱’,’20 次 ARM A( 4, 4) 的估计图’) ; hold on;
R8= zer os( 16, 8) ; %ARMA( 8, 8)的 R r8= zer os( 16, 1) ; %ARMA( 8, 8) 的 r for i= 1: 16
r8( i, 1) = - Ry ( 264+ i) ; for j= 1: 8
R8( i, j) = Ry( 264+ i- j) ; End End R8 r8
a8= inv ( R8 ’* R8) * R8 ’* r8%利用最小二乘法得到的 a的估计参数
2.3基于ARMA的功率谱
该算法原理是事先估计AR 参数,同时对已知数据x ( n)用FIR 滤波器, 将其输出近似 MA ( q) 过程, 并求解MA( q)参数的方法求出b,从而实现 ARMA 模型的参数估计。
下面的程序是给出对已知数据x ( n) ,用FIR 滤波器A = 1+ a( k) z ( - k)进行滤波,通过滤波器的输出模拟一个近似的MA 过程来估计MA 参数的方法: y= zeros( 1, 256) ; for i= 1: 256
y( i) = y 1( 200+ i) ; end
ny= [ 0: 255] ;
z= fliplr ( y ) ; nz= - fliplr( ny) ;
nb= ny( 1) + nz( 1) ; ne= ny( length( y) ) + nz( lengt h( z) ) ; n= [ nb: ne] ; Ry= conv( y, z) ;
R4= zero s( 8, 4) ; %ARM A( 4, 4) 的 R r4= zer os( 8, 1) ; %ARMA( 4, 4) 的 r
for i= 1: 8; r 4( i, 1) = - Ry ( 260+ i) ; fo r j= 1: 4; R4( i, j) =Ry( 260+ i- j) ; End End
R4 %R 矩阵 r4 %r矩阵
a4= inv( R4 ‘* R4) * R4 ‘* r4%利用最小二乘法得到的 a的估计参数
二、常见谱估计法的比较及适用情况
通过实验仿真可以直观地看出以下特性:( 1)功率谱估计中的相关函数法和周期图法所得到的结果是一致的, 其特点是离散性大, 曲线粗糙, 方差较大, 但是分辨率较高。 ( 2)平均周期图法和平滑平均周期图法的收敛性较好, 曲线平滑, 估计的结果方差较小, 但是功率谱主瓣较宽, 分辨率低。 这是由于对随机序列的分段处理引起了长度有限所带来的 Gibbs 现象而造成的。( 3)平滑平均周期图法与平均周期图法相比, 谱估值比较平滑, 但是分辨率较差。 其原因是给每一段序列用适当的窗口函数加权后, 在得到平滑的估计结果的同时, 使功率谱的主瓣变宽, 因此分辨率有所下降。
在自相关法中AR 功率谱曲线能基本区分间隔较远的信号,而对于间隔绩效的频率点分辨率和检出能力不好,必须通过提高阶次而达到效果从而提高对硬件的要求能力;而Burg 算法中通过改变, 其分辨率佳,在阶次低的情况下也能较好的分辨出间隔小的频率点;在MA 的功率谱曲线中,其分辨能力较差, 但由于它的平滑效果很好而得到利用。总之在 ARMA 模型中结 合AR, MA 算法的优点从而有较精确的谱估计合脚优良的谱分辨率性能,在功率谱曲线中可以很好地得到分辨和平滑效果,但它的实际参数估计比较麻烦, 在没有复杂的应用时,多数用AR 模型来估计功率谱参数。
经典功率谱估计的分辨率反比于有效信号的长度, 但现代谱估计的分辨率可以不受此限制。这是因为对于给定的 N点有限长序列 x(n), 虽然其估计出的自相关函数也是有限长的, 但是现代谱估计的一些隐含着数据和自相关函数的外推, 使其可能的长度超过给定的长度, 不象经典谱估计那样受窗函数的影响。因而现代谱的分别率比较高, 而且现代谱线要平滑得多, 从上图可以清楚看出。