現(xiàn)代數(shù)字信號(hào)處理-第三章-3-2016

1、1v譜分解定理譜分解定理 vAR模型法模型法vLevision-Durbin算法算法vAR模型的穩(wěn)定性及其階的確定模型的穩(wěn)定性及其階的確定vAR譜估計(jì)的性質(zhì)譜估計(jì)的性質(zhì)vAR模型參數(shù)提取方法模型參數(shù)提取方法vAR譜估計(jì)的異常現(xiàn)象及其補(bǔ)救措施譜估計(jì)的異常現(xiàn)象及其補(bǔ)救措施3.2 AR模型模型參數(shù)譜估計(jì)法步驟參數(shù)譜估計(jì)法步驟l為被估計(jì)的隨機(jī)過程選定一個(gè)合理的模型，這有賴于對(duì)隨機(jī)過程進(jìn)行的理論分析和實(shí)驗(yàn)l根據(jù)已知觀測(cè)數(shù)據(jù)估計(jì)模型的參數(shù)，這涉及各種算法的研究l用估計(jì)得到的模型參數(shù)計(jì)算功率譜。2322221*1*2*22)()()()()()()()()1()()(ARMA)()()(0)(jjjjxxx

2、xkkkkkkeAeBeHeSzAzAzBzBzHzHzSabzazbzAzBzHnu在下列關(guān)系：輸出功率和輸入功率存）系數(shù)。為后饋（）系數(shù)，為前饋（其線性系統(tǒng)函數(shù)：的白噪聲序列，方差為均值為輸入為4譜分解定理譜分解定理l任何實(shí)平穩(wěn)隨機(jī)信號(hào)任何實(shí)平穩(wěn)隨機(jī)信號(hào)ynyn的有理功率譜的有理功率譜S Syyyy(z)(z)都可唯一都可唯一地表示成下列最小相位形式地表示成下列最小相位形式 式中，式中， 為常系數(shù)，為常系數(shù)，B(z)B(z)是有理函數(shù)，即是有理函數(shù)，即 其中，其中，N(z)N(z)和和D(z)D(z)都是最小相位多項(xiàng)式。都是最小相位多項(xiàng)式。l譜分解定理保證了平穩(wěn)隨機(jī)信號(hào)模型的存在。任何平譜

3、分解定理保證了平穩(wěn)隨機(jī)信號(hào)模型的存在。任何平穩(wěn)隨機(jī)信號(hào)穩(wěn)隨機(jī)信號(hào)y yn n可以看成是由白噪聲序列可以看成是由白噪聲序列 激勵(lì)一個(gè)激勵(lì)一個(gè)因果穩(wěn)定的線性時(shí)不變系統(tǒng)因果穩(wěn)定的線性時(shí)不變系統(tǒng)B(z)B(z)產(chǎn)生的輸出。產(chǎn)生的輸出。21( )( ) ()yySzB z B z2( )( )( )N zB zD z n3.2.1 3.2.1 AR模型法模型法lP階AR模型差分方程5的逆濾波器。模型就是預(yù)測(cè)誤差濾波器）（)(AR)(A)(),()(.)1 (1)(A00)()()(1)(.Xmin221p21111nA zHzneEkaaZpaZazmmkmrakmramrzAzHXaXapkppkAk

4、pkAkAnpAnpAn6v基本公式基本公式 基于基于ARAR模型的譜估計(jì)由下式計(jì)算：模型的譜估計(jì)由下式計(jì)算：)()()(12zAzAzSxx212221)()(pkkjkjjxxeaeAeS故要求知道：模型的階數(shù)故要求知道：模型的階數(shù)p和和p個(gè)參數(shù)以及激勵(lì)源方差個(gè)參數(shù)以及激勵(lì)源方差2v基本思路基本思路 把這些參數(shù)與已知或估計(jì)的自相關(guān)函數(shù)聯(lián)系起來，把這些參數(shù)與已知或估計(jì)的自相關(guān)函數(shù)聯(lián)系起來，構(gòu)成著名的構(gòu)成著名的Yule-Walker方程，迭代求解該方程得到新方程，迭代求解該方程得到新的參數(shù)的參數(shù)7AR模型法模型法v 推導(dǎo)推導(dǎo) 推導(dǎo)方法推導(dǎo)方法 - 通過對(duì)通過對(duì)Sxx(z)求求z反變換來獲得反變

5、換來獲得Y-Wa方程方程 - 直接根據(jù)模型的差分方程導(dǎo)出直接根據(jù)模型的差分方程導(dǎo)出 具體推導(dǎo)具體推導(dǎo) 考慮考慮) 1 (/1)(0pkkkzazH其差分方程為1)()()(01anuknxanxpkk從而1( )( ) ()() ( ) ()xxpk xxkrmE x n x n ma rm kE x n u n m8AR模型法模型法l具體推導(dǎo)（續(xù)）具體推導(dǎo)（續(xù)） 考慮上式第二項(xiàng)的計(jì)算。設(shè)考慮上式第二項(xiàng)的計(jì)算。設(shè)ARAR模型的脈沖響應(yīng)為模型的脈沖響應(yīng)為h(n),h(n),在在方差為方差為 白噪聲作用下產(chǎn)生白噪聲作用下產(chǎn)生輸出輸出x(n)x(n)，故有，故有 20,00,)()()(22mmmh

6、mnunxE于是有于是有)2(0, )(0,)()(121mkmramkmramrxxpkkxxpkkxx這就是著名的這就是著名的Y-W方程。方程。0)m(h0:時(shí)注m9AR模型法模型法l具體推導(dǎo)（續(xù)）具體推導(dǎo)（續(xù)） 為求為求AR模型參數(shù)，應(yīng)先由模型參數(shù)，應(yīng)先由(2)式第二式選擇式第二式選擇m0的的p個(gè)方程求個(gè)方程求出出p個(gè)模型參數(shù)，然后代入第一個(gè)方程求出個(gè)模型參數(shù)，然后代入第一個(gè)方程求出 。 設(shè)已知自相關(guān)函數(shù)的頭設(shè)已知自相關(guān)函數(shù)的頭p+1個(gè)值為個(gè)值為r(0),r(1),r(p)，則，則(2) 式可表示為式可表示為2) 3(001)0()2() 1()() 1() 1 ()0() 1 ()()

7、2() 1 ()0(21paarprprprprrrrprrrr10Levision-Durbin算法算法l用線性方程組的常用解法(例如高斯消元法)求解式(3)，運(yùn)算量較大，利用系數(shù)矩陣的特征可構(gòu)成一些高效算法，Levision-Durbin算法是其中最著名、應(yīng)用最廣泛的一種。它是一種按階次進(jìn)行遞推的算法。lLevision-Durbin算法的關(guān)鍵是要推導(dǎo)出由AR(k)模型的參數(shù)計(jì)算AR(k+1)模型的參數(shù)的迭代計(jì)算公式。有三種方法： 1.根據(jù)AR(1)、AR(2)、AR(3)各階模型的Y-W方程的求解結(jié)果歸納出一般的迭代計(jì)算公式。 2.引入缺口函數(shù)的概念。 3.矩陣擴(kuò)充法。下面將介紹用這種方

8、法推導(dǎo)Levision算法。11Levision-Durbin算法算法l設(shè)已求得k階Y-W方程 的參數(shù) ，現(xiàn)求解k+1階Y-W方程v Levision算法的推導(dǎo)算法的推導(dǎo)2,1,1(0)(1)(2)( )(1)(0)(1)(1)0( )(1)(2)(0)0kkk krrrr karrrr kar kr kr kr 2,1,2,kkk kkaaa211,11,1,11(0)(1)( )(1)(1)(0)(1)( )0( )(1)(0)(1)0(1)( )(1)(0)0kkkkkkrrr kr karrr kr kar kr krrar kr krr 12Levision-Durbin算法算法l為

9、此，將k階方程的系數(shù)矩陣增加一列和增加一行，成為下列形式的“擴(kuò)大方程” 擴(kuò)大方程中的Dk由下式定義v Levision算法的推導(dǎo)算法的推導(dǎo)2,1,(0)(1)( )(1)1(1)(0)(1)( )0( )(1)(0)(1)0(1)( )(1)(0)0kkk kkrrr kr krrr kr kar kr krrar kr krrD ,00(1),1kkk ikiDa r kia BACK13Levision-Durbin算法算法l利用系數(shù)矩陣的Toeplitz性質(zhì)，將擴(kuò)大方程的行倒序，同時(shí)列也倒序，得到下列“預(yù)備方程”l將待求解的k+1階Y-W方程的解表示成擴(kuò)大方程的解和預(yù)備方程的解的線性組合

10、形式v LevisionLevision算法的推導(dǎo)算法的推導(dǎo),12(0)(1)( )(1)00(1)(0)(1)( )0( )(1)(0)(1)(1)( )(1)(0)1kk kkkDrrr kr krrr kr kar kr krrar kr krr 1,1,1,11,11,111001kkk kkkkk kkkkaaaaaaa14Levision-Durbin算法算法l或 式中， 是待定系數(shù)，稱為反射系數(shù)。上式各項(xiàng)都右乘以k+1階系數(shù)矩陣，得到 由該式可求出v Levision算法的推導(dǎo)算法的推導(dǎo)1,1,1,1,2,kik ikk kiaaaik 1k221120000kkkkkkDD12

11、kkkD2222111(1)kkkkkkD15Levision-Durbin算法算法l由擴(kuò)大方程的第一個(gè)方程可求出l從上面的推導(dǎo)中可歸納出如下由k階模型參數(shù)求k+1階模型參數(shù)的計(jì)算公式：l對(duì)于AR(p)模型，遞推計(jì)算直到k+1=p為止。v Levision算法的推導(dǎo)算法的推導(dǎo)2,1(0)( )kkk iira r i2,1,0012222111,1,11,11(0)( )(1),1(1),1,2, ;kkk iikkk ikikkkkkkkik ikk kikkkra r iDa r kiaDaaaika 功率譜求解步驟功率譜求解步驟l對(duì) x(n) 求系數(shù) a(1).a(p)l計(jì)算A， 得到H

12、=1/A(z)l將一個(gè)方差為的白噪聲通過該系統(tǒng)，所得即為待估計(jì)信號(hào)X(n)的功率譜。l注意： lA是白化濾波器， 但是輸出不總是,只有當(dāng)X(n)是一個(gè)AR過程時(shí)，輸出白噪聲才成立lp,p階FIR預(yù)測(cè)誤差濾波器退化為IIR預(yù)測(cè)誤差濾波器。 16174.2.4 AR模型的穩(wěn)定性及其階的確定模型的穩(wěn)定性及其階的確定AR模型的穩(wěn)定性及其階的確定lAR(p)模型穩(wěn)定的充分必要條件是H(z)的極點(diǎn)(即A(z)的根)都在單位圓內(nèi)。穩(wěn)定的AR(p)模型將具有以下性質(zhì)：l(1)H(z)的全部極點(diǎn)或A(z)的所有根都在單位圓內(nèi)。l(2)自相關(guān)矩陣是正定的。l(3)激勵(lì)信號(hào)的方差(能量)隨階次增加而遞減，即 。l(

13、4)反射系數(shù)的模恒小于1，即 。222120p1,1,2,.,kkp18AR模型的穩(wěn)定性及其階的確定模型的穩(wěn)定性及其階的確定l通常事先并不知道AR模型的階。階選得太低，功率譜受到的平滑太厲害，平滑后的譜可能會(huì)分辨不出真實(shí)譜中的兩個(gè)峰。階選得太高，固然會(huì)提高譜估計(jì)的分辨率，但同時(shí)會(huì)產(chǎn)生虛假的譜峰或譜的細(xì)節(jié)。l因此，要估計(jì)AR(p)過程，就應(yīng)該把AR(k)模型的階選得等于或大于p,即 ，但k不能太大。當(dāng)選擇 時(shí)，如果自相關(guān)函數(shù)的估計(jì)是精確的，那么AR(k)模型參數(shù)的估計(jì)為l式中， 是模型參數(shù)的精確值。這樣，用AR(k)模型能夠得到AR(p)過程的精確譜估計(jì)( )。v討論如何確定討論如何確定AR模型

14、的階模型的階kpkp,1,2,0,1,2,p ik iaipaippk,p iakp19AR模型的穩(wěn)定性及其階的確定模型的穩(wěn)定性及其階的確定l實(shí)際上自相關(guān)函數(shù)估計(jì)是有誤差的，那么AR(k)模型的階究竟選擇得偏高好還是偏低好？這主要應(yīng)從譜估計(jì)的質(zhì)量來考慮。lAR模型譜估計(jì)方法，既要估計(jì)AR模型參數(shù)，又要估計(jì)模型的階。l一種簡(jiǎn)單而直觀的確定AR模型的階的方法，是不斷增加模型的階，同時(shí)觀察預(yù)測(cè)誤差功率，當(dāng)其下降到最小時(shí)，對(duì)應(yīng)的階便可選定為模型的階。l另一種簡(jiǎn)單方法是觀察各階模型預(yù)測(cè)誤差序列的周期圖，當(dāng)它最接近于平坦(白色譜)時(shí)即對(duì)應(yīng)于最佳的階。另外還可用其他幾種不同的誤差準(zhǔn)則來作為確定模型階的依據(jù),

15、見圖p130。v討論如何確定討論如何確定AR模型的階模型的階l準(zhǔn)則l最終預(yù)測(cè)誤差準(zhǔn)則lAIC準(zhǔn)則l自回歸傳遞函數(shù)準(zhǔn)則(CAT)l經(jīng)驗(yàn)法： N/3 N/2 之間20214.2.3 AR譜估計(jì)的性質(zhì)譜估計(jì)的性質(zhì)l能夠?qū)ψ韵嚓P(guān)函數(shù)進(jìn)行外推，是AR譜估計(jì)分辨率高的根本原因。l設(shè)要估計(jì)一個(gè)AR(p)過程的譜，已知它的自相關(guān)函數(shù)的p+1個(gè)取樣值的估計(jì)值為 ，將它們代入Y-W方程，用Levision算法求解，得到AR(p)模型參數(shù)的估計(jì)值，然后將其代入譜計(jì)算公式，便得到AR(p)過程的譜估計(jì)，即l另一方面，譜與自相關(guān)序列之間存在傅立葉變換關(guān)系，可以得到vAR譜估計(jì)隱含著自相關(guān)函數(shù)的外推譜估計(jì)隱含著自相關(guān)函數(shù)

16、的外推(0), (1), ( )rrr p12*( )( )()ARSzA z A z( )( )mARmSzr m z(1)(2)22l比較以上(1)和(2)二式，可得到l假定濾波器 是因果的，且有h(0)=1，則上式左端 。l上式右端 ， 這里 是 的系數(shù)。l因此得到l對(duì)于 ，上式變?yōu)閘或?qū)懗?這里假設(shè)vAR譜估計(jì)隱含著自相關(guān)函數(shù)的外推譜估計(jì)隱含著自相關(guān)函數(shù)的外推1211*( )( )()mmA zr m zA z 1( )( )H zA z2*2()( ),0hmm m 0()()( ) (),0pla mr ma l r mlm( )a m( )A z20()( ) (),0plma

17、l r mlm mp0( ) ()0,pla l r mlmp1( )( ) (),plr ma l r ml mp (0)1a(3) 見(2)(4)參見(2)23AR譜估計(jì)的性質(zhì)譜估計(jì)的性質(zhì)l對(duì)于m=0,1,2,p，式(3)為l或?qū)懗蒷可以看出，式(5)與前面的Y-W方程相同，只是這里用 代替了Y-W方程中的 。式(4)說明對(duì)于在mp范圍內(nèi)的 值，現(xiàn)在我們并沒有認(rèn)為它們等于零，而認(rèn)為它們的值應(yīng)按該式進(jìn)行外推。vAR譜估計(jì)隱含著自相關(guān)函數(shù)的外推譜估計(jì)隱含著自相關(guān)函數(shù)的外推20 ,0( ) ()0,1,2,plma l r mlmp211( ) (),0()( ) (),1,2,plpla l

18、r mlmr ma l r mlmp( )r m( )r m( )r m(5)24AR譜估計(jì)的性質(zhì)譜估計(jì)的性質(zhì)l最大熵譜估計(jì)是基于將一段已知的自相關(guān)序列進(jìn)行明顯地外推，以得到未知的自相關(guān)取樣值，從而去除因?qū)ψ韵嚓P(guān)序列加窗而使譜估計(jì)特性變壞的弊端。l若已知 如何外推求得r(p+1),r(p+2),?最合理的方法：外推后的自相關(guān)序列所對(duì)應(yīng)的時(shí)間序列應(yīng)當(dāng)具有最大熵。這意味著，在具有已知的p+1個(gè)自相關(guān)取樣值的所有時(shí)間序列中，該時(shí)間序列將是最隨機(jī)或不可預(yù)測(cè)的，或者說它的譜是最平坦或最白的。通過由這樣的外推得到的自相關(guān)序列求出的譜稱為最大熵譜估計(jì)(MESE)。l選擇最大熵準(zhǔn)則的合理性在于：對(duì)未知自相關(guān)值

19、所加的約束最少，因而對(duì)應(yīng)的時(shí)間序列的隨機(jī)性越大，故可得到一個(gè)具有最小偏差的解。vAR譜估計(jì)與最大熵譜估計(jì)譜估計(jì)與最大熵譜估計(jì)(MESE)等效等效(0), (1), ( )rrr p25AR譜估計(jì)的性質(zhì)譜估計(jì)的性質(zhì)l高斯隨機(jī)過程最大熵譜估計(jì)可由下式表示 式中，a(m)可根據(jù)Y-W方程由已知的p1個(gè)自相關(guān)函數(shù)取樣值求取。l由上式可以看出，在已知 的情況下，對(duì)于高斯隨機(jī)過程，MESE與AR(p)是等效的。vAR譜估計(jì)與最大熵譜估計(jì)等效譜估計(jì)與最大熵譜估計(jì)等效2221( )1( )MESEpjfmmSfa m e(0), (1), ( )rrr p26AR譜估計(jì)的性質(zhì)譜估計(jì)的性質(zhì) AR(p)參數(shù)可以作

20、為參數(shù)可以作為p階線性預(yù)測(cè)系數(shù)來求取，準(zhǔn)則是階線性預(yù)測(cè)系數(shù)來求取，準(zhǔn)則是使預(yù)測(cè)誤差功率最小。使預(yù)測(cè)誤差功率最小。x(n)的線性預(yù)測(cè)值為的線性預(yù)測(cè)值為預(yù)測(cè)誤差為預(yù)測(cè)誤差為預(yù)測(cè)誤差功率為預(yù)測(cè)誤差功率為根據(jù)梯度為根據(jù)梯度為0，可以得到與，可以得到與Y-W方程類似的方程。方程類似的方程。 vAR譜估計(jì)等效于預(yù)測(cè)濾波譜估計(jì)等效于預(yù)測(cè)濾波1( )()pkkx na x nk ( )( )( )e nx nx n2( )minE en11(0)(1)(2)( )(1)(0)(1)(1)0( )(1)(2)(0)0prrrr parrrr par pr pr pr 27AR譜估計(jì)的性質(zhì)譜估計(jì)的性質(zhì)l容易知道，

21、當(dāng) 時(shí)， 剛好對(duì)應(yīng)于最佳預(yù)測(cè)誤差濾波器。lAR模型法和預(yù)測(cè)誤差濾波法互為逆濾波，AR模型是把一個(gè)白噪聲序列作為系統(tǒng)的激勵(lì)得到x(n)，而預(yù)測(cè)誤差濾波可以看作是把x(n)作為系統(tǒng)的激勵(lì)而得到一個(gè)白噪聲序列的輸出，所以預(yù)測(cè)誤差濾波器也是白化濾波器，它去掉了AR過程的相關(guān)性，從而在輸出端得到白噪聲。vAR譜估計(jì)等效于預(yù)測(cè)濾波譜估計(jì)等效于預(yù)測(cè)濾波kkaa228等同于線性預(yù)測(cè)等同于線性預(yù)測(cè) 102, 01minpkpkx nkx nke nx nx nk x nkEenk ARAR模型參數(shù)與線性預(yù)測(cè)器參數(shù)相同模型參數(shù)與線性預(yù)測(cè)器參數(shù)相同29等同于最優(yōu)白化濾波等同于最優(yōu)白化濾波1 21 21 21 2ex

22、plnthe geometric mean of , 01the arithmetic mean of xxxxxxxxxSfdfSfdfSfSfAR模型參數(shù)也可以通過最大化預(yù)測(cè)誤差濾波器模型參數(shù)也可以通過最大化預(yù)測(cè)誤差濾波器Prediction Error Filter (PEF)輸出信號(hào)的譜平坦度輸出信號(hào)的譜平坦度spectral flatness來獲得。來獲得。PEF(0)max min(0)(0)xxexeeeeRRR預(yù)測(cè)誤差譜平坦度預(yù)測(cè)誤差譜平坦度最大等效于最小相最大等效于最小相位濾波器位濾波器30AR 譜估計(jì)可以被看作是最優(yōu)白化操作的結(jié)譜估計(jì)可以被看作是最優(yōu)白化操作的結(jié)果果( )e

23、 n( )x n111/ ( ) 1/(1( )pkA zk zAR modelPrediction Error Filter31格形濾波器格形濾波器l將AR(p)模型參數(shù)ak,i看成一個(gè)序列(這里i=0,1,k)，并用多項(xiàng)式表示為 它的倒序多項(xiàng)式為 于是，Levinson算法中，由AR(k)的參數(shù)ak,i計(jì)算AR(k+1)的參數(shù)ak+1,i的遞推公式 可用多項(xiàng)式表示為l或,00( ),1kikk ikiA za za1( )()RkkkAzzA z1,1,11,11,1,2, ;kik ikk kikkkaaaika 111( )( )( )RkkkkAzA zz Az111( )( )(

24、)RRkkkkAzz AzA z(1)(2)32格形濾波器格形濾波器l將式(1)和式(2)寫成矩陣形式 上式由低階到高階進(jìn)行遞推運(yùn)算，稱為前向遞推運(yùn)算。lx(n)前向預(yù)測(cè)值為 前向預(yù)測(cè)誤差為 或lx(n-k)的后向預(yù)測(cè)值為 后向預(yù)測(cè)誤差為 或111111( )( )1( )( )kkkRRkkkAzA zzAzAzz,1( )()kk iix na x n i,1()()kk iix n ka x n ki (4)(3)(6)(5)(7),00( )( )( )(),1kkk ikie nx nx na x n ia( )( )( )kkEzA z X z,00 0( )()()()(),1k

25、kkk ik k ikiie nx n kx n ka x n k iax n ia ( )( ) ( )RkkE zA z X zBACK33格型濾波器的性質(zhì)格型濾波器的性質(zhì)l各級(jí)參數(shù)(反射系數(shù))的模小于1，一般情況下可以保證濾波器穩(wěn)定l級(jí)間是去耦的，因此當(dāng)各級(jí)分布調(diào)到最佳時(shí)可以使濾波器達(dá)到全局最佳。343.2.5 AR模型參數(shù)提取方法模型參數(shù)提取方法l在實(shí)際應(yīng)用中，常需根據(jù)信號(hào)的有限個(gè)取樣值來估計(jì)AR模型的參數(shù)，應(yīng)用較多的有三種方法：Y-W法或自相關(guān)法；協(xié)方差法；Burg法。這三種方法都可以用由時(shí)間平均代替集合平均的最小平方準(zhǔn)則推到得到。l理論上，AR模型參數(shù)是根據(jù)預(yù)測(cè)誤差功率最小的準(zhǔn)則來

26、確定的，該準(zhǔn)則表示為 或l即有 式中， 是x(n)的p1階自相關(guān)矩陣，R,1,1Tpp paaa(9)BACK2( )minpEen 2( )minpEen 22( )( )TppEenEena Ra35從有限數(shù)據(jù)集合提取模型參數(shù)從有限數(shù)據(jù)集合提取模型參數(shù)(1) 前向預(yù)測(cè)誤差和后向預(yù)測(cè)誤差前向預(yù)測(cè)誤差和后向預(yù)測(cè)誤差.(2) Yule-Walker法法(3) 協(xié)方差法協(xié)方差法(4) Burg法法理論上，理論上，ARAR模型參數(shù)是根據(jù)模型參數(shù)是根據(jù)預(yù)測(cè)誤差功率預(yù)測(cè)誤差功率最小的準(zhǔn)則確定最小的準(zhǔn)則確定 2minminnen e n e n 120minNpnen 12minNnpen 122minN

27、npenen數(shù)據(jù)加窗，效率高數(shù)據(jù)加窗，效率高,保證保證PEF最小相位最小相位不加窗，效率高不加窗，效率高, ,潛在不潛在不穩(wěn)定因素穩(wěn)定因素不加窗，更多數(shù)據(jù)不加窗，更多數(shù)據(jù)- -更好的估計(jì)和更低誤差；更好的估計(jì)和更低誤差；最小化復(fù)合全局誤差。最小化復(fù)合全局誤差。（自相關(guān)法）（自相關(guān)法）3611xnxnxnpxnp,1,1pppppaaa ,1,1,pppppaaa ,1pp iix nax ni ,1pp iix npax npi ,00( )( ), 1ppp ipienx nx nax nia ,00 ,00()(), 1, 1ppp ipippi replace ip pipienx np

28、x npax npiaax nia 前向預(yù)測(cè)誤差前向預(yù)測(cè)誤差：后向預(yù)測(cè)誤差后向預(yù)測(cè)誤差37AR模型參數(shù)提取方法模型參數(shù)提取方法l用最小平方的時(shí)間平均準(zhǔn)則代替集合平均準(zhǔn)則，有 或 上式中， 可由式(4)得到l 是由長(zhǎng)度為p1的預(yù)測(cè)誤差濾波器沖激響應(yīng)序列1,ap,1,ap,p與長(zhǎng)度為N的數(shù)據(jù)序列x(0),x(1),x(N-1)進(jìn)行卷積得到的，所以它的長(zhǎng)度為N+p，這決定了上式中的求和項(xiàng)數(shù)。vY-W法（自相關(guān)法）法（自相關(guān)法） ,00( )(),1ppp ipienax nia(11)(10)BACK1201( )minNppnenN 120( )minNppnen ( )pen( )pen38AR

29、模型參數(shù)提取方法模型參數(shù)提取方法l將式(11)代入式(10)，得 式中， 是由取樣自相關(guān)序列 構(gòu)成的N階取樣自相關(guān)矩陣。式(12)和式(9)等效，只是用取樣自相關(guān)矩陣 取代了自相關(guān)矩陣的 。l因此，用時(shí)間平均最小化準(zhǔn)則同樣可以導(dǎo)出Y-W方程組，只是方程組中的 要用 取代。vY-W法法R(12)RRRR10( )( ) (),01Nkkr kx n x nkkN 12,0,0( )()NppTpp ip jni jena r ij a a Ra39 0 0 0 (1) 1 1 (0) 0 0 0 x Nx n x nx n px n px ,1,1,1 pppppaaa ,1,1,1 pp pp

30、 paaa ,1,1,1 pp pp paaa (1 ) (2) (1 ) 1 1 ( ) (1 ) (1 ) (0)xNxNxN px n x nx n px n pxpxpx x ,1,1,1 pp pp paaa ,1,1,1 pp pp paaa ,1,1,1 pp pp paaa 自相關(guān)法自相關(guān)法(Yule-Walker法法)計(jì)算前向預(yù)測(cè)誤差原理計(jì)算前向預(yù)測(cè)誤差原理: :協(xié)方差法計(jì)算前向預(yù)測(cè)誤差原理協(xié)方差法計(jì)算前向預(yù)測(cè)誤差原理: :(0)e( )e n(1 )e N p ( )e p( )e n(1 )e N兩端加零兩端加零兩端不加零兩端不加零兩端加零兩端加零注意點(diǎn)：注意點(diǎn)： 40l

31、41AR模型參數(shù)提取方法模型參數(shù)提取方法l用下列時(shí)間平均的最小平方準(zhǔn)則代替集合平均的最小平方準(zhǔn)則： 該式與自相關(guān)法的主要區(qū)別是求和范圍不同?，F(xiàn)在的求和范圍是p(N-1)。這意味著沒有“加數(shù)據(jù)窗”的不合理假設(shè)。l與式(12)類似，可推導(dǎo)出 其中，自相關(guān)矩陣的估計(jì)為 這里，自相關(guān)序列的估計(jì)為l一般情況下， 不是Toeplitz的， 半正定，這是與自相關(guān)法不同的, 無法保證最小相位，無法保證極點(diǎn)都在單位圓內(nèi)。v協(xié)方差法協(xié)方差法R( , )r i jR1( , )() ()Nn pr i jx ni x nj(13)BACK12( )minNpn pen12( )NTpn pena Ra42AR模型參

32、數(shù)提取方法模型參數(shù)提取方法l自相關(guān)法的計(jì)算效率高，且能保證預(yù)測(cè)誤差濾波器是最小相位的，但數(shù)據(jù)兩端要附加零取樣值，實(shí)際上是數(shù)據(jù)加窗，這將使參數(shù)估計(jì)的精度下降。特別當(dāng)數(shù)據(jù)段很短時(shí)，加窗效應(yīng)就更加嚴(yán)重。協(xié)方差法計(jì)算效率也高，但潛在著不穩(wěn)定因素。自相關(guān)法和協(xié)方差法都是直接估計(jì)AR參數(shù)。lBurg法則一方面希望利用已知數(shù)據(jù)段兩端以外的未知數(shù)據(jù)(但它對(duì)這些未知數(shù)據(jù)不作主觀臆測(cè))，另一方面又總是設(shè)法保證使預(yù)測(cè)誤差濾波器是最小相位的。Burg法與自相關(guān)法和協(xié)方差法不同，它不直接估計(jì)AR參數(shù)，而是先估計(jì)反射系數(shù)，然后利用Levinson遞推算法由反射系數(shù)來求得AR參數(shù)。vBurg法法43AR模型參數(shù)提取方法模型

33、參數(shù)提取方法lBurg法首先要估計(jì)反射系數(shù)，所使用的準(zhǔn)則是前向和后向預(yù)測(cè)誤差功率估計(jì)的平均值最小準(zhǔn)則。在這里，預(yù)測(cè)誤差功率估計(jì)仍然用時(shí)間平均來代替集合平均。Burg法估計(jì)反射系數(shù)的準(zhǔn)則表示為l該式的求和范圍與協(xié)方差法相同。前向和后向預(yù)測(cè)誤差濾波器的工作都是在數(shù)據(jù)段上進(jìn)行的(數(shù)據(jù)段兩端不需要補(bǔ)充零)。l由式(14)，求 對(duì) 的偏導(dǎo)數(shù)并令其等于0，即vBurg法法(修正協(xié)方差法修正協(xié)方差法)(14)p(15)122( )( )minNppn penen1( )( )2( )( )0Nppppn ppppe ne ne ne n44AR模型參數(shù)提取方法模型參數(shù)提取方法l由于濾波器運(yùn)算未超出已知數(shù)據(jù)段

34、的范圍，因此，式(8b)對(duì)pnN-1是成立的，由式(15)得到 或者寫成 由上式解出 可以證明 ，這就保證了預(yù)測(cè)誤差濾波器具有最小相位性質(zhì)。vBurgBurg法法111( )(1)( )( )0Nppppn pen enen en1111111( )(1)(1)(1)( )( )0Nppppppppn penenenenen en111122112( )(1)( )(1)Nppn ppNppn pen enenen1p45AR模型參數(shù)提取方法模型參數(shù)提取方法l綜上，Burg法估計(jì)AR(p)模型參數(shù)的具體計(jì)算步驟如下：l(1)確定初始條件l(2)確定k-1階AR參數(shù)(迭代計(jì)算時(shí)，k的值從1開始選

35、取) 、 ( )l(3)利用下式計(jì)算vBurg法法00( )( )( ),01e ne nxnn N 122001( )NnxnN1( )kAz21k11knN 111122112( )(1)( )(1)Nkkn kkNkkn ken enenenk46AR模型參數(shù)提取方法模型參數(shù)提取方法l(4)用下式計(jì)算l(5)用下式計(jì)算 和 ( )l(6)計(jì)算k階均方誤差l(7)回到步驟(2)，進(jìn)行下一次迭代。vBurg法法( )kA z,1,1,1,2,1;k ikikkk ik kkaaaika 111( )( )1( )(1)kkkkkkenenenen( )ken( )ken1k n N 2221

36、(1)kkk47AR譜估計(jì)的異?，F(xiàn)象及其補(bǔ)救措施譜估計(jì)的異?，F(xiàn)象及其補(bǔ)救措施l虛假譜峰 補(bǔ)救措施：模型的階不宜選得太高，最高不應(yīng)超過數(shù)據(jù)記錄長(zhǎng)度的一半。l譜線分裂 補(bǔ)救措施：同時(shí)調(diào)整所有反射系數(shù)，使預(yù)測(cè)誤差功率真正達(dá)到最小。l噪聲對(duì)AR譜估計(jì)的影響 AR譜估計(jì)對(duì)觀測(cè)噪聲比較敏感，噪聲會(huì)使譜峰展寬，從而導(dǎo)致分辨率下降，而且會(huì)使譜峰偏離正確的位置。 為減小噪聲對(duì)AR譜估計(jì)的惡化影響，一般可使用下列四種方法:(1)采用ARMA譜估計(jì)方法；(2)對(duì)數(shù)據(jù)進(jìn)行濾波，減小噪聲；(3)采用高階AR模型；(4)補(bǔ)償自相關(guān)函數(shù)或反射系數(shù)估計(jì)中噪聲的影響。3.3 3.3 最大熵譜估計(jì)法最大熵譜估計(jì)法l48*(k)(

37、e )(x)110 ,(x)2(e )(k)(e )11ee0,(k)2(e )為最大化，令， 則上式變?yōu)閑jwxjwexejwjkwjkwxjwexrPHdwkprPrPdwkprPl49l502221112212211221(e )()(u)e(p+1)P() 1e e(p+1)PppjwmemHpHpwwppHwwPe aPeP3.4 3.4 最大似然估計(jì)法最大似然估計(jì)法l51Maximum Likelihoodl523.5 互協(xié)方差估計(jì)與互譜估計(jì)互協(xié)方差估計(jì)與互譜估計(jì)l533.6 特征分解法譜估計(jì)特征分解法譜估計(jì)l54l55求解過程求解過程l56二階二階l57l58211(e)jwMH

38、iiipPa e v3.6 Pisarenko 諧波分解諧波分解l59l用Pisarenko 求解噪聲中估計(jì)兩個(gè)復(fù)指數(shù)分量頻率和功率，假設(shè)兩個(gè)復(fù)指數(shù)加噪聲的隨機(jī)過程自相關(guān)序列的前三個(gè)值為60(0)6xr(1)1.927054.58522xrj(2)3.427053.49541xrj l6161.927054.585223.427053.495411.927054.5852261.927054.585223.427053.49541 1.927054.585226xjjRjjjj1=15.89512=1.10493=1.00001230.57630.00000.27400.65180.27850

39、.30060.22440.53420.00010.01000.32090.74920.40340.41160.28300.64800.40970.0058jjjVv vvjjjjjjl最小的特征值為1， 其對(duì)應(yīng)的特征矢量l求根得到62min0.27850.30060.32090.74920.40970.0058jvjj /310.50.8660jzje2 /520.30990.9511jzje13w225wlV1,V2 對(duì)應(yīng)的功率功率方程：其中 ： 得到 63121()2.9685jwV e221()2.9861jwV e222()0.0139jwV e122()0.0315jwV e2112

40、222.96852.98610.03150.0139wwPP 2min1w12P 23P l不實(shí)用原因l已知復(fù)指數(shù)個(gè)數(shù)l假設(shè)噪聲為加性白噪聲l需要求解自相關(guān)矩陣的最小特征值和特征矢量，計(jì)算量高l改進(jìn)： MUSIC算法6465帶噪聲經(jīng)過MMSE處理66四、四、MA模型譜估計(jì)模型譜估計(jì)67RMSE68/htakeda/kernelreg/kernelreg.htm697071MA模型譜估計(jì)模型譜估計(jì)lMA譜估計(jì)以全零點(diǎn)模型為基礎(chǔ)，將其用于估計(jì)窄帶譜時(shí)得不到高分辨率，但用于MA隨機(jī)過程時(shí)，由于MA隨機(jī)過程的功率譜本身具有寬峰窄谷的特點(diǎn)，故能得到精確估計(jì)

41、。l對(duì)于ARMA模型差分方程 其中，輸入激勵(lì)u(n)是均值為0、方差為 的白噪聲序列。輸出功率譜和輸入功率譜之間存在下列關(guān)系010( )()(),1pqkkkkx na x nkb u nka 設(shè)2*2*2*1( )()1( )( )()1( )()xxB z BzSzH z HzA z Az(1)72MA模型譜估計(jì)模型譜估計(jì)l由式(1)可得 對(duì)上式兩端取逆Z變換，分別得到 這里，假設(shè)h(n)是實(shí)序列。由上兩式得到lh(n)是因果序列，即n0時(shí)h(n)0，故上式右端有*22*1()1( ) ( )( )() ( )1()xxBzSz A zB zHB zzAz10( ) ( )( )*()pxxxxmk xxkZSz A zrmaa rmk12*22*01() ( )()*()qmkkZHB zhmbb h kmz200()()pqk xxkkka rmkb h km(2)00()( ),0,1,()0,1qq mqkkmk mkkkb h kmbh kmqb h kmmq或(3)73MA模型譜估計(jì)模型譜估計(jì)l將式(3)代入式(2)，得 這就是ARMA模型參數(shù)與自相關(guān)函數(shù)之間的關(guān)系式。l當(dāng)a0=1且ak=0(k=1,2,p)時(shí)，由式(4)可得出MA模型參數(shù)與信號(hào)的自相關(guān)函數(shù)之間的關(guān)系式。注意此時(shí)h(k)=bk,故有2101()( ),0,1,( )(),1pq m