現(xiàn)代數(shù)字信號處理課件：功率譜估計

功率譜估計4.1經(jīng)典功率譜估計4.2AR模型功率譜估計的方法和性質4.3最大熵譜估計方法4.4最大似然譜估計4.5互協(xié)方差估計與互譜估計4.6特征分解法譜估計4.1經(jīng)典功率譜估計

4.1.1BT法

根據(jù)維納-辛欽定理，1958年Blackman和Tukey給出了這一方法的具體實現(xiàn)，即先由N個觀察值xN(n)估計出自相關函數(shù)r

x(m)，求其傅立葉變換，以此變換結果作為對功率譜

Px(ω)的估計。

如果我們得到的是x(n)的N個觀察值x(0)，x(1)，…，x(N－1)，令

xN(n)=a(n)·x(n)

(4.1.1)

其中a(n)是數(shù)據(jù)窗，對于矩形窗

計算rx(m)的估計值的一種方法是(4.1.2)

的均值為若a(n)是矩形窗，則所以，偏差(4.1.4)由此可以看出:

(1)這種自相關函數(shù)的估計是一個有偏估計，且估計的偏差是，當N→+∞時,。因此是rx(m)的漸近無偏估計。

(2)對于一個固定的N，當|m|越接近于N時，估計的偏差越大。

(3)由式(4.1.3)可看出，是真值rx(m)和三角窗函數(shù)

,,的乘積，q(m)的長度是2N－1，它是由矩形數(shù)據(jù)窗ar(n)的自相關所產(chǎn)生的。的方差是(4.1.5)而(4.1.6)假定x(n)是零均值的高斯隨機信號，有

E[x(n)x(k)x(n+m)]=r2x(n－k)+rx(n－k－m)rx(n－k+m)

所以

將上式和式(4.1.3)代入式(4.1.5)得(4.1.7)令n－k=i，式(4.1.7)可寫成

(4.1.8)在大多數(shù)情況下，rx(m)是平方可求和的，所以當N→+∞時，var[rx(m)]→0，又因為limbias[rx(m)]=0，所以對于固定的延遲|m|，rx(m)是rx(m)的一致估計。對由式(4.1.2)得到的自相關函數(shù)估計rx(m)進行傅立葉變換:^^^^(4.1.9)其中，v(m)是平滑窗，其寬度為2M+1，以此作為功率譜估計，即為BT譜估計。因為用這種方法求出的功率譜是通過自相關函數(shù)的估計間接得到的，所以此法也稱為間接法。4.1.2周期圖法

周期圖法是把隨機信號的N個觀察值xN(n)直接進行傅立葉變換，得到XN(ejω)，然后取其幅值的平方，再除以N，作為對x(n)真實功率譜Px(ω)的估計。以Pper(ω)表示周期圖法估計的功率譜，則^(4.1.10)其中(4.1.11)a(n)為所加的數(shù)據(jù)窗，若a(n)為矩形窗，則因為這種功率譜估計的方法是直接通過觀察數(shù)據(jù)的傅立葉變換求得的，所以習慣上又稱之為直接法。周期圖法功率譜估計的均值為(4.1.12)(4.1.13)其中4.1.3周期圖法與BT法的關系

式(4.1.9)中取M為其最大值N－1，且平滑窗v(m)為矩形窗，則

令n+m=l，上式可變成所以(4.1.14)由此可見，周期圖法功率譜估計是BT法功率譜估計的一個特例，當間接法中使用的自相關函數(shù)延遲M=N－1時，二者是相同的。

1.Welch法

假定觀察數(shù)據(jù)是x(n)，n=0，1，…，N－1，現(xiàn)將其分段，每段長度為M，段與段之間的重疊為M－…k。如圖4.1所示，第i個數(shù)據(jù)段經(jīng)加窗后可表示為

xiM(n)=a(n)x(n+ik)，i=0,1,…,L－1,n=0,1,…,M－1

其中，k為一整數(shù)，L為分段數(shù)，它們之間滿足如下關系:

(L－1)k+M≤N

該數(shù)據(jù)段的周期圖為

(4.1.15)其中(4.1.16)U為歸一化因子，使用它是為了保證所得到的譜是真正功率譜的漸進無偏估計。由此得到平均周期圖圖4.1數(shù)據(jù)分段方法(4.1.17)如果x(n)是一個平穩(wěn)隨機過程，每個獨立的周期圖的期望值是相等的，根據(jù)式(4.1.15)和式(4.1.16)有(4.1.18)其中(4.1.19)A(ω)是對應M個點數(shù)據(jù)窗a(n)的傅立葉變換，若M值較大，則Q(ω)主瓣寬度較窄，如果Px(ω)是一慢變的譜，那么認為Px(ω)在Q(ω)的主瓣內(nèi)為常數(shù)，這樣式(4.1.17)可以寫成(4.1.20)為了保證Welch法估計的譜是漸進無偏的，必須保證

(4.1.21)或(4.1.22)根據(jù)Parseval定理，式(4.1.22)可寫成所以歸一化因子U應取成(4.1.23)的方差表達式為

(4.1.24)如果x(n)是一個平穩(wěn)隨機過程，上式的協(xié)方差僅僅取決于i－l=r，令(4.1.25)式(4.1.24)可寫成單求和表示式(4.1.26)其中表示某一數(shù)據(jù)段的周期圖方差，即而是與的歸一化協(xié)方差，如果

各個數(shù)據(jù)段的周期圖之間的相關性很小，那么式(4.1.26)可近似寫成(4.1.27)這也就是說,平均周期圖的方差減小為單數(shù)據(jù)段圖方差的1/L。

2.Bartlett

對應Welch法，如果段與段之間互不重疊，且數(shù)據(jù)窗選用的是矩形窗，此時得到的周期圖求平均的方法即為Bartlett法。可以從上面討論的Welch法得到Bartlett法有關計算公式，第i個數(shù)據(jù)段可表示為

xiM(M)=x(n+iM)，i=0，1,…,L－1,n=0,1,…,M－1

(4.1.28)

其中，LM≤N。該數(shù)據(jù)段的周期圖為

(4.1.29)其中(4.1.30)平均周期圖為(4.1.31)其數(shù)學期望為(4.1.32)其中(4.1.33)將式(4.1.33)與式(4.1.13)相比，取平均情況下A(ω)的主瓣寬度是不取平均情況下A(ω)的主瓣寬度的N/M。由此可知，取平均以后，由式(4.1.32)與式(4.1.33)計算的平均周期圖偏差要比式(4.1.12)與式(4.1.13)計算的平均周期圖偏差大，同時分辨率也下降。而平均周期圖的方差仍可應用式(4.1.26)計算，由于數(shù)據(jù)段非重疊，各數(shù)據(jù)段的相關性比Welch法各數(shù)據(jù)段的相關性要小，因此平均周期圖的方差更趨向于式(4.1.27)的理論結果，但要注意，在N一定的情況下，此時所能分的段數(shù)比Welch法有重疊情況下所能分的段數(shù)L小，因此總的來說，Welch法的計算結果要比Bartlett法好。4.1.5經(jīng)典功率譜估計性能比較

為了對幾種經(jīng)典功率譜估計方法的性能進行比較，我們采用參考文獻27中的數(shù)據(jù)。信號表示為

其中包含有四個復正弦，其歸一化頻率分別是f1=0.15，f2=0.16，f3=0.252，f4=－0.16,Ak對應不同的系數(shù)，可得到不同的信噪比，本數(shù)據(jù)在f1處的信噪比為64dB，在f2處的信噪比為54dB，在f3處的信噪比為2dB，在f4處的信噪比為30dB。y(n)是一個復值的噪聲序列，其功率譜為其中，σ2=0.01，b(k)是模型系數(shù)。

x(n)的真實功率譜曲線如圖4.2(a)所示，注意其頻率范圍是從－0.5~0.5，即－π~π。令歸一化頻率f1和f2相差0.01，目的是檢驗算法的分辨能力；f3的信噪比很低，目的是檢驗算法對弱信號的檢測能力。

現(xiàn)取N＝128，圖4.2(b)示出了該數(shù)據(jù)段直接求出的周期圖，所用數(shù)據(jù)窗為矩形窗，由于主瓣過零點寬度B=2/128=0.015625>0.01，所以f1和f2不能完全分開，只是在波形的頂部能看出兩個頻率分量。

圖4.2(c)是利用Welch平均法求出的周期圖，共分四段，每段32點，沒有重疊，使用Hamming窗，這時譜變得較平滑，但分辨率降低。圖4.2功率譜估計方法比較圖4.2(d)也是用Welch平均法求出的周期圖，共分七段，每段32點，重疊16點，使用Hamming窗，譜變得更加平滑，分辨能力和圖4.2(c)大體一致。

圖4.2(e)是用BT法求出的功率譜曲線，M=32，沒有加窗；圖4.2(f)也是用BT法求出的歐尼功率譜曲線，M=16，使用了Hamming窗。

4.2AR模型功率譜估計的方法和性質

4.2.1AR模型功率譜估計的引出

由第2章討論可知，p階AR模型滿足如下差分方程:

xAn+a1xAn－1+…+apxAn－p=εn

其中，a1，a2，…，ap為實常數(shù)，且ap≠0；εn是均值為0、方差為σ2ε的白噪聲序列，也就是說，隨機信號xAn可以看成是白噪聲εn通過一個系統(tǒng)的輸出，如圖4.3所示。圖4.3AR模型信號圖4.3中

(4.2.1)而A(z)=1+a1z－1+…+apz－p

(4.2.2)已經(jīng)證明(4.2.3)其中，rA(m)是AR模型的自相關函數(shù)，尤其對于0≤m≤p，由式(4.2.3)可寫出矩陣方程為

(4.2.4)這就是AR模型的正則方程，又稱Yule-Walker方程。對于一個p－1階預測器，預測值為(4.2.5)其中，h(k)=－a(k)，預測誤差為(4.2.6)其中，a(0)=1。p階預測誤差濾波器Ap(z)如圖4.4所示。圖4.4p階預測誤差濾波器圖中，

Ap(z)=1+a(1)z－1+…+a(p)z－p

當E[e2(n)]達到其最小值E[e2(n)]min時，必滿足Yule-Walker方程

(4.2.7)當x(n)就是圖4.3所產(chǎn)生的p階AR過程xAn，也即xn=xAn或rx(m)=rA(m)時，m=0，1，…，p，必滿足關系式(4.2.8)或

A(z)=Ap(z)

(4.2.9)

此時，預測誤差濾波器Ap(z)就是AR模型H(z)的逆濾波器，實際上也就是一個白化濾波器，而且它的輸出en得到了完全的白化，也即en是一個方差為σ2ε的白噪聲εn。

通過上面的分析可以看出，對于一個p階的AR過程xn，如果首先建立階數(shù)等于p或大于p的預測誤差濾波器Ap(z)，然后以1/Ap(z)構成一個AR模型，那么以方差為E[e2(n)]min的白噪聲εn通過此線性系統(tǒng)，其輸出功率譜必定與待估計的隨機信號的功率譜完全相同，因此，模型H(z)可以完全表示出AR(p)過程的xn功率譜，它們的關系即是

(4.2.10)當然，實際待估計的隨機信號xn可能是一個階數(shù)大于p的AR過程，也可能根本就不是一個AR過程，但仍可以采用上述方法建立一個p階的AR模型，作為對隨機信號xn的功率譜估計，此功率譜估計可作為xn真實功率譜的一個近似，其步驟是:

(1)對此隨機信號xn建立p階的線性預測誤差濾波器，求得系數(shù)a(1)，a(2)，…，a(p)和E[e2(n)]min；

(2)令A(z)=1+a1z－1+…+apz－p，其中a1=a(1)，…，ap=a(p)，并構成一線性系統(tǒng)(4.2.11)那么將一方差為σ2ε的白噪聲εn通過該系統(tǒng)，其輸出的功率譜可作為待估計隨機信號xn的功率譜估計(4.2.12)其中，σ2ε=E[e2(n)]min。有下述幾點需要注意:

(1)預測誤差濾波器Ap(z)實際上就是一個白化濾波器，但一般情況下，它不能將隨機信號xn進行完全白化，所以圖4.4中預測誤差濾波器的輸出en并不是白噪聲εn,因此E[e2(n)]min表示的也是非白噪聲的方差；只有當被估計的隨機信號xn本身就是一個AR(p)過程,且當預測誤差濾波器的階數(shù)大于或等于p時，其輸出en才是一個白噪聲εn，E[e2(n)]min

即表示此白噪聲的方差。

(2)對比式(4.2.4)和式(4.2.7)，當ak=a(k)，k=1，2，…，p和E[e2(n)]min時，用這樣的方法所建立的p階AR模型，其自相關函數(shù)rA(m)和待估計隨機信號xn的自相關函數(shù)rx(m)必有如下關系:

(4.2.13)由式(4.2.13)可見，用AR(p)模型來對任意的隨機信號建模，一定具有自相關函數(shù)的匹配性質，隨著模型階數(shù)p的增大，匹配的程度會越來越好；當p→+∞時，AR模型的自相關函數(shù)和隨機信號xn的自相關函數(shù)可達到完全的匹配，即(4.2.14)而當p→+∞時，p階FIR預測誤差濾波器退化為一個IIR預測誤差濾波器，它一定可將隨機信號xn進行完全的白化。為了說明這個問題，我們又重畫出預測誤差濾波器(如圖4.5所示)，圖中G(z)是隨機信號xn的建立模型。圖4.5AR模型與預測誤差濾波器當p→+∞時因此，預測誤差濾波器(4.2.15)這證明了上述預測誤差濾波器可將隨機信號完全白化的結論。4.2.2AR模型譜估計的性質

1.隱含的自相關函數(shù)延拓的特性

在前面討論的經(jīng)典BT法功率譜估計中，假定由給定的數(shù)據(jù)

xN(n)，n=0，1，…，N－1，可估計出自相關函數(shù)rx(m)，m=－(N－1)~(N－1)，在這個區(qū)間以外，用補零的方法將其外推，對此求其傅立葉變換

^(4.2.16)就可得到BT法的功率譜估計PBT(ω)，此PBT(ω)的分辨率顯然是隨著信號長度N的增加而提高的。^^而在AR模型譜估計中，上述限制不再存在。雖然給定的數(shù)據(jù)xN(n)，n=0，1，…，N－1,是有限長度，但現(xiàn)代譜估計的一些方法，包括AR模型譜估計法，隱含著數(shù)據(jù)和自相關函數(shù)的外推，使其可能的長度超過給定的長度。前面討論的AR模型的建立，用到了單步預測的概念，預測值為這樣x(n)可能達到的長度是1~(N－1+p)，如果在遞推的過程中，用x(n)代替x(n)，那么還可繼續(xù)不斷地外推。同樣，從AR模型的建立過程看，AR過程的自相關函數(shù)必滿足^^由此式可見，AR模型的自相關函數(shù)在0~p范圍內(nèi)與rx(m)完全匹配，而在這區(qū)間外，可用遞推的方法求得。自相關函數(shù)

rA(m)實際上就是被估計信號xn的自相關函數(shù)的估計

將其進行傅立葉變換，就可得到隨機信號xn的功率譜(PSD)估計，即(4.2.17)比較式(4.2.16)和(4.2.17)可見，AR模型法避免了窗函數(shù)的影響，因此它可得到高的譜分辨率，同時它所得出的功率譜估計Px(ω)與真實的功率譜Px(ω)偏差較小。圖4.6示出了AR模型譜估計和BT譜估計法的比較。^^圖4.6AR模型譜估計和BT譜估計性能比較2.AR模型的穩(wěn)定性

AR模型穩(wěn)定的充分必要條件是H(z)=1/A(z)的極點(即A(z)的零點)必須在單位圓內(nèi)，而且這一條件也是保證x(n)是一個廣義平穩(wěn)過程所必需的。這很容易證明，如果H(z)有一個極點在單位圓外，那么x(n)的方差將趨于無窮，則x(n)是非平穩(wěn)的。

因為AR模型的系數(shù)是由Yule-Walker方程解得的，可以證明如果(p+1)×(p－1)的全相關矩陣

是正定的，AR模型的系數(shù)具有非零解，此時預測誤差濾波器A(z)一定具有最小相位的特性；換句話說，AR模型H(z)一定是一個穩(wěn)定的全極點濾波器。使用AR模型對純正弦信號建模是不合適的，因為此時Rp可能出現(xiàn)奇異或非正定的情況，但在信號處理中經(jīng)常要用正弦信號作為試驗信號以檢驗某個算法或系統(tǒng)的性能。為克服自相關矩陣奇異的情況，最常用的方法是加上白噪聲εn，這樣det(Rp)就不會等于零。

3.譜的平坦度

前面的討論已經(jīng)指出，AR(p)的系數(shù)ak就是預測誤差功率最小時的p階線性預測誤差濾波器的系數(shù)，由于預測誤差濾波器是一個白化濾波器，它的作用是去掉隨機信號x(n)的相關性，在自己的輸出端得到白噪聲εn，因此在這一節(jié)中，把白化的概念加以推廣，表明AR參數(shù)也可以使用預測誤差濾波器Ap(z)的輸出過程具有最大的譜平坦度的方法得到。利用譜平坦度的概念可以把AR譜估計得到的結果看成是最佳白化處理的結果。

功率譜密度的譜平坦度可定義為

(4.2.18)它是Px(ω)的幾何均值與算術均值之比，可以證明

0≤εx≤1

如果Px(ω)有很多峰(也就是它的動態(tài)范圍很大)，例如，在由p個復正弦所組成的隨機信號的式子中，Ak、ωk為常量，jk是在(－π，+π)范圍內(nèi)均勻分布的隨機變量，則此種信號的功率譜具有最大的動態(tài)范圍。將上式代入式(4.2.18)，分子顯見為零，因此εx=0。但如果Px(ω)是一個常數(shù)(也就是它的動態(tài)范圍為零)，也即相當于xn是一個白噪聲，則由式(4.2.18)顯見εx=1，由此可見，譜平坦度εx直接度量了譜的平坦程度。

現(xiàn)設預測誤差濾波器

為最小相位，輸入時間序列x(n)是任意的(不一定是AR過程)，按照使輸出誤差序列e(n)的譜平坦度最大的準則來確定預測系數(shù)，為此，引入下述結果:如果Ap(z)是最小相位，則

(4.2.19)計算預測誤差濾波器輸出過程e(n)的平坦度。因為對上式兩端取指數(shù)并除以，得由于對于隨機信號x(n)來說，ξx和rx(0)均是固定的，因此要使ξe最大，必須使re(0)最小，因為re(0)=E[e2(n)]，因此使預測誤差e(n)的功率譜平坦度最大和使p階預測誤差濾波器輸出的誤差功率最小是等效的，亦即條件maxξe和條件a完全等效。如果x(n)本身就是一個AR(p)過程，也即其中，A(ω)=1+a1e－jω+…+ape－jpω?，F(xiàn)使其通過一個p階預測誤差濾波器

Ap(ω)=1+a(1)e－jω+…+a(p)e－jpω

在滿足maxξe的條件下，一定有A(ω)=Ap(ω)，也即aak=a(k),

k=1，2，…，p此時預測誤差濾波器輸出的誤差序列e(n)一定是一個白噪聲序列。反之，如果將AR(p)通過一個k階預測誤差濾波器k＜p，同樣在滿足maxξe的條件下，誤差序列e(n)不可能是一個白噪聲序列，這一結果與前面討論的AR模型譜估計的引出中所得的結果是完全一致的。a在這里有一個重要的概念需強調，即預測誤差濾波器的輸入、輸出功率譜總滿足關系式在滿足E[e2(n)]min的條件下，根據(jù)求得的Ap(ω)建立AR模型:其中，σ2ε=E[e2(n)]min，是一個常數(shù)。比較這兩式可見，用PA(ω)作為Px(ω)的一個估計，其估計的好壞完全取決于Pe(ω)與一個常量相逼近的程度，換句話說，在建立AR模型時，正是由于用σ2ε代替了Pe(ω)，才使得建立的AR模型功率譜PA(ω)中喪失了很多Pe(ω)的重要細節(jié)，而只有當誤差序列e(n)是一個白噪聲序列，Pe(ω)是一個常量，并且就等于σ2ε時，才能得到Px(ω)=PA(ω)=Px(ω)。^4.2.3AR模型參數(shù)提取方法

由前面所討論的AR模型建立過程可知，AR模型參數(shù)必滿足Yule-Walker方程:對于某隨機過程來說，如果已知rx(0)、rx(1)、…、rx(p)，就可由上式求解出(σ2ε，a1，…，ap)，但實際上，對于一個所要估計的隨機過程{xn}來說，通常已知的僅僅是有限的觀察數(shù)據(jù)，如x(0)，x(1)，…，x(N－1)，而其自相關函數(shù)值通常是未知的，因此通常是首先求得自相關函數(shù)的估計值，然后再通過某種算法，得到AR模型參數(shù)的估計值。下面討論常用的幾種AR模型參數(shù)的提取方法。

1.自相關法

假定觀察到的數(shù)據(jù)是x(0)，x(1)，…，x(N－1)，而對于無法觀察到的區(qū)間(即n<0和n>N－1)，x(n)的樣本假定為零，預測誤差功率的表示式為

(4.2.20)式(4.2.20)的求和限又可以寫為[0，N－1+p]，這是因為在此區(qū)間外誤差表示式總為0的緣故。為了得到預測誤差功率的極小值，將式(4.2.20)對a(l)求微分并令其等于0，得

(4.2.21)令(4.2.22)代入式(4.2.21)得(4.2.23)寫成矩陣形式，這組方程為(4.2.24)求出白噪聲方差σ2ε的估計值σ2ε，它是

^將式(4.2.21)代入上式，得將式(4.2.25)與式(4.2.23)合并，得(4.2.26)(4.2.25)或寫成(4.2.27)例4.1

利用有限長加窗數(shù)據(jù)其中，a(n)為矩形窗，x∞(n)=rnu(n)，|r|<1，采用自相關法對其用一階的AR模型建模。解:

首先計算rx(0)和rx(1):^^

由此得而a(1)的真正解應該是－r，這說明當觀察數(shù)據(jù)為有限值N時，其解總是有偏的，直到N→+∞時，a(1)才收斂于－r。^^

2.協(xié)方差法

假定觀察到的數(shù)據(jù)是x(0)，x(1)，…，x(N－1)，寫出預測誤差功率的表示式為(4.2.28)比較式(4.2.28)和式(4.2.20)可見，協(xié)方差法與自相關法的區(qū)別主要在于預測誤差功率求和式的上下限取的不同。由于協(xié)方差法對于觀察區(qū)間[0，N－1]外的x(n)樣本并未假定為零，這就要求式(4.2.28)中x(n－k)總是落在觀察區(qū)間[0，N－1]中，為此預測誤差功率的求和上下限必須取在[p，N－1]之間。將式(4.2.28)對a(l)求微分，井令其等于0，以得到預測誤差功率的極小值，得

(4.2.29)令將其代入式(4.2.29)，得寫成矩陣形式，這組方程為

(4.2.30)求出白噪聲方程σ2ε的估計值σ2ε:^將式(4.2.29)代入上式，得將上式與式(4.2.30)合并，得

或寫成(4.2.31)比較式(4.2.31)和式(4.2.27)，表面上看，協(xié)方差法和自相關法最終所得到的正則方程具有相同的形式，但實際上兩者具有完全不同的內(nèi)容，主要表現(xiàn)在以下兩個方面。式(4.2.31)中自相關矩陣Rp是對稱的半正定矩陣，且不具有Toeplitz性質，故式(4.2.31)的求解不能采用Levinson遞歸算法來進行，雖然現(xiàn)在已經(jīng)提出一系列協(xié)方差法的求解算法，但總的來說，其算法求解的復雜性仍遠遠超過自相關法。另外由于Rp的非正定性，它也不能保證所得到的預測誤差濾波器具有最小相位特性，換句話說，利用協(xié)方差法估計出的AR模型極點不能保證在單位圓內(nèi)。舉個例子，當p=1和N＝2時，有^^而因此

顯然a(1)的幅度也可能大于或等于1。采用協(xié)方差法對信號進行建模，能夠較好地反映出信號真正的模型。例4.2

采用與自相關法討論中相同的例子，根據(jù)x(n)，n=0，1，…，N－1的觀察數(shù)據(jù)采用協(xié)方差法對其用一階的AR模型建模。首先計算rx(1，0)和rx(1，1)，即^^^由此得

由此可見，只要N≥2，所建立的模型就與信號的真正模型嚴格相同。同理可算出由此求得這說明一階AR模型和數(shù)據(jù)的真正模型完全一致。

3.修正的協(xié)方差法

假定觀察到的數(shù)據(jù)是x(0)，x(1)，…，x(N－1)，分別寫出前向預測與后向預測的表示式

和其中a(k)是AR模型系數(shù)，然后分別寫出前向預測誤差功率和后向預測誤差功率的表示式則前向預測誤差功率ρf

和后向預測誤差功率ρb

的平均值ρ為

^^^(4.2.32)修正的協(xié)方差法采用前向和后向預測誤差功率的平均值ρ為極小的方法估計AR參數(shù)。為使式(4.2.32)極小化，可以將ρ

相對于a(l)求微分，并令其等于0，得^^經(jīng)整理得

令式(4.2.34)可寫成(4.2.35)(4.2.36)l=1，2，…，p將式(4.2.36)寫成矩陣形式，即求出白噪聲方差σ2ε的估計值σ2ε:^其中利用了式(4.2.36)的結果，最后得

將上式與式(4.2.36)合并得(4.2.37)或寫成修正的協(xié)方差法的一個優(yōu)點是，它的誤差功率的計算(見式(4.2.32))是在相對于協(xié)方差法多一倍的數(shù)據(jù)點上進行，這在觀察數(shù)據(jù)長度很短的情況下，是非常有利的，但是這要求信號在正反兩個方向上呈現(xiàn)相同的特性。例如，噪聲中疊加上多個正弦信號就具有這種特性，因為正弦信號在正反兩個方向上看起來是相同的，而噪聲的自相關函數(shù)也是與方向無關的。另一方面，若對于一個在兩個方向呈現(xiàn)不同特性的信號，如一個指數(shù)衰減的信號，從相反的方向看，就是一個指數(shù)增長的信號，對這種信號如仍采用修正的協(xié)方差法來進行AR模型的估計，可能就會得到不好的結果。例4.3

采用與自相關法和協(xié)方差法討論中相同的例子。根據(jù)x(n)，n=0，1，…，N－1的觀察數(shù)據(jù)，采用修正的協(xié)方差法對其用一階的AR模型建模。

解:

首先計算rx(1，0)和rx(1，1)，即

^^和由此得因為模型的真正解是a(1)<－r，這說明此模型參量的估計是有偏的，且與數(shù)據(jù)的長度N無關。因為r<1，所以|a(1)|>|a(1)|，因此一階AR模型的極點比真正的極點更接近于單位圓。換句話說，用修正協(xié)方差法對一隨機信號進行建模，可以得到一“銳化”的譜分析結果，因此用修正協(xié)方差法對隨機信號建模，可以很好地用來進行高分辨率譜估計。

式(4.2.37)中自相關矩陣Rp是一個非Toeplitz矩陣，另外，雖然Rp是正定的(純正弦信號情況除外)，但它也不能保證所得到的預測誤差濾波器具有最小相位特性，即利用修正協(xié)方差法所估計出的極點與協(xié)方差法一樣，不能保證在單位圓內(nèi)。但是有一個例外，在一階的情況下，^^^p=1時用修正協(xié)方差法所建立的AR模型又總是穩(wěn)定的，換句話說，所估計出的極點總是在單位圓內(nèi)。這是因為在p=1時，由式(4.2.35)可知和由此得因為在觀察數(shù)據(jù)x(n)不恒為0的情況下，上式的分母總是大于分子，因此|a(1)|<1。^4.2.4AR模型階次的選擇

在AR譜估計中模型階次的選擇是一個關鍵問題，正如前面所討論的：階次太低將會導致平滑的譜估計結果，如圖4.7所示；而階次太高，將會產(chǎn)生虛假譜峰，并且估計的方差也會增大。

一種簡單而直觀的模型階次估計方法是基于預測誤差功率，因為對于所有討論過的AR模型參數(shù)估計方法，預測誤差功率都是隨模型階次的增加而減小，但是我們不能簡單地只把監(jiān)測預測誤差功率的減小作為確定模型階次的方法，還必須考慮模型階次的增加會導致譜估計方差的增大。圖4.7AR譜估計值與階次的關系

1.最終預測誤差準則

最終預測誤差(FPE)準則是通過使下式達到最小來估計模型階次的，即(4.2.39)其中，ρk是k階AR模型的預測誤差功率估計值，N是觀察數(shù)據(jù)的長度?？梢钥闯?，盡管ρk隨k的增大而減小，但(N+k)/(N－k)卻隨k的增大而增大，它表示了由于AR模型參數(shù)估計的不精確而作出的對預測誤差功率的一個修正，所以FPE(k)將有一個最小值，這個FPE(k)的最小值所對應的階便是最后確定的階。該準則的實際應用表明，雖然對于AR過程來說效果很好，但在處理地球物理數(shù)據(jù)時，一般都認為這一準則確定的階次偏低。^^

2.阿凱克信息準則

阿凱克(Akaike)信息準則(AIC)是通過使下式達到最小來估計模型階次的，即

AIC(k)=Nlnρk+2k

(4.2.40)

可以證明AIC表示的是AR模型估計的PDF和數(shù)據(jù)的真實PDF之間的庫爾貝克-利布勒(Kullback-Leibler)距離的估計值，這個準則不僅可用來確定AR模型的階次，還可用于MA模型和ARMA模型階次的確定。

AIC和FPE的性能是相近的，建議對于短觀察數(shù)據(jù)使用AIC。對于長觀察數(shù)據(jù)(N→+∞),這兩種估計方法將得到相同的模型階次估計。^

3.自回歸傳遞函數(shù)準則

自回歸傳遞函數(shù)準則(CAT)是通過使下式達到最小來估計模型階次的，即(4.2.41)這里其中，ρi是i階AR模型的預測誤差功率估計值。CAT準則選擇AR模型階次，使得由該模型估計出的預測誤差濾波器最接近于最佳無限長濾波器。^

4.3最大熵譜估計方法

給定一個寬平穩(wěn)過程的自相關rx(k)在|k|≤p時的值，要解決的問題是如何外推出k>p

時的rx(k)值。用re(k)表示外推的值，顯然應對re(k)加一些約束。例如，若(4.3.1)則Px(ejw)應對應于一個合法的功率譜，即Px(ejw)對所有w都是實值的且非負的。但一般來說,只約束Px(ejw)是實值和非負的還不能保證獲得唯一的外推結果，因此必須對可允許的外推再加上一些約束。Burg就提出一個這樣的約束，就是使隨機過程的熵達到最大化。由于熵是隨機性或不確定性的一個測度，最大熵外推等價于要找出使x(n)盡可能“白”(隨機)的自相關序列re(k)。在某種意義上，這種外推對x(n)加了盡可能少的約束或最少量的結構要求。對功率譜而言，這對應于約束Px(ejw)“盡可能平坦”。如圖4.8所示，不同的外推方法會得到不同的功率譜形狀，其中最上面的外推方法有最“平坦”的譜。圖4.8自相關序列的不同外推方法及其對應功率譜對一個功率譜為Px(ejw)的高斯隨機過程，其熵是

(4.3.2)因此對于高斯過程，若已知部分自相關序列rx(k)(|k|≤p)，其最大熵功率譜是使式(4.3.2)最大化,同時約束Px(ejw)的逆DTFT在|k|≤p時等于給定的自相關值，即(4.3.3)為求使熵最大化的re(k)值，取H(x)對r＊e(k)的導數(shù)并使之等于0，即(4.3.4)由式(4.3.1)可得，再代入上式就有

(4.3.5)定義Qx(ejw)=1/Px(ejw)，則式(4.3.5)表明Qx(ejw)的逆DTFT應是一個有限長序列qx(k),在|k|>p時值為0，即因此且高斯過程的最大熵功率譜Pmem(ejw)是一個全極點功率譜，為

^(4.3.6)利用譜因子化定理，式(4.3.6)可以表達為若用矢量ap=[1，ap(1)，…，ap(1)]和e=[1，ejw，…，ejpw]T來表示，則MEM譜可寫為(4.3.7)確定了MEM譜的形式后，剩下的是求系數(shù)ap(k)和b(0)。由于式(4.3.3)給定的約束，必須選擇這些系數(shù)使得Pmem(ejw)的逆DTFT產(chǎn)生的自相關序列與給定的rx(k)在|k|≤p時的值相匹配。在4.2節(jié)已看到，若系數(shù)ap(k)是如下自相關正則方程的解:^(4.3.8)且若則所獲得的自相關將滿足式(4.3.6)的約束。因此，MEM譜估計是(4.3.9)其中，ap是方程(4.3.8)的解。顯然，這時的MEM譜估計就是AR模型譜估計中自相關法的結果，只是表達為上式的形式。當然，若對信號的類型或熵的表述采用不同的定義，則最大熵譜估計還有其他的結果，這里不再深入討論。例4.4(噪聲中復指數(shù)的MEM估計)考慮一個白噪聲中含有復指數(shù)的隨機過程

x(n)=A1ejnw1+w(n)

式中A1=|A1|ejj，j是在[－π，π]區(qū)間均勻分布的隨機變量，白噪聲w(n)的方差為σ2w，求其p階MEM譜。

解:先解自相關正則方程以獲得AR系數(shù)ap(k)。x(n)的(p+1)×(p+1)階自相關陣是

Rx=P1e1eH1+σ2wI

其中

e1=[1，ejw1，…，ejpw1]T

P1=|A1|2Rx的逆陣是

因此(4.3.11)其中，u1=[1，0，…，0]T。

MEM譜估計應為由于

其中，WR(ejw)是矩形窗的DTFT，因此MEM譜變成為確定εp值，注意由于ap(0)=1，因此由式(4.3.11)得解上式得MEM譜又變成

與最小方差法的結果一樣，MEM譜的峰值發(fā)生在頻率w=w1處，且若信噪比很大，即P1>>σ2w，則MEM譜估計在w=w1處的值近似為因此MEM譜的峰值正比于復指數(shù)分量的平方。

4.4最大似然譜估計

當信號形式為y(t)=Aejω0t時，第k個傳感器的接收信號為

zk(t)=yk(t+τk)+vk(t)，1≤k≤p

(4.4.1)

其中，τk是在第k個傳感器上的傳輸延遲，vk(t)為在第k個傳感器上的白噪聲過程。傳輸延遲為其中，·表示向量的點積，zk為第k個傳感器的位置，k0是被傳輸?shù)男盘柕姆较?，c是光速。對于一個均勻線性傳感器陣列(見圖4.9)，有

zk=kdlx式中，lx為沿x軸的單位向量，d是傳感器之間的距離。圖4.9傳感器陣列由圖4.9可知

zk·k0=kdsinθ0

其中，θ0為接收信號的到達角。該傳感器的輸出延遲τm秒后，再乘以權系數(shù)am，將乘積項求和，可得如下形式的輸出:

(4.4.2)第m個傳感器的輸出為zm(t－τm)=y(t)+vm(t－τm)

(4.4.3)且有若vm(t)，m=(1，…，p)是互不相關的均值為0、方差為σ2v的隨機過程，則均方功率為(σ2y為信號功率)

E[x2(t)]=p2σ2y+p2σ2v

(4.4.4)

因此在理想情況下，輸出信噪比是任一傳感器上的信噪比的p倍。波束形成的目的之一是通過調節(jié)延遲來收集來自該信號源的信號能量，而拒絕所有其他信號(包括噪聲在內(nèi))。

若定義A=[a1，…，ap]T，Z=[zn，…，zn－p]T，則其采樣輸出為

其方差為式中，R=E[ZZT]，且Δt為采樣間隔。

在高斯環(huán)境下最大似然估計與最小方差估計一致，所以對x(t)的估計即為

(4.4.5)

由于為常數(shù)，故應使，即:

使方差在下列約束條件下為最小:

式中，ED=[1，D1，D2，…，DP]，D=ejωΔt。

最小化可減小來自k0以外的其他方向的能量，約束條件使接收器的增益為1。對于正弦輸入，頻率ω0處的正弦輸出信號的增益為1。引入拉格朗日乘子，構造目標函數(shù):

(4.4.6)其中，λ是拉格朗日乘子。令則可求出其中，EHD=(E*D)T。代入方差表達式得當向量ED中的頻率ω等于信號頻率時，上式就代表信號功率。當頻率ω的值遍歷定義域時，上式就可以看做是過程x(t)的功率譜密度的估值。這樣就得到了最大似然譜估計的值，為(4.4.7)可見，自相關函數(shù)矩陣R→ΦML(ω)。

可以證明:最大似然譜估計ΦML(ω)和最大熵譜估計

ΦME(ω)之間有密切聯(lián)系，即:最大似然譜估計的倒數(shù)等于1~p階的最大熵譜估計的倒數(shù)和。^^^式中，為m階AR模型(最大熵)譜估計。通常，先用自回歸AR模型法計算，m=1，2，…，p,再用上式計算最大似然譜估計。一般來說，最大似然譜估計的分辨率要比最大熵譜估計的分辨率低。

4.5互協(xié)方差估計與互譜估計

假設x(n)和y(n)都是零均值的隨機信號，因此γxy(m)=

φxy(m)。與自協(xié)方差的估計相同，互協(xié)方差估計為

(4.5.1)當y(n)=x(n)時，上式就完全轉化為自協(xié)方差估計式。式(4.5.1)的期望值為其中，jxy(m)是隨機信號x(n)和y(n)的互相關。同理

合并上述二式，得(4.5.2)可見，γxy(m)是協(xié)方差jxy(m)的一個漸近無偏估計，并且估計的方差與N成反比。為了估計互功率譜，取γxy(m)的傅立葉變換得到譜估計^^(4.5.3)若y(n)=x(n)，上式顯然轉化為周期圖。因為γxy(m)一般沒有對稱特性，所以Γxy(ω)一般是復函數(shù)。容易證明^^顯然從上式可知，Γxy(ω)是互功率譜的漸近無偏估計，但是完全與周期圖中的情況一樣，隨著N的增加，Γxy(ω)的方差并不趨向于零。為了把方差減小并把估計平滑，必須對根據(jù)較短的記錄段做出的估計進行平均或加窗處理。^^

4.6特征分解法譜估計

前面我們介紹了基于模型的寬平穩(wěn)隨機過程的功率譜估計問題，信號可以建模為白噪聲激勵的線性移不變?yōu)V波器的輸出。另一個重要的模型是x(n)為白噪聲與多個復指數(shù)分量之和，即

(4.6.1)這里假設幅度Ai是復數(shù)，有其中,fi是互不相關的且在區(qū)間[-π,π]內(nèi)均勻分布的隨機變量。4.6.1自相關陣的特征分解

為理解自相關陣的特征分解可作為頻率估計的一種方法，先來看一個一階諧波過程:

x(n)=A1ejnw1+w(n)

它表示白噪聲中有一個復指數(shù)分量，復指數(shù)的幅度是Ai=|Ai|ejf1，f1是均勻分布的隨機變量，w(n)是方差為σ2w的白噪聲。前面已介紹過，該x(n)的自相關序列是

rx(k)=P1ejkw1+σ2wδ(k)

其中，P1=|A1|2是復指數(shù)分量的功率。因此，x(n)的M×M階自相關陣是信號的自相關陣Rs和噪聲的自相關陣Rn的和，即

Rx=Rs+Rn

(4.6.2)其中信號自相關陣的具體形式是(4.6.3)其秩為1。噪聲自相關陣是個對角陣Rn=σ2wI，它是滿秩的。注意，若我們定義

e1=[1，ejw1，ej2w1，…，ej(M－1)w1]T

(4.6.4)則Rs可用e1表示為Rs=P1e1eH1由于Rs的秩為1，因此Rs只有一個非零特征值。又由

Rse1=P1(e1eH1)e1=P1e1(eH1e1)=MP1

得Rs的非零特征值應等于MP1，且e1就是對應的特征矢量。另外，由于Rs是共軛對稱陣，因此剩下的特征矢量v2,v3，…，vM應與e1正交(共軛對稱陣的不同特征值對應的特征矢量是相互正交的)，即

e1Hvi=0，i=2，3，…，M

(4.6.5)

最后，若用λsi表示Rs的特征值，則有

Rsvi=(Rs+σ2wI)vi=λsivi+σ2wvi=(λsi+σ2w)vi

(4.6.6)

這說明Rx的特征矢量和Rs的特征矢量是相同的，只是Rx的特征值是

λi=λsi+σ2w其中，Rx的最大特征值是λmax=MP1+σ2w，剩下的Rx的特征值都等于σ2w。這樣，就可以由Rx的特征值和特征矢量獲得所有的有關x(n)的參數(shù)，具體方法如下:

(1)對自相關陣Rx進行特征分解，其最大特征值應等于MP1+σ2w，剩下的特征值都是σ2w。

(2)利用Rx的特征值求信號功率P1和噪聲方差σ2w得

(3)用對應于最大特征值的特征矢量vmax確定頻率w1，例如，利用vmax的第二個系數(shù)vmax(1)，應有

例4.5

設x(n)是一個一階諧波過程，為白噪聲加一個復指數(shù)信號x(n)=A1ejnw1+w(n),設其自相關陣為,求x(n)中復指數(shù)信號的頻率。

解:

由于x(n)的自相關陣Rx的特征值為

特征矢量為因此白噪聲的方差是復指數(shù)分量的功率是

最后得復指數(shù)分量的頻率為這里要指出的是，對一階諧波過程，實際上可以更容易地求出該過程的參數(shù)。由于

，所以立即可得，w1=π/4，而有了P1，白噪聲的方差就可以確定如下:實際上，以上介紹的方法在頻率估計中的實用價值是有限的，因為它要求精確地已知自相關矩陣。雖然也可以用估計的自相關陣代替Rx，但這樣做后，其最大特征值只是近似等于MP1+σ2w，對應的特征矢量也只是近似為e1，而特征值和特征矢量可能對rx(k)的很小的估計誤差就非常敏感。因此一般不用一個特征矢量來估計復指數(shù)的頻率，而利用的是它們的加權平均。設vi表示Rx的噪聲特征矢量，即特征值為

σ2w的特征矢量，vi(k)表示vi的第k個分量。如果計算vi各系數(shù)的DTFT為(4.6.7)則式(4.6.5)所給定的正交性條件就意味著Vi(ejw)在w=w1處(即復指數(shù)的頻率處)將等于0，因此，若我們構造的頻率估計函數(shù)如下:

(4.6.8)則Pi(ejw)在w=w1處的值將很大(理論上是無窮大)。這樣，通過對頻率估計函數(shù)的峰值進行定位，就可以估計出復指數(shù)的頻率。但是，式(4.6.8)只利用了一個特征矢量，因此可能對Rx的估計誤差很敏感，所以可以考慮利用所有噪聲特征矢量的加權平均，即^(4.6.9)其中，ai是適當選擇的常數(shù)。

現(xiàn)在來考慮白噪聲中有兩個復指數(shù)分量的情況，即

x(n)=A1ejnw1+A2ejnw2+w(n)

其中，Ai=|Ai|ejj1(i=1，2)，是復指數(shù)分量的幅值，w1和w2是其頻率，且w1≠w2。若w(n)是方差為σ2w的白噪聲，則x(n)的自相關序列是

其中，P1=|A1|2，P2=|A2|2。顯然其自相關陣可寫成如下幾項的求和:

其中，Rs=P1e1eH1+P2e2eH2代表的是Rx中的信號分量，且Rs的秩為2；而Rn=σ2wI是Rx中的噪聲分量，Rn是一個對角陣。Rx的更簡潔的分解形式是

Rx=EPEH+σ2wI

(4.6.10)

其中，E=[e1，e2]是M×2的矩陣，且由兩個信號矢量e1和e2組成，P=diag[P1，P2]是信號功率值組成的對角陣。

除了像式(4.6.10)那樣將Rx分解為兩個自相關陣之和外，也可以對Rx進行特征分解。設vi和λi分別表示Rx的特征矢量和特征值，且特征值按降序排列為

λ1≥λ2≥…≥λM

由于Rx=Rs+σ2wI，因而有

λi=λsi+σ2w

其中，λsi是Rs的特征值。由于Rs的秩等于2，所以Rs只有兩個非零特征值，且都大于零(Rs是非負定的)，因此Rx的前兩個特征值都大于σ2w，其他剩下的特征值等于σ2w。這樣Rx的特征值和特征矢量可以分為兩組:第一組是兩個特征值大于σ2w的特征矢量，稱為信號特征矢量，張成一個二維的子空間，稱為信號子空間；第二組是那些特征值等于σ2w的特征矢量，稱為噪聲特征矢量，張成一個M－2維的子空間，稱為噪聲子空間。注意Rx是共軛對稱矩陣，所以其特征矢量vi形成一個正交集合，這樣，信號和噪聲子空間就應是相互正交的，也就是說，對信號子空間中的任意矢量u和噪聲子空間的任意矢量v，有uHv=0。與單個復指數(shù)分量時的情況不同，對兩個復指數(shù)分量加上白噪聲的情況，信號特征矢量一般不等于e1和e2(因為信號矢量一般不是正交的，而特征矢量總是相互正交的)，但e1和e2將位于信號特征矢量v1和v2所張成的子空間中，且由于信號子空間和噪聲子空間相互正交，所以e1和e2與噪聲特征矢量vi是正交的，即

(4.6.11)這樣，與單個復指數(shù)分量一樣，復指數(shù)的頻率w

1和w

2可以用如下形式的頻率估計函數(shù)來進行估計:最后來看更一般的情況，即p個不同頻率的復指數(shù)加上白噪聲構成一個寬平穩(wěn)過程，其自相關序列應是其中，Pi=|Ai|2是第i個分量的功率。它的M×M的自相關陣是(4.6.12)其中，ei=[1，ejwi，ej2wi，…，ej(M－1)wi]T(i=1，2，…，p)是p個線性獨立的矢量。就像前面一樣，式(4.6.12)可以更簡潔地表達為

Rx=EPEH+σ2wI

(4.6.13)其中，E=[e1…ep]是p個信號矢量ei組成的M×p的矩陣，P=diag[P1，…，Pp]是信號功率值組成的對角陣。由于Rx=

Rs+σ2wI(λsi是Rs的特征值)，且Rs的秩是p，因此Rx的前p個特征值將大于σ2w，其他剩下的特征值等于σ2w。這樣Rx的特征值和特征矢量可以分為兩組:第一組是特征值大于σ2w的信號特征矢量v1，…，vp，另一組是特征值等于σ2w的噪聲特征矢量vp+1，…，vM。假設特征矢量都歸一化為單位模，則Rx可以分解為若寫成矩陣形式，則有

Rx=VsVssVHs+VnVnnVHn

(4.6.14)

其中，Vs=[v1…vp]是信號特征矢量組成的M×p矩陣，Vn=

[vp+1…vM]是噪聲特征矢量組成的M×(M－p)矩陣，Vss和

Vnn分別是特征值λi=λsi+σ2w和λi=σ2w組成的對角陣。

后面我們將要把一個矢量投影到信號子空間或噪聲子空間，完成這種投影的投影矩陣Ps和Pn應是

Ps=VsVHs，Pn=VnVHn

(4.6.15)

與前面一樣，信號子空間和噪聲子空間的正交性可用于估計復指數(shù)的頻率。具體來說，由于每個信號矢量e1…ep都屬于信號子空間，所有的ei與每個噪聲特征矢量都是正交的，即

eHivk=0，i=1，2，…，p，k=p+1，p+2，…，M

結果頻率可用如下形式的頻率估計函數(shù)進行估計:(4.6.16)下面各節(jié)將基于上式給出幾種不同類型的頻率估計算法。首先介紹Pisarenko諧波分解，它就采用式(4.6.16)的頻率估計函數(shù)進行估計，并取M=p+1，aM=1。4.6.2Pisarenko諧波分解

Pisarenko分解中，假定x(n)是由p個復指數(shù)加上白噪聲組成，且復指數(shù)的個數(shù)p是已知的。又假設已知p+1個自相關序列值，或它們已被估計出來，因此有一個(p+1)×(p+1)的自相關陣，噪聲子空間的維數(shù)是1，它由最小特征值

λmin=σ2w所對應的特征矢量張成。用vmin表示該噪聲特征矢量，則vmin將與每個信號矢量ei正交，即

(4.6.17)因此，

在每個復指數(shù)頻率wi(i=1，2，…，p)處都等于零。這樣，噪聲特征矢量的Z變換(稱為特征濾波器)就在單位圓上有p個零點，即

(4.6.18)而且由特征濾波器的根就可以獲得各復指數(shù)的頻率。若不想求Vmin(z)的根，也可以構造如下的頻率估計函數(shù):它是式(4.6.16)在M=p+1，ap+1=1時的特例。由于PPHD(ejw)在復指數(shù)的頻率處其值很大(理論上是無窮大)，對PPHD(ejw)進行峰值定位就可得到頻率的估計。雖然PPHD(ejw)是功率譜的形式，但它稱為偽譜(有時稱為特征譜)，因為它并沒有包含復指數(shù)的任何功率信息，也不包含噪聲的功率信息。^^^一旦確定了復指數(shù)的頻率以后，就可以由Rx的特征值求功率譜Pi了。假設信號子空間的特征矢量v1，…，vp已被歸一化，即vivHi=1，由于

Rxvi=λivi，i=1，2…，p

將上式兩邊左乘vHi，則得

vHiRxvi=λivHivi=λi，

i=1，2…，p

(4.6.19)將式(4.6.12)的Rx表達式代入上式，有簡化后，得(4.6.20)注意，求和項中的|eHkvi|2是頻率wk處信號子空間特征矢量vi的DTFT的平方幅度，即有

|eHkvi|2=|Vi(ejwk)|2

其中，。因此式(4.6.20)可以寫成

上式代表p個線性方程，它有p個未知量Pk，寫成矩陣形式為(4.6.21)求解該方程組就得到功率Pk。

例4.6(用Pisarenko法在噪聲中估計兩個復指數(shù)分量)假設有兩個復指數(shù)加白噪聲的隨機過程，其自相關序列的前三個值是

rx(0)=6

rx(1)=1.92705+j4.58522

rx(2)=－3.42705+j3.49541

試用Pisarenko法求復指數(shù)分量的頻率和功率。

解:

由于p=2，必須對3×3的自相關陣Rx進行諧波分解。首先有

其特征矢量為

特征值為λ1=15.8951，λ2=1.1049，λ3=1.0000因此，最小特征值是λ3=1.0000，對應的特征矢量是于是得特征濾波器Vmin(z)，求得其根為z1=0.5+j0.8660=ejπ/3，z2=0.3099+j0.9511=ej2π/5因此兩個復指數(shù)的頻率是為求復指數(shù)的功率，必須計算信號特征矢量v1和v2的DTFT的平方幅值在w1和w2處的值，有|V1(ejw1)|2=2.9685，|V1(ejw2)|2=2.9861|V2(ejw1)|2=0.0315，|V2(ejw2)|2=0.0139因此式(4.6.21)就變成其中，σ2w=λmin=1。求解上式得P1=2，P2=3。例4.7(用Pisarenko法估計單個正弦信號)設x(n)是隨機相位正弦波加白噪聲形成的隨機過程:

x(n)=Asin(nw0+φ)+w(n)

已知rx(0)=2.2，rx(1)=1.3，rx(2)=0.8，則其3×3的自相關矩陣是

Rx=Toep{2.2，1.3，0.8}

可計算出其特征值為

λ1=4.4815，λ2=1.4，λ3=0.7185

因此白噪聲功率是σ2w=λmin=0.7185。Rx的特征矢量是

注意，各特征矢量都有對稱性。為估計正弦波的頻率，應求特征濾波器Vmin(z)的根。Vmin(z)是v3的Z變換,即

Vmin(z)=0.4437(1－1.755z－1+z－2)

可求出Vmin(z)的根為z=e±jw0，其中w0滿足

2cosw0=1.755，即w0=0.159π

最后用式(4.6.21)估計正弦波的功率。由于x(n)只包含一個正弦波，實際上可以簡化計算，不需求解方程。由于該隨機過程的自相關序列是因此，，而已知rx(0)=2.2，σ2w=0.7185，因此有A2=2.963。雖然Pisarenko諧波分解在數(shù)學上很精致，但它并不實用，原因之一是必須已知復指數(shù)的個數(shù)。在自相關值精確已知的理想情況下，可以先估計復指數(shù)的個數(shù)，然后檢驗最小特征值的重復性，就可以確定復指數(shù)的個數(shù)。但自相關是估計值時，就不便采用重復性檢驗準則了，因為很少會出現(xiàn)兩個最小特征值相等的情況。Pisarenko分解的另一個局限性是它假定加性噪聲是白噪聲，但實際中一般都不是這樣，因此會導致頻率估計的偏差。雖然可以修正Pisarenko方法使其適用于非白噪聲的情況，但要求已知加性噪聲的功率譜。從計算量的角度來看，Pisarenko諧波分解要求出信號自相關矩陣的最小特征值和特征矢量，對高階問題，這需要很大的計算量，但可以用基于Levinson-Durbin遞歸的迭代算法來求最小特征值和特征矢量，以提高計算效率。4.6.3MUSIC算法

假設隨機過程x(n)有p個復指數(shù)分量和方差為σ2w的白噪聲，Rx是其M×M的自相關矩陣，且M