【導(dǎo)讀】礎(chǔ)上的功率譜估計(jì)。在實(shí)際應(yīng)用中往往不能獲得具體信號(hào)的表達(dá)式,需要根據(jù)有。的各種方法，且對(duì)每種方法的估計(jì)質(zhì)量做了數(shù)學(xué)推導(dǎo)，并給出仿真程序及仿真圖。論了AR模型參數(shù)法[2]。同時(shí)論文將通過(guò)對(duì)經(jīng)典譜估計(jì)和現(xiàn)代譜估計(jì)的實(shí)現(xiàn)方法。果，為此應(yīng)用更為廣泛[3]。

　　

【正文】 從相關(guān)函數(shù)估計(jì)值經(jīng)傅立葉變換后得到功率譜的估計(jì) 值。從原理上說(shuō)這種方法是直截了當(dāng)?shù)?，可是由此得出的譜估計(jì)值不僅有偏，而且方差很大，它還不隨數(shù)據(jù)長(zhǎng)度的增大而減小到零。因此必須采取其他措施改善它的估計(jì)性能，才能使得這種方法有實(shí)用價(jià)值。常用的兩種辦法是相關(guān)函數(shù)加窗和分段平均周期圖法，這兩種辦法都有實(shí)用價(jià)值，它們都采用了 FFT 算法，使得計(jì)算量大大減小。 歸結(jié)起來(lái)，自相關(guān)法和周期圖法的主要優(yōu)點(diǎn)是； （ 1） 計(jì)算量小 ： （ 2） 功率譜估計(jì)值正比于正弦波信號(hào)的功率。 （ 3）對(duì)于長(zhǎng)記錄數(shù)據(jù)，是一種良好的實(shí)用模型。 26 它們的主要缺點(diǎn)有： 1）頻率分辨率（區(qū)分兩個(gè)鄰近頻率分量的能力）不高。 這是因?yàn)樗鼈兊念l率分辨率約為數(shù)據(jù)長(zhǎng)度的倒數(shù)，且與數(shù)據(jù)的特征或其信噪比無(wú)關(guān)，而實(shí)際應(yīng)用中一般不可能獲得很長(zhǎng)的數(shù)據(jù)記錄； 2）經(jīng)典譜估計(jì)方法在工程中都是以離散傅立葉變換為基礎(chǔ)的，它隱含著對(duì)無(wú)限長(zhǎng)數(shù)據(jù)序列進(jìn)行加窗處理（加了一個(gè)有限寬的矩形窗）。矩形窗的頻譜主瓣不是無(wú)限窄的，且有旁瓣存在，這將導(dǎo)致能量向旁瓣中“泄漏”，主瓣變得模糊不清。嚴(yán)重時(shí)，會(huì)使主瓣產(chǎn)生很大失真，甚至主瓣中的弱分量被旁瓣中的強(qiáng)泄漏所淹蓋。 3）弱信號(hào)被強(qiáng)信號(hào)的旁瓣淹沒(méi)； 4）需采用某種平滑或平均措施以改善譜估計(jì)的統(tǒng)計(jì)特性； 5）某些加窗的相關(guān)函數(shù) 會(huì)使功率譜估計(jì)值出現(xiàn)負(fù)值 [22]。 27 4 現(xiàn)代譜估計(jì) 平穩(wěn)隨機(jī)信號(hào)的參數(shù)模型 由上一章討論可知，經(jīng)典功率譜估計(jì)方法的方差性較差，分辨率較低。方差性能差的原因是無(wú)法實(shí)現(xiàn)功率譜密度原始定義中的求均值和求極限的運(yùn)算。分辨率低的原因，對(duì)周期圖法是假定了數(shù)據(jù)窗以外的數(shù)據(jù)全為零，對(duì)自相關(guān)法是假定了在延遲窗以外的自相關(guān)函數(shù)全為零。當(dāng)然，這種假定是不符合實(shí)際的，正是由于這些不符合實(shí)際的假定產(chǎn)生了經(jīng)典譜估計(jì)較差的分辨率。 在第一章已經(jīng)簡(jiǎn)潔的介紹了現(xiàn)代譜估計(jì)的基本方法，這些方法技 術(shù)的目標(biāo)在于努力改善譜估計(jì)的分辨率。參數(shù)模型法是現(xiàn)代譜估計(jì)的主要內(nèi)容，也是本章討論的重點(diǎn)，參數(shù)模型法的思路如下： （ 1）假定所研究的過(guò)程 )(nx 是由一個(gè)輸入序列 ??nu 激勵(lì)一個(gè)線性系統(tǒng) ??zH的輸出。 （ 2）由已知的 )(nx ，或其自相關(guān)函數(shù) ??mrx 來(lái)估計(jì) ??zH 的參數(shù)。 （ 3）由 ??zH 的參數(shù)來(lái)估計(jì) )(nx 的功率譜。 ??Hz是一個(gè)因果的線性移不變離散時(shí)間系統(tǒng)，當(dāng)然，它應(yīng)該是穩(wěn)定的，其單位抽樣響應(yīng) ??hn是確定的。輸出序列 ??xn可以是平穩(wěn)的隨機(jī)序列，也可以是確定性的時(shí)間序列。若 ??xn是確定性的，那么 ??un是一個(gè)沖激序列，若 ??xn 是隨機(jī)序列，那么 ??un應(yīng)是一個(gè)白噪聲序列。 工程實(shí)際中所遇到的功率譜大體分為三種，一種是“平滑”，即白噪聲的譜，另一種是“線譜”，這是由一個(gè)或多個(gè)純正弦信號(hào)所組成的信號(hào)的功率譜，這兩種是極端的情況；介于二者之間的是既有峰值又有谷值，這種譜稱為 ARMA 譜。顯然，由于 ARMA 模型是一個(gè)零極點(diǎn)模型，它易于反映功率譜中的峰值和谷值。不難想象， AR 模型易于反映譜中的峰值，而 MA 模型易于反映譜中的谷值。 AR,MA 和 ARMA 是功率譜估計(jì)中最主要的參數(shù)模型。本章將會(huì)詳細(xì)的討論AR 模型參數(shù)的計(jì)算、譜的性能及其他算法（如線性預(yù)測(cè)、最大熵譜估計(jì)等）的關(guān)系，最后簡(jiǎn)要給出 MA 模型及 ARMA 模型譜估計(jì)算法 [23]。 28 AR 模型的正則方程與參數(shù)計(jì)算 正則方程的求導(dǎo) 參數(shù)模型法功率譜估計(jì)的主要思想是：將廣義平穩(wěn)的過(guò)程 ??xn表示成一個(gè)輸入序列 ??un激勵(lì)線性系統(tǒng) ??Hz的輸出；由已知的 ??xn或其自相關(guān)函數(shù) ??mrx來(lái)估計(jì) ??Hz的參數(shù)；由 ??Hz的參數(shù)估計(jì) ??xn的功率譜。 AR 模型又稱為自回歸模型，它是一個(gè)全極點(diǎn)模型，其當(dāng)前輸出是現(xiàn)在輸入和過(guò)去輸入的加權(quán)和，表示如下 (其中 ??un為白噪聲序列； p 為 AR 模型的階數(shù) ): ? ? ? ? ? ?nuknxanX pk k ???? ??1 (41) ? ? ? ?kpk kzazAzH??????11112 (42) 由隨機(jī)信號(hào)通過(guò)線性系統(tǒng)理論知輸出序列的功率譜 ? ? 2121 ????? pkkjkjxeaep?? ? (43) 其中 2? 為白噪聲序列的方差，因此進(jìn)行功率譜估計(jì)，必需求得 AR 模型的參數(shù) ka (k=l， 2? p)及 2? 。 假定 ??un、 ??un都是平穩(wěn)的隨機(jī)信號(hào)， ??un為白噪聲，方差為 2? ，現(xiàn)在，我們希望建立 AR 模型的參數(shù) ka 和 ??xn的自相關(guān)函數(shù)的關(guān)系，也即 AR 模型的正則方程 [24]。 將方程（ 41）兩邊同乘以 ? ?mnx ? ，并求得均值，最后得到 ? ? ? ? ? ?mrkmramr xuxpk kx ???? ?? 1 (44) 又因?yàn)? ? ? ? ?? ? ? ???? ???? 00 00 2 mh mmnxnuE ? (45) 29 由 Z 變換的定義， ? ? ? ?lim 0z H z h?? ?，在（ 42）式中，當(dāng) ??z 時(shí)，有 ?? 10?h ，綜合（ 44）和（ 45）兩式，有 ? ?? ?? ?? ? ? ???????????????????pkxkxpkkxmkramkmramr12101? (46) 在上面的推導(dǎo)中，應(yīng)用了自相關(guān)函數(shù)的偶對(duì)稱性。上式寫成矩陣形式，即 ? ? ? ? ? ? ? ?? ? ? ? ? ? ? ?? ? ? ? ? ? ? ?? ? ? ? ? ? ? ?2120 1 2 11 0 1 1 02 0 0 2 01 2 0 0 0x x x xx x x xx x x xx x x xr r r r pr r r r p ar r r r p ar p r p r p r??? ?????? ????????????????????????????? ?? ???? 上述兩式即是 AR模型的正則方程，又稱 YuleWalker 方程 [25]。 AR 模型參數(shù)求解的典型算法 用線性方程組的常用解法 (例如高斯消元法 )解 Yule— Walker 方程，需要的運(yùn)算量數(shù)量級(jí)為 3p ，但若利用系數(shù)矩陣的對(duì)稱性和 Toeplitz 性質(zhì)，則可構(gòu)成一些高效算法， Levinson— Durbin 算法是其中最著名、應(yīng)用最廣泛的一種，這種算法的運(yùn)算量數(shù)量級(jí)為 2p 。這是一種按階次進(jìn)行遞推的算法，即首先以 ? ?0AR 和 ??1AR模型參數(shù)作為初始條件，計(jì)算 ? ?2AR 模型參數(shù)；然后根據(jù)這些參數(shù)計(jì)算 ??3AR 模型參數(shù)等，一直到計(jì)算出 ? ?ARp 模型參數(shù)為止，當(dāng)整個(gè)迭代計(jì)算結(jié)束后，不僅求得了所需要的 P 階 AR 模型參數(shù)，而且還得到了所有各低階模型的參數(shù)。 根據(jù)線性預(yù)測(cè)理論知：一個(gè) P 階 AR 模型的 1p? 個(gè)參數(shù)同樣可用來(lái)構(gòu)成 P階的最佳線性預(yù)測(cè) 器，其預(yù)測(cè)的最小均方誤差等于 AR 模型激勵(lì)白噪聲的能量，即 AR 模型是在最小方差意義上對(duì)數(shù)據(jù) 的擬合。 “ 前向預(yù)測(cè)”是利用 n 之前的 P 個(gè)值對(duì) ??xn性預(yù)測(cè)，如公式 (48)、 (49)、(410)所示；與之對(duì)應(yīng)的“后向預(yù)測(cè)”公式為 (411)、 (412)、 (413)，其中 ??en為為預(yù)測(cè)誤差， P 預(yù)測(cè)誤差功率， f 表示前向預(yù)測(cè)， b 表示后向預(yù)測(cè)。模型參數(shù)算法就是基于 上述最小均方誤差時(shí)由模型參數(shù)估計(jì)信號(hào)功率的方法，主要有以下幾種經(jīng)典算法： 30 ? ? ? ? ? ?knxkanx pk f ???? ??1 ? ? ? ? ? ?fe n x n x n??? ? ?? ?2neEp ff ? ? ? ? ? ? ?knxkanx pk b ??? ??? 1 ? ? ? ? ? ?pnpnxne xb ???? ? ? ?? ?2neEp fbb ? 自相關(guān)法 (BT 法 )。用自相關(guān)法進(jìn)行功率譜估計(jì)，但估計(jì)時(shí)令前向預(yù)測(cè)誤差功率最小，即對(duì) ??nef 前后都加窗構(gòu)成， Wiener— Hopf 方程系數(shù)為 Toepli tz 矩陣，使用 Levinson— Durbin 算法可方便快速的求解 AR 系數(shù)。因此自相關(guān)法也是已知所有 AR 系數(shù)求解方法中簡(jiǎn)單的一種，但譜分辨率相對(duì)較差 [26]。 Burg算法。用 Burg算法進(jìn)行功率譜估計(jì)時(shí)令前后向預(yù)測(cè)誤差功率之和最小，即對(duì) ??neb 、 ??nef 前后都不加窗，使用 Levinson— Durbin 遞推可快速的求解 AR系數(shù)。 Burg 算法是建立在數(shù)據(jù)基礎(chǔ)之上的，避免了先計(jì)算自相關(guān)函數(shù)從而提高計(jì) 算速度；是較知為通用的方法，計(jì)算不太復(fù)雜，且分辨率優(yōu)于自相關(guān)個(gè)法，但對(duì)于白噪聲加正弦信號(hào)有時(shí)會(huì)出現(xiàn)譜線分裂現(xiàn)象。 計(jì)算步驟如下：①由初始條件 ? ? ? ? ? ?nxene nbb ?? 00 ，再由式? ? ? ?? ?? ? ? ?pmnenenenek NmnNmnbmfmNmnbmfm,.. .,2,11121 1 21211111 ??????? ?????????????? (47) 求出 ?1k ； ②由 ? ? ? ? 21010 ???? ? Nnx nxNr 得 1?m 時(shí)的參數(shù) 11,1 ?? ?ka , ? ?01 211??? ???????? ??xrk? ③求出 ??nef1 和 ??neb1 ，再 ?mk 由估計(jì) ?2k ； ④依照 Levinson 遞推關(guān)系，求 2?m 時(shí)的參數(shù) , 2,21,2 ?? aa 及 2?? ； ⑤重復(fù)上述過(guò)程，直到 pm? ，求出所有階次時(shí)的 AR 參數(shù)。 31 (1)改進(jìn)協(xié)方差算法。同 Burg 算法一樣，改進(jìn)協(xié)方差算法進(jìn)行功率譜估計(jì)時(shí)令前后向預(yù)測(cè)誤差功率之和最小，即對(duì) ??neb 、 ??nef 前后都不加窗，但得到的協(xié)方差矩陣不是 Toeplitz 矩陣，因此正則方程不能用 Levinson 遞推算法求解。Marple 于 1980 年提出實(shí)現(xiàn)協(xié)方差方程求解的快速算法，大大提高了譜估計(jì)的性能 [27]。 AR 模型譜估計(jì)的實(shí)現(xiàn)及性質(zhì) 譜估計(jì)的步驟 我們可以按如下步驟估計(jì) ??nx 的功率譜： ① 由 ??nxN 估計(jì) ??nx 的自相關(guān)函數(shù)，得 ?xr ， m=0,1,....p； ②用 ?xr 代替上述地推算法中的 ??mrx ，重新求解 YuleWalker 方程，這時(shí)求出的 AR 模型參數(shù)是真實(shí)參數(shù)的估計(jì)值，即 ???1a 、 ???2a … ???pa 、 ?p? ； 將這些參數(shù)代入（ 3）式，得到 ??nx 的功率譜 ? ??jxeP 的估計(jì)。 AR 模型譜估計(jì)的性質(zhì) 譜的平滑特性 AR 模型是一有理分式，估計(jì)出的譜平滑，不需要像周期圖那樣再做平滑或平均，因此，不需要為此去犧牲分辨率。 譜的分辨率 AR 譜估計(jì)的頻率分辨率，要優(yōu)于經(jīng)典譜估計(jì)方法。其原因在于求解 AR 模型 參 數(shù) 的 過(guò) 程 ， 實(shí) 際 上 意 味 著 將 根 據(jù) ? ? ? ? ? ?1,...,1,0 ?Nxxx 估計(jì)的? ? ?? ? ?Mrrr xxx ,...,1,0 按一定準(zhǔn)則進(jìn)行了外推。 模型譜的匹配性質(zhì) 由維納 辛欽定理及 AR 譜對(duì)應(yīng)的一個(gè)無(wú)限長(zhǎng)的自相關(guān)函數(shù)可得 ? ? ? ?? ?2???? j jjx eA ePeP ? ? ? ? ?2m in?? ? jjAR eAeP ? 若用 AR譜去匹配信號(hào)的譜，則誤差系列的譜應(yīng)由常數(shù)譜來(lái)匹配，體現(xiàn) ??zA 的 32 白化性質(zhì)， 從整體上看 ? ?jARPe?將均勻地和 ? ?jxPe?相跟隨。當(dāng)均值為 1 時(shí)， ? ?jARPe?將在 ? ?jxPe? 的上下波動(dòng)，即在有的區(qū)域， ? ? ? ?jjAR xP e P e???，而在另外的區(qū)域 ，? ? ? ??? jxjAR ePeP ? 。因?yàn)?? ?? ? 121 ???? ?? ?? ?? deP eP jARjx ，整個(gè)積分值為 1，所以又可得到譜匹配的局部性質(zhì)，也就是說(shuō) ? ??jAReP 是 ? ??jxeP 的包絡(luò)的一個(gè)好的近似。 譜的統(tǒng)計(jì)性質(zhì) 嚴(yán)格的分析 AR 譜的方法比較困難，目前尚未有一個(gè)解析表達(dá)式。粗略的講，AR譜的方差反比于數(shù)據(jù) )(nxN 的長(zhǎng)度 N和信噪比 SNR[28]。 AR 模型階次 p 的選擇 準(zhǔn)則（最終預(yù)測(cè)誤差準(zhǔn)則） ? ? ? ?? ?112 ?? ??? mN mNmFPE m? 隨著 m的增加，使 ? ?mFPE 達(dá)到最小值時(shí)的 pm? 。 2. AIC 準(zhǔn)則（信息論準(zhǔn)則