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

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

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

3、分解定理保證了平穩(wěn)隨機信號模型的存在。任何平穩(wěn)隨機信號穩(wěn)隨機信號y yn n可以看成是由白噪聲序列可以看成是由白噪聲序列 激勵一個激勵一個因果穩(wěn)定的線性時不變系統(tǒng)因果穩(wěn)定的線性時不變系統(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ù)測誤差濾波器）（)(AR)(A)(),()(.)1 (1)(A00)()()(1)(.Xmin221p21111nA zHzneEkaaZpaZazmmkmrakmramrzAzHXaXapkppkAk

4、pkAkAnpAnpAn6v基本公式基本公式 基于基于ARAR模型的譜估計由下式計算：模型的譜估計由下式計算：)()()(12zAzAzSxx212221)()(pkkjkjjxxeaeAeS故要求知道：模型的階數(shù)故要求知道：模型的階數(shù)p和和p個參數(shù)以及激勵源方差個參數(shù)以及激勵源方差2v基本思路基本思路 把這些參數(shù)與已知或估計的自相關(guān)函數(shù)聯(lián)系起來，把這些參數(shù)與已知或估計的自相關(guān)函數(shù)聯(lián)系起來，構(gòu)成著名的構(gòu)成著名的Yule-Walker方程，迭代求解該方程得到新方程，迭代求解該方程得到新的參數(shù)的參數(shù)7AR模型法模型法v 推導(dǎo)推導(dǎo) 推導(dǎo)方法推導(dǎo)方法 - 通過對通過對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ù)） 考慮上式第二項的計算。設(shè)考慮上式第二項的計算。設(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:時注m9AR模型法模型法l具體推導(dǎo)（續(xù)）具體推導(dǎo)（續(xù)） 為求為求AR模型參數(shù)，應(yīng)先由模型參數(shù)，應(yīng)先由(2)式第二式選擇式第二式選擇m0的的p個方程求個方程求出出p個模型參數(shù)，然后代入第一個方程求出個模型參數(shù)，然后代入第一個方程求出 。 設(shè)已知自相關(guān)函數(shù)的頭設(shè)已知自相關(guān)函數(shù)的頭p+1個值為個值為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)，運算量較大，利用系數(shù)矩陣的特征可構(gòu)成一些高效算法，Levision-Durbin算法是其中最著名、應(yīng)用最廣泛的一種。它是一種按階次進行遞推的算法。lLevision-Durbin算法的關(guān)鍵是要推導(dǎo)出由AR(k)模型的參數(shù)計算AR(k+1)模型的參數(shù)的迭代計算公式。有三種方法： 1.根據(jù)AR(1)、AR(2)、AR(3)各階模型的Y-W方程的求解結(jié)果歸納出一般的迭代計算公式。 2.引入缺口函數(shù)的概念。 3.矩陣擴充法。下面將介紹用這種方

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ù)矩陣增加一列和增加一行，成為下列形式的“擴大方程” 擴大方程中的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ì)，將擴大方程的行倒序，同時列也倒序，得到下列“預(yù)備方程”l將待求解的k+1階Y-W方程的解表示成擴大方程的解和預(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ù)。上式各項都右乘以k+1階系數(shù)矩陣，得到 由該式可求出v Levision算法的推導(dǎo)算法的推導(dǎo)1,1,1,1,2,kik ikk kiaaaik 1k221120000kkkkkkDD12

11、kkkD2222111(1)kkkkkkD15Levision-Durbin算法算法l由擴大方程的第一個方程可求出l從上面的推導(dǎo)中可歸納出如下由k階模型參數(shù)求k+1階模型參數(shù)的計算公式：l對于AR(p)模型，遞推計算直到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對 x(n) 求系數(shù) a(1).a(p)l計算A， 得到H

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

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

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

15、見圖p130。v討論如何確定討論如何確定AR模型的階模型的階l準(zhǔn)則l最終預(yù)測誤差準(zhǔn)則lAIC準(zhǔn)則l自回歸傳遞函數(shù)準(zhǔn)則(CAT)l經(jīng)驗法： N/3 N/2 之間20214.2.3 AR譜估計的性質(zhì)譜估計的性質(zhì)l能夠?qū)ψ韵嚓P(guān)函數(shù)進行外推，是AR譜估計分辨率高的根本原因。l設(shè)要估計一個AR(p)過程的譜，已知它的自相關(guān)函數(shù)的p+1個取樣值的估計值為 ，將它們代入Y-W方程，用Levision算法求解，得到AR(p)模型參數(shù)的估計值，然后將其代入譜計算公式，便得到AR(p)過程的譜估計，即l另一方面，譜與自相關(guān)序列之間存在傅立葉變換關(guān)系，可以得到vAR譜估計隱含著自相關(guān)函數(shù)的外推譜估計隱含著自相關(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對于 ，上式變?yōu)閘或?qū)懗?這里假設(shè)vAR譜估計隱含著自相關(guān)函數(shù)的外推譜估計隱含著自相關(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譜估計的性質(zhì)譜估計的性質(zhì)l對于m=0,1,2,p，式(3)為l或?qū)懗蒷可以看出，式(5)與前面的Y-W方程相同，只是這里用 代替了Y-W方程中的 。式(4)說明對于在mp范圍內(nèi)的 值，現(xiàn)在我們并沒有認為它們等于零，而認為它們的值應(yīng)按該式進行外推。vAR譜估計隱含著自相關(guān)函數(shù)的外推譜估計隱含著自相關(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譜估計的性質(zhì)譜估計的性質(zhì)l最大熵譜估計是基于將一段已知的自相關(guān)序列進行明顯地外推，以得到未知的自相關(guān)取樣值，從而去除因?qū)ψ韵嚓P(guān)序列加窗而使譜估計特性變壞的弊端。l若已知 如何外推求得r(p+1),r(p+2),?最合理的方法：外推后的自相關(guān)序列所對應(yīng)的時間序列應(yīng)當(dāng)具有最大熵。這意味著，在具有已知的p+1個自相關(guān)取樣值的所有時間序列中，該時間序列將是最隨機或不可預(yù)測的，或者說它的譜是最平坦或最白的。通過由這樣的外推得到的自相關(guān)序列求出的譜稱為最大熵譜估計(MESE)。l選擇最大熵準(zhǔn)則的合理性在于：對未知自相關(guān)值

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

20、為參數(shù)可以作為p階線性預(yù)測系數(shù)來求取，準(zhǔn)則是階線性預(yù)測系數(shù)來求取，準(zhǔn)則是使預(yù)測誤差功率最小。使預(yù)測誤差功率最小。x(n)的線性預(yù)測值為的線性預(yù)測值為預(yù)測誤差為預(yù)測誤差為預(yù)測誤差功率為預(yù)測誤差功率為根據(jù)梯度為根據(jù)梯度為0，可以得到與，可以得到與Y-W方程類似的方程。方程類似的方程。 vAR譜估計等效于預(yù)測濾波譜估計等效于預(yù)測濾波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譜估計的性質(zhì)譜估計的性質(zhì)l容易知道，

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

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

23、 n( )x n111/ ( ) 1/(1( )pkA zk zAR modelPrediction Error Filter31格形濾波器格形濾波器l將AR(p)模型參數(shù)ak,i看成一個序列(這里i=0,1,k)，并用多項式表示為 它的倒序多項式為 于是，Levinson算法中，由AR(k)的參數(shù)ak,i計算AR(k+1)的參數(shù)ak+1,i的遞推公式 可用多項式表示為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)寫成矩陣形式 上式由低階到高階進行遞推運算，稱為前向遞推運算。lx(n)前向預(yù)測值為 前向預(yù)測誤差為 或lx(n-k)的后向預(yù)測值為 后向預(yù)測誤差為 或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各級參數(shù)(反射系數(shù))的模小于1，一般情況下可以保證濾波器穩(wěn)定l級間是去耦的，因此當(dāng)各級分布調(diào)到最佳時可以使濾波器達到全局最佳。343.2.5 AR模型參數(shù)提取方法模型參數(shù)提取方法l在實際應(yīng)用中，常需根據(jù)信號的有限個取樣值來估計AR模型的參數(shù)，應(yīng)用較多的有三種方法：Y-W法或自相關(guān)法；協(xié)方差法；Burg法。這三種方法都可以用由時間平均代替集合平均的最小平方準(zhǔn)則推到得到。l理論上，AR模型參數(shù)是根據(jù)預(yù)測誤差功率最小的準(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ù)測誤差和后向預(yù)測誤差前向預(yù)測誤差和后向預(yù)測誤差.(2) Yule-Walker法法(3) 協(xié)方差法協(xié)方差法(4) Burg法法理論上，理論上，ARAR模型參數(shù)是根據(jù)模型參數(shù)是根據(jù)預(yù)測誤差功率預(yù)測誤差功率最小的準(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ù)- -更好的估計和更低誤差；更好的估計和更低誤差；最小化復(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ù)測誤差前向預(yù)測誤差：后向預(yù)測誤差后向預(yù)測誤差37AR模型參數(shù)提取方法模型參數(shù)提取方法l用最小平方的時間平均準(zhǔn)則代替集合平均準(zhǔn)則，有 或 上式中， 可由式(4)得到l 是由長度為p1的預(yù)測誤差濾波器沖激響應(yīng)序列1,ap,1,ap,p與長度為N的數(shù)據(jù)序列x(0),x(1),x(N-1)進行卷積得到的，所以它的長度為N+p，這決定了上式中的求和項數(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因此，用時間平均最小化準(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法法)計算前向預(yù)測誤差原理計算前向預(yù)測誤差原理: :協(xié)方差法計算前向預(yù)測誤差原理協(xié)方差法計算前向預(yù)測誤差原理: :(0)e( )e n(1 )e N p ( )e p( )e n(1 )e N兩端加零兩端加零兩端不加零兩端不加零兩端加零兩端加零注意點：注意點： 40l

31、41AR模型參數(shù)提取方法模型參數(shù)提取方法l用下列時間平均的最小平方準(zhǔn)則代替集合平均的最小平方準(zhǔn)則： 該式與自相關(guān)法的主要區(qū)別是求和范圍不同?，F(xiàn)在的求和范圍是p(N-1)。這意味著沒有“加數(shù)據(jù)窗”的不合理假設(shè)。l與式(12)類似，可推導(dǎo)出 其中，自相關(guān)矩陣的估計為 這里，自相關(guān)序列的估計為l一般情況下， 不是Toeplitz的， 半正定，這是與自相關(guā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)法的計算效率高，且能保證預(yù)測誤差濾波器是最小相位的，但數(shù)據(jù)兩端要附加零取樣值，實際上是數(shù)據(jù)加窗，這將使參數(shù)估計的精度下降。特別當(dāng)數(shù)據(jù)段很短時，加窗效應(yīng)就更加嚴(yán)重。協(xié)方差法計算效率也高，但潛在著不穩(wěn)定因素。自相關(guān)法和協(xié)方差法都是直接估計AR參數(shù)。lBurg法則一方面希望利用已知數(shù)據(jù)段兩端以外的未知數(shù)據(jù)(但它對這些未知數(shù)據(jù)不作主觀臆測)，另一方面又總是設(shè)法保證使預(yù)測誤差濾波器是最小相位的。Burg法與自相關(guān)法和協(xié)方差法不同，它不直接估計AR參數(shù)，而是先估計反射系數(shù)，然后利用Levinson遞推算法由反射系數(shù)來求得AR參數(shù)。vBurg法法43AR模型參數(shù)提取方法模型

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

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

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

36、(1)kkk47AR譜估計的異常現(xiàn)象及其補救措施譜估計的異?，F(xiàn)象及其補救措施l虛假譜峰 補救措施：模型的階不宜選得太高，最高不應(yīng)超過數(shù)據(jù)記錄長度的一半。l譜線分裂 補救措施：同時調(diào)整所有反射系數(shù)，使預(yù)測誤差功率真正達到最小。l噪聲對AR譜估計的影響 AR譜估計對觀測噪聲比較敏感，噪聲會使譜峰展寬，從而導(dǎo)致分辨率下降，而且會使譜峰偏離正確的位置。 為減小噪聲對AR譜估計的惡化影響，一般可使用下列四種方法:(1)采用ARMA譜估計方法；(2)對數(shù)據(jù)進行濾波，減小噪聲；(3)采用高階AR模型；(4)補償自相關(guān)函數(shù)或反射系數(shù)估計中噪聲的影響。3.3 3.3 最大熵譜估計法最大熵譜估計法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 最大似然估計法最大似然估計法l51Maximum Likelihoodl523.5 互協(xié)方差估計與互譜估計互協(xié)方差估計與互譜估計l533.6 特征分解法譜估計特征分解法譜估計l54l55求解過程求解過程l56二階二階l57l58211(e)jwMH

38、iiipPa e v3.6 Pisarenko 諧波分解諧波分解l59l用Pisarenko 求解噪聲中估計兩個復(fù)指數(shù)分量頻率和功率，假設(shè)兩個復(fù)指數(shù)加噪聲的隨機過程自相關(guān)序列的前三個值為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， 其對應(yīng)的特征矢量l求根得到62min0.27850.30060.32090.74920.40970.0058jvjj /310.50.8660jzje2 /520.30990.9511jzje13w225wlV1,V2 對應(yīng)的功率功率方程：其中 ： 得到 63121()2.9685jwV e221()2.9861jwV e222()0.0139jwV e122()0.0315jwV e2112

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

41、。l對于ARMA模型差分方程 其中，輸入激勵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模型譜估計模型譜估計l由式(1)可得 對上式兩端取逆Z變換，分別得到 這里，假設(shè)h(n)是實序列。由上兩式得到lh(n)是因果序列，即n0時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模型譜估計模型譜估計l將式(3)代入式(2)，得 這就是ARMA模型參數(shù)與自相關(guān)函數(shù)之間的關(guān)系式。l當(dāng)a0=1且ak=0(k=1,2,p)時，由式(4)可得出MA模型參數(shù)與信號的自相關(guān)函數(shù)之間的關(guān)系式。注意此時h(k)=bk,故有2101()( ),0,1,( )(),1pq m