ch6_4现代功率谱估计

1、数字信号处理(Digital Signal Processing)信号与系统系列课程组国家电工电子教学基地离散随机序列的特征描述平稳随机序列通过LTI系统经典功率谱估计现代功率谱估计现代谱估计简介问题提出平稳随机信号的参数模型AR模型参数与自相关函数的关系Y-W方程的L-D递推算法AR模型参数与线性预测滤波器的关系伯格(Burg)递推算法功率谱估计问题提出经典法存在的主要问题:方差性能不好，不是Px(W)的一致估计平滑周期图和平均周期图改善了周期图的方差 性能，但却降低了谱分辨率和增大了偏差。1.2.3.可能使短序列的功率谱估计出现错误的结果出现问题的主要原因:将观测数据以外的数据一律视为零，

2、与实际不符。功率谱估计参数模型法的基本思想根据所研究信号的先验知识，对观测数据以外的数据作出某种比较合理的假设。假设信号是白噪声通过LTI系统产生的。由观测数 据估计LTI系统模型的参数。最后由LTI系统模型的参数得出功率谱。hky kRn = s 2d n输入白噪声的自相关函数h2 sPy (W ) =H (ejW )2h平稳随机信号的参数模型AR模型(AutoRegressive model )11H ( z ) =pA( z )1 + an= 0z - nnMA模型(Moving Average model)ql =1H ( z ) = 1 +- lb zlql = 0- lb z= B

3、( z )lARMA模型(AutoRegressive- MovingAverage model )H ( z ) =pA( z )n= 01 +z - nan功率谱估计平稳随机信号的参数模型若输入白噪声的功率谱P(W ) = s2h则输出序列的功率谱为22P(W ) =P(W ) = sH (e jW )H (e jW2)hx功率谱估计若能确定模型中各参数an, bn和s2就可以求得功率谱Px(W)AR模型参数与自相关函数的关系描述p阶AR模型的差分方程为pyk + ap (n) yk - n =hk n=1将上式两端同乘以yk-m再求数学期望pEyk - myk + ap(n)Eyk -

4、myk - n = Eyk - mhkn=1Rym - nRymRyhmpRy m + n=1ap (n)Ry m - n = Ryh m功率谱估计AR模型参数与自相关函数的关系pRy m + a p (n)Ry m - n = Ryh mn=1由于输入是白噪声信号Rhn=s2dnRyh m = h- m * Rh m = s2 h-m对于因果系统, p阶AR模型的自相关函数与模型参数的关系pRm + (n) Rm - n = s 2d mm = 0,1,L, paypyn=1Yule-Walker(Y-W) 方程功率谱估计功率谱估计AR模型参数与自相关函数的关系Y-W 方程的矩阵表示pRm

5、+ m - n = s 2d mm = 0,1,L, pa(n) Rypyn=1例: p=3 时的Y-W 方程 Ry 0Ry 3sRy 1Ry 212 RR2 a(1) 1R0R1000 = yyyy3 Ry 2Ry 1a3 (2)Ry 1Ry 0 RR0a(3)3R2R1 yyyy3若已知Ryn ，由Y-W方程解出各参数a3(1), a3 (2),a3 (3), s2。则可由AR模型参数获得功率谱Py(W)的估计值。Y-W方程的L-D递推算法一阶Y-W方程的解 Ry 0 Ry 1Ry 1 s121 = Ry 0a1 (1)0a(1) = - Ry 1解此方程得1Ry 02 )s= R0 +

6、R1a(1) = R0(1- a(1)21yy1y1功率谱估计Y-W方程的L-D递推算法二阶Y-W方程的解 RR2s220R11yyyRy 1 a2 (1) = Ry 1Ry 000 RR0a(2)2R1yyy22 -2R2 + aR0RR1(1) R1a(2) = -yyy= -y1y2R0 - Rs22y121y1 - RR0R1R2yyyy(1) = -= a(1) + aa(2)a(1)2121R0 -22R1yy2 ss(2)= 1 -a= R0 + R1a(1) + R2122(2)2a功率谱y 估计2yy22Y-W方程的L-D递推算法p阶Y-W方程的递推解p -1Ry p + a

7、p -1 (n) Ry p - na( p) = - n=1ps2p -1ap (n) = ap-1 (n) + ap ( p)ap-1 ( p - n)(n = 1,2,L, p -1)2ss= 1 -2p2p -1ap ( p)功率谱估计Y-W方程的L-D递推算法L-D算法估计功率谱的步骤(1)计算自相关函数的估计值(2)求解一阶模型参数关函数的估计值Ry 1a1 (1) = - Rs2 )= R0(1- a21(1)0y1y功率谱估计Y-W方程的L-D递推算法L-D算法估计功率谱的步骤(3)由递推算法求解p阶模型参数Ry p + ap -1 (n) Ry p - np -1a( p) =

8、 - n=1ps2p -1ap (n) = ap-1 (n) + ap ( p)ap-1 ( p - n)(n = 1,2,L, p -1)2ss= 1 -2p2p -1ap( p)(4)求出功率谱估计s2(W ) =PAR2p1+ a(n)e- jWnp功n=率1谱估计利用L-D算法估计频谱的MATLAB函数Pxx,f = pyulear(x,p,NFFT,Fs)x：进行功率谱估计的输入有限长序列；p:模型的阶数NFFT：DFT的点数；Fs ：绘制功率谱曲线的抽样频率，默认值为1；P率谱估计值；f：Pxx值所对应的频率点例：利用L-D算法进行谱估计一序列含有白噪声和两个频率的余弦信号，xk

9、= cos(0.3k ) + cos(0.4k ) + hk 利用L-D算法估计该序列的功率谱。N=128;p=40; NFFT=2048;Fs=2;n=0:N-1;randn(state,0); x=cos(0.3*pi*n)+cos(0.4*pi*n)+randn(size(n); P,f=periodogram(x,NFFT,2) ;Py,fy=pyulear(x,p,NFFT,2); subplot(211);plot(f,10*log(P);grid;title( Periodogram ); axis(0 1 -60 60);subplot(212);plot(fy,10*log(

10、Py);grid;title( Yule); axis(0 1 -60 60);谱估计结果p=40,N=128Periodogram6040200-20-40-6000.10.20.30.40.50.60.70.80.91Yule6040200-20-40-6000.10.20.30.40.50.60.70.80.91AR模型参数与前向线性预测滤波器的关系前向线性预测滤波器例:前向线性预测(p=2)yk = -a2 (1) yk -1 - a2 (2) yk - 2前向预测误差yk k = yk - yk = yk + a(1) yk -1 + a(2) yk - 2ef222前向预测误差滤波

11、器系统函数1+ a2 (1) z-1+ a2 (2) z-2ykA(z)efk2功率谱估计yk-2yk-1AR模型参数与前向线性预测滤波器的关系均方误差最小的前向预测误差滤波器k = yk + a(1) yk -1 + a(2) yk - 2ef222= min Eek 2 E ff22正交准则Eek yk - n= 0,n = 1,2f若2则估计的均方误差达到最小。 且=min Eek = Eek yk E ff2f222功率谱估计AR模型参数与前向线性预测滤波器的关系预测系数及均方误差的确定k = yk + a(1) yk -1 + a(2) yk - 2ef222Eek yk -1=R1

12、 + a(1) R0 + a(2) R1 = 0fy2y2y2Eek yk - 2=Ry 2 + a2 (1) Ry 1 + a2 (2) Ry 0 = 0f2= R0 + a(1) R1 + a=E ff(2) R2Eek yk y2y2y22 RR2 E f0R11yyy2Ry 1 a2 (1) = 0 Ry 1Ry 0 RR0a(2)2R10yyy2功率谱估计AR模型参数与后向线性预测滤波器的关系后向线性预测滤波器例:后向线性预测(p=2)yk - 2 = -a2 (1) yk -1 - a2 (2) yk yk - 2后向预测误差= yk - 2 + a(1) yk -1 + aek

13、= yk - 2 -yk - 2b(2) yk 222后向预测误差滤波器系统函数E b ( z)(2) + a-+ z-=Ab ( z) = a12 2Y ( z)(1) z22功率谱估计yk-1ykAR模型参数与后向线性预测滤波器的关系预测系数及均方误差的确定E=ek yk - nn = 0,1b0,由2=min Eek =k yk - 2E bb2b 2Ee22可得 RR2 E b0R11yyy200Ry 1 a2 (1) = Ry 1Ry 0 RR0a(2)2R1yyy2比较前向预测系统和后向预测系统的方程后可得系数a2(1),a2(2)满足的方程相同= E= sE fb222功率A谱b

14、(估z)计=z -2A(z)在计算相关函数估计时，对N个观测数据以外的数据作零的假设，故谱估计误差较大。伯格(Burg)递推算法基本思想：直接从观测的数据利用线性预测器的前向和后向预测的总均方误差之和为最小的准则来估计反射系数，进而通过L-D算法的递推公式求出AR 模型的优化参数。功率谱估计伯格(Burg)递推算法前向预测误差的递推公式e k =yk + a (1)yk - 1 + a (2)yk - 2f2阶前向预测误差1阶后向预测误差L-D算法的递推公式222e k = yk - 1 + a (1)ykb11a2(1) = a1(1) + a2(2)a1(1)e k = yk +a (1)

15、 + a (2)a (1) yk - 1 + a (2)yk - 2f21212= yk + a (1)yk - 1 + a (2)(yk - 2 + a (1)yk - 1)12e k +1e k =K e k - 1ffb212 1K2=a2(2)2阶预测器的反射系数功率谱估计伯格(Burg)递推算法预测误差的递推公式一般地e k =ek +K ek - 1ffbp-1pp-1p同理可得后向预测误差的递推公式e k =ek - 1 +bbfK ekp-1pp-1pKp=ap(p)为p阶预测器的反射系数。功率谱估计伯格(Burg)递推算法预测误差滤波器的格形结构e k =ek +K ek -

16、 1ffbp-1pp-1pe k =- 1 + K ekbbfekp-1pp-1pe fe fk k e fe fykkkp012 K pK1K2K pK1K2z -1z -1z-1eb k eb k eb k eb k012p功率谱估计伯格(Burg)递推算法反射系数Kp 的确定前向和后向均方预测误差的总和为N -1E= ef k2 +b2e k pppk = p由EekfekbN -1=+ = 0ppp2ekf2ekbppKKfKk = ppppN -1- 2可得k - 1bp-1ekep-1K=k = ppN -1efk2 + ek - 12 bp-1功p-率1 谱估计k = p伯格(B

17、urg)递推算法伯格(Burg)递推算法步骤(1)确定初始条件1N -1s=2y 2 k k =ek = yk e fb000Nk =0(2)从p=1开始迭代计算:计算AR模型参数N -1- 2 k e bk - 1e fp-1p-1k = pK=pN -1k = pk 2 +k - 12 e fbp-1ep-1s 2= (1- K 2 )s 2递推p阶均方误差p-1pp功率谱估计伯格(Burg)递推算法伯格(Burg)递推算法步骤(3)递推高一阶前、后向预测误差ek =k+ Kk - 1fefebp-1pp-1pek =k - 1+ Kbebefkp-1pp-1p(4)若阶数小于p，则阶数加

18、1，回到步骤(2)进行下一次迭代，直到达到预定阶数p。s2(W ) =(5)估计功率谱PAR2p1 + ap (n)e- jWnn=1功率谱估计Burg算法估计频谱的MATLAB函数x：进行功率谱估计的输入有限长序列；p:模型的阶数NFFT：DFT的点数；Fs ：绘制功率谱曲线的抽样频率，默认值为1；P率谱估计值；f：Pxx值所对应的频率点功率谱估计Pxx,f = pburg (x,p,NFFT,Fs)利用Burg法进行谱估计程序N=512;NFFT=1024;Fs=2;p=40;n=0:N-1;randn(state,0);x=cos(0.3*pi*n)+cos(0.32*pi*n)+randn(size(n);P,f=pyulear(x,p,NFFT,2);Pw,f2=pburg(x,p,NFFT,2);subplot(211);plot(f,10*log(P);grid;titl