第二章 現(xiàn)代譜估計(jì)現(xiàn)代譜估計(jì)概述AR模型譜估計(jì)線性預(yù)測Burg算法ARMA模型譜估計(jì)擴(kuò)展Prony方法多重信號分類(MUSIC)法12.1 現(xiàn)代譜估計(jì)概述經(jīng)典譜估計(jì)的主要問題基于信號參數(shù)模型的譜估計(jì)方法2經(jīng)典譜估計(jì)：分辨率低（受窗函數(shù)長度的限制）;

方差性能不好；方差和分辨率之間的矛盾。對平穩(wěn)信號建模:

用于功率譜估計(jì):提高分辨率，減小方差；也可用于信號的特征提取，預(yù)測，編碼及數(shù)據(jù)壓縮等。32.1 現(xiàn)代譜估計(jì)概述

在第一章中介紹了用參數(shù)模型來描述隨機(jī)信號的方法，如果能確定信號x(n)的信號模型，根據(jù)信號觀測數(shù)據(jù)求出模型參數(shù)，系統(tǒng)函數(shù)用H(z)表示，模型輸入白噪聲方差為σw2，信號的功率譜用下式求出：按照這種思路，功率譜估計(jì)可分成三個(gè)步驟：（1）選擇合適的信號模型；（2）根據(jù)x(n)有限的觀測數(shù)據(jù)，或者它的有限個(gè)自相關(guān)函數(shù)估計(jì)值，估計(jì)模型的參數(shù);

（3）計(jì)算模型輸出功率譜。4參數(shù)模型與模型譜信號測量數(shù)據(jù)1、基本思路

模型（參數(shù)）（MA,AR,ARMA）確定隨機(jī)w(n)S(n)u(n)

逼近即：信號的當(dāng)前值是由其過去的值和輸入信號現(xiàn)在與過去的值按照模型參數(shù)線性加權(quán)組合得到。52.自回歸－滑動(dòng)平均模型(ARMA(p,q)）

（Pole-zeromodel）X(n)的功率譜6推導(dǎo)：于是：73.滑動(dòng)平均模型—MA(q)（Allzeromodel）84.自回歸模型—（Allpolesmodel）95.模型參數(shù)ModelParameters:10

已知有限長數(shù)據(jù)或有限長的相關(guān)函數(shù)，估計(jì)的步驟為：①建立參數(shù)模型(MA,AR,ARMA,階次p,q)②

由或估計(jì)的參數(shù)(解法)

－－參數(shù)估計(jì)③由的參數(shù)估計(jì)【注】估計(jì)ARMA或MA的參數(shù)一般需解一組非線性方程組，估計(jì)AR模型參數(shù)通常只需解一組線性方程組。6.參數(shù)估計(jì)與譜估計(jì)11例.(白噪聲中的AR過程)12故即：白噪聲中的AR(P)過程是一個(gè)ARMA(p,q)過程，其激勵(lì)噪聲是白的，方差為其中：132.2 AR模型譜估計(jì)AR模型的正則方程Levinson-Durbin算法AR譜估計(jì)的自相關(guān)法AR模型階次的選擇AR模型譜估計(jì)的性質(zhì)14一、AR模型的正則方程AR模型：差分方程：15目標(biāo)：找到已知參數(shù)和未知參數(shù)的關(guān)系， 以便求解未知參數(shù)：已知參數(shù)：求解方法：由下面的差分方程入手：兩邊同乘，求均值未知參數(shù)：?=+--=pkknuknxanx1)()()(16一、AR模型的正則方程假定、都是實(shí)平穩(wěn)的隨機(jī)信號，為白噪聲，方差為，為服從AR過程的因果信號。由AR模型的差分方程，有將上式兩邊同乘以，并求均值，得

17一、AR模型的正則方程

(a)式中，為AR模型的單位取樣響應(yīng)。由z變換的性質(zhì) ，當(dāng)時(shí)，有。將之代入上式，有

(b)18一、AR模型的正則方程綜合式(a)與式(b)，有

在上述推導(dǎo)中，應(yīng)用了實(shí)信號自相關(guān)函數(shù)的偶對稱性，即。由上式可得個(gè)方程，寫成矩陣形式為19一、AR模型的正則方程上述兩式即為AR模型的正則方程，又稱Yule-Walker方程。

或20上式不考慮對稱性表示為如下形式可簡單的表示為其中，

是

階AR模型的系數(shù)向量，

是

維全零列向量，

定義為(6)21AR模型的正則方程也可表示為22式(8)簡單地表示為其中，

因矩陣

是非奇異的，有將

代入式(7)中即得到

。23例

已知

滿足實(shí)

模型，即滿足如下差分方程其中，

是均值為零，方差為

的白噪聲。

模型階數(shù)

，得到二階的Yule-Walker方程

取AR24所以，可以解得同樣可以用

和

來表示

和

，即2526一、AR模型的正則方程需要指出的是，上式中的自相關(guān)矩陣為Toeplitz矩陣；若是復(fù)過程，那么，則其自相關(guān)矩陣是Hermitian對稱的Toeplitz矩陣。這類矩陣具有一系列好的性質(zhì)，利用這些性質(zhì)，可以找到快速求解AR模型參數(shù)的高效算法。

若取估計(jì)

則27二、Levinson-Durbin算法

Levinson-Durbin遞推算法是求解Yule-Walker方程的快速有效算法，這種算法利用了方程組系數(shù)矩陣（自相關(guān)矩陣）所具有的一系列好的性，使運(yùn)算量大大減少。其推導(dǎo)的方法有多種，這里只介紹一種較為簡便的推導(dǎo)方法。

在實(shí)際應(yīng)用中，為避免矩陣求逆運(yùn)算，降低計(jì)算量，通常并不直接求解正則方程，而是根據(jù)自相關(guān)矩陣的Toeplitz性質(zhì)，利用Levinson-Durbin迭代算法進(jìn)行求解。28二、Levinson-Durbin算法定義

為的第

個(gè)系數(shù)；

為

階AR模型輸入白噪聲的方差；

階AR模型思路：利用Toeplitz

矩陣特點(diǎn)，由低階高階反射系數(shù)29二、Levinson-Durbin算法計(jì)算

階AR模型的參數(shù)，由(6)得

對于

，若已知

階AR模型的參數(shù)

和

容易解得，模型的正則方程為30二、Levinson-Durbin算法(a)的參數(shù)，要求解的m階Yule-Walker方程為31二、Levinson-Durbin算法

(b)32二、Levinson-Durbin算法(c)為此，將式(a)的系數(shù)矩陣增加一行和增加一列，成為下式：33二、Levinson-Durbin算法式中利用前述的系數(shù)矩陣的特點(diǎn)，將式(c)的行倒序，同時(shí)列也倒序，得到34二、Levinson-Durbin算法

(d)35二、Levinson-Durbin算法將待求解的m階Yule-Walker方程表示成式(c)和式(d)的線性組合形式，即

(e)36二、Levinson-Durbin算法或式中，是待定系數(shù)，稱為反射系數(shù)。式(e)兩邊各右乘以m階系數(shù)矩陣，得到

(f)37二、Levinson-Durbin算法由式(f)可求出由式(c)的第一個(gè)方程可求出從上面的推導(dǎo)中可歸納出由m-1階模型參數(shù)求m階模型參數(shù)的計(jì)算公式如下：38二、Levinson-Durbin算法對于AR(p)模型，遞推計(jì)算直到p階為止。39三、AR譜估計(jì)的自相關(guān)法已知N點(diǎn)觀測數(shù)據(jù)和AR的階數(shù)p，則AR譜估計(jì)可按下述步驟來進(jìn)行：由已知的估計(jì)令

40三、AR譜估計(jì)的自相關(guān)法用代替L-D遞推算法式中的，對于，重新求解Yule-Walker方程，這時(shí)求出的AR模型參數(shù)是真實(shí)參數(shù)的估計(jì)值，即和將這些參數(shù)代入下式，得到的功率譜的估計(jì)，即41三、AR譜估計(jì)的自相關(guān)法若在（0，2π）內(nèi)對進(jìn)行N點(diǎn)均勻抽樣，則得到離散譜式中，。離散譜，用FFT計(jì)算42Levinson-Durbin算法流程圖43三、AR模型階次的選擇FPE準(zhǔn)則（最終預(yù)測誤差準(zhǔn)則）隨著m的增加，使達(dá)到最小值時(shí)的。AIC準(zhǔn)則（信息論準(zhǔn)則）前者表征將隨著m的增加而單調(diào)下降，后者表示計(jì)算誤差將隨著m的增加而增長。44五、AR模型譜估計(jì)的性質(zhì)AR模型估計(jì)出的功率譜具有一系列好的性質(zhì)，現(xiàn)分別討論如下：平滑特性AR模型是一有理分式，估計(jì)出的譜平滑，不需要像周期圖那樣再做平滑或平均，因此，不需要為此去犧牲分辨率。45五、AR模型譜估計(jì)的性質(zhì)

2AR模型功率譜的分辨率BT法的功率譜

的分辨率是隨著信號長度N

的增加而提高的。而AR模型法避免了窗函數(shù)的影響，因此它可得到高的譜分辨率，同時(shí)它所得出的功率譜估計(jì)與真實(shí)的功率譜偏差很小。AR譜估計(jì)的頻率分辨率，要優(yōu)于經(jīng)典譜估計(jì)方法。其原因在于求解AR模型參數(shù)的過程，實(shí)際上意味著將根據(jù)估計(jì)的按一定準(zhǔn)則進(jìn)行了外推。462.AR譜的分辨率經(jīng)典譜估計(jì)：假定：分辨率反比于N

，即對間接法：分辨率還要降低47AR模型包含了對的“預(yù)測”或“外推”。實(shí)際上，這包含著自相關(guān)函數(shù)的“外推”。令：AR譜AR譜對應(yīng)的自相關(guān)函數(shù)可以證明：AR模型自相關(guān)函數(shù)匹配性質(zhì)48五、AR模型譜估計(jì)的性質(zhì)3AR模型的穩(wěn)定性AR模型穩(wěn)定的充分必要條件是

的極點(diǎn)

必須在單位圓內(nèi)，

而且這一條件也是保證了

是一

個(gè)廣義平穩(wěn)過程。

4AR譜估計(jì)與最大熵譜估計(jì)的等效性對于高斯隨機(jī)信號，最大熵譜估計(jì)和AR模型譜估計(jì)是等效的。49五、AR模型譜估計(jì)的性質(zhì)5AR模型譜估計(jì)方法的不足在實(shí)際應(yīng)用時(shí)，發(fā)現(xiàn)AR模型在譜估計(jì)中存在一些缺點(diǎn)。有些缺點(diǎn)和模型本身有關(guān)，有些則和采用的求解模型參數(shù)的方法有關(guān)。AR譜估計(jì)的分辨率受到加性觀測噪聲的影響。求AR模型時(shí)所使用信號的信噪比越低，AR譜的分辨率越差。如果待估計(jì)信號是含有噪聲的正弦信號，譜峰的位

50置易受初相位的影響，且在有的算法中還可能出現(xiàn)“譜線分裂（spectrallinesplitting）”的現(xiàn)象。通過算法的改進(jìn)和其它一些措施可以較好的克服這些缺點(diǎn)。

譜估計(jì)的質(zhì)量受到階次p的影響。階次太低，譜太平滑，反映不出譜峰。階次選得過大可能會(huì)產(chǎn)生虛假的譜峰。512.3 線性預(yù)測前向線性預(yù)測后向線性預(yù)測格形濾波器52一、前向線性預(yù)測 已知n時(shí)刻以前的m個(gè)信號數(shù)據(jù)，用這m個(gè)數(shù)據(jù)來線性預(yù)測n時(shí)刻信號的值，如圖所示，預(yù)測值為式中，上標(biāo)f表示前向預(yù)測。53一、前向線性預(yù)測其預(yù)測誤差為

(a) ——稱此預(yù)測器為m階前向線性預(yù)測器。

令，由此解得

將式(a)代入上式，得

54一、前向線性預(yù)測

(b)由最小均方誤差的表達(dá)式及正交性原理可知

(c)聯(lián)立式(b)與式(c)得

55一、前向線性預(yù)測

(d)——前向線性預(yù)測的Wiener-Hopf方程

解此方程則得m階線性預(yù)測器的最佳參數(shù)及。56一、前向線性預(yù)測上式與AR模型參數(shù)的正則方程式極其相似，若令，則有，成立。這說明，對于同一個(gè)p階的AR隨機(jī)信號，其AR模型和同階的最佳線性預(yù)測器模型是等價(jià)的。所以有

(f)即p階線性預(yù)測器的輸出是一個(gè)白噪聲序列。57應(yīng)等于AR模型激勵(lì)白噪聲的功率。由LP的含意，因此AR模型也可以看作是在最小平方意義上對數(shù)據(jù)的擬合；上面等效的含意是：由于LP包含了對數(shù)據(jù)的外推，因此，對應(yīng)的譜估計(jì)所用數(shù)據(jù)的范圍比實(shí)際的應(yīng)有擴(kuò)展，因此可以提高分辨率。線性預(yù)測器的誤差序列等效于激勵(lì)A(yù)R模型的白噪聲序列；58一、前向線性預(yù)測結(jié)論：對于給定的隨機(jī)信號，若其最佳前向線性預(yù)測器的階次等于的AR模型階次時(shí)，其前向線性預(yù)測誤差為白噪聲序列。所以階次等于AR模型階次的最佳前向預(yù)測誤差濾波器實(shí)際上是AR模型的逆系統(tǒng)，即白化濾波器。

59AR模型白化濾波器線性預(yù)測器60二、后向線性預(yù)測與前向線性預(yù)測對應(yīng)的還有后向線性預(yù)測器。即由n時(shí)刻以后的p個(gè)數(shù)據(jù)來預(yù)測，即式中，上標(biāo)b代表后向預(yù)測。在實(shí)際工作中，總是用同一組數(shù)據(jù)來同時(shí)實(shí)現(xiàn)前向和后向預(yù)測，這樣上式可改寫為61二、后向線性預(yù)測預(yù)測誤差但習(xí)慣上常將寫成，即仿照前向預(yù)測器的推導(dǎo)方法，同樣可導(dǎo)出下列公式：最佳均方誤差及

62上述的關(guān)系還是集總平均。對實(shí)際的信號：單個(gè)樣本有限長，求均值要簡化，對取代

的范圍？6364N點(diǎn)數(shù)據(jù)，前向預(yù)測誤差序列范圍6566上三角+中間塊+下三角：上、下加窗；中間塊：上、下不加窗；中間塊+上三角：下不加窗、上加窗；中間塊+下三角：上不加窗、下加窗；67二、后向線性預(yù)測對于實(shí)序列有及 若為復(fù)數(shù)序列，則

68三、格形濾波器對于p階的前、后向預(yù)測誤差，有如下遞推公式成立式中，稱為反射系數(shù)，。且691.格形濾波器結(jié)構(gòu)702.格型濾波器的性質(zhì)（1）各階后向預(yù)測誤差相互正交。用公式表示如下：各階后向預(yù)測誤差相互正交的結(jié)果，使濾波器前后級互相解耦，對于系統(tǒng)最小化問題化為一系列獨(dú)立的對每一級局部最小化問題。用作自適應(yīng)濾波時(shí)，各級可選用不同的自適應(yīng)步長，使收斂速度提高。另外，為提高線性預(yù)測性能，需要增加一節(jié)或幾節(jié)，可以只對新增加的級進(jìn)行獨(dú)立的調(diào)節(jié)，達(dá)到輸出均方誤差最小，無需再調(diào)節(jié)前面的系數(shù)。712.格型濾波器的性質(zhì)

(2)平穩(wěn)隨機(jī)序列可由自相關(guān)函數(shù)或反射系數(shù)表征。按照Levinson-Durbin遞推公式，已知rxx(0)，k1，k2，…，kp，從一階開始，可以推出全部的預(yù)測系數(shù)ap,1，ap,2，…，ap,p和σ2p，把得到的這些數(shù)據(jù)代入Yule-walker方程，可求得信號的自相關(guān)函數(shù)rxx(0)，rxx(1)，rxx(2)，…，rxx(p)。以上說明平穩(wěn)隨機(jī)序列可由自相關(guān)函數(shù)表征，也可由rxx(0)，k1，k2，…，kp表征。

(3)前向預(yù)測誤差濾波器是最小相位濾波器，即它的全部零點(diǎn)在單位圓內(nèi)。722.4 Burg算法Burg算法的基本概念Burg算法存在的問題改進(jìn)的協(xié)方差算法73一、Burg算法的基本概念基本思想 自相關(guān)法進(jìn)行AR譜估計(jì)時(shí)，是遵循以下思路進(jìn)行的：由觀測的信號數(shù)據(jù)先估計(jì)自相關(guān)函數(shù)。根據(jù)估計(jì)的自相關(guān)函數(shù)，利用Levinson-Durbin遞推算法求解模型參數(shù)、。由得出的AR模型參數(shù)計(jì)算信號的功率譜。74一、Burg算法的基本概念1967年提出的Burg算法在一定程度上改善了這種狀況。它所遵循的計(jì)算思路是：由觀測的信號數(shù)據(jù)先估計(jì)反射系數(shù)。根據(jù)估計(jì)出的反射系數(shù)，利用Levinson-Durbin算法遞推出AR模型參數(shù)、。由得出的AR模型參數(shù)計(jì)算信號的功率譜。75一、Burg算法的基本概念由于在計(jì)算中避開了估計(jì)自相關(guān)函數(shù)，而直接從輸入數(shù)據(jù)計(jì)算AR模型參數(shù)所以減小了計(jì)算誤差，從而改善了的頻率分辨率。Burg算法的另一特點(diǎn)是：使用前向、后向預(yù)測誤差平方和最小的原則來估計(jì)，而不是象自相關(guān)法那樣只按前向預(yù)測誤差的方差最小的原則導(dǎo)出其正則方程。762.Burg算法的推導(dǎo)算法推導(dǎo)令應(yīng)滿足在上式中代入格形濾波器公式，可得式1772.Burg算法的推導(dǎo)估計(jì)出后，階次m時(shí)的AR模型參數(shù)系數(shù)仍然由Levinson-Durbin算法遞推求出，即有式(3)式(4)式(5)783.Burg算法的計(jì)算步驟計(jì)算步驟由初始條件，再由式1求出；由得時(shí)的參數(shù)由及式2求出和，再由式1

估計(jì)；依照式(3)~式(5)，求出時(shí)的參數(shù)、及；重復(fù)上述過程，直到，求出所有階次時(shí)的AR參數(shù)。793.Burg算法的計(jì)算步驟若定義式(1)的分母為那么可以證明，可以由和遞推計(jì)算，即這樣，可以有效地提高計(jì)算速度。Burg算法的優(yōu)點(diǎn)：頻率分辨率高；所得的AR模型穩(wěn)定；計(jì)算效率高。80Burg算法流程81二、Burg算法存在的問題

Burg算法由于避免了直接計(jì)算自相關(guān)函數(shù)，所以估計(jì)誤差、頻率分辨率要高于自相關(guān)法。但由于它仍然用Levinson-Durbin算法來求AR模型參數(shù)，因此，仍存在譜峰分裂與偏移問題。當(dāng)使用單頻正弦信號（周期為T）加白噪聲作試驗(yàn)數(shù)據(jù)時(shí)，有如下規(guī)律：信噪比高時(shí)容易發(fā)生。由于此時(shí)譜峰突出，因此已不再是AR過程，而是退化的ARMA過程。82二、Burg算法存在的問題當(dāng)信號長度（為采樣周期）時(shí)，對任何整數(shù)k，任何的初始相位，譜峰無分裂。當(dāng)信號長度，的初始相位為0或時(shí)，譜峰無分裂；當(dāng)信號長度，的初始相位為的奇數(shù)倍時(shí)，譜峰分裂嚴(yán)重；N充分大時(shí)，分裂現(xiàn)象逐漸減弱。83三、協(xié)方差法與修正協(xié)方差法

1.協(xié)方差法

這種方法和自相關(guān)法一樣，仍然利用使預(yù)測誤差功率最小的方法求模型參數(shù)，但由觀測數(shù)據(jù)求預(yù)測誤差功率的公式如下式：將該式對比自相關(guān)法中求預(yù)測誤差功率的公式，不同的是求和限不同。該公式中使用的觀測數(shù)據(jù)均已得到，不需要在數(shù)據(jù)兩端補(bǔ)充零點(diǎn)，因此比較自相關(guān)法去掉了加窗處理的不合理假設(shè)。為求得模型參數(shù)仍然應(yīng)用復(fù)梯度法使式達(dá)到最小，公式如下：841、協(xié)方差法式中

白噪聲的方差為

851、協(xié)方差法由觀測數(shù)據(jù)x(n)(n=0,1,2,…,N-1)，利用上面三個(gè)公式可以求出模型的參數(shù)：{api

(i=1,2,3,…,p)；σ2w}。按照定義，式中的cxx(j,k)可以稱作協(xié)方差函數(shù)，它有兩個(gè)變量，因此也適合于非平穩(wěn)隨機(jī)信號。式中的協(xié)方差矩陣是埃爾米特(Hermitian)矩陣，cxx(k,j)=c*xx(j,k)，是半正定的。這種方法近似于自相關(guān)法。一些實(shí)驗(yàn)結(jié)果說明它的分辨率優(yōu)于自相關(guān)法，另外對于純正弦信號數(shù)據(jù)，可以有效地估計(jì)正弦信號的頻率。862、改進(jìn)的協(xié)方差算法基本思想：由觀測的數(shù)據(jù)直接估計(jì)，而不需要估計(jì)，從而無需使用Levinson-Durbin算法。具體說來，其思路是首先定義：然后，令也就是說，在令為最小時(shí)，不是僅令其相對為最小，而是令其相對都為最?。╩由1至p）。

87

3.遞推算法：由求，由遞 推，還是直接由遞推各算法之間的主要區(qū)別：1.的取值范圍，即選擇那一個(gè)？

2.僅用前向預(yù)測，還是前后向都預(yù)測？即 令最小，還是最?。?？88

例

已知信號的4個(gè)觀察數(shù)據(jù)為x(n)={x(0),x(1),x(2),x(3)}={2,4,1,3}，分別用自相關(guān)法和協(xié)方差法估計(jì)AR(1)模型參數(shù)。解:（1）自相關(guān)法：89

(2)協(xié)方差法:902.5 ARMA模型譜估計(jì)噪聲對AR模型譜估計(jì)的影響MA模型譜估計(jì)ARMA模型譜估計(jì)91一、噪聲對AR譜估計(jì)的影響設(shè)是一個(gè)p階AR過程，它被測量噪聲污染后成為。即如果是方差為的白噪聲，且與不相關(guān)，則有92二、MA譜估計(jì)的計(jì)算

(2.5.3) (2.5.4) (2.5.5)由式(2.5.3)得

(2.5.6)

93二、MA譜估計(jì)的計(jì)算將上式兩邊同乘以，并求均值，得

(2.5.7)

式中，。94二、MA譜估計(jì)的計(jì)算將x(n)樣值響應(yīng)代入，并注意到是方差為的白噪聲，有

(2.5.8)對MA(q)模型，由式(2.5.4)

，得95二、MA譜估計(jì)的計(jì)算所以，可以求出MA(q)模型的正則方程，即有MA(q)的功率譜為

等效于經(jīng)典譜估計(jì)中的自相關(guān)法，即MA譜估計(jì)等效為信號長度為q+1的自相關(guān)法譜估計(jì)。96三、ARMA模型譜估計(jì)ARMA(p,q)模型的差分方程式中，。類似地，可導(dǎo)出其正則方程如下：97三、ARMA模型譜估計(jì)式中，是系數(shù)和的函數(shù)，前q+1個(gè)方程是高度非線性的。從第q+1個(gè)方程開始是線性的，可以解出AR部分的系數(shù)，將上式中的第二個(gè)方程寫成如下展開形式：(2.5.13)98三、ARMA模型譜估計(jì)上式雖然可解出AR部分的系數(shù)，但存在以下兩個(gè)問題：由于式中的真實(shí)自相關(guān)函數(shù)是未知的，因此只能使用估計(jì)值來代替，且要用到大延遲的估計(jì)值（最大延遲是p+q），而對于給定的信號長度，這將造成估計(jì)很不準(zhǔn)確。因而，也就不能得到AR部分系數(shù)的準(zhǔn)確估計(jì)。

99三、ARMA模型譜估計(jì)式中階次p和q都是未知的，需要事先指定。而p和q的不正確指定有可能導(dǎo)致式(2.5.13)的系數(shù)矩陣奇異。因此，在實(shí)際應(yīng)用中，對式(2.5.13)采用更一般的形式，即取L個(gè)方程，這里，即式中， 100三、ARMA模型譜估計(jì)由此得到的最小二乘解為

求得ARMA(p,q)模型中的AR參數(shù)，余下的任務(wù)就是求解MA部分的參數(shù)。101三、ARMA模型譜估計(jì)利用求得的AR系數(shù)先得到一個(gè)FIR系統(tǒng)為序列經(jīng)此FIR系統(tǒng)濾波，得到一個(gè)輸出序列ARMA(p,q)模型與FIR系統(tǒng)級聯(lián)，近似于模型。

102三、ARMA模型譜估計(jì)因此，可以利用輸出序列估計(jì)自相關(guān)序列并按MA(q)模型譜估計(jì)公式來得到MA譜，即得到MA譜估計(jì)后，利用下式即可求得ARMA譜估計(jì)1032.6最大熵譜估計(jì)與最大似然譜估計(jì)一、最大熵譜估計(jì)

1.利用最大熵的原則外推自相關(guān)函數(shù)

按照Shannon對熵的定義，

當(dāng)隨機(jī)變量X取離散值時(shí)，熵的定義為

式中，pi是出現(xiàn)狀態(tài)i的概率。當(dāng)X取連續(xù)值時(shí)，熵的定義為

（a）104一、最大熵譜估計(jì)式中，p(x)是X的概率密度函數(shù)，對于離散隨機(jī)序列，概率密度函數(shù)用聯(lián)合概率密度函數(shù)代替。顯然,熵代表一種不確定性，最大熵代表最大的不確定性，或者說最大的隨機(jī)性。下面我們研究對于有限的自相關(guān)函數(shù)值不作任何改變，對于未知自相關(guān)函數(shù)用最大熵原則外推，即不作任何附加條件的外推方法。假設(shè)x(n)是零均值正態(tài)分布的平穩(wěn)隨機(jī)序列，它的N維高斯概率密度函數(shù)為式中

105按照（a）式，x(n)信號的熵為

（b）

式中，det(Rxx(N))表示矩陣Rxx(N)的行列式。上式表明，為使熵最大，要求det(Rxx(N)最大。106若已知N+1個(gè)自相關(guān)函數(shù)值rxx(0)，rxx(1)，…，rxx(N)，下面用最大熵方法外推rxx(N+1)。設(shè)rxx(N+1)確實(shí)是信號自相關(guān)函數(shù)的第N+2個(gè)值，根據(jù)自相關(guān)函數(shù)的性質(zhì)，由N+2個(gè)自相關(guān)函數(shù)組成的矩陣為107它必須是非負(fù)定的矩陣，

即

將行列式展開，det(Rxx(N+1))是rxx(N+1)的二次函數(shù)，該二次函數(shù)系數(shù)的符號是：(-1)1+N+2(-1)1+N+1=-1,且det(Rxx(N+1))對rxx(N+1)的二次導(dǎo)數(shù)是-2det［Rxx(N-1)］，它是負(fù)值，負(fù)值表示det[Rxx(N+1)]對rxx(N+1)的一次導(dǎo)數(shù)是減函數(shù)，det[Rxx(N+1)]作為rxx(N+1)的函數(shù)，凹口向下，那么只有一個(gè)最大值。為選擇rxx(N+1)使det[Rxx(N+1)]最大，解下列方程：108用數(shù)學(xué)歸納法，得到

上式是rxx(N+1)的一次函數(shù)，可以解出rxx(N+1)。繼續(xù)再將rxx(N+1)代入Rxx(N+2)和det[Rxx(N+2)]中，求det[Rxx(N+2)]對rxx(N+2)的最大值，得到rxx(N+2)；以此類推，可推出任意多個(gè)其它自相關(guān)函數(shù)值，而不必假設(shè)它們?yōu)榱悖@就是最大熵譜估計(jì)的基本思想。(2.6.7)1092.最大熵譜估計(jì)與AR模型譜估計(jì)的等價(jià)性

我們已經(jīng)知道AR模型信號自相關(guān)函數(shù)與模型參數(shù)服從Yule-Walker方程，即m≥1m=0將m≥1的情況寫成矩陣形式：

（2.6.8）110式中，ai=hA(i)，i=1,2,…,N，ai是AR模型系數(shù)。解該方程，可以得到模型系數(shù)ai，即111112在（2.6.8）式中，令m=N+1，得到

將以上求出的系數(shù)a1,a2,…,aN代入上式，求出rxx(N+1)。而最大熵外推自相關(guān)函數(shù)的公式是(2.6.7)式，按照該公式的最后一行展開，得到（2.6.12）113（2.6.13）上式即是最大熵外推自相關(guān)函數(shù)的公式，對比（2.6.12）式，兩公式完全一樣，證明了AR模型功率譜估計(jì)和最大熵譜估計(jì)的等價(jià)性。這里最大熵外推自相關(guān)函數(shù)等價(jià)于已知N+1個(gè)自相關(guān)函數(shù)，匹配一個(gè)N階AR信號模型的系數(shù)。一旦通過解Yule-Walker方程，解出模型參數(shù)，最大熵譜估計(jì)用下式計(jì)算信號功率譜：（4.7.14）114

思想：最大似然譜估計(jì)是用一個(gè)FIR濾波器實(shí)現(xiàn)，該濾波器對所關(guān)心頻率的正弦信號，可以無失真地通過，而對于其它頻率的信號，讓其頻響盡可能地小，亦即將它們盡可能地濾除。此時(shí)，濾波器輸出的均方值，就作為信號的功率譜估計(jì)。

二、最大似然譜估計(jì)——最小方差譜估計(jì)115MVDR濾波器原理

圖M抽頭的FIR濾波器FIR濾波器的輸入為隨機(jī)過程

，輸出為

定義輸入信號向量和權(quán)向量分別為

116則輸出可表示為輸出的平均功率可以表示為

117其中，

假設(shè)濾波器輸入信號

為

其中，

是加性白噪聲，

為復(fù)幅度。

118設(shè)感興趣的期望信號是角頻率為

的復(fù)正弦信號，

則選擇濾波器權(quán)向量

的原則是

，使復(fù)正弦信號

無失真地通過濾波器，而盡量抑制其余頻率的信號和噪聲。

設(shè)信號

通過濾波器的響應(yīng)為

，則

119定義向量則所以，當(dāng)滿足

時(shí)，

無失真地通過

濾波器。同時(shí)考慮到要使其它復(fù)正弦信號和噪聲盡量被抑制，濾波器權(quán)向量應(yīng)滿足

(1)

約束

，這是為了使

無失真地通過濾波器；

120(2)

輸出平均功率

最小，達(dá)到抑制其它頻率信號和噪聲的目的。考慮更一般的情況，設(shè)期望無失真通過系統(tǒng)的信號頻率為

，且令

，此時(shí)，濾波器權(quán)向量

應(yīng)滿足

應(yīng)用拉格朗日乘子法，構(gòu)造代價(jià)函數(shù)為

121令梯度

，得

考慮到相關(guān)矩陣

是非奇異的，所以

將上式代入到

中

，得

122最優(yōu)權(quán)向量為

濾波器的最小輸出功率為

在工程實(shí)際中，常用

自相關(guān)矩陣的估計(jì)

替換

則最優(yōu)權(quán)向量的估計(jì)為

123最大似然譜估計(jì)稱為最小方差譜估計(jì)更為合適，但由于習(xí)慣也可以仍稱為最大似然譜估計(jì)。而MVDR譜估計(jì)為應(yīng)該指出，此時(shí)并不是真正意義上的信號功率譜，只是描述了信號功率譜的相對強(qiáng)度。

124算法（MVDR信號頻率估計(jì)算法）

步驟1

由

個(gè)觀測樣本

相關(guān)矩陣

估計(jì)樣本步驟2

在

內(nèi)改變

，畫出

，峰值位置

就是信號角頻率的估計(jì)值。

；

1252.7 擴(kuò)展Prony方法擴(kuò)展Prony方法對于N點(diǎn)復(fù)數(shù)序列采用的估計(jì)模型是一組p個(gè)具有任意幅值、相位、頻率與衰減因子的指數(shù)函數(shù)，即式中，為幅值；為相位（單位：弧度）；為衰減因子；表示振蕩頻率；代表采樣間隔。126為敘述方便起見，先將上式寫成如下形式

(a)式中，，。并定義特征多項(xiàng)式

(b)式中，。顯然，式(a)中的即為此多項(xiàng)式的根。設(shè)實(shí)信號用x(n)表示，F(xiàn)IR濾波器系統(tǒng)函數(shù)用A(z)表示：

1272.7 擴(kuò)展Prony方法由（a）式可得兩邊同乘，并求和，得上式等于零是因?yàn)榈诙€(gè)求和項(xiàng)恰好是式(b)位于根處的特征多項(xiàng)式。上式意味著，滿足遞推差分方程式1282.7 擴(kuò)展Prony方法令，則該式表明：此過程是一個(gè)特殊的ARMA(p,p)過程。1292.7 擴(kuò)展Prony方法不妨令可得這就將參數(shù)估計(jì)問題轉(zhuǎn)化為AR模型參數(shù)估計(jì)問題了。當(dāng)找到之后，利用特征方程式求出。然后，利用式(b)可構(gòu)成以下方程：1302.7 擴(kuò)展Prony方法式中，； ； 。這里Φ是一維Vandermonde矩陣。1312.6 擴(kuò)展Prony方法使為最小的解為從而，正弦信號的振幅、相位、衰減因子與頻率可由以下算法得出：振幅相位衰減因子頻率1322.8 多重信號分類(MUSIC)法相關(guān)矩陣的特征分解基于信號子空間的頻率估計(jì)基于噪聲子空間的頻率估計(jì)改進(jìn)的MUSIC法133一 相關(guān)矩陣的特征分解基本思想：將自相關(guān)矩陣中的信息空間分解成兩個(gè)子空間，即信號子空間和噪聲子空間。設(shè)式中

134一 相關(guān)矩陣的特征分解可得自相關(guān)函數(shù)：式中，是第i個(gè)復(fù)正弦的功率。

若定義矢量

是由復(fù)正弦波的N個(gè)取樣值構(gòu)成的矢量，則可寫為式中稱為信號矢量

135一 相關(guān)矩陣的特征分解令是由的N個(gè)取樣數(shù)據(jù)構(gòu)成的矢量，即是由白噪聲的N個(gè)取樣值構(gòu)成的矢量，即

由前述式可得的自相關(guān)矩陣為

136一 相關(guān)矩陣的特征分解若再定義信號自相關(guān)矩陣及噪聲自相關(guān)矩陣如下：則、都是N階方陣，其秩分別為M和N。而可寫為

若，顯然，是奇異的，但由于的存在而正定

137一 相關(guān)矩陣的特征分解現(xiàn)對方陣做特征分解，得：式中，為特征向量所對應(yīng)的特征值，即并且且特征向量之間是正交的，即

138一 相關(guān)矩陣的特征分解由于的秩為M，故其特征值中必有個(gè)0。因此，可寫成如下形式稱為主特征向量。它所張的空間稱為信號子空間，其維數(shù)為M；而由所張的空間稱為噪聲子空間，其維數(shù)為。顯然，這兩個(gè)子空間互為正交補(bǔ)空間。

139一 相關(guān)矩陣的特征分解對于單位陣I，也可表示為特征向量的外積，即因此，可表示為

140二基于信號子空間的頻率估計(jì)若舍去特征向量，僅保留信號子空間，那么我們將用秩為M的相關(guān)陣來近似，這樣可大大提高信號的信噪比?；诰仃?，再用前面所介紹的方法來估計(jì)的功率譜，將得到好的頻率估計(jì)和功率譜估計(jì)。

141三、基于噪聲子空間的頻率估計(jì)對于的非零特征值所對應(yīng)的特征矢量，有成立。從而可得因?yàn)闃?biāo)量，所以有

142三、基于噪聲子空間的頻率估計(jì)即信號子空間的基底可以表示為另一組正交基的線性組合。所以與