第四章功率譜估計4.1引言4.2經典譜估計4.3現代譜估計中的參數建模4.4

AR模型譜估計的性質4.5

AR譜估計的方法4.6最大熵譜估計與最大似然譜估計4.7特征分解法譜估計NO.24.1引言對信號和系統進行分析研究、處理有兩類方法：一類是在時域進行，我們前面學習的維納濾波、卡爾曼濾波和自適應濾波都屬于這種方法；本章則是在頻率域進行研究的一類方法。這兩類方法都是信號處理的重要方法。對確定性信號傅里葉變換是在頻率域分析研究的理論基礎，但對于隨機信號，其傅里葉變換不存在，因此轉向研究它的功率譜。譜，就是信號的某些特征在頻域隨頻率的分布。功率譜反映了隨機信號功率的分布特性，有著很廣泛的應用：在雷達信號處理中，回波信號的頻譜提供了運動目標的位置、強度、速度等信息；在聲納系統中，為了尋找艦艇或潛艇也要對混有噪聲的信號進行譜分析；在語音處理中，譜分析用來探測語音語調共振峰；在電子戰中，還利用頻譜對目標進行分類。按照Weiner-Khintchine定理，信號的功率譜和其自相關函數服從傅里葉變換關系：實際所能得到的隨機信號的長度總是有限的，用有限長度的信號計算得到的功率譜只是真實功率譜的一個估計。功率譜估計分為經典譜估計和現代譜估計兩大類。上式稱做功率譜的定義，對于平穩隨機信號，服從各態歷經定理，集合平均可以用時間平均代替，由上式還可以推出功率譜的另一個定義：將計算自相關函數中的集合平均用時間平均代替：代入功率譜定義式，得令l=n+m,則

上式中x(n)是觀測數據，Pxx(ejω)是隨機變量，必須對Pxx(ejω)取統計平均值，得到該式被認為是功率譜的另一定義，周期圖法譜估計

Weiner-khintchine定理表明功率譜是無限多個自相關函數的函數，但觀測數據只有有限個,只能得到有限個自相關函數。按照上面的定義公式求功率譜,也需要無限多個觀測數據。因此根據有限個樣本數據，計算隨機信號的功率譜，是一個功率譜的估計問題現代譜估計是以信號模型為基礎，下圖表示x(n)的信號模型，輸入白噪聲w(n)的均值為0，方差為σ2w，x(n)的功率譜為：如果由觀測數據能夠估計出信號模型的參數，則信號的功率譜可以按照上式計算出來，這樣,估計功率譜的問題就變成了估計信號模型參數的問題。信號模型有很多種，如AR模型、MA模型等等，針對不同的情況，需要選擇不同的模型。現代譜估計的質量比經典譜估計的質量有很大的提高。但遺憾的是，尚無任何理論能指導選擇一個合適的模型，只能根據功率譜的一些先驗知識，或者說一些重要的譜特性，來選擇模型。4.2經典譜估計BT法是先估計自相關函數,然后進行傅里葉變換得到功率譜。設對隨機信號x(n)，只觀測到一段樣本數據，n=0,1,2,…,N-1。根據這一段樣本數據估計自相關函數,有兩種估計方法，即有偏自相關函數估計和無偏自相關函數估計。有偏自相關函數估計的誤差相對較小，這種估計是一種漸近一致估計：4.2.1

BT法對上式進行傅里葉變換，得到BT法的功率估計:為了減少譜估計的方差，經常用窗函數w(m)對自相關函數進行加權，此時譜估計公式為式中-(M-1)≤m≤(M-1)其它為了采用FFT計算傅里葉變換，必須將求和域(-M+1,M-1)移到(0~L-1)，功率譜的計算公式為：上式也被稱為加權協方差譜估計。它要求加窗后的功率譜仍是非負的，這樣窗函數w(m)的選擇必須滿足一個原則，即它的傅里葉變換必須是非負的，例如巴特利特窗就滿足這一條件。

k=0,1,2,…,L-10≤m≤M-1M≤m≤L-ML-M+1≤m≤L-1按照有偏自相關函數公式估計自相關函數，已經證明這是漸近一致估計，但經過傅里葉變換后得到功率譜的估計，功率譜估計卻不一定仍是漸近一致估計，可以證明它是非一致估計，是一種不好的估計方法。BT法中用有偏自相關函數進行估計時，它和用周期圖法估計功率譜是等價的，因此BT法的估計質量和周期圖法的估計質量是一樣的。將功率譜的另一定義式重寫如下:如果忽略上式中求統計平均的運算，假設觀測數據為:x(n)

0≤n≤N-1，便得到周期圖法的定義:4.2.2周期圖法用周期圖法計算功率譜框圖由周期圖法功率譜估計公式推導它與BT法的等價關系令

m=k-n,即k=m+n1.周期圖與BT法的等價關系方括號中的部分是有偏自相關函數的計算公式：利用有偏自相關函數的BT法和周期圖法是等價的關系已知自相關函數的估計值，m=-(N-1),…-1,0,1,…,N-1，按照BT法求功率譜的統計平均值：有偏自相關函數的統計平均值在第一章中已確定，將結果代入上式,得到2.周期圖法譜估計質量分析1)周期圖的偏移式中兩序列乘積的傅里葉變換，其頻域服從卷積關系WB(ejω)稱為三角窗的譜函數。上式表明，周期圖的統計平均值等于它的真值與三角窗函數頻譜的卷積，因此周期圖是有偏估計，但當N→∞時，wB(m)→1，三角窗函數的頻譜趨近于δ函數，周期圖的統計平均值趨于它的真值，因此周期圖屬于漸近無偏估計。周期圖的方差的精確表達式很繁冗，為分析簡單：假設x(n)是實的零均值的正態白噪聲信號，方差是σx2，即功率譜是常數σx2，其周期圖用IN(ω)表示，N表示觀測數據的長度。按照周期圖的定義:推導得周期圖的方差：2)周期圖的方差當N趨于無限大時，周期圖的方差并不趨于0，而趨于功率譜真值的平方，即所以無論怎樣選擇N，周期圖的方差總是和σ4x同一個數量級。信號功率譜的真值是σ2x，這說明周期圖的方差很大。用這種方法估計的功率譜在σ2x附近起伏很大，所以周期圖是非一致估計，是一種很差的功率譜估計方法。圖4.2.2白噪聲的周期圖可以看到，隨著N的增大，功率譜曲線并沒有趨向真值，而是起伏變的越來越劇烈。周期圖法估計功率譜不是一致估計，均方誤差很大，使估計出的功率譜不可靠。其頻率分辨力低是根本缺點，原因在于觀測數據只有一段。由于BT法和周期圖法具有等效的功率譜估計質量，因此BT法也不是一致估計，分辨率低，估計誤差大。4.2.3經典譜估計方法的改進基本思想：對隨機變量進行觀測，得到L組獨立的記錄數據，用每一組數據求其周期圖，然后將L個周期圖加起來求平均。這樣得到的周期圖，其方差將是用一組數據得到的周期圖的方差的1/L。1.平均周期圖法假設隨機信號x(n)的觀測區間為：0≤n≤M-1，共進行了L次獨立觀測，得到L組記錄數據，某一組記錄數據用xi(n),i=1,2,…,L表示，第i組的周期圖為：將得到的L個周期圖進行平均，作為信號x(n)的功率譜估計：為了分析偏移，對上式求統計平均:其中，窗函數的頻譜為表明平均周期圖仍然是有偏估計，偏移和每一段的數據個數M有關；由于M<N,平均周期圖的偏移比周期圖的偏移大，表現在三角譜窗主瓣的寬度比周期圖主瓣的寬度寬。所以其分辨率更加降低，因此也可以說，偏移的大小反映分辨率的低與高。對平均周期圖求方差，由于是L次獨立觀測，L個周期圖相互獨立，因此平均周期圖的方差為即平均周期圖的估計方差是周期圖的方差的1/L。即：是以分辨率的降低為代價換取了估計方差的減少平均周期圖法這種方法是用一適當的窗函數W(ejω)與周期圖進行卷積，來達到使周期圖平滑的目的。式中是有偏自相關函數-(M-1)≤n≤M-12.窗函數法周期圖的窗函數法和前面提到的BT法的加權協方差譜估計是類似的。窗函數法中，周期圖和窗函數的頻譜卷積得到功率譜，等效于在頻域對周期圖進行修正，使周期圖通過一個線性系統，濾除掉周期圖中的快變成分，譜窗函數需具有低通特性。對上式求統計平均,得到有偏自相關函數的期望等效于真值與三角函數的乘積上式表明,周期圖的窗函數法仍然是有偏估計,其偏移和wB(m)、w(m)兩個窗函數有關，如果w(m)窗的寬度比較窄，當M比N小得多時，|M|<<N，則wB(m)~1,上式可以近似寫成：由于窗函數w(n)比wB(n)寬度窄，所以其頻譜的主瓣更寬，利用窗函數法可以平滑周期圖，減少估計誤差，但是偏移加大了，使分辨率降低。這種方法和平均周期圖法一樣，首先把數據長度為N的信號x(n)分成L段，每一段數據長度為M，N=LM。然后把窗函數w(n)加到每一個數據段上，求出每一段的周期圖，形成修正的周期圖，再對每一個修正的周期圖進行平均。第i段的修正周期圖為3.修正的周期圖求平均法同樣,將每一段的修正的周期圖之間近似看成互不相關,最后功率譜估計為這種在計算周期圖之前，先對各數據段加窗函數的方法，使平均周期圖的估計方差減少，當然分辨率同樣變低。但這種方法對窗函數沒有限制；此外，分段時，相鄰的兩段可以有重疊，進一步使方差減少，可以重疊50%。總之，傳統的功率譜估計方法無論采取哪一種改進方法，總是以減少分辨率為代價來換取估計方差的減少，提高分辨率的問題無法根本解決。對于由白噪聲和正弦信號或者窄帶信號組成的隨機信號，在計算周期圖之前，一般應該給數據加窗，如果不加數據窗，相近的低電平信號可能被高電平信號的旁瓣淹沒掉。在下圖a)中，ω/π=0.12處的正弦信號的旁瓣幾乎掩蓋了在ω/π=0的信號，經過給數據加窗處理后，大大壓低了旁瓣，使低電平信號清晰可見，如圖b)所示。但由于主瓣加寬，功率譜波峰變寬了，降低了信號的分辨率。一般來說，兩個等幅的正弦信號的頻率相隔很近，可以不加數據窗，頻率間隔應該大于2π/N才能分辨。利用數據窗減少窄帶過程周期圖的旁瓣a)沒有數據窗b)加哈明數據窗4.3現代譜估計中的參數建模在第一章中學習了用參數模型來描述隨機信號的方法。如果能根據觀測數據求出信號的模型(即H(z))，假設模型的系統函數為H(z)，輸入白噪聲方差為σw2，則信號的功率譜就表示為：功率譜估計可分成三個步驟：1)選擇合適的信號模型；2)根據x(n)有限的觀測數據，或者它的有限個自相關函數估計值，估計模型的參數;3)由模型的參數計算功率譜。選擇模型主要考慮的是模型能夠表示譜峰、譜谷和滾降的能力。對于具有尖峰的譜，應該選用具有極點的模型，如AR和ARMA模型；對于具有平坦的譜峰和深谷的信號，可以選用MA模型；既有極點也有零點的譜應選用ARMA模型，相對地說,ARMA模型適用范圍較寬。對于滾降太快的譜，沒有一種模型可以準確地表示功率譜，可以選用高階的AR模型近似表示。如果選擇不合適，例如，選用MA模型去估計具有尖峰的功率譜，估計效果會很差。4.3.1模型選擇下圖表示的是用MA模型估計二階AR信號功率譜的例子，圖中a)、b)、c)中的AR信號譜峰較平坦，用二階MA模型功率譜擬和真實譜時，差別較大，隨著階數的提高，估計的譜愈來愈近似于真實的譜。但是對于d)、e)、f)中AR信號的譜峰很窄，在用MA信號模型擬和時，直到MA模型階數提高到10階，其效果仍很差對AR(2)信號的模型選擇對AR(2)信號的模型選擇實際信號中一般都含有和信號不相關的噪聲，對帶有噪聲的信號，如果信號是AR模型，由于噪聲的存在需要用ARMA模型；如果用AR模型，則需要階數更高。選擇模型的另一個考慮是盡量減少模型參數。當然這和選擇模型是否合適有關系，雖然三種信號模型均有普遍應用價值，當模型選擇合適時，估計的功率譜和實際的譜擬合得好，如果不合適，只有提高階數才能得到較近似的譜，這樣需要估計的參數增多，同樣也會降低譜估計的質量。因此應該在選擇模型合適的基礎上，盡量減少模型的參數。go根據觀測數據確定模型的參數對各種功率譜估計方法不盡相同。介紹模型參數和信號自相關函數之間的關系，這些關系在功率譜估計中起著很重要的作用假設模型的差分方程和系統函數分別如下式所示：4.3.2模型參數和自相關函數之間的關系在第一章已推導出系統輸出功率譜與輸入功率譜之間的關系:式中,Pxx(z)和Pww(z)分別表示模型輸出信號x(n)和輸入信號w(n)功率譜的z變換形式;H(z)表示模型的系統函數。將上式兩邊同乘以A(z)，得1.ARMA模型的系數和信號自相關函數之間的關系go將上式進行z反變換，公式左邊的反變換為公式右邊的反變換為式中由于模型H(z)是因果的，h(n)=0，n＜0,可以得到m=0,1,2,3,…,q

m≥q+1令ARMA模型輸出自相關函數與模型參數之間的關系式ARMA模型的系數和自相關函數之間的關系如果能由信號的觀測數據估計出信號的自相關函數，可以按照上式求出ARMA模型的參數，再由模型參數求出信號的功率譜。但是在公式中由于項的存在，模型輸出自相關函數和模型參數之間的關系是非線性的，從而增加了估計功率譜的困難。但是當m>q時，上式是一個線性方程go上式共有p個方程。可以用該方程首先計算出AR部分的p個系數hA(i),

i=1,2,3,…,p;然后代入關系式，設法求出MA部分的系數。將上式以矩陣形式表示：AR模型的系統函數為：H(z)=1/A(z)，相當于ARMA模型中B(z)=1的情況，這樣在公式中m≥1m=0

2.AR模型的系數和信號自相關函數之間的關系m=0,1,2,3,…,q

m≥q+1h(0)=?整個關系式用矩陣表示為可以用模型參數表示為：為自相關矩陣，它滿足 ，且沿任一對角線的元素相等，它也是正定矩陣上面推導出的關系式或模型系數關系式確定了AR模型參數(包括模型輸入噪聲方差)和信號自相關函數之間的關系。這是一個線性方程，如果能夠由信號的觀測數據求出信號的自相關函數，可以按照該公式，通過解線性方程得到模型參數.MA模型的系統函數H(z)=B(z)，相當于ARMA模型中A(z)=1,hA(n)=δ(n)的情況，此時h(n)=hB(n)，由ARMA模型可得到MA模型的系數和自相關函數的關系m=0,1,…,q

m≥q+1MA模型的參數和自相關函數之間也是非線性關系3.MA模型的系數和信號自相關函數之間的關系m=0,1,2,3,…,q

m≥q+1三種信號模型的參數和信號自相關函數之間的關系。這些關系式提供了一種估計功率譜的方法，即首先根據信號觀測數據估計信號自相關函數，然后按照所選擇的信號模型，解上面相應的方程，求出模型參數，最后按照下式求出信號的功率譜:參數模型信號的功率譜傳統的功率譜估計方法無論采取哪一種改進方法，總是以減少分辨率為代價來換取估計方差的減少，其最大缺點是：頻譜的分辨率低MA,ARMA模型：自相關函數和模型參數之間呈非線性的關系，增加了估計功率譜的困難用n時刻之前的p個數據為：x(n-1)，x(n-2),…,x(n-p)，預測n時刻的數據x(n)：預測誤差為：4.4.1AR模型的線性預測4.4

AR模型譜估計的性質其中，為預測值，為預測系數e(n)表示線性一步預測誤差api(i=1,2,3,…，p)表示預測器的系數，它和線性預測器單位脈沖響應h(n)差一負號對差分方程進行z變換輸入是觀測數據，輸出的是預測誤差令He(z)=E(z)/X(z)，得到預測誤差濾波器的系統函數稱He(z)為一步線性預測誤差濾波器，其作用是將信號x(n)轉換成預測誤差e(n)。一般認為e(n)具有白噪聲的性質，因此He(z)也稱為白化濾波器。預測誤差濾波器/白化濾波器AR模型的系統函數為:當api=ai(i=1,2,3,…,p)時，He(z)和H(z)互為逆濾波器，

He(z)=1/H(z)=A(z)白化濾波器的系統函數:預測誤差濾波器AR模型當AR模型的階數與線性預測誤差濾波器的階數相同時，二者互為逆濾波器。AR模型和預測誤差濾波器的級聯相當于一個全通網絡。w(n)是AR模型輸入白噪聲，故預測誤差e(n)具有白噪聲的性質，所以，預測誤差濾波器也叫白化濾波器。并且，最小預測誤差功率就是激勵源白噪聲的方差σw2。由于AR模型具有這種特性，因而AR模型法也稱為線性預測AR模型法。AR模型必須是因果穩定的，即極點均在單位圓內，才能保證信號x(n)是平穩隨機信號，由于AR模型H(z)和預測誤差濾波器He(z)互為逆濾波器，所以He(z)應為最小相位系統。4.4.2預測誤差濾波器的最小相位特性自相關函數和AR模型的模型參數之間的關系服從尤勒-沃克(Yule-Walker)方程：m≥1m=0上式中，對于m≥1的情況，公式本身就是一個遞推方程，如果已由觀測數據計算出p+1個自相關函數，用,m

=

0,1,2,,…,p表示，對于m＞p的情況,可以用該公式外推得到：4.4.3AR模型隱含自相關函數延拓特性0≤m≤p

m

>p

式中，系數hA(l)是用前p+1個自相關函數求出的參數因此AR模型隱含著自相關函數外推的特性。在經典譜估計BT法中，自相關函數只能限于由觀測數據計算出的有限個自相關函數，其它的認為是0，造成了譜估計分辨率低、模糊。正是AR模型具有自相關函數外推的特性，使AR估計的譜具有高分辨率的優點。信號頻域的分辨率與時域長度的關系：4.5

AR譜估計的方法自相關法的出發點是選擇AR模型的參數使預測誤差功率最小，預測誤差功率為假設信號x(n)的數據區在0≤n≤N-1范圍，有p個預測系數，N個數據經過沖激響應為api(i=0,1,2,…,p)的濾波器，輸出預測誤差e(n)的長度為N+P4.5.1自相關法——列文森(Levenson)遞推法由信號的觀測數據，估計信號的功率譜顯然,e(n)的長度長于數據的長度，上式中數據x(n)的兩端需補充零點，這相當于無窮長的信號經過加窗處理，得到長度為N的數據。用上式對系數api求微分的方法使預測誤差功率最小，得到上面的矩陣就是Yule-Walker方程式中的自相關函數采用有偏自相關估計m=0,1,2,…,pm=-p,-p+2,…,-1因此自相關法是基于解Yule-Walker方程的一種方法。首先由信號的觀測數據估計出自相關函數，再解該方程，得到模型參數，便可求出信號的功率譜。因此該方法也稱為Yule-Walker法。但是直接解該方程，需要計算逆矩陣，很不方便。利用Yule-Walker方程中自相關矩陣的性質，可以導出Levenson-Durbin遞推法，這是一種高效的解方程方法。簡稱為列文森遞推法。i=1,2,3,…,k-1由k=1開始遞推，遞推到k=p，依次得到{a11,σ21}，{a21,a22,σ22}，…，{ap1,ap2,…,app,σ2p}。AR模型的系數以及模型輸入白噪聲方差求出后，信號功率譜為：上式表明： ，說明隨著階數增加，預測誤差功率將減少或者不變，為此要求|akk|≤1,akk稱為反射系數。遞推公式提供了一種確定模型階數的方法，如模型的階數未知，由低階開始遞推，當遞推到M階時，預測誤差滿足要求，則停止遞推，選AR模型的階數為M。遞推法效率高，當階數變化時，無需重新計算如果已知N個觀測數據（x(n),0≤n≤N-1），利用列文森遞推法計算功率譜的計算流程圖如下圖所示利用列文森遞推法計算功率譜的流程圖k=1,2,3,…,p-1該方法和自相關法一樣，仍利用預測誤差功率最小方法求模型參數，但求預測誤差功率的公式不同：對比自相關法求預測誤差功率的公式，不同的是求和限不同。該公式使用的觀測數據是已知的，不需要在數據兩端補充零點，因此與自相關法相比，去掉了加窗處理的不合理假設。4.5.2協方差法與修正協方差法1.協方差法白噪聲的方差為仍然使用梯度最小的方法求模型參數：協方差函數觀測數據x(n)(n=0,1,2,…,N-1)，利用上面公式可以求出模型的參數：{api(i=1,2,3,…,p);σ2w}。式中的協方差函數cxx(j,k)，有兩個變量，因此也適合于非平穩隨機信號。式中的協方差矩陣是埃爾米特（Hermitian）矩陣，是半正定的。這種方法近似于自相關法。一些實驗結果說明它的分辨率優于自相關法，另外對于純正弦信號數據，可以有效地估計正弦信號的頻率。修正協方差法使用前向和后向預測誤差平均值最小的方法，估計AR模型的參數，從而得到信號的功率譜。信號的前向和后向預測分別如下：2.修正協方差法式中apk是AR模型的參數最小預測誤差平均功率是模型輸入白噪聲的方差，即ρp=σ2w，前、后向預測誤差平均功率為前向和后向預測誤差功率ρpe、ρpb分別用下式表示和協方差法一樣，上式僅對用到的觀測數據的預測誤差求和。為了使預測誤差平均功率最小，求ρp對apk(k=1,2,3,…,p)的微分,或者用復梯度法求，得到化簡并寫成矩陣形式為：白噪聲的方差估計值為例:已知信號的四個觀察數據為x(n)={x(0),x(1),x(2),x(3)}={2,4,1,3}，分別用自相關法和協方差法估計AR(1)模型參數。1)自相關法：解2)協方差法:x(n)={2,4,1,3}4.5.3伯格(Burg)遞推法設信號x(n)觀測數據區間為：0≤n≤N-1，前向、后向預測誤差功率分別用ρp,e和ρp,b表示，預測誤差平均功率用ρp表示，公式分別為前向、后向預測誤差遞推公式如下：將上式帶入平均誤差公式中，得到求預測誤差平均功率ρp最小時的反射系數kp，上式就是利用伯格遞推法求第p個反射系數的公式將伯格遞推法求AR模型參數的遞推公式總結如下：利用伯格遞推法求AR模型參數的流程圖如下圖所示伯格遞推法流程圖書159頁n=p…N-1p=0

計算AR模型參數應用較多的為以下三種方法：1、自相關法---（Levinson法）2、協方差法/改進協方差法3、Burg法自相關法（Levinson法）計算簡單，但由于假設數據外為零使得分辨率相對較差協方差法與自相關法的主要區別是：求和范圍不同。數據段兩端不需要添加任何零取樣值，沒有假設數據段N以外的數據等于零。譜估計性能較好，但潛在不穩定因素改進協方差法譜估計性能最好，但計算過于復雜Burg法不直接估計AR參數，而是先估計反射系數，然后利用Levinson算法由反射系數計算AR參數計算反射系數時，使用前向、后向預測誤差功率的平均值最小準則。計算不太復雜，且給出了較好的譜估計質量.是較為通用的方法在AR模型譜估計中，模型階數的選擇是一個關鍵問題。一般模型的最好選擇是先驗未知的，實際中需預先選定模型階次。如果是純p階AR信號，選擇模型階次k＜p時，將產生對譜的平滑作用，降低譜的分辨率，如下圖所示。圖中，AR信號p=4，選擇模型的階次k=2，產生的平滑作用使兩個峰變成一個峰，分辨率明顯降低。如果選擇k≥p，且假定觀測的數據沒有誤差（沒有干擾），估計的參數應是：4.5.4關于AR模型階次的選擇AR模型階次太小時的平滑作用因此,對于純p階AR信號，應選擇階次k≥P。如果是白噪聲中的AR信號(觀測數據有誤差或者信號中含有白噪聲)，此時選擇ARMA模型合適，如選擇了AR模型，其階次應加大，較低的階次會使譜估計產生偏移，降低分辨率。當然,這也和信噪比有關，信噪比愈低，平滑作用愈嚴重，愈需要高的階次，因此信噪比低應選高的階次。一般來說，階次愈高，分辨率愈高；但階次太高，會因參數過多而使估計誤差加大，譜峰分裂，因此，對于白噪中的AR信號，其階次的選擇應折衷考慮4.6最大熵譜估計與最大似然譜估計4.6.1最大熵譜估計按照Shannon對熵的定義，當隨機變量取離散值時:式中pi是i的概率。當X取連續值時，熵的定義為1.利用最大熵的原則外推自相關函數式中,p(x)是X的概率密度函數，對于離散隨機序列,概率密度函數用聯合概率密度函數代替。顯然,熵代表一種不確定性，最大熵代表最大的不確定性，或者說最大的隨機性。下面我們研究對于有限的自相關函數值不作任何改變，對于未知自相關函數用最大熵原則外推，即不作任何附加條件的外推方法。假設x(n)是零均值正態分布的平穩隨機序列，它的N維高斯概率密度函數為式中

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

（4.6.3）

式中det(Rxx(N))表示矩陣Rxx(N)的行列式，由上式表明為使熵最大，要求det(Rxx(N)最大。若已知N+1個自相關函數值rxx(0),rxx(1),…,rxx(N)，下面用最大熵方法外推rxx(N+1)。設rxx(N+1)確實是信號自相關函數的第N+2個值，根據自相關函數的性質，由N+2個自相關函數組成的矩陣為（4.6.4）

它必須是非負定的矩陣，

即

（4.6.5）

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

（4.6.7）

上式是rxx(N+1)的一次函數，可以解出rxx(N+1)。繼續再將rxx(N+1)代入Rxx(N+2)和det(Rxx(N+2))中，求det(Rxx(N+2))對rxx(N+2)的最大值，得到rxx(N+2);以此類推，可推出任意多個其它自相關函數值，而不必假設它們為零，這就是最大熵譜估計的基本思想。2.最大熵譜估計與AR模型譜估計的等價性我們已經知道AR模型信號自相關函數與模型參數服從Yule-Walker方程，即

m≥1m=0將m≥1的情況寫成矩陣形式：

式中,ai=hA(i),i=1,2,3,…,N,ai是AR模型系數。解該方程，可以得到模型系數ai,即（4.6.8）

（4.6.9）

（4.6.10）

(4.6.11)在（4.3.6）式中，令m=N+1，得到

(4.6.12)將以上求出的系數a1,a2,…,aN代入上式，求出rxx(N+1)。而最大熵外推自相關函數的公式是(4.6.7)式，按照該公式的最后一行展開，得到（4.6.13）

上式即是最大熵外推自相關函數的公式，對比（4.6.12）式，兩公式完全一樣，證明了AR模型功率譜估計和最大熵譜估計的等價性。這里最大熵外推自相關函數等價于已知N+1個自相關函數，匹配一個N階AR信號模型的系數。一旦通過解Yule-Walker方程，解出模型參數，最大熵譜估計用下式計算信號功率譜：

（4.6.14）

最大似然譜估計是用一個FIR濾波器實現，該濾波器對所關心頻率的正弦信號，可以無失真地通過，而對于其它頻率的信號，讓其頻響盡可能地小，亦即將它們盡可能地濾除。此時,濾波器輸出的均方值，就作為信號的功率譜估計。設實信號用x(n)表示，FIR濾波器系統函數用A(z)表示：

輸出y(n)為

(4.6.15)4.6.2最大似然譜估計——最小方差譜估計式中

輸出信號的均方值為

(4.6.16)上式中T表示轉置,H表示共軛轉置,Rp=E［XXT］是Toeplith自相關矩陣,為求,必須先求FIR濾波器的系數。求這些系數的原則是：在所關心頻率ωi處，信號x(n)無失真地通過，即在ωi處的傳輸函數為1：

式中

(4.6.17)另外一個原則是在ωi附近的頻率分量盡量衰減掉，即ω≠ωi處,濾波器輸出y(n)的均方差最小,即(4.6.16)式最小,此時作為信號x(n)的功率譜估計 。因此,最大似然譜估計稱為最小方差譜估計更為合適，但由于習慣也可以仍稱為最大似然譜估計。在以上原則下，使方差最小的濾波器系數和分別為［30］、［31］

應該指出,此時并不是真正意義上的信號功率譜,只是描述了信號功率譜的相對強度。

下面分析最小方差譜估計與AR模型譜估計之間的關系：

（4.6.20）

(4.6.21)4.7特征分解法譜估計無論是實正弦波還是復正弦波，都可以用一個退化AR模型表示，設P個實正弦波組成的信號用下式表示：

（4.7.1）

式中，初相位θi是在區間（-π，π）均勻分布的隨機變量，首先分析下面的三角恒等式：

-π＜ω＜π

4.7.1正弦波用退化AR模型表示令x(n)=sin(ωn+θ),則上式變為

（4.7.2）

將上式進行Z變換，得到

（4.7.3）

這樣(4.7.2)式的特征多項式為

（4.7.4）

上式的兩個根分別是:z1=ejω,z2=e-jω,它們共軛成對，且模為1。由這兩個根可以確定正弦波的頻率。對比AR模型的系統函數，可以把正弦波信號用一個特殊的AR（2）模型表示，括弧中的2表示模型是二階的。該AR模型的激勵白噪聲方差趨于0，極點趨于單位圓。通常稱為退化的AR模型。這一模型系數有兩個，即2cosω和1，（4.7.2）式是模型的差分方程。

對于P個實正弦波,特征多項式是

（4.7.5）

上式是z-1的2P階多項式，可以表示為

（4.7.6）

注意上式中的系數ak(k=1,2,3,…,2P),必須保證它的根共軛成對。考慮到根共軛成對，也可表示為

（4.7.7）

這樣由(4.7.6)式，P個正弦波組合的模型用下面2P階差分方程描述

（4.7.8）

對于復正弦波情況，P個復正弦波組成的信號是

（4.7.9）

用一個退化的AR（p）模型表示的差分方程為

（4.7.10）

其特征多項式為

(4.7.11)其根為

1≤i≤P

注意這里的根不是共軛成對出現的。總結以上P個正弦波組合是一個退化的AR(2P)過程，獨立參量個數為P個；P個復正弦波的組合是退化的AR（P）過程，獨立參量個數仍為P個。實正弦過程相應的退化AR過程的階數比復正弦情況的階數高1倍。

白噪聲中正弦波組合的信號為

（4.7.12）

式中，w(n)為白噪聲,且

將(4.7.12)式中x(n)的用AR(2P)表示，即將(4.7.8)式帶入(4.7.12)式中，

得到

（4.7.13）

4.7.2白噪聲中正弦波組合用一特殊的ARMA模型表示將（4.7.12）式中的n用n-i代替,x(n-i)=y(n-i)-w(n-i)再將上式帶入(4.7.13)式,得到

(4.7.14)上式可以看成一個特殊的ARMA(2P,2P)模型，括弧中的兩個2P分別表示ARMA模型系統函數分子和分母的階次。它與一般的ARMA模型比較，有三方面不同：

(1)它的AR部分和MA部分具有相同的參數，它們存在共同的因子；(2)由于特征多項式(4.7.6)式的根的模為1，故AR部分特征多項式不滿足平穩性條件，MA部分特征多項式也不滿足可逆性條件；（3）AR部分的y(n)=x(n)+w(n)，y(n)是含白噪聲的觀測值，而通常為信號的x(n)不含白噪聲。這種特殊的ARMA模型結構，不能用一般的ARMA模型結構求解。下面介紹特征分解技術。將（4.7.14）式寫成矩陣形式：

YTA=WTA

(4.7.15)式中

4.7.3特征分解法譜估計用向量Y左乘(4.7.15)式,并取數學期望,得到

E［YYT］A=E［YW

T］A

(4.7.16)式中

將上面關系式帶入(4