功率譜估計(jì)2012-1

1、*1五、功率譜估計(jì)*2內(nèi)容n經(jīng)典譜估計(jì)n基于參數(shù)模型的功率譜估計(jì)n基于特征分解法的譜估計(jì)*35.1 概述n意義信號(hào)分析工具（分析信號(hào)頻域特征）q確定性信號(hào)：Fourier變換頻譜特性n絕對(duì)可積q周期信號(hào)nFourier級(jí)數(shù)q隨機(jī)信號(hào)：不滿足絕對(duì)可積條件n能量無(wú)限，平均功率有限n頻域分析：功率譜密度*45.1 概述n功率譜定義TTjwtTdtetxX)()(TXESTTx2)(lim)(2NNnnjjNenxeX)()(NeXEeSjNNjx2)(lim)(2連續(xù)隨機(jī)信號(hào)離散隨機(jī)信號(hào)*55.1 概述n廣義平穩(wěn)隨機(jī)信號(hào)的功率譜與其自相關(guān)函數(shù)的關(guān)系 ：維納辛欽（Wiener-Khinchin）定理

2、mnjxemRS)()(*65.1 概述n定義q譜分析：利用有限的樣本數(shù)據(jù)估計(jì)一個(gè)平穩(wěn)隨機(jī)信號(hào)的功率譜密度。q從觀察樣本估計(jì)信號(hào)的功率分布情況n應(yīng)用q尋找數(shù)據(jù)的“隱周期性”q語(yǔ)音處理n基音提取n估計(jì)聲音語(yǔ)調(diào)共振峰n基于譜減法的語(yǔ)音增強(qiáng)算法 *75.1 概述n應(yīng)用q生物工程n功率譜密度的峰形和波形巓癇病的發(fā)作周期q無(wú)源聲納信號(hào)處理n功率譜密度位置目標(biāo)的方向q雷達(dá)信號(hào)處理：回波信號(hào)的功率譜密度分析n功率譜密度峰值的寬度運(yùn)動(dòng)目標(biāo)的位置n功率譜密度峰值的高度運(yùn)動(dòng)目標(biāo)的強(qiáng)度n功率譜密度峰值的位置運(yùn)動(dòng)目標(biāo)的速度*8應(yīng)用例子應(yīng)用例子1，腦電波分析腦電波信號(hào)1（正常的腦電波）腦電波信號(hào)2（癲癇病人腦電波）*9

3、信號(hào)1兩種不同方法的功率譜估計(jì)周期圖AR模型*10信號(hào)2兩種不同方法的功率譜估計(jì)周期圖AR模型*112,人體磁場(chǎng)探測(cè)周期圖BT不同方法得到的功率譜估計(jì)：不同方法得到的功率譜估計(jì)：*12AR模型MUSIC*133,地震序列*144,信道估計(jì)123456-0.6-0.4-0.200.20.40.60.81123456-0.6-0.4-0.200.20.40.60.8輸入H(z)輸出噪聲)(zH真實(shí)值估計(jì)值0dB噪聲10dB噪聲)()()(2XHY*155.1 概述n發(fā)展歷史q19世紀(jì)末（1899），舒斯特（Schuster）周期圖q1930年，維納：維納辛欽定理q20世紀(jì)50年代，BT法q1967

4、年，Burg最大熵譜法q1968年，ParzenAR譜估計(jì)法q1971年，Van Den Bos證明最大熵譜法與AR譜估計(jì)等效q1969年，Capon提出最大似然譜估計(jì)法*165.1 概述n分類q經(jīng)典譜估計(jì)法、現(xiàn)代譜估計(jì)法q線性譜估計(jì)法、非線性譜估計(jì)法：線性、非線性運(yùn)算q非參量方法、參量方法*17譜估計(jì)算法一覽譜估計(jì)算法現(xiàn)代譜估計(jì)經(jīng)典譜估計(jì)周期圖PeriodogamBTAR,ARMA模型最大熵最大似然基于特征值分解Pisarenko諧波分解MUSIC*185.1 概述n譜估計(jì)的質(zhì)量q偏差q方差q譜分辨率*195.1 概述n功率譜估計(jì)的問(wèn)題n影響估計(jì)質(zhì)量的因素q估計(jì)方法 q信號(hào)模型 q外推方法

5、*205.2經(jīng)典譜估計(jì)的基本方法n周期圖法（直接法）n相關(guān)圖法（間接法）*21周期圖法n概念q定義（Periodogram）NeXEeSjNNjx2)(lim)(2對(duì)遍歷信號(hào)：21)02)(1)(1)(limlimNnnjNjNNjwxenxNeXNeS22)()2() 1 (1)(NjjjeNxexexNS*22周期圖法n有限的觀察值2102)(1)(1)(NnnjjNxenxNeXNS1)1()()(NNmmjxemRS*23 )1()1(1)1(101010)(10102)()()()(1)()(1)()(1)(1)(NNmjmxNNmNkjmNNNiNkkijNNNkjkNNijiNj

6、xemRkimemkxkxNekikxkxNekxeixNeXNS10)()(1)(NnNNxmnxnxNmR其中)()();()(xxxSSmRmRN證明：*24估計(jì)的質(zhì)量n偏差q功率譜估計(jì)1)1()()(NNmmjxxemRESE1)1()()(NNmjmxxemRS估計(jì)的均值：)()(mREFSExx*25估計(jì)的偏差由3.2.2：)()()(|)(mRmwmRNmNmRExxx其中：NmNmNmNmw|01|0|)(*26估計(jì)的偏差dWSWSSExxx)()(21)()()()()()(mRmwFSExx所以：可見，功率譜估計(jì)的均值是真實(shí)功率譜與Bartlett窗的Fourier變換的卷

7、積1)1(2)()2sin()2sin(1)(NNmmjemwNNW*27估計(jì)的偏差：W(),N=800.20.40.60.81-100-50050Magnitude (dB)Normalized Frequency*28估計(jì)的偏差n結(jié)論)()();()(ExxxSSEmRmRN漸近無(wú)偏估計(jì)實(shí)際中數(shù)據(jù)長(zhǎng)度N有限，偏差總是存在：)()()(21)(xxxSdWSSbia例子)()*2sin()(nunfnx*29估計(jì)的質(zhì)量n估計(jì)的方差q功率譜估計(jì)在兩個(gè)頻率處的協(xié)方差，定義為)()()()()(),(221121xxxxxSESSESESSCov220102201021)()()(21)()()(

8、21)(),(dDDSNdDDSNSSCov*30估計(jì)的方差)(, 01|, 1)(00DotherNnndF其中：212200)()()()(21)(SEdDDSNSVar0)()(ar2xxSESVN結(jié)論：不是一致估計(jì)，估計(jì)的方差比要大)(xS2| )(|xS*31窗函數(shù)的影響1，對(duì)估計(jì)方差的影響2200)()()()(21)(xxSEdDDSNSVar假設(shè)窗函數(shù)旁瓣為0，主瓣寬度為B)2(2BB當(dāng)頻率落在以下區(qū)間：有：0)()(00DD2)()(xxSESVar結(jié)論：方差在特定頻率范圍達(dá)到最小*32窗函數(shù)的影響2，對(duì)估計(jì)協(xié)方差的影響當(dāng)頻率滿足以下關(guān)系：有：結(jié)論：主瓣寬度越小，譜估計(jì)起伏越

9、大（相關(guān)度為0意味估計(jì)值差異）220102201021)()()(21)()()(21)(),(dDDSNdDDSNSSCovB|210)()(0)()(20102010DDDD0)(),(21SSCov*3321)(10D)(20D(d)(c)(20D02B)(0D2B2B(b)(0D)(0D(a)21)(10D02B02B2B02B2B2B2B*34估計(jì)的方差n隨著N的增大，周期圖的譜估計(jì)起伏增大*35估計(jì)的分辨率n定義：q真實(shí)譜中兩個(gè)靠得很近的譜峰能被分辨的能力dWSWSSExxx)()(21)()()(分辨率取決于W()的主瓣寬度例子*36相關(guān)圖法n1958年由Blackman與Tuk

10、ey提出n改善相關(guān)函數(shù)的估計(jì)方法，對(duì)周期圖平滑處理改善方差性能n周期圖中自相關(guān)函數(shù)的估計(jì)nm接近N時(shí)，估計(jì)出的自相關(guān)函數(shù)的偏差較大n根據(jù)維納辛欽定理mnjemRS)()(mNnmnxnxNmR10|)|()(1)(*37相關(guān)圖法n步驟q估計(jì)自相關(guān)函數(shù)) 11(, |)|()(1)(|10NmNmnxnxNmRmNnMMmmjemvmRS)()()(q對(duì)自相關(guān)函數(shù)進(jìn)行Fourier變換|M|=N-1窗函數(shù)*38相關(guān)圖法n周期圖法利用觀測(cè)序列直接給出譜估計(jì)，又稱直接法。與之比較，BT譜估計(jì)先估計(jì)出自相關(guān)函數(shù)，再估計(jì)功率譜，又稱間接法。BT譜估計(jì)基于自相關(guān)函數(shù)，也稱相關(guān)圖法)()()(VSSPERB

11、T*39周期圖法和相關(guān)圖法的比較n當(dāng)M=N-1時(shí)，BT法與周期圖法估計(jì)出的功率譜是相同 n當(dāng)M快速算法nToeplitz矩陣的性質(zhì) qHermitian對(duì)稱 q任一條對(duì)角線上的元素都相等 n利用Toeplitz矩陣性質(zhì)，快速的遞推算法) 0() 1()() 1() 0() 1 ()() 1 () 0(xxxxxxxxxRpRpRpRRRpRRRR*685.3.2 AR參數(shù)譜估計(jì)與最佳線性預(yù)測(cè)器的關(guān)系n問(wèn)題：q已知k階的參數(shù)，怎樣用k階的參數(shù)求出k+1階的參數(shù)?nLevinson-Durbin算法q階數(shù)遞推 : 0p書上使用了線性預(yù)測(cè)模型的概念推導(dǎo)*691 線性預(yù)測(cè)的基本概念利用x(n-p),x

12、(n-p+1),x(n-1)的線性組合預(yù)測(cè)x(n)pkkfknxnx1)()(*70最佳線性預(yù)測(cè)器)()()(nxnxnef預(yù)測(cè)誤差：212)()()(pkkknxnxEneE預(yù)測(cè)均方誤差：*71正交性原理xe01xy2x投影最佳線性預(yù)測(cè)器，最小0mpmnxnxmnxEf, 2, 1, 0)()()(pkxkxkmRmR1)()()() 0()()()(1minkRRnxnxnxExpkkxf最小預(yù)測(cè)均方誤差：0)()(nemnxE*72最佳線性預(yù)測(cè)器n線性預(yù)測(cè)WienerHopf方程：ak和min001)0()1()()1()0()1()()1()0(1pxxxxxxxxxRpRpRpRRR

13、pRRR與AR模型的YuleWalker方程相似*73線性預(yù)測(cè)與AR模型關(guān)系n對(duì)于同一個(gè)隨機(jī)信號(hào)：min2, 2, 1,pkkk結(jié)論：p階的AR模型可以構(gòu)造一個(gè)p階的最佳線性預(yù)測(cè)器對(duì)于p階的AR過(guò)程信號(hào)x(n)，應(yīng)用p階線性預(yù)測(cè)器，預(yù)測(cè)誤差為：)()()()()()(1nuknxnxnxnxnepkkf預(yù)測(cè)誤差是白噪聲，預(yù)測(cè)過(guò)程又稱為白化*74線性預(yù)測(cè)與AR模型關(guān)系)(2th)(nx)(1zA)(nu(a)(nx)(ne)(zA(b)(1zA)(ne)( nx)( nx)(nx(c)ppzaza1111ppzaza111ppzaza11*75前、后向線性預(yù)測(cè)的關(guān)系前向預(yù)測(cè)：利用p個(gè)n時(shí)刻之前的

14、數(shù)據(jù)預(yù)測(cè)x(n)；后向預(yù)測(cè)：利用某時(shí)刻n-p之后p個(gè)數(shù)據(jù)預(yù)測(cè)x(n-p)；*76前、后向線性預(yù)測(cè)的關(guān)系21)()()()()()()(neEnxnxneknxkanxffffpkffForwardPrediction21)()()()()()()(neEpnxpnxnekpnxkapnxbbbbpkbbBackwardPrediction*77前、后向線性預(yù)測(cè)的關(guān)系利用正交性原理，可以得到后向預(yù)測(cè)的WeinerHopf方程：0bmbpmkmRkamRkRkaRpkxbxpkxbxb,2,1, )()()()()()0(11min比較前、后向預(yù)測(cè)的WeinerHopf方程，利用Toeplitz矩

15、陣fbminmin)()(kakabf*785.3.3 Levinson-Durbin算法n遞推算法：q低階AR參數(shù) 高階AR模型參數(shù)nAR模型與線性預(yù)測(cè)模型的參數(shù)一致q遞推線性預(yù)測(cè)系數(shù) AR模型參數(shù)nm階的線性預(yù)測(cè)系數(shù)：am1,am2,ammfbmminmin預(yù)測(cè)誤差功率（最小預(yù)測(cè)均方誤差）：*795.3.3 Levinson-Durbin算法后向零階遞推)(nx) 1( nx)(nx1階遞推前向) 2( nx) 1( nx)(nx前向前向后向后向002階遞推*805.3.3 Levinson-Durbin算法p=1時(shí)，011aR01)0() 1 () 1 ()0(111aRRRR21102

16、1111111)0()0() 1 (aaRRRa解方程組)0()()()()(22000RnxEneEnxneff*815.3.3 Levinson-Durbin算法定義前向預(yù)測(cè)誤差和后向預(yù)測(cè)誤差的相關(guān)系數(shù)m+1為：)1()(1neneEbmfmmmimimimRamR11)1() 1(根據(jù)最佳預(yù)測(cè)的正交原理：*82證明：m階=m+1階001)0() 1()() 1()0() 1 ()() 1 ()0(1mmmmaaRmRmRmRRRmRRR0mmaRm階要求m+1階：0001)0() 1 ()() 1() 1 ()0() 1()()() 1()0() 1 () 1()() 1 ()0(11,

17、 1, 11 , 1mmmmmmaaaRRmRmRRRmRmRmRmRRRmRmRRRM+1階*83110001)0() 1 ()() 1() 1 ()0() 1()()() 1()0() 1 () 1()() 1 ()0(mmmmmaaRRmRmRRRmRmRmRmRRRmRmRRR對(duì)自相關(guān)矩陣進(jìn)行行和列的擴(kuò)展mimimimRamR11)1() 1(M+1階*84顛倒方程組的行次序：mmmmmaamRmRRRmRmRRRRRmRmRRRmRmR0001) 1()() 1 ()0()() 1()0() 1 () 1 ()0() 1()()0() 1 ()() 1(11顛倒方程組的列次序：mmm

18、mmaaRRmRmRRRmRmRmRmRRRmRmRRR0010)0() 1 ()() 1() 1 ()0() 1()()() 1()0() 1 () 1()() 1 ()0(11*85要求m+1階：0001)0() 1 ()() 1() 1 ()0() 1()()() 1()0() 1 () 1()() 1 ()0(11, 1, 11 , 1mmmmmmaaaRRmRmRRRmRmRmRmRRRmRmRRR引入反射系數(shù)km+1mmimkmmmimRamRK111)1() 1(M+1階*86mmmmmmmmmmmmKaaKaaRRmRmRRRmRmRmRmRRRmRmRRR00001001)0

19、() 1 ()() 1() 1 ()0() 1()()() 1()0() 1 () 1()() 1 ()0(111111對(duì)比m+1階：0001)0() 1 ()() 1() 1 ()0() 1()()() 1()0() 1 () 1()() 1 ()0(11, 1, 11 , 1mmmmmmaaaRRmRmRRRmRmRmRmRRRmRmRRRmmkmkmmmmmmkmmmmkkmmmmkmRamRkkmkakaaka1112111,1, 111, 1)1() 1()1 (, 1,*875.4.1 AR譜估計(jì)的相關(guān)函數(shù)法)0()()(0)()(220000RnxEneEKnxneff*885.

20、4.1 AR譜估計(jì)的相關(guān)函數(shù)法)0()1 ()1 ()0() 1 () 1()()(210211111111RKKRRKanxanxnef*895.4.1 AR譜估計(jì)的相關(guān)函數(shù)法1222111222112112112122212)1 () 1 ()2()1 ()()2() 1()()(KRaRKaaKaKananxanxanxnef*90Levinson-Durbin算法步驟由觀察數(shù)據(jù)x(n)，選擇合適的模型階數(shù)p，估計(jì)自相關(guān)函數(shù)R(m);利用Levinson-Durbin算法，遞推求出ap1,ap2,app, p22121)(pkkjpkpAReaSmmkmkmmmmmmmmmkmmmmkk

21、mkmRamRKKKamkaKaamKnxER11121111,11,1,1020)1()1()1(1p1,2,.,0)()0(*91結(jié)論n運(yùn)算量q遞推過(guò)程的每一步運(yùn)算量為O(m)q所有遞推過(guò)程的計(jì)算量為O(m2)n反射系數(shù)的性質(zhì)：qm=m20q1 1-km20q| km | 1q當(dāng)| km |=0時(shí)，遞推算法停止qm2 m-12 12 02q隨著模型階數(shù)的增加，預(yù)測(cè)誤差功率逐漸減少q從最小均方誤差準(zhǔn)則的意義，算法收斂*92線性預(yù)測(cè)的格型濾波器結(jié)構(gòu)前、后向預(yù)測(cè)誤差的關(guān)系：初始條件為)()()(00nxnenebf)() 1()() 1()()(1111nekneneneknenefmmbmbm

22、bmmfmfmpmneEneneEneEneneEKbmbmfmfmbmfmm,.,2 , 1,) 1() 1()()() 1()(21112111*93線性預(yù)測(cè)的格型濾波器結(jié)構(gòu)濾波器的系數(shù)僅是反射系數(shù)優(yōu)點(diǎn)：模塊化較小的舍入誤差保證穩(wěn)定性1k)(nebp)( nx1z)(nx2kpk1z1z)(nefp)(2nef)(1nef)(0nef)(0neb)(1neb)(2neb*945.4.1 AR譜估計(jì)的相關(guān)函數(shù)法nMATLAB函數(shù)qLevinson-Durbin遞歸算法 na,e,k=levinson(r,n) *955.3.4 Burg算法nBurg算法q自相關(guān)函數(shù)估計(jì)不準(zhǔn)確不估計(jì)q通過(guò)前、

23、后向線性預(yù)測(cè)系數(shù)之間的關(guān)系，由觀察數(shù)據(jù)直接求反射系數(shù)q分辨率較自相關(guān)函數(shù)法好q缺點(diǎn)：譜線分裂（本來(lái)只有一個(gè)譜線的位置分裂為2條譜線），對(duì)于高階模型可能產(chǎn)生虛假峰值*96相關(guān)函數(shù)法的基本思路：0m212)()()(pkkknxnxEneE預(yù)測(cè)均方誤差最小x(n)(mR)(ARS*970mk12)(1NmnbmbmnemNbmfmm2112)(1NmnfmfmnemNBurg算法的基本思路：預(yù)測(cè)均方誤差最小前、后向線性預(yù)測(cè)均方誤差之和最小*98)()1()()1()()(1111nekneneneknenefmmbmbmbmmfmfm0mkpmnenenenekNmnbmNmnfmNmnbmfmm

24、, 2, 1,) 1()() 1()(21211211112121121| ) 1(| ) 1(| |1 Nemekbmfmmmmmkmm,21m*99)|1(1112, 1, 1,mmmmmmkmmmkmkmkkmkk步驟：0階初始條件，求出和1階AR參數(shù)遞推計(jì)算m階的線性預(yù)測(cè)誤差和反射系數(shù)；計(jì)算m階時(shí)的AR模型系數(shù)重復(fù)上述過(guò)程，直到m=p，求出所有階次的AR參數(shù)計(jì)算功率譜1k*100n特點(diǎn)q不需估計(jì)自相關(guān)函數(shù)q分辨率比自相關(guān)函數(shù)法更好，但會(huì)出現(xiàn)“譜線分裂”，高階產(chǎn)生虛假峰值qN較小時(shí)，性能不亞于Levinson算法q滿足|ki|1，AR模型總是穩(wěn)定的q對(duì)于有噪聲的正弦信號(hào)，Burg算法存在

25、對(duì)正弦初相位的敏感問(wèn)題*101修正協(xié)方差法譜估計(jì))(21bffb121)()()(1NpnpkffknxkanxpNpNnpkbbknxkanxpN1021)()()(1)()(kakafb0)(iafbpiiaiaf, 2, 1),()( 不同*1020)()()()()()()()(110111pNnpkNpnpkinxknxkanxinxknxkanxpNpNnNpnpkpNnNpninxnxinxnxinxknxinxknxka1011101)()()()()()()()()( pNnNpnxinxknxknxinxpNkic101)()()()()(21),(*103)0 ,()0

26、, 2()0 , 1 ()( )2( ) 1 ( ),()2 ,() 1 ,(), 2()2 , 2() 1 , 2(), 1 ()2 , 1 () 1 , 1 (pcccpaaappcpcpcpcccpcccxxxxxxxxxxxxpNnpkNpnpknxknxkanxnxknxkanxpN10111min)()()( )()()()( )()(21pkxxkckac1min), 0()( )0, 0(不是Toeplitz矩陣協(xié)方差法的正則方程*1045.3.5 AR譜估計(jì)的性質(zhì)nAR譜比經(jīng)典譜平滑nAR譜分辨率比經(jīng)典譜要高 n隨著階數(shù)的增加，譜與真實(shí)譜就越接近 n方差反比于N和信噪比 n如

27、果自相關(guān)矩陣正定，則由Yule-Walker方程求出的AR模型是穩(wěn)定的*105AR譜估計(jì)的不足n信號(hào)的信噪比關(guān)系較大，信噪比低，則方差大，分辨率低 n如果信號(hào)是含噪聲的正弦信號(hào)，其譜峰易受信號(hào)初相位的影響，并且可能出現(xiàn)“譜線分裂”的現(xiàn)象 n譜的質(zhì)量受階數(shù)p的影響大，p取值小，則過(guò)于平滑，精度不夠，p太大，則可能會(huì)產(chǎn)生虛假的譜峰 *106AR模型階數(shù)的選擇最終預(yù)測(cè)誤差準(zhǔn)則（FPE）kkNkNkFPE)(AIC信息論準(zhǔn)則kNkAICx2ln)(CAT準(zhǔn)則kkiiNkCAT111)(1階數(shù)太小：過(guò)于平滑，分辨率不足階數(shù)太大：產(chǎn)生虛假譜峰，譜線分裂懲罰項(xiàng)*1075.3.6 MA模型及其正則方程qkkn

28、ukbnunx1)()()()(212)(1)(qkkjjekbeSqmqmmkbkbmkbkbmrmqkqmkx, 0, 1, 0, )()()()()(022qqmmjxjjMAjemreBePeS)()()()(22qqmmjxjBTemreS)()(非線性方程，很難直接求解*108MA模型nMA模型譜估計(jì)等效于經(jīng)典譜估計(jì)的相關(guān)圖法，分辨率較低n單純?yōu)榱俗V估計(jì)，沒(méi)有必要求解MA模型nMA模型參數(shù)的求解方法q譜分解法q高階的AR模型近似MA模型q最大似然估計(jì)法（最小二乘法）*1095.7.2 MA模型參數(shù)求解方法1)(11)(1)(kkzkazAzH)()(1)()(1zBzkbzHzHqkkq1)()(zBzA)()()()(1mkmakbmaqkz反變換由柯爾莫可洛夫定理式中a(0)=1，m0實(shí)際上，e(m)不為0 )( map*111MA模型參數(shù)求解步驟n由觀測(cè)數(shù)據(jù)建立一個(gè)p階的AR模型，pqn求p階AR模型系數(shù)n利用AR模型參數(shù)建立一個(gè)q階的線性預(yù)測(cè)，再次利用AR系數(shù)求解方法，求解b(k)n計(jì)算功率譜212)(1)(qkkjjekbeSAIC準(zhǔn)則用于階次的判斷：qNqAICMA2)ln()(*1125.7.3 ARMA譜估計(jì)qmkmrkaqkk