1、第二章 高斯噪声背景下谐波恢复数学模型2.1 高斯过程高斯过程又称正态随机过程,它是一种普遍存在和重要的随机过程.所谓高斯过程 ,)(t即指它的任意 n 维概率密度函数由下式表示的过程,即 njk kjjn axBtxf 12/121/121exp)(;,( 式中 ;)(;)(22kkk atEta|B|-归一化协方差矩阵的行列式,即 1121212 nnbB-行列式 |B|中元素 的代数余因子;jkBjk-归一化协方差函数:jkb kjkjjk attEb)()(由上式可以看出,正态随机过程的 n 维分布仅由各随机变量的数学期望,方差和协方差决定.一维高斯正态分布的概率密度函数可写为: 2)
2、(exp21)(af式中, 及 是两个常量(均值和方差 ).a高斯噪声一般分为白噪声和有色噪声。功率谱密度在整个频域内都是均匀分布的噪声，被称之为白噪声，即 2)(0nP白噪声的自相关函数为 )()(0R显见，白噪声的自相关函数仅在 时才不为零；这说明，白噪声只有在零点才相0关，而它在任意两个时刻上的随机变量都是不相关的。有色噪声与白噪声不同，它的功率谱在整个频段上不是均匀分布的。22 谐波恢复的数学模型高斯噪声背景下的谐波恢复，主要是利用特征子空间分析的方法，对观测值进行处理从而估计出原始信号的频率等参量特征，即完成了在噪声背景下对信号的恢复。我们首先对特征子空间进行分析。从几何意义上说协方
3、差几何空间=信号子空间+ 噪声子空间我们所要做的就是从大的空间抽取低秩子空间，对信号进行分解处理。下面将要介绍的Pisarenko,MUSIC,Prony,ESPRIT 等方法其核心思想都是由此而来的。设观测信号模型：y(n)=x(n)+w(n)pinwjienx1)()(其中， 为信号幅度i为谐波信号频率iw为相位在 均匀分布的参量i,p 为谐波个数w(n)为零均值方差为 的高斯白噪声2构造 协方差矩阵m )0()()1( 20)1()()(yyy yyy RMRR 设 R=S+W ：S 表示信号协方差矩阵，W 表示噪声协方差矩阵。 HipiisS21其中： )1(expiiiMjjs222
4、00IW下面我们对协方差矩阵 R 进行分析：R=S+W HipiisS21首先我们来说明 R 为什么等于 S+W.y(n)的自相关函数)()exp( )()( exp)(exp )()(exp)(exp()(221 21* 2*11* 11 *mj miljjnEjj mnwEnjnjEjjmnyripiil lilpi ill lili pl lli ii l lli iiy 对应矩 2212 212122001)2(exp)1(exp )(exp1)( 1)exp(1 )(exp)(exp)(1 ii iii iipiipi iiiiipiHii Mjj Mjjj jj IjjjjIsWS
5、阵各项可知 R=S+W.现在对矩阵 R 进行分析：(1)设 分别是矩阵 S 的特征值和特征向量，则iv,i=1MiivHipiS1讨论：因为 Mp，信号阵的秩必为 p。所以 S 阵有 p 个非零特征值（ i=1p）,M-p 个零特征值（i=p+1M） 。(2)W= I 由于 所以2MiHiv1MiHivW12pi MpiHiHiiiiii vvR11221)(讨论：（1）R 阵与 S 阵具有相同的特征向量 ，且ijivji0（2） 是信号特征对pii1),(是噪声特征对Mpivi1),(分别记： ,1MpsvGV由信号向量张成的子空间叫做信号子空间，而有噪声向量张成的子空间叫做噪声子空间。（3
6、）矩阵 R 的特征值 piii 12 由以上结论我们可以看出，信号向量和噪声子空间中的所有向量（包括它们的线性组合）是正交的。即 MpkkHivs10)(由于此结论十分重要我们在此进行证明。 22111200,pppiHiipiHiiiBsEvSvs)exp(1iiijjs piHii HpHpHssssEBS12 21221100,0)()(1111pHpHppvEBvS01p是 正 定 的展开上式即有 0s1Hipv证毕；利用信号向量和噪声空间的正交性进行信号恢复的方法称之为噪声子空间方法。第三章 算法概述与分析定义观测信号的空间协方差矩阵为R=Ey(n)yH(n) (3.1)设噪声方差为
7、 ，由假设得2R=EAs(n)+w(n)As(n)+w(n)H=E As(n)+w(n) As(n) H + w H (n)=E As(n)+w(n) s H (n) A H+ w H (n)=E As(n) s H (n) A H + w(n) s H (n) A H +As(n) w H (n)+ w(n) w H (n)=EAs(n) s H (n) A H+Ew(n) s H (n) A H+EAs(n) w H (n) +Ew(n) w H (n)=ASA H+ IM (3.2)2其中，S=Es(n) s H (n)，I M 为单位阵。注意到(3.1)式表达的是观测信号的“统计平均”
8、 。我们认为通信中的随机过程是平稳随机过程，而平稳随机过程具有各态历经性，所以在这里，我们可以用观测信号的“时间平均”来代替“统计平均” ，即 )()()( )()()(21 22 1211 nynynyEnyRMMMH )()()( )()()(21 22 1211 nyEnyEnyE MMM )0()0()(2122121MMRR R 中各元素的估计为： NnHjiij y1)(故可得 R 的估计为(3.3)n1)(Ry3.1 Pisarenko 谐波分解法在 Pisarenko 谐波分解法中，考虑的是由 p 个实正弦组成的确定性过程 (3.1.1) )2sin()(1iipifAx我们假
9、定，初始相位 是在 均匀分布的独立随机变量，它在一次实现中为常量。i,先来推导过程x(n)满足的差分方程。为此，先来考虑单个正弦波的情况，即.)2sin()fx将三角恒等式 )1(2sin)co(2)(i)2sin( ffff中的正弦函数换乘 x(n)后，得到二阶差分方程式 0)2(nxx上式两边去 z 变换，得 0)()2cos(121zXzf这样，我们就得到特征多项式 )(21f它有一对共轭负根即 fjefjfz2)sin()co(共轭根的模为 1，即 ，由它们可决定正弦频率1|2z/)Re(/arctnImiiif通常我们取正的频率。显然，如果 p 个实正弦波没有重复频率的话，则着 p
10、个频率也应由特征多项式 0)(201*kpkpi ii zazz或 (3.1.2)20kpka的根决定。易知， ，且系数 是对称的，即1|izka(3.1.3)0(2pakpk与式(3.1.2)对应的差分方程为(3.1.4)()(21inxpi进一步的，我们来考虑白噪声中的正弦波过程：y(n)=x(n)+w(n) (3.1.5)其中w(n)是一零均值、方差为 的高斯白噪声，即2ww(n)WN(0, )。 将式(3.1.5)代入式(3.1.4) ，立即可得（3.1.6）)()()( 2121 inwainyapipi 这表明白噪声中的正弦波过程y(n)是一个特殊的 ARMA(2p,2p)过程，其
11、 AR 参数与 MA参数相同。 为了推导 AR 参数满足的方程，令(3.1.7)Tpnwnwayy)2(,)1(, 于是式(3.1.6)可以写成一下矩阵形式：(3.1.8)ayT用向量 y 左乘式(3.1.8)两边，并取数学期望得到(3.1.9)aywETT若 ，则显知，)()(lnyElRy IwExEywRpRpyTTT yyyy yyyT 2)( )0()12(120)1)()( 将上述关系代入式(3.1.9),得特征方程(3.1.10)aRwy2式(3.1.10)表明， 是自相关矩阵 的特征值，而 a 是对应特征值 的特征向量。式2y 2w(3.1.10)组成了由 Pisarenko
12、发展的谐波分解方法的基础。这样一来，谐波恢复问题就转化成了自相关矩阵 的特征分解。y归纳起来，利用 Pisarenko 分解法确定阶数 p 和正弦波频率 f 的步骤如下。(1)计算矩阵(m2p)0()12()( 120)()()(yyy yyy RpR 的特征分解。（2）求矩阵 R 最小特征值，即为 ，并用 表示 的多重度。然后将 R 将结为 m+1- ,00即取前 m+1- 行和列组成新矩阵 R1,重复这一步直至从某一阶降至下一阶，最小特征值不变。这样，即可定出 p，与最小特征值对应的特征向量各分量用 表示。pa210,（3）求多项式的根，有 0210pzaza（4）将上式的根记为 ，其中
13、，由此定出 （取正*2,p )ex(tfjiiif值）无噪声的正弦波过程(2.11) 式的自相关函数为 piixfPR1)2cos()(其中， 。因此，白噪声中的正弦波 y(n)=x(n)+w(n)的自相关函数2/iiAP,故有)()(wxyRpi iiypiwiyfPR112)0(cos)(0rDPrDjfpijiij Tyy11)2cos(,)写成一式为 )()2cos()()( 212 xwpiiwy RfR因此接收信号 y(n)的功率谱由上式的傅里叶变换可得 pi iiwy ffPfP12 )()()( 下面进行仿真实验：信号： )25.01()5.0()tjtjeenx噪声：高斯白噪
14、声，方差为 1，均值（E）为 00 024 0 05 15 2 0 024 0 0 0 05 15 2 0 0 0 0 0 0 0 04 0 0 05 15 0 0 0 0 0 0 05 0 0 05 15 E =0.4715 1.4349Pisarenko 谐波分解法在理论上首次揭示出了相关矩阵对应于最小特征值的特征向量的系数即为正弦波过程的 ARMA 模型的 AR 参数。但是作为一种算法，pisarenko 算法并不有效，实用。3.3 MUSIC 法从原理上讲，Pisarenko 分解时将自相关矩阵分解成信号子空间（对应与大特征值的特征向量）和噪声子空间（对应于噪声方差的特征向量） 。由于
15、特征多项式的根决定正弦频率，而多项式的系数 又是由对应于噪声方差012pzaz ia的特征向量决定的，所以 Pisarenko 分解实质是一种噪声子空间方法。Pisarenko 方法不是一种十分有效的实际算法，但 MUSIC 和下面将要介绍的 ESPRIT 方法却十分有效。MUSIC 是多重信号分类（Multiple Signal Classification Characterrization）的英文缩写。信号处理中的几个重要问题都可以归结为估计下列模型的参数：(3.3.1),1)()(NtetxAty式中， 是带噪声的数据向量， 是信号幅值的向量， 代表nmC1)(nCtx 1)(mCte
16、加性噪声，而矩阵 具有下面的特殊结构：nm)(3.3.2),1naA这里， 是实的参数， 是第 I 个信号与 y(t)之间的传输向量，且i1(miCTn,1对以上模型，我们做以下假定。A1:mn,且对应于不同 值的向量 是线性独的。)(aA2: 0)(,)(,0)(2teEIteEte TH这里 H 代表向量的共轭转置， T 代表转置。A3：矩阵 是非奇异的，且 Nm.当满足假设条件时，观测向量 y(t)的协)(txP方差矩阵由（3.3.3）IPAtyERHH2)()( 给定。为方便计，我们将 简记为 A。注意到 R 是一对称阵，其特征值分解具有下列形式： ，其中，HUR2对角线元素 叫做 R
17、 的特征值。,212mdiag2ii在 A1 条件下，矩阵 A 显然是非奇异的，即 rank(A)=n 从而，在 A3 条件下 的HAP秩也为 n。因此 R 的特征值必须满足下列关系： ),1(,2mnii 不妨令 分别对应的特征向量为 ，而 ,对应的特征向量为n,1 s,1 n，定义mg, ,1nmgGsSS 和 G 分别叫做信号子空间和噪声子空间。于是，特征矩阵 U 分为两个子矩阵 U=S|G。我们的兴趣是研究噪声子空间 G 和矩阵 A 之间的关系。一方面，由于 和 G 分别是2R 的特征值和对应的特征向量，故有特征方程(3.3.4)R2另一方面，用 G 右乘式(3.3.3)，又有(3.3
18、.5)GAPH2综合上面两式给出 0从而有 。)()(GAPAPGHH由于 P 是非奇异的，所以 (3.3.6)0上式又可写作(3.3.7),(0)(1nHaG由于 U 是酉阵，故 IGSSUH|于是式(3.3.7)也可写作(3.3.8),(0)(1nHHaSIaMUSIC 算法的基本思想是对真实相关阵 R 使用性质(3.3.7)或(3.3.8) 。在实际应用中，R 是未知的，但它可从观测数据一致估计，即用下式计算 NtHty1)(类似于 R 的特征分解，令 表示 R 的归一化特征向量，且按特征值降,1nmngs序排列，定义（3.3.9）)()()(aGfH= (3.3.10)SI一般地说，当
19、 nm-n 时使用（ 3.3.9） ；反之，使用(3.3.10)。使 最小的 n 个 值即为)(f所求。下面叙述 MUSIC 算法的步骤：（1） 利用观测值作出样本自相关矩阵 R。（2） 对 R 进行特征值分解，求出 n 以及 。nns,11 （3） 利用式（3.3.9）或(3.3.10)计算出 n 个 值。下面进行仿真实验：信号： )25.01()5.0()tjtjeenx噪声：高斯白噪声，方差为 1，均值（E）为 0İ 5 1 3 5 4 5 5 5 60 0 024 0 5 1 5 2 5 3 5 4 5 5 5 601234 4 0 0 0 0 0 0 0 024 İ 5 1 5 2 5 5 4 5 5 5 60 0 0 0 0 0 0 5 0 5 1 5 2 5 3 5 4 5 5 5 60 0 0 0 0 0 0 0 0 0 0 05
