《谱估计的分类及应用》由会员分享，可在线阅读，更多相关《谱估计的分类及应用（8页珍藏版）》请在金锄头文库上搜索。
1、本文格式为Word版，下载可任意编辑谱估计的分类及应用 一、谱估计的分类 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
2、=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;%产生含有
3、噪声的序列 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
4、 太大时, 谱曲线起伏加剧, 若 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,
5、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), 并在周期图计算
6、前直接加进去, 加窗的优点是无论什么样的窗函数均可使谱估计非负。二是在分段时, 可使各段之间有重叠,这样会使方差减小。 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
7、=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_Px
8、x2); 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)
9、 +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 estim
10、ate ) ; 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 )
11、) ) ; 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 %
12、 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=
13、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)举行滤波
14、,通过滤波器的输出模拟一个近似的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) ; %
15、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)平均周期图法和平滑平均周期图法的收敛性较好, 曲线平滑, 估计的结果方差较小, 但是功率谱主瓣较宽, 辨识率低。 这是由于对随机序列的
