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

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

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

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

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

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

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

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

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

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

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

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

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

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

=

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

m

>p

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

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

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

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

（4.6.3）

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

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

即

（4.6.5）

將行列式展開(kāi)，det(Rxx(N+1))是rxx(N+1)的二次函數(shù)，該二次函數(shù)系數(shù)的符號(hào)是：(-1)1+N+2(-1)1+N+1=-1，且det(Rxx(N+1))對(duì)rxx(N+1)的二次導(dǎo)數(shù)是-2det［Rxx(N-1)］，它是負(fù)值，負(fù)值表示det(Rxx(N+1))對(duì)rxx(N+1)的一次導(dǎo)數(shù)是減函數(shù)，det(Rxx·(N+1))作為rxx(N+1)的函數(shù)，凹口向下，那么只有一個(gè)最大值。為選擇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))對(duì)rxx(N+2)的最大值，得到rxx(N+2);以此類(lèi)推，可推出任意多個(gè)其它自相關(guān)函數(shù)值，而不必假設(shè)它們?yōu)榱悖@就是最大熵譜估計(jì)的基本思想。2.最大熵譜估計(jì)與AR模型譜估計(jì)的等價(jià)性我們已經(jīng)知道AR模型信號(hào)自相關(guān)函數(shù)與模型參數(shù)服從Yule-Walker方程，即

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

式中,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)式，按照該公式的最后一行展開(kāi)，得到（4.6.13）

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

（4.6.14）

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

輸出y(n)為

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

輸出信號(hào)的均方值為

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

式中

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

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

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

（4.6.20）

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

（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）

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

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

（4.7.5）

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

（4.7.6）

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

（4.7.7）

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

（4.7.8）

對(duì)于復(fù)正弦波情況，P個(gè)復(fù)正弦波組成的信號(hào)是

（4.7.9）

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

（4.7.10）

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

(4.7.11)其根為

1≤i≤P

注意這里的根不是共軛成對(duì)出現(xiàn)的。總結(jié)以上P個(gè)正弦波組合是一個(gè)退化的AR(2P)過(guò)程，獨(dú)立參量個(gè)數(shù)為P個(gè)；P個(gè)復(fù)正弦波的組合是退化的AR（P）過(guò)程，獨(dú)立參量個(gè)數(shù)仍為P個(gè)。實(shí)正弦過(guò)程相應(yīng)的退化AR過(guò)程的階數(shù)比復(fù)正弦情況的階數(shù)高1倍。

白噪聲中正弦波組合的信號(hào)為

（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)上式可以看成一個(gè)特殊的ARMA(2P,2P)模型，括弧中的兩個(gè)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)是含白噪聲的觀測(cè)值，而通常為信號(hào)的x(n)不含白噪聲。這種特殊的ARMA模型結(jié)構(gòu)，不能用一般的ARMA模型結(jié)構(gòu)求解。下面介紹特征分解技術(shù)。將（4.7.14）式寫(xiě)成矩陣形式：

YTA=WTA

(4.7.15)式中

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

E［YYT］A=E［YW

T］A

(4.7.16)式中

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