课程材料数字信号处理ii chapter

1、第五章现代谱估计2013/11/17MMVCLAB15.10 白噪声中正弦波频率的估计2013/11/17MMVCLAB2 问题的提出：Ø 估计淹没在噪声中的正弦波的频率与幅度，是信号 处理中最有实际应用价值的技术之一。Ø 估计淹没在噪声中的正弦波的频率，也是评价谱估 计方法性能的经常采用的谱估计问题。2013/11/17MMVCLAB3二 问题的定义：已知长度为N的离散随机时间序列模型为：Måi=1j (w n+j )x(n) =+ v(n) ,n = 0, 1 " N-1Aeiii其中：Ai、wi是待估计的未知常数。ji是0, 2p)内均匀分布的v(

2、n)是均值为0，方差为s2的随量。白噪声。试根据这N个随机采样样本估计Ai、wi，i=1, 2 . M。2013/11/17MMVCLAB4三 白噪声中正弦波序列的性质：统计性质（时间平均和集合平均）：1）一阶与MEx(n) = Eå Ai ej (w n+j )+ v(n) = 0iii=1R (n, m) = Ex(n + m)x*(n)xMMååk =1j(n+m)w +jw +j )= E+ v(n + m)+ v(n)- j (nAe A eiikkiki=1Må iw+ s 2d (m)=jmPeii=1= R (m) ,P =| A |2x

3、ii结论：¾ 白噪声中的正弦波序列为平稳的、遍历的随机过程。¾ 白噪声中的正弦波序列的频率与其功率谱密度函数中谱峰的位置相对应。2013/11/17MMVCLAB5+¥xåxS (w ) =R (m)e- jmwm=-¥M= å Pd (w -w ) + s 2ii i=1Mx(n) = å Aej(win+ji ) + v(n)ii=12）白噪声中正弦波序列与ARMA模型的关系：以只有一个复正弦波序列的情况为例：pq目标：x(n) = å a x(n-i) + åb u(n-j)ej(wAn+ s(+n

4、)-j)令：s(n) =ij=x(n)1)v(n)i=1j=0s(=n+1)x(n-v(n-1)= Ajew (n=s(ne)-jw= e- jw-+1)v+j(n-1)1)AR和MA参数相同的、特殊的- AR1)MA(1,1)模型+v(n v(n) +x(n)v(nx(nv(n-1)特殊性：1- ejw1- ejwz-11.H (z) =；z-12. 极点ejw 在3. v(n) = u(n)。¾可以证明：具有M个复正弦波的序列也可以用类似的圆上，系统不稳定，不适合用ARMA模型法求解；A2R01M3/A1(1M/1,7M) 模型来表示。MMVCLAB6x(n)-= e1jw)+v

5、(n)e jw-四 问题的解法：n 一般谱估计的方法；n 根据问题特性寻找新的方法。2013/11/17MMVCLAB7（一）一般谱估计的方法：1） 周期图方法n 简便。n 两个正弦波频率相互很接近时无法分辨。2） AR模型法：n 分辨率高。n 对噪声和正弦波的相位敏感，容易引起谱峰移动。n 采用修正协方差法估计AR模型参数对噪声和正弦波相位相对来说敏感性较低，能提供相对稳定的高分 辨率的估计结果。2013/11/17MMVCLAB83）修正协方差AR谱估计法：（1） 基本思路：n 利用AR模型法的性质，特别是AR模型法与线性预 测谱估计等效。n AR模型的参数与线性E(z) = A(z)&#

6、215;X(z)。n 对于平稳的随机时间序列，用时间平均代替集合平 均，即：AR模型参数的提取化为以下优化问题：滤波器的参数相同：Ee+ (n)2Minak , k=1, 2 " pMinak , k=1, 2 " ppEe+ (n)2 + e- (n)2pp2013/11/17MMVCLAB9（2）估计准则：Minak , k=1, 2 " pN-1åee (n) +e- (n)2+=2ppn=p(1)(2)(4)(3)n =0n =N-1与Burg法相同。2013/11/17MMVCLAB10（3）求解方法：与Burg法不同之处在于直接对AR参数求解

7、，即： ¶e¶api= 0 ,i = 1, 2 " p2013/11/17MMVCLAB11（4）求解结果：éC(1, 1)C(1, p)ù éap1 ù""%"éC(1, 0)ùC(1, 2)C(2, 2)#C(p, 2)C(2, p)ú êaúêC(2, 1)êC(2, 0)úú êp2 ú = - êêúú êúê

8、;êú#C(p, p)ú êaúêC(p, 1)êC(p, 0)úû ëêpp ûúës 2ëûp= C(0, 0) + å a C(0, i)pi其中：i=1ìN-1üN-11åx (n - j)x(n - i) +åx (n - p + i)x(n - p +=*C(i, j)j)ýþíN - pîn=pn=p2013/11/17MMVCLAB

9、12（5） 估计性能：n 当信噪比很高时，得到的频率估计是无偏估计且方 差接近于Cramer - Rao极限。n 当信噪比很低时，估计性能却很差，出现偏倚同时 有大的方差。（6） 信噪比低时的改进方法：n 适当增加AR模型的阶数。2013/11/17MMVCLAB13（二）根据问题特性寻找新的方法：I)最大似然法噪声是白噪声。II)特征分解频率估计法通过研究自相关矩阵中正弦波信号和白噪声相互关系 的性质提出的。2013/11/17MMVCLAB14I）最大似然法：nA(I)符号s表示：wn+i j=)e j,jje wi=x(1)= 1,1)1)1)j(ni-(n)AeA" 2iii

10、iiTx =x(0)(0) v(0)A"x(Ns =sTs(1)""s(Niv =iiiTv(Nv(1)jjjwiw "- wi1) TA =j2j(Ne1eeeeiiciis = AcieiiiMåAci ex=+vii = 12013/11/17MMVCLAB15Mx( )e= åjv(win+ji ) +( )nii=1（2）基本思想：Mx = åAciei + vi=1Mx - åAi=1ÞN(0, s2I)eciiv N(0,s2I)2013/11/17MMVCLAB16（3）考虑x(n)中只有

11、一个复正弦信号的情况：x = Ac1e1 + v = s1+ v目标函数：11-( x-s )*T ( x-s )11p(x - s ) =s2e1pNdet(s2I)L(A, w ) = (x - s )*T (x - s )MinAc1 , w1c11112013/11/17MMVCLAB17j( w n+j )jjjw nsi (n) = Ai eii= Ai ei ei,i = 1, 2 " M（a）求解方法：¶ L(Ac1, w1 ) = 0¶Ac1¶ L(Ac1, w1 ) = 0¶w12013/11/17MMVCLAB18（b）求

12、解结果：2N-11Nåx(n) e-jnw1Maxw1n=0= 1NN-1å-jnwA1x(n) e1n=0éùN-1ê Imå x(n) e-jnw1új1 = ar ctg ên=0úN-1êêëúúûån=0-jnwRex(n) e1（c）结论：在白噪声中含有单个复正弦信号，其频率的最大似然估计，可根据数据的周期图的最大值所在的频率位置求出来。2013/11/17MMVCLAB19(4)考虑x(n)中有M个复正弦信号的情况：Mx

13、= åAciei + vi=1目标函数：L(Ac1 , Ac2 "AcM ; w1, w2 "wM )MinAci , wi ,i=1, 2 " MMM= (x - åAi=1åi=1(x -*Te )Ae )ciicii2013/11/17MMVCLAB20MinL(A, w ) = (x - s )*T (x - s )Ac1 , w1c1111(a)求解方法：¶ L(Ac1 , Ac2 "AcM ; w1, w2 "wM ) = 0 ,i = 1, 2 "M¶Aci¶

14、L(Ac1 , Ac2 "AcM ; w1, w2 "wM ) = 0 ,i = 1, 2 "M¶wi2013/11/17MMVCLAB21(b)求解结果：x*TET (E*ET )-1E*x ,E = e eTMaxw1 , w2 "wM(c)结论：" e12M当白噪声中有多个复正弦信号时，最大似然估计是非线性方程的最大化求解。最大似然估计容易陷入局部极值点，难以得到好的估计结果。如果各正弦波的频率能用周期图进行分辨，那么 最大似然估计结果将对应于周期图诸最大值所在的频率。2013/11/17MMVCLAB22II) 特征分解频率估

15、计法：(1) 问题的提出：n 在低信噪比情况下，AR谱估计结果不理想，为了提 高谱估计结果的精度，要降低噪声的影响。n 通过研究信号和噪声的特性及其关系，可能找到估 计白噪声中正弦波频率的其他思路。Mx = åAcieii=1+ v2013/11/17MMVCLAB23(2)符号s表示：wn+i j=)e j,jje wi=x(1)= 1,1)1)1)j(ni-(n)AeA" 2iiiiiTx =x(0)(0) v(0)"x(Ns =sTs(1)""s(Niv =iiiTv(Nv(1)jjjwiw "- wi1) TA =j2j(NA

16、e1eeeeiiciis = AcieiiiMåAci ei = 1x=+vi2013/11/17MMVCLAB24(3)基本思路：自相关矩阵是功率谱密度估计中很重要的一个量，几乎所有AR谱估计的方法都要用到自相关矩阵。因此通过研究自相关矩阵中正弦波信号和白噪声的相 互关系来达到以上目的。*TRR = Ess*Tvv*TsR = Ev2013/11/17MMVCLAB25Mx=åAci e i+v i = 1= s+v(4)自相关矩阵的性质（理论上）：性质1：Rx= Rs +Rv，且Rs的秩为M， Rv的秩为N。x = åAciei + vi=1M证明：M

17、9;i=1*TR =Exx =+ s2I ,=| A*T|2P e ePxiiiiiMå*TR =P e eR = s2Isiiivi=1Rx = Rs+ Rv2013/11/17MMVCLAB26jwj2wj(N-1)wTei = 1eiei"ei Mx = åAciei + vi=1= s + v性质2： Rx和 Rs有相同的前M个主特征向量。证明：R =R+R=R+2Isxsvs-( delt d(etl2-I)s=RR-s0l=0sI)-(R-=0(l2I- )us=0sxxR-ss xilIs()u s= lsi + l ,2= usi= ui ,i =

18、 1, 2 " Mi = 1, 2 "Muxilxis 2 +,+M1, M=2 " =Ni2013/11/17MMVCLAB27N -1åi=0TqR性质3：主特征向量u1, u2 . uM的子空间和e1,e2 . eM证明：的子空间相同。MMM= åuu= ål= ål*eTu*T*TRRP esiiisiisissiiisi = 1i=1i=1e1, e2 . eMRs列空间的基。Rs列空间的基。u1, u2. uM也iN(N-1)1)1)2013/11/17MMVCLAB28x =x(1)"T x(s =s

19、(1)"sT iiiiv =v(1)"Tv(Nis = AcieMx(0)x=åAcsi e i+v(0 i = 1= s+v v(0)性质4：定义u1, u2 . uM为信号子空间， uM+1,uM+2. uN为噪声子空间，则ei, i=1, 2 . M与噪声子空间正交。证明：Mei = åckuk k=1 uk"k > Mei,2013/11/17MMVCLAB29R = ULU*T ,其中：U = u , u" u性质5：x12NMN= å(li=1å+ s2 )u u*T +s u u2*Tsiiii

20、ii=M+1证明：RU = ULÞR = ULU*Ti = 1, 2 " Ml= l + s 2 ,xisi= s 2 ,lxiuxiuxii = M +1, M + 2 "N= usi =uii = 1, 2 " M i = M +1, " N30,= ui,2013/11/17MMVCLAB(5)信号子空间频率估计：基本原理：利用性质5忽略噪声子空间的影响，而只保留信号子空间中特征向量的信息，这样就能有效地提高信噪比，然后用AR模型法中的Yule - Walker 方程组估计AR模型的参数，进而用估计出的功率谱密度求正弦波的频率，就能较好地

21、保证估计精度。2013/11/17MMVCLAB31MNR= å(l + s 2 )u u*T + å s 2u u*Txsiiiiii=1i=M+1MMR = ål u u *T = ål u u *T ssisisisiiii=1i=1(a) 实现方法：§ 用随机采样样本估计p×p阶自相关矩阵Rx；§ 求Rx的前M个主特征值li和主特征向量ui，并用下式代替Rx：MR =å ul*u Txiiii = 1§ 用Yule - Walker方程组估计AR(p)模型的参数： Rx A p = -rxp

22、29; Rx (i) = s= 1 2ai,a0ri =0 (2) "R A = a=T ,(p)Ta" a12pxx2013/11/17MMVCLAB32=L(b)结果： =u - -1 A-1M= å 1*T irRulipxii=1圆的M个零点所对应的角度即为正弦最靠近波频率的估计结果。(c) 分析：不管模型阶数p增加到多少，用M个特征值和特征向量来估计模型参数。这样就既能增加模型的阶，而又不易产生虚假谱峰（因为提高了信噪比）。2013/11/17MMVCLAB33MR =å ul*u T xiii i = 1Rx A p = -rxURLU

23、22;R U=U*T(6) 噪声子空间频率估计：(a) 基本原理：e jwie j2wiej(N-1)wiei = 1T ,"i = 1, 2 " M,"k > M信号子空间ei uk< e , u*> = 0 ,"k > Mik2013/11/17MMVCLAB34(b) Pisarenko谐波分解（PHD）方法：Ø 问题：考虑N=M+1的特殊情况，即采样样本长度比 正弦波的数目多1。Ø 原理：Rx= Rs +Rv，且Rs的秩为M， Rv的秩为M+1。将估计出的(M+1)×(M+1)阶自相关矩阵Rx

24、对角化后由前M个主特征向量u1, u2 . uM子空间，而最后一个特征向量uM+1空间，即有：信号噪声子< uM+1, ei >=0,*i=1, 2 . M2013/11/17MMVCLAB35Ø 实现：N-1åun=0(n +1) e-jnwi= 0,i = 1, 2 "MM+1Ø 结果：将uM+1看成一个时间序列，其z变换位于上零点的角度就是正弦波的频率。圆2013/11/17MMVCLAB36Ø 注意事项：9自相关矩阵元素的估计值必须是自相关函数的有偏估计才能保持理论上的自相关函数应该满足的 性质(如正定性)。9当正弦波的数目小于M时， uM+1的z变