高阶谱讲义第3章

1、第3章高斯有色噪声中的谐波恢复3.1模型假设我们关心这样一类有噪观测值y(n) = x(n) + w(n) = a (n)S (w ) + w(n)(3.1)i=1其中，S ()为信号波形(通常已知)，w为未知常数，a (n)为非高斯随机过程 或随机变量。附加噪声w(n)假设为能谱密度未知的高斯有色噪声，且w(n)与x(n) 相互独立。在噪声中谐波恢复RHN (Retrieve Harmonic in Noise)问题中，如果信号x(n) 包含p个复数谐波分量，且各个谐波频率各不相同。这时S. (wi) = a i exp( jwn)， a = exp( jw )，其中a和w为未知常数，平为独

2、立地服从同一分布的(i.i.d.)随机 iii ii变量，且平.在-兀,兀)上服从均匀分布，则x(n) = 2a expj(w n + 甲)(3.2)i=1称作复数谐波信号。如果x(n)包含p个实数谐波分量，则x(n) = Xa cos(w n + 甲)(3.3)i=1称作实数谐波信号。我们感兴趣的问题是怎样由有噪观测值y(n)估计谐波数目p，谐波频率w i及谐波幅度a J。3.2谐波信号的高阶累积量特性1.谐波信号的高阶累积量特性由第1章可知，零均值实数高斯随机过程y/勺二阶、三阶和四阶累积量分 别为cum( y , y ) = E y y i ji jcum( y , y , y ) =

3、E y y y i j ki j kcum(y , y , y , y ) = Ey y y y Ey y Ey y i j k li j k li j k l-Eyy Ey y - Eyy Ey y (3.4)i k j li l j k由于式(3.1)中过程y(n)可能为复数，因此需要定义复数过程的高阶累积量。复数过程的k阶累积量可以有2k种不同的定义方式(k项中的任何一项都可以取共轭或不取共轭)。对于不同的随机过程，需要根据其特征选取不同的定义方式。 这可以从下面的引理加以说明。引理3.1设随机变量中在-兀,兀)上服从均匀分布，且s = ejP，则s的各阶 累积量分别为E s = 0cu

4、m(s,s) = 0，cum(s*,s) = E|s|2 = 1cum( s, s, s) = cum( s *, s, s) = 0cum( s, s, s, s) = cum( s *, s, s, s) = 0cum(s*,s*,s,s) = E|s|4 -E(s2)2 一 2E|s|2E|s|2 = -1该引理的证明由累积量的定义式(3.4)容易得到。由此可知，单个谐波s = ejcp的三阶累积量恒等于零，对于一般的谐波信号， 我们有如下定理定理3.1谐波信号式(3.2)和式(3.3)的三阶累积量恒等于零，即cum(x(n),x(n + m ),x(n + m ) = cum(x*(n

5、),x(n + m ),x(n + m ) = 0(3.5)对于观测过程y(n)，我们有下列推论：推论3.1设w(n)为高斯(白色或有色)噪声，x(n)为谐波信号式(3.2)或式(3.3)，且w(n)与x(n)相互独立，则y(n) = x(n) + w(n)的三阶累积量恒等于零，即(3.6)cum(y(n), y(n + m ), y(n + m ) = cum(y* (n), y(n + m ), y(n + m ) = 0该推论说明高斯噪声中的谐波恢复问题不能利用三阶累积量信息来处理。下面讨论复数谐波信号的四阶累积量。由于式(3.2)中信号x(n)是平稳的， 所以(x(n),x(n + m

6、)x(n + m2)x(n + m3)的四阶累积量必须与变量n无关。由于RHN问题的复数谐波信号x(n)可以表示成复数指数和的形式，故四项中的两项 必须取共轭。这样，我们给出四阶累积量的定义如下：(m , m , m ) = cum(y* (n), y*(n + m ), y(n + m ), y(n + m )123123(3.7) TOC o 1-5 h z 其中cum(x ,x ,x ,x )的表达式见式(3.4)。 0123与传统的定义方式一样，我们定义零均值复数平稳随机过程y(n)的自相关函 数为r (m) = cum(y*(n), y(n + m)(3.8)现在讨论谐波信号的四阶累

7、积量特性。定理3.2对于模型(3.1)，设气(n)为零均值独立非高斯随机过程，其四阶累积量为c4 .(mm2,m3),w(n)为高斯噪声，则y(n)的四阶累积量为c(m,m ,m ) = s*()S*()S ()S()c(m,m ,m )(3.9)4, y123n .n+min+min+mi4,a.123.=11231若气(n)为零均值独立的非高斯随机变量气，且气的二阶和四阶累积量分别为Y(i)和Y(i)，则2,a4,ac (m ,m ,m ) = 2S*( )S*( )S( )S ( )Y (i)(3.10)4, y 123n i n+mi n+min+mi 4,ai=1123且r (m)

8、= 2LS*(3 )S()Y(i) + r (m)(3.11)yn i n+mi 2,awi =1定理3.3对于式(3.2)中的复数谐波信号x(n)，其四阶累积量为|a |4 exp j(-m + m + m )其自相关函数为4 e x j!D (m) k定理3.4对于式(3.3)中的实数谐波信号x(n)，其四阶累积量为c (m , m , m ) = -8a 4 cos (m - m - m ) + cos (m - m - m )4, x1238 kk 123k 231= 1+ cos w (m - m - m )其自相关函数为r (m)=8 a 4 cos (m)= 1由定理3.3和定理

9、3.4可以得到下列推论：推论3.2式(3.2)中复数谐波信号的四阶累积量的一维对角切片为c (m) = c (m, m, m) = -8 |a,= 1推论3.3式(3.3)中实数谐波信号的四阶累积量的一维对角切片为c (m) = c (m,m,m) = -28a4 c o c( m)，k=1卜面给出广义复数谐波信号和广义实数谐波信号的定义。定义3.2如果(3.12)(n) = 8 a 2 exp j(% n + 甲)k=1其中，ak,气和%如式(3.2)定义3.3如果则称(n)为广义复数谐波信号。(3.13) (n) = 8 a 2 c o(gb n + 甲)k=1其中，ak, %和%k如式(

10、3.3)，则称(n)为广义实数谐波信号。容易证明，谐波信号的四阶累积量的一维对角切片与广义谐波信号的自相 关函数之间有如下关系。定理3.5复数谐波信号(式(3.2)的四阶累积量的一维对角切片(m)与广 义复数谐波信号(式(3.12)的自相关函数七(m)的关系为 TOC o 1-5 h z c (m) = -r (m) 4,尤定理3.6实数谐波信号(式(3.3)的四阶累积量的一维对角切片(m)与广义实数谐波信号(式(3.13)的自相关函数r (m)的关系为/、3/、c (m) = - r (m) 4/4定理3.5和定理3.6表明，谐波信号的四阶累积量的一维对角切片正好与具 有相同频率及相应幅度的

11、广义谐波信号的自相关函数相等(若不考虑因子-1或 -3/4)。2.线性预测方程我们知道，复数谐波信号满足以下具有零输入的AR (p)模型X a (m)x(n - m) = 0m=0其中a(0) = 1，且多项式A(z) = a(m)z-m的根为z = e (i = 1,2, ,p)。而实数m=0谐波信号则满足以下AR(2p)模型艺 a (m)x(n - m) = 0m=0其中多项式A(z) = a(m)z-m的根为z = e+汽(i = 1,2, , p)。于是有下列定理。m=0定理3.7对于有噪观测值y(n) = x(n) + w(n)，其中w(n)为白噪声，且rw(m) = b26(m)，

12、如果x(n)为复数谐波信号，则基于自相关的线性预测方程La(m)r (k - m) =b 2a(k) ywm=0成立；如果x(n)为实数谐波信号，则基于自相关的线性预测方程La(m)r (k - m) =b 2a(k) ywm=0成立。定理3.8对于有噪观测值y(n) = x(n) + w(n)，其中w(n)为i.i.d.非高斯白噪声，且c (m ,m ,m ) = y 6(m ,m ,m )。如果x(n)为复数谐波信号，则基于四4, w 1 234, w1 2 3阶累积量的线性预测方程 a(m)c (m , m , k 一 m) = y a(k)m=0成立；如果x(n)为实数谐波信号，则方程

13、La(m)c (m , m , k 一 m) = y a(k)m=0成立。定理3.9对于有噪观测值y(n) = x(n) + w(n)，其中w(n)为高斯噪声(白色或有色)，如果x(n)为复数谐波信号，则基于四阶累积量的线性预测方程a(m)c (m , m , k - m) =0m=0(3.14)成立；如果x(n)为实数谐波信号，则方程a(m)c (m , m , k - m) =0m=0(3.15)成立。定理3.73.9表明，基于自相关的线性预测方程仅适用于白噪声情形，而基 于四阶累积量的线性预测方程不仅适用于白噪声情形，而且还适用于高斯有色噪 声的情形。正是基于线性预测(3.14)和(3.

14、15)估计高斯有色噪声中的谐波信号模型 参数a(m)，进而估计谐波频率气。我们知道，在所有的自相关谱估计法中，SVD-TLS方法和ESPRIT方法分 别是线性预测法和特征结构法中性能最为优越的方法，而自相关SVD-TLS方法 已经推广到四阶累积量情形。因此，下面将ESPRIT方法推广到四阶累积量情形， 给出了四阶累积量ESPRIT方法(FOC-ESPRIT)。3.3高斯白噪声中的谐波恢复 设有噪观测值为y (n) = x(n) + w( n)其中，x(n)为谐波信号，w(n)高斯白噪声1 .自相关函数Pisarenko法自相关线性预测方程实数情况:a(m)r (k - m) =b 2a(k)

15、ywm=0复数情况:a(m)r (k - m) =b 2a(k) ywm=0其中：r (m) =02 d (m)， r ()是y (n)的自相关函数。证明:,/ y(n) = x(n) + w(n),x(n) = y (n) - w(n)(实数情况)又 a(m)x(n - m) =0m=0Epa(m) y (n - m) =J a(m)w(n - m)两边同乘以j*(n-k)再取期望，得a(m) E y *(n k) y(n m)a(m)E y *(n k) w(n m)m=0m=0其中Ey*(n k)y(n m) = r (n m)yE y *(n k) w(n m) = E x * (n

16、k) + w* (n m)w(n m)=Ew* (n k)w(n m) = r (n m) = b 26 (k m)wwa(m)b 2 6 (k m) = b 2 a(k) w预测方程变为：a(m)r (k m) =Zym=0m=0讨论：若上式中 k = 2p +1,2p + 2,4p +1其矩阵形式为:.暑a(m)r (k m) =0m=0p + 1) ry (2p + 2):ry (2 p) r (2p +1):U) ry (2)a(0) 1:匕(4 p +1)r (4p)ry (2 p +1)_a (2 p)y=0(文)若相关函数已知，则可计算出a(m),m = 0,2p,/ y(n)

17、= x(n) + w(n) = a. c o c n(+ 甲.)+ w(n)m=0相关函数为：ry (0) =b w+ Orq.=1吐2对上式按k = 1,2,.,p展开，得矩阵其中：FS = rc ocsc oosc oos12pc o2cc o2oc o2c1:2:pc o psoc o psoc o pc12pF =r =,(1),r (2)，r (p)Ty yyS = SS 2，.,Sp TS =a 2/2, i = 1,2,p恢复方法：根据观测值y ()计算自相关函数r (m)(2p +1 m 4p +1)根据大矩阵，求 a (1),a (2 p)。 A(z)=艺a(m)z-m 求.

18、(i = 1,2,.,p)。m=0(4)根据(FS = r)S = F-1r 求 a . p。2.阶统计量方法(1)三阶累积量,/ y(n) = x(n) + w(n):.c(m , m ) = c(m , m ) + c (m , m ) = 03, y123,x 123,w12(2)四阶累计量c(m , m , m ) = c(m , m , m )丰 04, y1234, x1233.4高斯有色噪声中的谐波恢复 设有噪观测值为r (m)=b 26 (m)r (m) = f (m)wa(m)x(n - m) = 0y (n) = x(n) + w( n)其中，x(n)是谐波信号，w(n)为

19、高斯有色W?:高有:四阶累积量线性预测方程：复数情况：Xa(m)c (m , m , k 一 m) = 0m=0实数情况：艺 a(m)c (m , m , k 一 m) = 0m=0证明：.x(n -m) = y(n -m) - w(n -m),（实情况）:a(m)y(n - m)=La(m)w(n - m) TOC o 1-5 h z m=0m=0定义四阶累积量：c4,v=c u m*(n - k ), y *(n - k ), y *(n - k ), a(m)y (n - m) 123Lm=0=La(m)cumy*(n - k ), y*(n - k ), y*(n - k ), y(n

20、 - m)123m=0=Za(m)c(k - k ,k - k ,k - m)4, y 12131m=0, a(m)(n - m)3m=0另外c = cum y*(n - k ), y*(n - k ), y*(n - k ) 4,v=a(m)c u it*(n - k ) + w* (n - k )-j L*(n - k ) + w*(n - k ) 1122m=0m)*(n - k ) + w*(n - k=a(m)cumL*(n - k ), x*(n - k ), x(n - k ),0-1123m=0+ Pa(m)cumx*(n - k ), x*(n - k ), x(n - k

21、), w(n - 123m=0=0所以 a(m)c(k - k , k - k , k - m) = 04, y12131m=0令：k - k = m , k - k = m , k = k1211321a(m)c (m , m , k - m) = 0m=0令： m =m =0a(m)c(k - m) = 04, ym=0证毕上式按k = 0,1,2, ,2p展开:a 4, y (0)c (1)气.1)c,y (0)4, y c4 (-p)c (1 - p)4, y:c 4 (-2 p). c 1-2p):4,y.a(0) :c (2 p)1- 4,yc(2p -1)4, yc(p)4, y

22、c (0)4, y_a(2 p)_=0频率估计方法：(1)由观测值y()计算四阶累积量c(m)- 2p m 2p4, y(2)根据式估计 a(0),a(1),a(2p)。(3)根据 A( z) = a (m) z -m 求 w (i = 1,2,., p)。 m=0幅度估计方法：,/ y(n) = x(n) + w(n):.c(m) = c(m) + c(m) = c (m)= 一一乙a4 cos(m) TOC o 1-5 h z 4, y4, x4,w4, x8 iii=1将上式按m = 0,1,2,.,p -1展开，得FS = r-111一F =cos w妃 cos2w-1:cos w2c

23、os2w2cos wc pcos 2wp:cos(p 一 1)w 1cos(p 一 1)w 2cos(p 一 1)w pS = a4, a4,.,a 4t TOC o 1-5 h z 88(p -1)c(0),.,一 c4, y3 4, y3:.S = F-1r步骤：由y(n)计算四阶累积量c4 y (m)0 m p ,构造mxl维向量y ()= D(),( +1), ,川 + 质-l)ry i (乃)=顷 +1),+ 2), ,)( + m)rw(n) = w(n), w(n +1)，w(n + m- l)ry()= *() +w()= Y a expj(co n +(p ) + w(n)*

24、 i大1 _e沁p们叫1=1用矩阵形式表达：y(n) = AS(n) + w(n)y.) = AC 知+ w( + l)其中-1 1 .e j!e j2 A =幻2气e,2%e j(w-i)ro e j(-i)2A 的列为o(co )J = 1,2,ia(a) = 1,erg/(mT)co0 = d i ajgi,e俱，,e叫 旋转因子矩阵S(n) L ej(叩+q),以 ej(2n+%),以 ej(%+?p)求y(n)自相关函数% Ey (n)y *(n)把大式代入上式，得R EAS(n)S * (n)A * + E&(n)&* (n) = AS A * + E&(n)&* (n)其中=ES

25、(n)S *(n) = EL ej(气n+甲)1L * e - j (气 n+中),.,l * e -j ( pn+中 p)diag 以 |2,L |2,L y(n)与y 1(n)互相关阵R Ey(n)y * (n) AS 0* A * + E0(n 知 * (n +1) TOC o 1-5 h z 冯ii设c与c分别为R 一 Em(n)&* (n) AS A *yic R一 E0(n知* (n +1) AS 0A*yy1 yy1i若 r (m) = Ex(n)x* (n + m)，则c AS A * c AS 0*Ar (0)r *(1)3)r (0)xx:r (m -1) r (m - 2) xx-r *(1)r *(2)xxr (0)r* =x.x:r (m - 2) r (m - 3) xxr*(m 一 1)r *(m 一 2)xr (0)xr *(m)r *(m -1)x:r*x求 人 i仅适用于白噪声。c x 人c xESPRIT方法