高阶谱第5章 高斯有色噪声中的谐波恢复.docx

高阶谱 第5章 高斯有色噪声中的谐波恢复第5章 高斯有色噪声中的谐波恢复 5.1 模型假设 我们关心这样一类有噪观测值 py(n)=x(n)+w(n)=åai=1i(n)Sn(wi)+w(n) (5.1) 其中，Sn(·)为信号波形，wi为未知常数，ai(n)为非高斯随机过程或随机变量。附加噪声w(n)假设为能谱密度未知的高斯有色噪声，且w(n)与x(n)相互独立。 在噪声中的谐波恢复RHN问题中，如果信号x(n)包含p个复数谐波分量，而且各个谐波频率各不相同。这时， ai=exp(jji)， 其中ai和wi为未知常数，ji为独立地Sn(wi)=aiexp(jwin)服从同一分布的(i.i.d.)随机变量，且ji在-p,p)上服从均匀分布，则 p x(n)=åai=1iexpjwi(+nji ) (5.2) 称为复数谐波信号。如果x(n)包含p个实数谐波分量，则 px(n)=åai=1icos(win+ji) (5.3) 称为实数谐波信号。 我们感兴趣的问题是如何由有噪观测值y(n)估计谐波数目p，谐波频率wi及谐波幅度ai。 5.2 谐波信号的高阶累积量特性 5.2.1 谐波信号的高阶累积量特性 由第1章可知，零均值实数高斯随机过程yi的二阶、三阶和四阶累积量分别(yi,yj)=Eyiyj, cum(yi,yj,yk)=Eyiyjyk 为 cumcum(yi,yj,yk,yl)=Eyiyjykyl-EyiyjEykyl -EyiykEyjyl-EyiylEyjyk (5.4) 1 由于式中的过程y(n)可能为复数，因此需要定义复数过程的高阶累积量。复数过程的k阶累积量可以有2k种不同的定义方式。对于不同的随机过程，需要根据其特征选取不同的定义方式。这可以从下面的引理加以说明。 引理5.1 设随机变量j在-p,p)上服从均匀分布，且s=ejj，则s的各阶累积量分别为 Es=0 cum(s,s)=0，cum(s,s)=Es=1 *2cum(s,s,s)=cum(s,s,s)=0 *cum(s,s,s,s)=cum(s,s,s,s)=0 *cum(s,s,s,s)=Es-E(s)-2EsEs=-1 *42222该引理的证明由累积量的定义式容易得到。 由此可知，单个谐波s=ejj的三阶累积量恒等于零，对于一般的谐波信号，我们有如下定理 定理5.1 谐波信号式和式的三阶累积量恒等于零，即 cum(x(n),x(n+m1),x(n+m2)=cum(x(n),x(n+m1),x(n+m2)=0 * (5.5) 对于观测过程y(n)，我们有下列推论： 推论5.1 设w(n)为高斯噪声，x(n)为谐波信号式(5.2)或式(5.3)，且w(n)与x(n)相互独立，则y(n)=x(n)+w(n)的三阶累积量恒等于零，即 cum(y(n),y(n+m1),y(n+m2)=cum(y*(n),y(n+m1),y(n+m2)=0 (5.6) 该推论说明高斯噪声中的谐波恢复问题不能利用三阶累积量信息来处理。 下面讨论复数谐波信号的四阶累积量。由于式(5.2)中信号x(n)是平稳的，所以x(n),x(n+m1),x(n+m2),x(n+m3)的四阶累积量必然与变量n无关。由于RHN问题的复数谐波信号x(n)可以表示成复数指数和的形式，故四项中的两项必须取共轭。这样，我们给出四阶累积量的定义如下： c4,y(m1,m2,m3)=cum(y(n),y(n+m1),y(n+m2),y(n+m3) (5.7) * 2 其中cum(x1,x2,x3,x4)的表达式见式(5.4)。 与传统的定义方式一样，我们定义零均值复数平稳随机过程y(n)的自相关函数为 * ry(m)=cum (5.8) (y(n),+y(n m) ) 现在讨论谐波信号的四阶累积量特性。 定理5.2 对于模型(5.1)，设ai(n)为零均值独立非高斯随机过程，其四阶累积量为c4,a(m1,m2,m3)，w(n)为高斯噪声，则y(n)的四阶累积量为 ipc4,y(m1,m2,m3)=åSi=1*n(wi)Sn+m1(wi)Sn+m2(wi)Sn+m3(wi)c4,ai(m1,m2,m3) (5.9) *若ai(n)为零均值独立的非高斯随机变量ai，且ai的二阶和四阶累积量分别为r2,a(i)和g4,a(i)，则 pc4,y(m1,m2,m3)=åSi=1*n(wi)Sn+m1(wi)Sn+m2(wi)Sn+m3(wi)g4,a(i) (5.10) *且 pry(m)=åSi=1*n(wi)Sn+m(wi)r2,a(i)+rw(m) (5.11) 定理5.3 对于式(5.2)中的复数谐波信号x(n)，其四阶累积量为 pc4,x(m1,m2,m3)=-åakk=14expjwk(-m1+m2+m3) (5.12) 其自相关函数为 p rx(m)=åk=1ak2expjw(km) (5.13) 定理3.4 对于式(5.3)中的实数谐波信号x(n)，其四阶累积量为 c4,x(m1,m2,m3)=-1p4kaå8k=1coswk(m1-m2-m3)+coswk(m2-m3-m1) +coswk(m3-m1-m2) (5.14) 其自相关函数为 rx(m)=1p2kaå2k=1coswk(m) 3 (5.15) 由定理5.3和定理5.4可以得到下列推论： 推论5.2 式(5.2)中复数谐波信号的四阶累积量的一维对角切片为 p c4,x(m)=c4,x(m,m,m)=-åakexp(jwkm) (5.16) k=14推论5.3 式(5.3)中实数谐波信号的四阶累积量的一维对角切片为 c4,x(m)=c4,x(m,m,m)=-åa8k=13p4kcosw(km) (5.17) 下面给出广义复数谐波信号和广义实数谐波信号的定义。 定义5.2 如果 x(n)=påk=1ak2expj(wkn+jk) (5.18) 其中，ak，wk和jk如式(5.2)，则称x(n)为广义复数谐波信号。 定义5.3 如果 x(n)=påak=12kco(wskn+jk) (5.19) 其中，ak，wk和jk如式(5.3)，则称x(n)为广义实数谐波信号。 容易证明，谐波信号的四阶累积量的一维对角切片与广义谐波信号的自相关函数之间有如下关系： 定理5.5 复数谐波信号(式(5.2)的四阶累积量的一维对角切片c4,x(m)与广(m)的关系为 义复数谐波信号(式(5.18)的自相关函数rx(m) (5.20) c4,x(m)=-rx定理5.6 实数谐波信号(式(5.3)的四阶累积量的一维对角切片c4,x(m)与广(m)的关系为 义实数谐波信号(式(5.19)的自相关函数rx c4,x(m)=-34r(m) x (5.21) 定理5.5和定理5.6表明，谐波信号的四阶累积量的一维对角切片正好与具有相同频率及相应幅度的广义谐波信号的自相关函数相等(若不考虑因子-1或-3/4)。 5.2.2 线性预测方程 (in+ji)+coswi(n-2)+ji=2cosw(i)coswi(n-1)+ji 例： cosw(in+ji) 设： x(n)=cosw4 则： x(n)=2cosw(i)x(n-1)-x(n-2) AR(2)过程 上式两边取z变换，得 1-2z-1cosw-2i+zX(z)=0 1-2z-1coswi+z-2=0 特征多项式 zwijwi1=ej， z2=e-， 特征根 特点：z1=z2=1,w-1i=tanIm(zi)/Re(zi) p 若： x(n)=åaicosw(in+ji) i=1p 则特征多项式为：Õ1-2z-1cosw2i+z-=0 i=12p 也可写成： åa(i)z-i=0i=0 其中，a(0)=1,a(i),i=1,2,L,2p，受限于根的模等于1 上式对应的差分方程为 2påa(m)x(n-m)=0，实谐波信号线性预测方程 m=02p 其中 A(z)=åa(m)z-m， 根z=e±jwim=0·复数谐波信号的线性预测方程满足以下具有零输入的AR(p)模型 påa(m)x(n-m)=0 m=0p其中a(0)=1，多项式A(z)=åa(m)z-m的根为z=ejwi(i=1,2,L,p)。 m=0·实数谐波信号则满足以下AR(2p)模型 2påa(m)x(n-m)=0 m=02p其中a(0)=1，多项式A(z)=åa(m)z-m的根为z=e±jwi(i=1,2,L,p)。 m=05 (5.22) (5.23) 若已知a(m),于是有下列定理。 m=0,1,2,L，可计算出wi。 定理5.7 对于有噪观测值y(n)=x(n)+w(n)，其中w(n)为白噪声，且rw(m)=swd(m)，如果x(n)2为复数谐波信号，则基于自相关的线性预测方程 (k-m)=swa(k) (5.24) 2påa(m)rm=0y成立；如果x(n)为实数谐波信号，则基于自相关的线性预测方程 2p 成立。 åa(m)rm=0y(k-m)=swa(k) (5.25) 2定理 5.8 对于有噪观测值y(n)=x(n)+w(n)，其中w(n)为i.i.d.非高斯白噪声，且c4,w(m1,m2,m3)=g4,wd(m1,m2,m3)。如果x(n)为复数谐波信号，则基于四阶累积量的线性预测方程 påa(m)cm=04,y(m1,m2,k-m)=g4,wa(k) (5.26) 成立；如果x(n)为实数谐波信号，则方程 2p 成立。 åa(m)cm=04,y(m1,m2,k-m)=g4,wa(k) (5.27) 定理5.9 对于有噪观测值y(n)=x(n)+w(n)，其中w(n)为高斯噪声(白色或有色)，如果x(n)为复数谐波信号，则基于四阶累积量的线性预测方程 påa(m)cm=04,y(m1,m2,k-m)=0 (5.28) 成立；如果x(n)为实数谐波信号，则方程 2påa(m)cm=04,y(m1,m2,k-m)=0 (5.29) 成立。 定理5.75.9表明，基于自相关的线性预测方程仅适用于白噪声情形，而基6 于四阶累积量的线性预测方程不仅适用于白噪声情形，而且还适用于高斯有色噪声的情形。正是基于线性预测(5.28)和(5.29)，来估计高斯有色噪声中的谐波信号模型参数a(m)，进而估计谐波频率wi。 我们知道，在所有的自相关谱估计法中，SVD-TLS方法和ESPRIT方法分别是线性预测法和特征结构法中性能最为优越的方法，而自相关SVD-TLS方法已经推广到四阶累积量情形。因此，下面将ESPRIT方法推广到四阶累积量情形，给出了四阶累积量ESPRIT方法(FOC-ESPRIT)。 5.3 高斯白噪声中的谐波恢复 设有噪观测值为 y(n)=x(n)+w(n) 其中，x(n)为谐波信号，w(n)高斯白噪声 5.3.1 自相关函数法 ·自相关线性预测方程 2p实数情况： åa(m)r2y(k-m)=swa(k) m=0p复数情况： åa(m)r2y(k-m)=swa(k) m=0其中： r2w(m)=swd(m)，ry(·)是y(n)的自相关函数。 证明： Qy(n)=x(n)+w(n)， x(n)=y(n)-w(n) 2p (实数情况) 又Qåa(m)x(n-m)=0 m=02p2p åa(m)y(n-m)=åa(m)w(n-m) m=0m=0两边同乘以y*(n-k)再取期望，得 2p2påa(m)Ey*(n-k)y(n-m)=åa(m)Ey*(n-k)w(n-m) m=0m=0其中 Ey*(n-k)y(n-m)=ry(k-m) Ey*(n-k)w(n-m)=Ex*(n-k)+w*(n-k)w(n-m) =Ew*(n-k)w(n-m)=r2w(k-m)=swd(k-m) 7 (5.30) (5.31) (5.32) (5.33) 将上述结果代入式(5.33)，得预测方程变为： 2p2pyåa(m)rm=0(k-m)=åa(m)swd(k-m)=swa(k) (5.34) 22m=0讨论：若上式中k=2p+1,2p+2,L,4p+1 2p则 åa(m)ry(k-m)=0 (5.35) m=0其矩阵形式为： éry(2p+1)rrêy(2p)Ly(1)ùéa(0)êry(2p+2)ry(2p+1)Lr)úùúêúy(2êMúêMMú= êêëry(4p+1)ry(4p)Lr(2p+1)úêMúêúyûúëa(2p)û若相关函数已知，则可计算出a(m),m=0,L,2p p Qy(n)=x(n)+w(n)=åaicowsin(+ji)+w(n) i=1相关函数为： ìp(0)=sa2ïr2iyw+ ïåíi=12ïpa2 ïry(k)=åicowsikîi=12对上式按k=1,2,L,p展开，得矩阵 FS=r écowscowsê12Lcowspù其中： F=êco2sw1co2sw2Lco2swúpúêMMú êêëcopsw1copsw2Lcopswúpûú r=ry(1),ry(2),L,ry(p)T S=S1,S2,L,SpT Si=a2i/2， i=1,2,L,p ·谐波恢复恢复方法： (1) 根据观测值y(n)计算自相关函数 ry(m) (2p+1£m£4p+1) 8 (5.36) (5.37) (2) 根据式(5.36)的矩阵，求a(1),L,a(2p)。 2p(3) 由A(z)=åa(m)z-m，求wi(i=1,2,L,p)。 m=0(4) 根据(FS=r)S=F-1r求a2i。 5.3.2 高阶统计量方法 1. 三阶累积量 Qy(n)=x(n)+w(n) c3,y(m1,m2)=c3,x(m1,m2)+c3,w(m1,m2)=0 2. 四阶累计量 c4,y(m1,m2,m3)=c4,x(m1,m2,m3)¹0 5.4 高斯有色噪声中的谐波恢复 设有噪观测值为 y(n)=x(n)+w(n) 其中，x(n)是谐波信号，w(n)为高斯噪声ì白色：ír2w(m)swd(m)î有色：rw(m)=f(m)·四阶累积量线性预测方程： p 复数情况： åa(m)c4,y(m1,m2,k-m)=0 m=02p 实数情况： åa(m)c4,y(m1,m2,k-m)=0 m=02p证明：Qx(n-m)=y(n-m)-w(n-m),åa(m)x(n-m)=0 m=02p2p (实情况)： åa(m)y(n-m)=åa(m)w(n-m) m=0m=0定义四阶累积量： c=cumé2pêy*(n-k1),y*(n-k*ù4,v2),y(n-k3),åa(m)y(n-m)ú ëm=0û 9 (5.38) (5.39) =åa(m)cumym=02p2p*(n-k1),y(n-k2),y(n-k3),y(n-m) * =åa(m)cm=04,y(k1-k2,k1-k3,k1-m) 另外 c4,v2pé*ù*=cumêy(n-k1),y(n-k2),y(n-k3),åa(m)w(n-m)ú m=0ëû2p =åa(m)cm=0umx(n-k1)+w(n-k1),x(n-k2)+w(n-k2), * x*(n-k3)+w*(n-k3),w(n-m) =åa(m)cumxm=02p2p*(n-k1),x(n-k2),x(n-k3),0 * +=0 åa(m)cumwm=0*(n-k1),w(n-k2),w(n-k3),w(n-m) 2p所以 åa(m)cm=04,y(k1-k2,k1-k3,k1-m)=0 令： k1-k2=m1,k1-k3=m2,k1=k 2påa(m)cm=04,y(m1,m2,k-m)=0 #证毕 令： m1=m2=0 2p åa(m)cm=04,y(k-m)=0 (5.40) 上式按k=0,1,2,L,2p展开： éc4,y(0)êc(1)ê4,yêMêêc4,y(2p)ëc4,y(-1)c4,y(0)c4,y(2p-1)LLMLc4,y(p)c4,y(-p)c4,y(1-p)LLMLc4,y(-2p)ùéa(0)ùúc4,y(1-2p)êMúúêú=0 (5.41) úêMúúêúc4,y(0)ûa(2p)úëû·频率估计方法： (1) 由观测值y(n)计算四阶累积量 10 4,y(m) -2p£m£2p c(2) 根据式(5.41)估计 a(0),a(1),L,a(2p)。 2p(3) 根据 A(z)=·幅度估计方法： (m)zåam=0-m，求wi(i=1,2,L,p)。 Qy(n)=x(n)+w(n) aå8i=1c4,y(m)=c4,x(m)+c4,w(m)=c4,x(m)=-3p4icos(mwi) (5.42) 将式(5.42)按m=0,1,2,L,p-1展开，得 FS=r (5.43) 1éêcosw1ê F=êcos2w1êMêêcos(p-1)w1ë1cosw2cos2w2cos(p-1)w2LLLMLùúcoswpúú cos2wpúúcos(p-1)wpúû14T S=a14,a24,L,ap 步骤： 8é8ùr=ê-c4,y(0),L,-c4,y(p-1)ú 3ë3ûS=Fr-1T (5.44) (1) 由y(n)计算四阶累积量 c4,y(m) 0£m£p-1 (2) 根据S=F-1r估计ai4(i=1,2,L,p) 5.5 四阶累积量子空间旋转不变法(FOC-ESPRIT) Roy等最先把ESPRIT方法应用于白噪声中的谐波恢复问题，该方法通过求矩阵对的广义特征值和广义特征向量来确定谐波参数，因而矩阵对的构造和矩阵对的广义特征分解是ESPRIT方法的两大步骤，同样地FOC-ESPRIT方法也包括这两大步骤。 11 数学基础：矩阵A,B,n´n方阵 Ax=lBx l为广义特征值，x为广义特征向量。 广义特征方程： F(l)=A-lB=0 由F(l)可解出l1Lln广义特征值。 (1) 构造矩阵对(A,B)。 (2) 求li，li=ejw®求wi。 i(3) 求幅度ai。 1. 四阶累积量矩阵对的构造 设mfp，构造m´1维向量 y(n)=y(n),y(n+1),L,y(n+m-1)T y1(n)=y(n+1),y(n+2),L,y(n+m)T w(n)=w(n),w(n+1),L,w(n+m-1)T p y(n)=x(n)+w(n)=用矩阵形式表达： åai=1iexpj(win+ji)+w(n) (5.45) y(n)=AS(n)+w(n) (5.46) y1(n)=AFS(n)+w(n+1) (5.47) 1éêjwe1ê其中 A=êej2w1êMêj(m-1)w1êëeA1eejw2j2w2LLLMee1jwpj2wpej(m-1)w2Lej(m-1)wpùúúú m´p维范德蒙矩阵 úúúûii的列为a(wi),i=1,2,L,p， a(wi)=1,ejw,L,ej(m-1)wT jw1F=diage,e1jw2,L,ejwp ) 旋转因子矩阵 ,L,apej(wpn+jp)Tj(wn+j)j(w,a2e S(n)=a1e12n+j2 12 ·求y(n)自相关函数 Ry=Ey(n)y(n) *把式(5.46)代入上式，得 Ry=EAS(n)S*(n)A*+Ew(n)w*(n)=AS1A*+Ew(n)w*(n) 其中 ìéa1ej(w1n+j1)ùïêú*-j(w1n+j1)*-j(wpn+jp)S1=ES(n)S(n)=EíêMae,L,ae1púïêaej(wpn+jp)úûîëpüïýïþ=diaga1,a2,La222p * ·y(n)与y1(n)互相关阵 Ryy1=Ey(n)y1(n)=AS1A+Ew(n)w(n+1) *设cy与cyy分别为 1 cy=Ry-Ew(n)w*(n)=AS1A* cyy=Ryy-Ew(n)w*(n+1)=AS1A* 11自相关ESPRIT方法就是通过求解矩阵对cy,cyy的广义特征值来确定旋转1因子矩阵F，进而得到谐波频率wk的。然而cy和cyy的获取必须以Ew(n)w*(n)1和Ew(n)w*(n+1)预知为前提。如果w(n)为白噪声，则可以通过对Ry进行特征分解来确定噪声方差sw，从而确定Ew(n)w(n)和Ew(n)w(n+1)；如果w(n)为高斯有*2色噪声，则无法确定Ew(n)w*(n)和Ew(n)w*(n+1)，这就是说，自相关ESPRIT方法仅适用于白噪声情形，而不适用于有色噪声情形。 实际上，由复数过程自相关的定义，即rx(m)=Ex(n)x*(n+m)可知 érx(0)êrx(1)=êêMêêrx(m-1)ërx(1)rx(0)rx(m-2)*LLML* cy=AS1Arx(m-1)ùú*rx(m-2)ú (5.48) úúrx(0)ûú* 13 cyy1érx*(1)êrx(0)*=AS1A=êêMêêërx(m-2)求rx(2)rx(1)rx(m-3)*LLML*rx(m)ùú*rx(m-1)ú (5.49) úú*rx(1)úûcyx=lcyy1x ¾¾®l¾¾®wi (Ax=lBx) ·构造四阶累积量矩阵的一维对角切片矩阵对 设 éc4,y(0)êc4,y(1)ê=êMêêëc4,y(m-1)*éc4(1),yêc4,y(0)=êêMêêëc4,y(m-2)c4,y(1)c4,y(0)c4,y(m-2)c4,y(2)c4,y(1)c4,y(m-3)*LLMLLLMLc4,yc4,y(m-1)ùú*c4,y(m-2)úúú*c4,y(0)úû* (5.50) c4,yy1*c4,y(m)ùú*c4,y(m-1)úúú*c4,y(1)úû (5.51) 考虑到w(n)为高斯有色噪声，由(5.16) pc4,y(m)=c4,x(m)=-åaii=14exp(jwim) 4 (5.52) 设 S=diag-a1,-a2,L,-ap则有 2. 恢复频率wi c4,y=ASAc4,yy1*44 üï (5.53) ý四阶累积量矩阵*=ASAïþ设l和e分别是矩阵对c4,y和c4,yy的广义特征值和广义特征向量，则 1 c4,ye=lc4,yye 1广义特征方程： c4,y-lc4,yy=0 1¾®wi 解出l¾说明：(l中含有wi) QASA*-lASF*A*=0 14 AS(I-lF*)A*=0 即 I-lF*=0 é1êêêëjwùée1úê-úêê1úûëOOùée-jw1úêúêjwpúêeûëOùúú=0 -jwepúû l=ejw， i=1,2,L,p 13. 恢复幅度ai Qc4,yei=lic4,yyei (5.54) 14c4,yei-lic4,yy1ei=0* 正交 AS(I-liF)Aei=0 因为矩阵AS(I-liF*)A*行空间与向量a(wj),i¹j描述的子空间相同。 ì0,i¹jï结论： ea(wj)=í *ïeia(wj)î*iì0,i¹jï a*(wj)ei=í *ïa(wj)eiî用式(5.54)左乘ei*，得 eic4,yei=lieiASFAei *é-a1ê*=eia(w1),a(w2),L,a(wp)êê0ë=-ai440O-apùéle-jw1iúêúê4úêûëOùéa*(w1)ùúêúMúêúei-jwêa*(w)úliepúpûûëea(wi)*i2*则 ai15 4=-eic4,yeiea(wi)*i2 i=1,2,L,p (5.55) ·FOC-ESPRIT方法步骤： (1) 由y(n)®c4,y(m),1mfp(2) 构造c4,y,c4,yy矩阵对 (3) 求c4,y,c4,yy的广义特征值li=ejw得wi,i=1,2,L,p及广义特征向量i1ei,i=1,2,L,p*(4) a4ic4,yeii=-e， e*a(w2ii)i=1,2,L,p 16