數(shù)字信號處理-第四章

第四章功率譜估計4.1引言4.2經(jīng)典譜估計4.3現(xiàn)代譜估計中的參數(shù)建模4.4

AR模型譜估計的性質(zhì)4.5

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

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

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

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

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

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

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

m≥q+1令A(yù)RMA模型輸出自相關(guān)函數(shù)與模型參數(shù)之間的關(guān)系式ARMA模型的系數(shù)和自相關(guān)函數(shù)之間的關(guān)系如果能由信號的觀測數(shù)據(jù)估計出信號的自相關(guān)函數(shù)，可以按照上式求出ARMA模型的參數(shù)，再由模型參數(shù)求出信號的功率譜。但是在公式中由于項(xiàng)的存在，模型輸出自相關(guān)函數(shù)和模型參數(shù)之間的關(guān)系是非線性的，從而增加了估計功率譜的困難。但是當(dāng)m>q時，上式是一個線性方程go上式共有p個方程?？梢杂迷摲匠淌紫扔嬎愠鯝R部分的p個系數(shù)hA(i),

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

2.AR模型的系數(shù)和信號自相關(guān)函數(shù)之間的關(guān)系m=0,1,2,3,…,q

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

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

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

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

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

=

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

m

>p

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

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

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

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

（4.6.3）

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

它必須是非負(fù)定的矩陣，

即

（4.6.5）

將行列式展開，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ù)，凹口向下，那么只有一個最大值。為選擇rxx(N+1)使det(Rxx(N+1)最大，解下列方程：(4.6.6)用數(shù)學(xué)歸納法，得到

（4.6.7）

上式是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);以此類推，可推出任意多個其它自相關(guān)函數(shù)值，而不必假設(shè)它們?yōu)榱?，這就是最大熵譜估計的基本思想。2.最大熵譜估計與AR模型譜估計的等價性我們已經(jīng)知道AR模型信號自相關(guān)函數(shù)與模型參數(shù)服從Yule-Walker方程，即

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

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

（4.6.9）

（4.6.10）

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

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

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

（4.6.14）

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

輸出y(n)為

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

輸出信號的均方值為

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

式中

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

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

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

（4.6.20）

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

（4.7.1）

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

-π＜ω＜π

4.7.1正弦波用退化AR模型表示令x(n)=sin(ωn+θ),則上式變?yōu)?/p>

（4.7.2）

將上式進(jìn)行Z變換，得到

（4.7.3）

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

（4.7.4）

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

對于P個實(shí)正弦波,特征多項(xiàng)式是

（4.7.5）

上式是z-1的2P階多項(xiàng)式，可以表示為

（4.7.6）

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

（4.7.7）

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

（4.7.8）

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

（4.7.9）

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

（4.7.10）

其特征多項(xiàng)式為

(4.7.11)其根為

1≤i≤P

注意這里的根不是共軛成對出現(xiàn)的?？偨Y(jié)以上P個正弦波組合是一個退化的AR(2P)過程，獨(dú)立參量個數(shù)為P個；P個復(fù)正弦波的組合是退化的AR（P）過程，獨(dú)立參量個數(shù)仍為P個。實(shí)正弦過程相應(yīng)的退化AR過程的階數(shù)比復(fù)正弦情況的階數(shù)高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模型系統(tǒng)函數(shù)分子和分母的階次。它與一般的ARMA模型比較，有三方面不同：

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

YTA=WTA

(4.7.15)式中

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

E［YYT］A=E［YW

T］A

(4.7.16)式中

將上面關(guān)系式帶入(4