基于廣義柯西分布的最大后驗準則頻譜估計方法

摘 要:經典的頻譜估計方法和現代的頻譜估計方法在低信噪比及小數據量的情況下，譜估計的分辨率和方差性能不能滿足實際應用需要。因此，提出一種高分辨率、高精度DFT變換的新方法，此方法特別適用于線性頻譜的估計。該方法基于最大后驗概率準則，建立廣義柯西-高斯分布模型，克服了短數據情況下的DFT變換分辨率低的缺點，具有收斂快、頻率分辨率高、頻率精度高的優點。計算機仿真結果證實了新方法的有效性。

關鍵詞:最大后驗概率; 離散傅里葉變換; 頻譜估計; 廣義柯西分布

中圖分類號:TN911.6 文獻標識碼:A

文章編號:1004-373X(2010)07-0017-04

New Method for Spectrum Estimation Based on Generalized Cauchy Distribution and MAP

SONG Jun-cai， ZHANG Shu

(College of Information and Communication Engineering， Harbin Engineering University， Harbin 150001， China)

Abstract: In the low SNR and small amount of data， the resolution and variance performance of spectral estimation can not meet the actual requirement by using classical or modern spectrum estimation methods. Therefore， a new high-resolution and high-precision method of DFT transform is proposed. It is suitable for estimation of linear spectra. Based on maximum a posteriori probability criterion， a generalized Cauchy-Gaussian distribution model to overcome the low resolution of DFT in the case of short data is established. The proposed method has advantages of fast convergence， high efficiency and high accuracy.The results of computer simulation show that the novel method is effective.

Key words: maximum a posterior probability; discrete Fourier transform; spectrum estimation; generalized Cauchy distribution

0 引 言

信號的頻譜分析是研究信號特征的重要手段之一，該技術在雷達、通信、震動、地震信號處理及電子監測領域有著廣泛的應用。理論上的傅里葉變換是對整個時域信號的變換，但現實中卻只能對有限長度的信號進行變換。經典的譜估計方法主要是周期圖法及改進的周期圖估計方法，但這些方法都存在著頻率分辨率低、頻譜泄露和方差性能不好等缺點。原因是周期圖頻譜估計是對數據加窗截斷，用有限數據樣本或自相關函數來估計無限個數據的頻譜?，F代的頻譜估計技術有基于ARMA模型的估計技術，基于預測外推的最大熵法(MEM)，最大似然Capon估計等，這些技術較經典頻譜估計有更高的頻譜分辨率。但帶來的問題是對模型的參數的估計，或者要求較高的信號噪聲比(SNR)等[1，2]。同時這些方法除了計算量大之外，也可能出現數值不穩定現象。Sacchi和Ulrych等人提出的基于最大后驗概率柯西-高斯模型的頻譜估計方法[3]已解決經典頻率估計和現代頻率估計技術中存在的缺點，該估計方法在小數據量的前提下有效地提高了頻譜估計的精度。

本文推廣了柯西-高斯模型，提出基于廣義柯西分布的頻譜估計方法，仿真結果顯示本方法較基于柯西分布的算法性能更好。

1 基于最大后驗準則的DFT

考慮一個包含N點的采樣時間序列，則該離散序列x0，x1，x2，…，xN-1的DFT變換可表示為:

Xk=∑N-1n=0xne-i2πnk/N，k=0，1，…，N-1

(1)

同樣，該離散序列的IDFT可表示為:

xn=1N∑N-1k=0Xkei2πnk/N， n=0，1，…，N-1

(2)

以下為了敘述方便，分別將x稱為數據樣本，X稱為頻域樣本。

現在希望由N個數據樣本得到M>N點的頻域樣本值。求解該問題的最直接方法為在N個數據樣本后補M-N個零，然后進行DFT變換。該方法雖然獲得了M個頻率點值，但這種基于原有N個樣本的插值算法，由于沒有增加原始數據信息，所以沒有減小時域窗而導致的頻譜旁瓣過大、頻譜稀疏性差的缺點。

當M>N時，將數據補M-N個零，進行DFT變換得到X。而對X進行IDFT時，以矩陣形式表示式(2)時，要改寫成:

x=FX

(3)

其中:x∈RN和X∈CM分別表示已知數據向量和其相應的DFT，F為N×M矩陣，Fn，k=1Mei2πnk/M。

已知當M>N時，由式(3)確定的線性方程組求解問題為非確定性病態問題，為了使式(3)有惟一解，必須增加相應的約束條件。根據文獻[5]，通過最小化下列表達式，可以獲得在某種條件限制下的惟一解:

J(X)=Φ(X)+‖x-FX‖22

(4)

式中:‖#8226;‖22表示Frobenius范數，‖X‖22=XHX=∑kX*kXk。約束項Φ(X)代表了增加的限制條件。從而在M>N時，惟一解的問題轉化為如何尋找到相應的Φ(X)。

假設數據包含噪聲，噪聲服從N(0，σ2n)分布，數據樣本x的條件分布概率表示為p(xX，σn)。設頻域樣本X的先驗分布概率為p(XσX)，則根據貝葉斯原理得到后驗分布概率:

p(Xx，σX，σn)=p(XσX)p(xX，σn)p(xX，σX，σn)

(5)

即式(5)是已知σX和σn的最大后驗概率估計器。由于式(5)的分母可以認為是常數，從而最大后驗等價于使式(5)的分子最大，因此不僅與數據中的噪聲分布有關，同時和X的先驗分布有關。下面將看到不同的先驗概率分布對X估計的影響。

(1) 當X的先驗概率分布為復高斯分布時，即:

p(XσX)=12πσ2XM-1e-(1/2σ2x)‖X‖22

(6)

式中:σ2X為信號功率。使式(5)最大，即等價于使下式最小:

Jgg=12σ2X‖X‖2+12σ2n‖x-FX‖2

(7)

其中:下角標gg表示數據樣本中的噪聲及頻域樣本都服從高斯分布。這種模型簡稱為GG模型。

將式(7)和式(4)相比，可以知道兩者具有相同的形式，即當假設X的先驗概率為高斯分布時，Φ(X)=‖X‖22。微分式(7)并使其為零，可得到:

=1M+λ-1FHx

(8)

其中因子λ=σ2n/σ2X為噪聲信號比。對比式(3)和式(8)，可知此時的DFT即為帶常數加權因子的補零DFT變換。當λ=0時，式(8)則退化為補零的DFT變換。

(2) 當X的先驗概率分布為柯西分布時，即:

p(XkσX)∝11+XkX*k2σ2X

(9)

可得到:Φ(X)=∑M-1k=0ln1+XkX*k2σ2X，從而:

Jcg=∑M-1k=0ln1+XkX*k2σ2X+12σ2n‖x-FX‖2

(10)

其中:下角標cg表示頻域樣本及數據樣本噪聲分別服從柯西分布和高斯分布。同樣這種模型簡記為CG模型。

經過如前的推導，得到的DFT變換為:

=[λQ-1+FHF]-1FHx

(11)

其中:Q是M×M對角矩陣:

Qii=1+XiX*i2σ2X， i=0，1，…，M-1

(12)

對比式(8)和(11)，可知此時的DFT仍為加權的補零DFT變換，但與式(8)不同的是式(11)中的λ=σ2n/σ2X被Q-1加權。當λ=0時，式(11)同樣退化為補零的DFT變換。

后面的仿真結果證明由式(11)得到的DFT優于式(8)的結果。

2 基于廣義柯西分布的MAP的DFT

分析上面的GG及CG模型所得到的結果可知，當X的先驗概率分布很窄時，如GG模型采用的高斯分布，由于不能充分利用數據樣本的統計特性，導致GG模型得到的結果與數據樣本補零后DFT的結果相同。而X的先驗概率展寬后，如CG模型采用的柯西分布，數據樣本的統計特性得到充分的利用，從而得到的結果要明顯優于GG模型的結果。因此為了得到性能優良的頻譜估計就要尋找具有良好拖尾特性的先驗分布[1]，以進一步改進基于MAP的頻譜估計特性。

2.1 廣義柯西分布

對于實隨機變量Xk，如果其概率密度函數為:

pm(XkσX)=am1+X2k2σ2Xm，

m>0.5(13)

則定義Xk服從參數為m的廣義柯西分布。參數m是大于0.5的實數，am是歸一化常數。

可以證明廣義柯西分布具有以下性質:

(1) 當0.5

(2) 當m=1時，廣義柯西分布即是常規定義下的柯西分布。

(3) 廣義柯西分布不同于下列分布fm(XkσX):

fm(XkσX)=bm1+Xk2σX2m， m>0.5

(14)

但當m為整數時兩者相等。

(4) 在所有的分布中，廣義柯西分布具有最大的散布特性，即廣義柯西分布具有最大的拖尾概率。

在圖1中給出了高斯、柯西及廣義柯西分布(m=0.6，m=0.9)的四條分布曲線。為了便于比較它們的展布特性，對密度函數的原點值作了歸一化處理。

由圖1可以清楚地看到在這些曲線中，高斯分布具有最差的拖尾特性?？挛鞣植纪衔蔡匦员雀咚狗植家忻黠@展寬。而對廣義柯西分布，隨著m值的減少，曲線的拖尾概率比柯西分布進一步增大。

圖1 三種分布的歸一化概率密度

2.2 基于廣義柯西分布的MAP的DFT

假設X的先驗概率p(XσX)服從廣義柯西分布:

pm(XkσX)∝11+‖Xk‖222σ2Xm， m>0.5

(15)

對于廣義柯西分布先驗概率，要使后驗概率最大等價于最小化Jgcg:

Jgcg=S(X)+12σ2n(x-FX)H(x-FX)

(16)

其中:下角標gcg表示廣義柯西-高斯分布模型(記為GCG模型)。同樣可以求得約束項S(X)為:

S(X)=∑M-1k=0ln1+‖Xk‖222σ2Xm

(17)

對式(16)求導并令導數為零，得到GCG模型下的DFT。

=[mλ(G1-m+G)-1+FHF]-1FHx

(18)

其中:G為M×M為對角矩陣，Gkk=‖Xk‖222σ2X。

對比式(11)和(18)，可知GCG模型的DFT仍為加權的補零DFT變換。與式(11)不同的是加權形式更加復雜，但當λ=0時，式(18)也退化為補零的DFT變換。

顯然當m=1時，式(18)轉化為式(11)，即柯西-高斯模型的解，只是廣義柯西-高斯分布模型解的一個特例。這個解的性能將在下一部分通過仿真來說明。

式(18)給出了基于廣義柯西分布的最大后驗概率準則下的頻域樣本估計值的表達式。表達式(18)具有簡潔的表達形式，但是在實際計算時并不方便，因為該等式的右邊含有與X有關的G。為此將式(18)變換為遞推迭代形式:

b(μ-1)=[λIN+FHG1-m+Gμ-1F]-1FHx

X(μ)=(G1-m+G)μ-1FHbμ-1

(19)

其中:μ表示迭代的次數，一般經過小于10次迭代后，即能得到滿意的結果。

3 仿真結果

下列仿真中，數據為幅度為1的0.2 Hz和0.21 Hz的正弦波及噪聲之和，數據樣本點數N=50。噪聲為σn=0.2的高斯白噪聲。用下式計算信噪比SNR=(P1+P2)/σ2n，得到SNR=14 dB。選用1 024點DFT，迭代次數μ=9。圖2表示GG，CG和GCG三種模型下的頻譜圖。圖標中m=1即為CG模型，GCG模型中選擇m=0.9，m=0.6兩種不同數值作為比較。

圖2 不同m值下GCG模型下頻譜圖對比

從圖2可以看出，GG模型給出的頻譜估計完全不能區分0.2 Hz和0.21 Hz兩個頻率分量，而CG模型和GCG模型均具有很好的頻率分辨力，可清晰地區分出0.2 Hz和0.21 Hz兩個頻率分量。但是GCG模型具有較CG模型要低50 dB以上的旁瓣抑制特性和更好的頻譜稀疏性。

圖3和圖4是取不同信噪比及不同m值時，GCG模型的DFT變換得到的0.2 Hz及0.21 Hz的頻率估計的相對誤差。圖3和圖4顯示CG和GCG兩種算法在高信噪比的情況下均具有良好的頻率估計精度。而在信噪比較低的情況(SNR<8 dB)，GCG算法(m=0.9，m=0.6)的頻率估計的相對誤差對信噪比的敏感性明顯優于CG模型(即m=1)的結果。

圖3 0.2 Hz處的估計偏差

圖4 0.21 Hz處的估計偏差

4 結 語

提出了一種高分辨DFT變換的新方法。根據最大后驗概率準則，在信號頻譜幅度的先驗概率為廣義的柯西分布，噪聲為高斯分布的假設下，通過最小化代價函數求解出GCG模型下的DFT變換的迭代公式。該方法是CG模型下DFT變換的推廣， 具有頻率分辨率高的優點，適用于短數據樣本下的頻率估計。理論分析和仿真實驗證明了新方法的正確性和有效性。

參考文獻

[1]KAY S M， MARPLE L M. Spectrum analysis: a modern perspective[J]. Proc. IEEE， 1981， 69(18): 216-223.

[2]CHEN C H. Nonlinear Maximum Entropy Spectral Analysis Method for Signal Recognition[M]. [S.l.]: Research Studies Press， 1982.

[3]SACCHI M D， ULRYCH T J， WALKER C J. Interpolation and extrapolation using a high-resolution discrete Fourier transform[J]. IEEE Trans. on Signal Processing， 1998， 46(9): 106-112.

[4]PAPOULIS A. Signal analysis[M]. [S.l.]: Mcgraw-Hill， 1977.

[5]TIKHONOV A H， GONCHARSKY V. ll- posed problems in the Natural Sciences[M].Moscow: MIR， 1987.

[6]JOHNSON N L， KOTZ S. Distributions in statistics: 1 and 2[M]. New York: Wiley，1970.

[7]RAO C R. Linear statistical inference and its applications[M]. New York: Wiley， 1973.

[8]程文波，王華軍.信號稀疏表示的研究及應用[J].西南石油大學學報:自然科學版，2008，30(5):148-151.

[9]張賢達.現代信號處理[M].北京:清華大學出版社，1995.

[10]徐盛.現代數字信號處理[M].北京:機械工業出版社，2005.