功率譜估計性能分析及Matlab仿真

1、功率譜估計性能分析及Matlab仿真1 引言隨機信號在時域上是無限長的，在測量樣本上也是無窮多的，因此隨機信號的能量是無限的，應該用功率信號來描述。然而，功率信號不滿足傅里葉變換的狄里克雷絕對可積的條件，因此嚴格意義上隨機信號的傅里葉變換是不存在的。因此，要實現(xiàn)隨機信號的頻域分析，不能簡單從頻譜的概念出發(fā)進行研究，而是功率譜1。信號的功率譜密度描述隨機信號的功率在頻域隨頻率的分布。利用給定的個樣本數(shù)據(jù)估計一個平穩(wěn)隨機信號的功率譜密度叫做譜估計。譜估計方法分為兩大類：經(jīng)典譜估計和現(xiàn)代譜估計。經(jīng)典功率譜估計如周期圖法、自相關(guān)法等，其主要缺陷是描述功率譜波動的數(shù)字特征方差性能較差，頻率分辨率低。方差

2、性能差的原因是無法獲得按功率譜密度定義中求均值和求極限的運算2。分辨率低的原因是在周期圖法中，假定延遲窗以外的自相關(guān)函數(shù)全為0。這是不符合實際情況的，因而產(chǎn)生了較差的頻率分辨率。而現(xiàn)代譜估計的目標都是旨在改善譜估計的分辨率，如自相關(guān)法和Burg法等。2 經(jīng)典功率譜估計經(jīng)典功率譜估計是截取較長的數(shù)據(jù)鏈中的一段作為工作區(qū)，而工作區(qū)之外的數(shù)據(jù)假設(shè)為0，這樣就相當將數(shù)據(jù)加一窗函數(shù)，根據(jù)截取的個樣本數(shù)據(jù)估計出其功率譜1。2.1 周期圖法( Periodogram )Schuster首先提出周期圖法。周期圖法是根據(jù)各態(tài)歷經(jīng)的隨機過程功率譜的定義進行的譜估計。取平穩(wěn)隨機信號的有限個觀察值，求出其傅里葉變換然

3、后進行譜估計周期圖法應用比較廣泛，主要是由于它與序列的頻譜有直接的對應關(guān)系，并且可以采用FFT快速算法來計算。但是，這種方法需要對無限長的平穩(wěn)隨機序列進行截斷，相當于對其加矩形窗，使之成為有限長數(shù)據(jù)。同時，這也意味著對自相關(guān)函數(shù)加三角窗，使功率譜與窗函數(shù)卷積，從而產(chǎn)生頻譜泄露，容易使弱信號的主瓣被強信號的旁瓣所淹沒，造成頻譜的模糊和失真，使得譜分辨率較低1。該方法基于Matlab實現(xiàn)的程序：clear all;load test x;N=4096;Fn=-0.5:1/N:0.5-1/N;px=fft(x,N);pmax=max(px);%歸一化px=px/pmax;px=10*log10(px

4、+0.000001);plot(Fn,fftshift(px);grid on;圖1 周期圖法 圖2 周期圖法 說明：(1) 本報告仿真中所采用的用于功率譜估計的數(shù)據(jù)文件來自參考文獻3的test.dat。該數(shù)據(jù)為128點復序列（圖3），由復數(shù)噪聲加上四個復正弦組成。其歸一化頻率分別是：。圖3 復序列 (2) 從仿真圖可以清晰看到，和不能完全分開，僅在波形的頂部能看出是兩個頻率分量；此外，當數(shù)據(jù)長度太大時（圖1），譜曲線呈現(xiàn)較大的起伏；當數(shù)據(jù)長度太小時（圖2），譜的分辨率又不好。據(jù)此，周期圖法不滿足一致性估計條件。2.2 自相關(guān)法( BT法)自相關(guān)法的理論基礎(chǔ)是維納辛欽定理。1958年Black

5、man和Tukey給出了這一方法的具體實現(xiàn)。對于平穩(wěn)隨機信號來說，其自相關(guān)函數(shù)是確定性函數(shù)，故其功率譜也是確定的。這樣可由平穩(wěn)隨機離散信號的有限個離散值求出自相關(guān)函數(shù)然后在內(nèi)對做傅里葉變換，得到功率譜該方法基于Matlab實現(xiàn)的程序：clear all;load test x;N=4096;Fn=-0.5:1/N:0.5-1/N;Mlag=64;rx=xcorr(x,Mlag,'unbiased');px=fft(rx,N);pmax=max(px);%歸一化px=px/pmax;px=10*log10(px+0.000001);plot(Fn,fftshift(px);gri

6、d on;圖4 自相關(guān)法不加窗 圖5 自相關(guān)法不加窗 圖6 自相關(guān)法使用漢明窗( Hamming )說明：(1) 該方法先由序列估計出自相關(guān)函數(shù)，然后對進行傅里葉變換，便得到的功率譜估計。當延遲與數(shù)據(jù)長度之比很小時，可以有良好的估計精度。(2) 圖4是用自相關(guān)法（BT法）求出的功率譜，沒有加窗；圖5也是用自相關(guān)法（BT法）求出的功率譜，沒有加窗；圖6同樣是采用自相關(guān)法求出的功率譜，使用了漢明窗。顯然，自相關(guān)函數(shù)的延遲越小，譜變得越平滑。2.3 Welch法該方法的基本原理是在對隨機序列分段時，使每一段有部分重疊，然后對每一段數(shù)據(jù)用一個合適的窗函數(shù)進行平滑處理，最后對各段譜求平均。這樣可得功率譜

7、其中 這里為窗函數(shù)。該方法基于Matlab實現(xiàn)的程序：clear all;load test x;N=4096;Fn=-0.5:1/N:0.5-1/N;xpsd=pwelch(x,hamming(33),16,N,'whole');mmax=max(xpsd);%歸一化xpsd=xpsd/mmax;xpsd=10*log10(xpsd+0.000001);plot(Fn,fftshift(xpsd);grid on;圖7 Welch法 不疊合 使用漢明窗( Hamming )圖8 Welch法 疊合16點 使用漢明窗( Hamming )圖9 Welch法 疊合16點 使用矩形

8、窗( Boxcar )圖10 Welch法 疊合16點 使用布萊克曼窗( Blackman )說明：(1) 因為Welch法允許各段數(shù)據(jù)交疊，所以數(shù)據(jù)段數(shù)會增加，使方差得到更大的改善，但是數(shù)據(jù)的交疊減小了每一段數(shù)據(jù)的不相關(guān)性，使方差的減小不會達到理論程度。另外，采用合適的窗函數(shù)可以減少信號的頻譜泄露，同時也可以增加譜峰的寬度，從而提高分辨率。(2) 圖7是利用Welch法求出的周期圖，共分四段，每段32點，沒有疊合，使用了漢明窗；圖8也是利用Welch法求出的周期圖，共分四段，每段32點，使用了漢明窗；圖9是利用Welch法求出的周期圖，共分四段，每段32點，交疊數(shù)為16，且使用了矩形窗；圖1

9、0是利用Welch法求出的周期圖，共分四段，每段32點，交疊數(shù)為16，使用了布萊克曼窗。從圖中可以看出，由矩形窗處理的譜估計的主瓣寬度最窄，分辨率最好，但是其起伏性較大，所以其方差特性最差。由漢明窗和布萊克曼窗得到的譜估計的主瓣寬度最寬，因此其分辨率相對較差，但其旁瓣較小，大大改善了由矩形窗處理的譜估計旁瓣較大所產(chǎn)生的譜失真。因此，選擇不同的窗函數(shù)其主瓣寬度不一樣，造成譜估計的分辨率也不相同。2.4 經(jīng)典功率譜估計的性能比較由以上的Matlab仿真圖形和相關(guān)結(jié)果分析，我們得到了經(jīng)典譜估計算法性能的直觀比較：(1) 周期圖法得到的功率譜分辨率最高，但是方差性能最差，功率譜起伏劇烈，容易出現(xiàn)虛假譜

10、峰。(2) 自相關(guān)法（BT法）由于使用了平滑窗對周期圖法估計的功率譜進行了平滑，因此方差性能較好，功率譜比周期圖法估計的要平滑，但其分辨率比周期圖法低。(3) Welch平均周期圖法是三種經(jīng)典功率譜估計方法中方差性能最好的，估計的功率譜也最為平滑，但這是以分辨率的下降及偏差的增大為代價的。綜合上述討論，我們可以對經(jīng)典譜估計的算法作大致的總結(jié)3：(1) 功率譜估計，不論是直接法還是間接法都可以用FFT快速計算，且物理概念明確，因而仍是目前較常用的譜估計方法。(2) 譜的分辨率較低，它正比于，是所使用的數(shù)據(jù)長度。(3) 方差性能不好，不是真實功率譜的一致估計，且增大時，功率譜起伏加劇。(4) 周期

11、圖的平滑和平均是和窗函數(shù)的使用緊密關(guān)聯(lián)的，平滑和平均主要是用來改善周期圖的方差性能，但往往又減小了分辨率且增加了偏差，沒有一個窗函數(shù)能使估計的功率譜在方差、偏差和分辨率各個方面都得到改善，因此使用窗函數(shù)只是改進估計質(zhì)量的一個技巧問題，并不能從根本上解決問題。3 現(xiàn)代功率譜估計由前一章的討論我們可知，經(jīng)典功率譜估計方法的方差性能較差，分辨率較低。而現(xiàn)代譜估計技術(shù)的目標都是旨在努力改善譜估計的分辨率。參數(shù)模型法是現(xiàn)代譜估計的主要內(nèi)容，參數(shù)模型主要分為AR模型、MA模型和ARMA模型。由于AR模型具有一系列好的性能，因此是被研究最多并獲得廣泛應用的一種模型。本報告中現(xiàn)代功率譜估計的仿真基于的是AR模

12、型。3.1 自相關(guān)法假定觀察到得數(shù)據(jù)為，而對于無法觀察到得區(qū)間（即），的樣本假定為0，觀測數(shù)據(jù)區(qū)間之外的數(shù)據(jù)為0，在均方誤差意義下使得數(shù)據(jù)的預測誤差功率最小。由于自相關(guān)矩陣是Toeplitz矩陣，而且又為正定的，故可利用Levinson-Durbin遞歸算法高效求解，得到AR模型參數(shù)。該方法基于Matlab實現(xiàn)的程序：clear all;load test x;N=4096;fn=-0.5:1/N:0.5-1/N;xpsd=pyulear(x,20,N);pmax=max(xpsd);xpsd=xpsd/pmax;xpsd=10*log10(xpsd+0.000001);plot(fn,fft

13、shift(xpsd);grid on;圖11 自相關(guān)法 圖12 自相關(guān)法 圖13 自相關(guān)法 說明：(1) 圖11、12和13是用自相關(guān)法求出的AR譜曲線，階次分別等于10，20和30?？梢钥闯觯陔A次較低時（圖11），分辨率和檢測能力均不好。當時，和處的兩個正弦剛剛可以分開，在和處的兩個正弦也可以檢出。因此必須通過提高階次來達到分辨出間隔較小的頻率點的效果。(2) AR模型的自相關(guān)法等效于對前向預測的誤差序列前后加窗，加窗的結(jié)果是使得自相關(guān)法的分辨率降低。數(shù)據(jù)越短，分辨率越不好。3.2 協(xié)方差法協(xié)方差法與自相關(guān)法的區(qū)別主要在于預測誤差功率求和式的上下限取得不同。由于協(xié)方差法對于觀察區(qū)間外樣本

14、并未假定為0，故預測誤差功率表達式中的總是落在觀察區(qū)間中，為此預測誤差功率的求和上下限必須在之間。但由此得到的自相關(guān)矩陣是對此的半正定矩陣，且不具有Toeplitz性質(zhì)，故不能采用Levinson-Durbin遞歸算法求解，因此得到的AR模型可能不穩(wěn)定。該方法基于Matlab實現(xiàn)的程序：clear all;load test x;N=4096;tn=-0.5:1/N:0.5-1/N;xpsd=pcov(x,10,N);pmax=max(xpsd);xpsd=xpsd/pmax;xpsd=10*log10(xpsd+0.000001);plot(tn,fftshift(xpsd);grid on

15、;圖14 協(xié)方差法 說明：(1) 可以看到譜圖在信號源頻率處：譜線狹窄突出，其他處譜線起伏較為劇烈。(2) 采用協(xié)方差法對信號進行建模，能夠較好地反映出信號真正的模型。3.3修正的協(xié)方差法AR譜估計的協(xié)方差算法基于的是最小化前向預測誤差。而修正的協(xié)方差算法基于的是最小化前向和后向預測誤差。這樣使得它的誤差功率的計算是在相對于協(xié)方差法多一倍的數(shù)據(jù)點上進行，這在觀察數(shù)據(jù)長度很短的情況下，是非常有利的，但這要求信號在正反兩個方向上呈現(xiàn)相同的特性。此外，由此得到的自相關(guān)矩陣不具有Toeplitz性質(zhì)，故其正則方程不能采用Levinson-Durbin遞歸算法求解。該方法基于Matlab實現(xiàn)的程序：cl

16、ear all;load test x;N=4096;fn=-0.5:1/N:0.5-1/N;xpsd=pmcov(x,10,N);pmax=max(xpsd);xpsd=xpsd/pmax;xpsd=10*log10(xpsd+0.000001);plot(fn,fftshift(xpsd);grid on;圖15 修正的協(xié)方差法 圖16 修正的協(xié)方差法說明：(1) 修正的協(xié)方差法較協(xié)方差法而言，譜估計圖大致相同，但前者在信號源頻率處的譜峰更加突出、尖銳，易于辨別。(2) 通過圖15和16的對比，我們可發(fā)現(xiàn)，在階次較高的情況下能得到非常滿意的結(jié)果。3.4 Burg法Burg算法是較早提出的建

17、立在數(shù)據(jù)基礎(chǔ)上的AR系數(shù)求解的有效算法。它基于最小化前向后向預測誤差的同時滿足Levinson-Durbin遞歸。對比其它的AR估計方法，Burg法避免了對自相關(guān)函數(shù)的計算，改而直接估計反射系數(shù)。對于短數(shù)據(jù)的估計，Burg法求出的AR功率譜密度估計非常逼近于真值。另外，它能確保產(chǎn)生一個穩(wěn)定的AR模型，并且能高效計算。Burg法由于具有上述優(yōu)點，所以分辨率比自相關(guān)法高，但對于混有白噪聲的正弦信號，有時可能會出現(xiàn)譜線分裂現(xiàn)象。該方法基于Matlab實現(xiàn)的程序：clear all;N=4096;fn=-0.5:1/N:0.5-1/N;xpsd=pburg(x,10,N);pmax=max(xpsd)

18、;xpsd=xpsd/pmax;xpsd=10*log10(xpsd+0.000001);plot(fn,fftshift(xpsd);grid on;圖17 Burg法 圖18 Burg法 說明：(1) 圖17和18是采用Burg算法對數(shù)據(jù)文件所做的功率譜估計，在階次的情況下，得到了非常滿意的效果。(2) 通過圖13和17的對比，我們可明顯看出Burg法的分辨率比自相關(guān)法高，在階次較低的情況下（圖17）也能較好的分辨出間隔小的頻率點。3.5 MUSIC法剛才所討論的自相關(guān)法、協(xié)方差法和Burg法都是基于參數(shù)建模的功率譜估計，而基于非參數(shù)建模的功率譜估計也是現(xiàn)代功率譜估計的重要內(nèi)容。該方法是基于自相關(guān)矩陣的特征分析或者特征值分解的功率譜估計，它將相關(guān)矩陣的特征向量空間分解為信號子空間和噪聲子空間，由此衍生出EV( Eigenvector )算法與MUSIC( Multiple Signal Classification )算法的信號功率譜估計4。其中EV譜估計與MUSIC算法譜估計都是基于噪聲子空間的功率譜估計。這類方法對線譜（正弦信號的譜）最合適，對檢測混有白噪聲的正弦信號很有效，特別是低信噪比的情況。MUSIC估計由下面方程給出此處是復正弦信號向量。該方法基于Matlab實現(xiàn)的程序：