隨機(jī)信號的功率譜估計(jì)

1、1第三章 隨機(jī)信號的功率譜估計(jì)鄭寶玉2內(nèi)內(nèi) 容容3v Levinson算法算法v Berg算法算法4) 1 (001)0()2() 1()() 1() 1 ()0() 1 ()()2() 1 ()0(1paarprprprprrrrprrrrMEM的核心是求解如下方程：的核心是求解如下方程：這個(gè)方程實(shí)際上是聯(lián)合這個(gè)方程實(shí)際上是聯(lián)合AR模型法和預(yù)測濾波法得出的。模型法和預(yù)測濾波法得出的。我們發(fā)現(xiàn)，方程我們發(fā)現(xiàn)，方程(1)有如下特點(diǎn)：有如下特點(diǎn)： 系數(shù)矩陣是一個(gè)系數(shù)矩陣是一個(gè)Toplitz矩陣，利用矩陣，利用Toplitz矩陣的性質(zhì)矩陣的性質(zhì) 可簡化方程求解?？珊喕匠糖蠼?。 實(shí)際問題中，一般只知

2、道信號的某些觀測值，而不知道實(shí)際問題中，一般只知道信號的某些觀測值，而不知道 其其AR模型階數(shù)，該階數(shù)也需要在方程求解過程中找到。模型階數(shù)，該階數(shù)也需要在方程求解過程中找到。 下面介紹兩種算法。下面介紹兩種算法。v引言引言5v 原理原理 假設(shè)已得到假設(shè)已得到k階線性預(yù)測系數(shù)階線性預(yù)測系數(shù)(預(yù)測濾波參數(shù)預(yù)測濾波參數(shù))，我們，我們來考慮求來考慮求k1階濾波參數(shù)。階濾波參數(shù)。k階濾波參數(shù)的矩陣方程為階濾波參數(shù)的矩陣方程為)2(001)0()2() 1()() 1() 1 ()0() 1 ()()2() 1 ()0()()()(1kkkkaarkrkrkrkrrrrkrrrr由于系數(shù)矩陣的由于系數(shù)矩陣

3、的Toplitz性質(zhì)，上式又有如下形式性質(zhì)，上式又有如下形式)2(001)0()2() 1()() 1() 1 ()0() 1 ()()2() 1 ()0()()(1)(kkkkaarkrkrkrkrrrrkrrrr6現(xiàn)考慮模型階數(shù)增加現(xiàn)考慮模型階數(shù)增加1, 即從即從k變?yōu)樽優(yōu)閗+1的情況。的情況。對于對于k+1模型模型, 有有)3(0001)0() 1 ()() 1() 1 ()0() 1()()() 1()0() 1 () 1()() 1 ()0()1()1(1)1()1(1)1(kkkkkkkaaarrkrkrrrkrkrkrkkrrkrkrrrR由于系數(shù)矩陣的由于系數(shù)矩陣的Toplit

4、zToplitz性質(zhì)，性質(zhì)，k+1k+1階系數(shù)矩陣階系數(shù)矩陣 可有兩種分塊形式。可有兩種分塊形式。 )1( kR7 利用這個(gè)性質(zhì)，可設(shè)利用這個(gè)性質(zhì)，可設(shè))4(001001)()()(1)()()(1)1(kkkkkkkkkaaaaR式中式中)5()1() 1(1) 1 ()() 1(1)()()(1kjkjTkkkjkrakraarkrkr8 比較比較(3)和和(4)，可知，當(dāng)，可知，當(dāng)(3)與與(4)等效；且有下列兩個(gè)遞推關(guān)系式：等效；且有下列兩個(gè)遞推關(guān)系式：)6(0)(k即當(dāng)下式成立時(shí)即當(dāng)下式成立時(shí))7(/ )1() 1(/)(1)()(kkjkjkjkrakr)8(10011)(1)()

5、()(1)1(1)1()1(1kkkkkkkkkkkaaaaaaa和)9()()1(kk由由(8)末式還可得末式還可得：)10()1(1kka(5)(10)構(gòu)成構(gòu)成Levinson算法基礎(chǔ)算法基礎(chǔ)。9 現(xiàn)用現(xiàn)用i表示遞推過程的階數(shù)，令表示遞推過程的階數(shù)，令i=k+1, 并設(shè)信號模型的并設(shè)信號模型的最大階數(shù)為最大階數(shù)為N，則有如下則有如下i=k+10，得，得 2) 置置i=k+11; 3) 由由(8)、(10)式計(jì)算式計(jì)算 4) 由由(7)、(10)式計(jì)算式計(jì)算 5)由由(6)、(7)和和(10)式計(jì)算式計(jì)算 6) 置置i =i+1; 7) 判別：若判別：若 轉(zhuǎn)轉(zhuǎn)3)；否則，結(jié)束程序。；否則，結(jié)

6、束程序。)0()0(r)11(/ )()()1(11)1()(iijijiijiraira)12() 1,.,1()1()()1()(ijaaaaijiiiijij)13()(1 2)()1()(iiiiaNi v 算法算法10 v 討論討論Levinson算法第算法第4步利用了一個(gè)重要遞推關(guān)系步利用了一個(gè)重要遞推關(guān)系(12)， 通常稱為通常稱為Levinson關(guān)系式關(guān)系式遞推過程產(chǎn)生一個(gè)濾波參數(shù)序列遞推過程產(chǎn)生一個(gè)濾波參數(shù)序列 通常稱為偏相關(guān)系數(shù)通常稱為偏相關(guān)系數(shù)遞推過程產(chǎn)生的遞推過程產(chǎn)生的 可用來監(jiān)視可用來監(jiān)視i階信號模型的均方階信號模型的均方 誤差估值。誤差估值。遞推結(jié)果的最終解為遞推結(jié)

7、果的最終解為 和和),.,1()(Niaii),.,1()(NjaNj)(N)(i遞推過程及結(jié)果遞推過程及結(jié)果11 v 討論討論 優(yōu)點(diǎn)優(yōu)點(diǎn)：計(jì)算簡單：計(jì)算簡單 缺點(diǎn)缺點(diǎn)：需根據(jù)有限觀測數(shù)據(jù)估計(jì)自相關(guān)序列：需根據(jù)有限觀測數(shù)據(jù)估計(jì)自相關(guān)序列r(n) 短數(shù)據(jù)序列時(shí)短數(shù)據(jù)序列時(shí)，自相關(guān)估計(jì)值誤差很大，引起預(yù)測，自相關(guān)估計(jì)值誤差很大，引起預(yù)測 濾波參數(shù)誤差，導(dǎo)致濾波參數(shù)誤差，導(dǎo)致“譜峰飄移譜峰飄移”和和“譜線分譜線分裂裂”( (即出即出 現(xiàn)虛假譜線現(xiàn)虛假譜線) ) 長數(shù)據(jù)序列時(shí)長數(shù)據(jù)序列時(shí)，自相關(guān)估計(jì)值雖精確，但計(jì)算量大。，自相關(guān)估計(jì)值雖精確，但計(jì)算量大。優(yōu)缺點(diǎn)優(yōu)缺點(diǎn)12 v前向預(yù)測與后向預(yù)測前向預(yù)測與

8、后向預(yù)測 考慮信號序列值：考慮信號序列值：)(, ) 1(,),2(),1(),(pnxpnxnxnxnx前向預(yù)測前向預(yù)測后向預(yù)測后向預(yù)測前向預(yù)測誤差：前向預(yù)測誤差：后向預(yù)測誤差：后向預(yù)測誤差：pkppkpkpkfaknxaknxanxnxnxne0011, )()()()( )()(pkpppkpkpkbbknxbknxbpnxnxnxne0101, )()()()( )()(pkpkba ,其中其中 分別為分別為p p階前、后向預(yù)測系數(shù)。階前、后向預(yù)測系數(shù)。 13 v Berg算法原理算法原理根據(jù)前面的基本概念，可知根據(jù)前面的基本概念，可知m m階前向預(yù)測誤差為階前向預(yù)測誤差為)14(.,

9、2 , 1, )()(0apminxanfmimim類似地，類似地， m m階后向預(yù)測誤差為階后向預(yù)測誤差為)14(.,2 , 1, )()(0*bpminxangmimimm再利用再利用LevinsonLevinson關(guān)系式：關(guān)系式：)15() 1,.,1( ,)()1()1()(mjaaaammmmjmmmjmj有有)16()() 1()()16() 1()()(1*111bnfngngangnfnfmmmmmmmm其中其中pmnxngnf,.,2 , 1; )()()(0014 v Berg Berg算法原理（續(xù)）算法原理（續(xù)）定義定義m m階前、后向預(yù)測誤差的功率為階前、后向預(yù)測誤差的

10、功率為)17()()(21P22pmnmmmngnf將將(16)代入代入(17)，并令，并令Pm對對 的偏導(dǎo)數(shù)為零的偏導(dǎo)數(shù)為零, ,得最佳得最佳mm)18(,.,2 , 1,) 1()() 1()(12*121*1pmngnfngnfpmnmmpmnmmm15 v Berg Berg算法算法 設(shè)已知有限數(shù)據(jù)序列設(shè)已知有限數(shù)據(jù)序列x(n),n=0,1,N，則可按下列步，則可按下列步驟計(jì)算預(yù)測濾波器系數(shù)，并在此基礎(chǔ)上計(jì)算功率譜。驟計(jì)算預(yù)測濾波器系數(shù)，并在此基礎(chǔ)上計(jì)算功率譜。1. 置置m=0, 計(jì)算初值計(jì)算初值)()()()(00200nxngnfnxPNn2. m=m+1,并按并按(19)計(jì)算反射

11、系數(shù)計(jì)算反射系數(shù)m3.3.計(jì)算濾波器系數(shù)：計(jì)算濾波器系數(shù)：) 1,.,1( ,)1()1()(mjaaamjmmmjmj4. 計(jì)算預(yù)測誤差功率計(jì)算預(yù)測誤差功率Pm：12)1 (mmmPP5.5.按按(16)(16)式計(jì)算濾波器輸出式計(jì)算濾波器輸出)(),(ngnfmm6. 置置m=m+1, 并重復(fù)步驟并重復(fù)步驟(2)-(5), 直到直到m=p。16 v Berg Berg算法（續(xù)）算法（續(xù)）最后，由最后，由Berg算法估計(jì)的濾波器系數(shù)算法估計(jì)的濾波器系數(shù)ppjPpja及),.,2 , 1( 計(jì)算功率譜密度：計(jì)算功率譜密度：21)(1)(pkjkpkpBurgeaPP17內(nèi)內(nèi) 容容18v 基本原

12、理基本原理v MV譜與譜與ME譜或譜或AR譜的關(guān)系譜的關(guān)系19v MVSE基本原理基本原理 三點(diǎn)說明三點(diǎn)說明 最小方差功率譜估計(jì)最小方差功率譜估計(jì)(MVSE)，又稱最大似然譜估，又稱最大似然譜估 計(jì)，但實(shí)際上它并不是最大似然譜估計(jì)；計(jì)，但實(shí)際上它并不是最大似然譜估計(jì)； 提出者提出者Capon,1969也把這個(gè)方法叫做高分辨率譜估也把這個(gè)方法叫做高分辨率譜估 計(jì)方法，但實(shí)際上其分辨率并不高于計(jì)方法，但實(shí)際上其分辨率并不高于AR模型法；模型法； 盡管這樣，但由于其思路獨(dú)特，仍有了解的必要。盡管這樣，但由于其思路獨(dú)特，仍有了解的必要。 下面，討論該方法的導(dǎo)出過程。下面，討論該方法的導(dǎo)出過程。20v

13、MVSE基本原理基本原理 算法推導(dǎo)算法推導(dǎo) 將隨機(jī)信號將隨機(jī)信號x(n)通過通過FIR濾波器濾波器A(z):) 1 ()()(0pkkzkazA則其輸出為則其輸出為)2()()()(*)()(0aXTpkknxkananxny其中其中TTpaaapnxnxnx)(),.,1 (),0()(),.,2(),1(aXy(n)的均方值，也就是的均方值，也就是y(n)的功率，由下式給出的功率，由下式給出：) 3()()()(*2aRaaXXapHHHHEnynyEnyE 式中 為r(0),r(1),r(p)構(gòu)成的Toeplitz矩陣。 若y(n)的均值為零，則 也是y(n)的方差。TpEXXR*21v

14、 MVSE基本原理基本原理 算法推導(dǎo)（續(xù)）算法推導(dǎo)（續(xù)） 求濾波器的系數(shù)，有兩個(gè)原則：求濾波器的系數(shù)，有兩個(gè)原則： 在某一給定頻率在某一給定頻率 處處,x(n),x(n)無失真通過無失真通過, ,這等效于要求：這等效于要求：)4(1)()exp()()(0aeiHpkiezkjkazAiji式中式中Tiiipjj)exp(),.,exp(, 1 )(e 在在 附近的頻率分量被拒絕附近的頻率分量被拒絕, ,即在即在 附近使附近使 為最小。為最小。 為同時(shí)滿足這兩個(gè)原則，必須滿足下式：為同時(shí)滿足這兩個(gè)原則，必須滿足下式：ii)5(1)(minaeaRaiHpH這就是這就是“最小方差最小方差”譜估計(jì)

15、譜估計(jì)的來歷。的來歷。22v MVSE基本原理基本原理 算法推導(dǎo)（續(xù)）算法推導(dǎo)（續(xù)） 利用利用Lagrange乘子法求解乘子法求解(5)式式,得最小方差濾波器系數(shù)為得最小方差濾波器系數(shù)為)()()(11ipiHipMVeReeRa)()(11ipiHMVeRe)()(1)(1eRepHMVP相應(yīng)的最小方差為相應(yīng)的最小方差為從而，估計(jì)的最小方差譜為從而，估計(jì)的最小方差譜為應(yīng)該注意應(yīng)該注意， 并不是真正意義上的功率譜，因?yàn)椴⒉皇钦嬲饬x上的功率譜，因?yàn)?對對 的積分并不等于信號的功率，但它描述了真正譜的的積分并不等于信號的功率，但它描述了真正譜的相對強(qiáng)度。相對強(qiáng)度。)(MVP)(MVP23v MV

16、譜與譜與AR譜的關(guān)系譜的關(guān)系對自相關(guān)矩陣的逆矩陣對自相關(guān)矩陣的逆矩陣 作作Cholesky分解，有分解，有)8(11HppppAPAR1pR其中其中 分別是分別是0階階p階階AR模型系數(shù)和激勵功率模型系數(shù)和激勵功率(即方差即方差)組成的矩陣，即組成的矩陣，即ppPA 、1) 1 ()2() 1()(1)2()3(1) 1 ()2(1) 1 (1ppppppppppaapapaaaaaaA2222120,ppP24v MV譜與譜與AR譜的關(guān)系譜的關(guān)系（續(xù)）（續(xù)）將將(8)(8)代入代入(7)(7)，得，得pkkkmkHpppHpHMVmjmaP020211/)exp()()()()()()(1e

17、APAeeRe于是，我們得到于是，我們得到MV譜與譜與AR譜之間的一個(gè)重要關(guān)系：譜之間的一個(gè)重要關(guān)系：pkARMVPP0)(1)(1其中kmkkARmjmaP022)exp()()(25內(nèi)內(nèi) 容容26v 自相關(guān)矩陣的特征分解自相關(guān)矩陣的特征分解v 基于子空間的頻率估計(jì)與信號估計(jì)基于子空間的頻率估計(jì)與信號估計(jì)27自相關(guān)矩陣的特征分解自相關(guān)矩陣的特征分解v 基本原理基本原理 設(shè)信號設(shè)信號x(n)由由M個(gè)復(fù)正弦加白噪聲組成個(gè)復(fù)正弦加白噪聲組成, 分別是第分別是第 i 個(gè)復(fù)正弦的功率和頻率個(gè)復(fù)正弦的功率和頻率, 則則x(n)的自相關(guān)函數(shù)為的自相關(guān)函數(shù)為iiA、) 1 ()(, )()exp()(21w

18、wwMiiixxkkjAkr式中正弦信號的幅度為式中正弦信號的幅度為 為白噪聲的功率。為白噪聲的功率。wiA, 如果由如果由( (p p+1)+1)個(gè)個(gè)rxx(n)組成自相關(guān)矩陣：組成自相關(guān)矩陣：)2()0() 1 () 1()() 1 ()0()2() 1() 1()2()0() 1 ()() 1() 1 ()0(*xxxxxxxxxxxxxxxxxxxxprrprprrrprprprprrrprprrrxxxxxxxxxxxxR28自相關(guān)矩陣的特征分解自相關(guān)矩陣的特征分解v基本原理（續(xù)）基本原理（續(xù)）且定義信號向量：且定義信號向量：)3(,.,2 , 1,)exp(),.,exp(, 1

19、MipjjTiiie則由則由(1)-(3), 有有)4(1MiwHiiipIeeAR 式中第一、二項(xiàng)分別為信號陣和噪聲陣式中第一、二項(xiàng)分別為信號陣和噪聲陣,前者最大秩為前者最大秩為M. 設(shè)設(shè) ，對，對Rp進(jìn)行特征分解得：進(jìn)行特征分解得：Mp )5()(111111pMiHiiwMiHiiwipiHiiwMiHiiipVVVVVVVVR式中式中Vi 是對應(yīng)于特征值是對應(yīng)于特征值 的特征向量的特征向量,特征向量相互正交，特征向量相互正交，即即i)(0)(1jijijHiVV29自相關(guān)矩陣的特征分解自相關(guān)矩陣的特征分解v 基本結(jié)論基本結(jié)論 從上面討論可以看出：從上面討論可以看出： Rp的所有特征向量

20、的所有特征向量V1, , Vp+1形成形成p+1維向量空間維向量空間(信息空間信息空間), 且且V1, , Vp+1相互正交。相互正交。 利用相關(guān)矩陣的特征值，可將利用相關(guān)矩陣的特征值，可將信息空間信息空間分成分成兩個(gè)子空間兩個(gè)子空間： 由主特征向量由主特征向量 V1, , VM 張成張成信號子空間信號子空間(主分量主分量) 其特征值分別為其特征值分別為 由特征向量由特征向量 VM+1, Vp+1 張成張成噪聲子空間噪聲子空間 其特征值均為其特征值均為 信號向量信號向量e1,eM和主特征向量和主特征向量V1, , VM張成相同的子空間張成相同的子空間 信號子空間。信號子空間。 結(jié)論結(jié)論:可在信

21、號子空間或噪聲子空間進(jìn)行譜估計(jì)和頻率估計(jì)可在信號子空間或噪聲子空間進(jìn)行譜估計(jì)和頻率估計(jì) 應(yīng)用應(yīng)用:借助噪聲子空間噪聲特性借助噪聲子空間噪聲特性,從信號子空間估計(jì)有用信號從信號子空間估計(jì)有用信號 下面考慮下面考慮pM和和pM兩種情況下基于噪聲子空間的估計(jì)問題兩種情況下基于噪聲子空間的估計(jì)問題wwMw,.,130v 自相關(guān)矩陣的特征分解自相關(guān)矩陣的特征分解v 基于子空間的頻率估計(jì)與信號估計(jì)基于子空間的頻率估計(jì)與信號估計(jì)- PHD方法方法- MUSIC方法方法31 v 理論基礎(chǔ)理論基礎(chǔ) 若若p=M, 則則(5)式中式中Rp僅有僅有一個(gè)噪聲向量一個(gè)噪聲向量VM+1,它所它所 對應(yīng)的對應(yīng)的 特征值就是噪

22、聲方差特征值就是噪聲方差 ，該特征值也是，該特征值也是Rp的最小特征的最小特征 值（因?yàn)榇藭r(shí)值（因?yàn)榇藭r(shí) ）?？梢宰C明，這時(shí)有）?？梢宰C明，這時(shí)有 定理定理 1 1 噪聲向量噪聲向量VM 與信號向量與信號向量ei(i1,M)都正交都正交,即即 w01p)6(,.,2 , 1,01MiMHiVe令)exp(ijz則根據(jù)定理則根據(jù)定理1 1有有)7(0)(v)V(11MkkMzkz其中其中TMMMM)(v),.,1 (v111V32 PHD算法的步驟算法的步驟 1)求求x(n)的自相關(guān)函數(shù)并構(gòu)成自相關(guān)矩陣的自相關(guān)函數(shù)并構(gòu)成自相關(guān)矩陣Rp,且設(shè)且設(shè) p=M2)對對Rp進(jìn)行特征分解進(jìn)行特征分解,得特征

23、值得特征值 及特征向量及特征向量 將其排序并找出最小的特征值將其排序并找出最小的特征值 及相應(yīng)的及相應(yīng)的3)將將 代入代入(7)，形成，形成 M 階多項(xiàng)式并求該多項(xiàng)式的根，階多項(xiàng)式并求該多項(xiàng)式的根， 得到得到x(n)的的M 個(gè)頻率個(gè)頻率4)由由(1)有有11,.,M1M1MV1MVM,.,1)8()()2() 1 ()exp()exp()exp()2exp()2exp()2exp()exp()exp()exp(21212121MrrrAAAjMjMjMjjjjjjxxxxxxMMMM5)再由(1)得： 故可求得故可求得)9()0(1MiwiAr)10()0(1MixxiwrA11,.,MVV3

24、3v理論基礎(chǔ)理論基礎(chǔ) 若噪聲子空間的向量不止一個(gè)，用類似的方法可以證明有若噪聲子空間的向量不止一個(gè)，用類似的方法可以證明有 定理2 信號向量ei與噪聲子空間的向量Vk都正交，即)10(1,.,1;,.,10,apMkMikiVe由于自相關(guān)矩陣Rp的特征向量 構(gòu)成一組正交基，因此有11,.,pVV)11(,.,.,11bspaneespanMMVV注意注意: (11)對應(yīng)于對應(yīng)于Mp的情況，在這種情況下，若再使用的情況，在這種情況下，若再使用(8), 則求出的則求出的V(z)將有將有p-M個(gè)多余零點(diǎn)。故不宜再使用個(gè)多余零點(diǎn)。故不宜再使用(8)計(jì)算。計(jì)算。34v MUSIC算法算法 基本思路基本思路 由于信號向量與噪聲子空間的各個(gè)向量都正交由于信號向量與噪聲子空間的各個(gè)向量都正交, 因此因此 信號向量與噪聲空間各向