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

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

3、觀察值，求出其傅里葉變換然后進(jìn)行譜估計(jì)周期圖法應(yīng)用比較廣泛，主要是由于它與序列的頻譜有直接的對(duì)應(yīng)關(guān)系，并且可以采用FFT快速算法來計(jì)算。但是，這種方法需要對(duì)無限長的平穩(wěn)隨機(jī)序列進(jìn)行截?cái)?，相?dāng)于對(duì)其加矩形窗，使之成為有限長數(shù)據(jù)。同時(shí)，這也意味著對(duì)自相關(guān)函數(shù)加三角窗，使功率譜與窗函數(shù)卷積，從而產(chǎn)生頻譜泄露，容易使弱信號(hào)的主瓣被強(qiáng)信號(hào)的旁瓣所淹沒，造成頻譜的模糊和失真，使得譜分辨率較低1。該方法基于Matlab實(shí)現(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;p

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

5、8年Blackman和Tukey給出了這一方法的具體實(shí)現(xiàn)。對(duì)于平穩(wěn)隨機(jī)信號(hào)來說，其自相關(guān)函數(shù)是確定性函數(shù)，故其功率譜也是確定的。這樣可由平穩(wěn)隨機(jī)離散信號(hào)的有限個(gè)離散值求出自相關(guān)函數(shù)然后在內(nèi)對(duì)做傅里葉變換，得到功率譜該方法基于Matlab實(shí)現(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.);plot(Fn,fftshift(px);gr

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

7、譜其中 這里為窗函數(shù)。該方法基于Matlab實(shí)現(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.);plot(Fn,fftshift(xpsd);grid on;圖7 Welch法 不疊合 使用漢明窗( Hamming )圖8 Welch法 疊合16點(diǎn) 使用漢明窗( Hamming )圖9 Welch法 疊合16點(diǎn) 使用矩形窗( Bo

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

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

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

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

12、 自相關(guān)法假定觀察到得數(shù)據(jù)為，而對(duì)于無法觀察到得區(qū)間（即），的樣本假定為0，觀測(cè)數(shù)據(jù)區(qū)間之外的數(shù)據(jù)為0，在均方誤差意義下使得數(shù)據(jù)的預(yù)測(cè)誤差功率最小。由于自相關(guān)矩陣是Toeplitz矩陣，而且又為正定的，故可利用Levinson-Durbin遞歸算法高效求解，得到AR模型參數(shù)。該方法基于Matlab實(shí)現(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.);plot(fn,fftshift(xpsd)

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

14、差功率表達(dá)式中的總是落在觀察區(qū)間中，為此預(yù)測(cè)誤差功率的求和上下限必須在之間。但由此得到的自相關(guān)矩陣是對(duì)此的半正定矩陣，且不具有Toeplitz性質(zhì)，故不能采用Levinson-Durbin遞歸算法求解，因此得到的AR模型可能不穩(wěn)定。該方法基于Matlab實(shí)現(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.);plot(tn,fftshift(xpsd);grid on;圖14 協(xié)方差法 說明：(1)

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

16、 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.);plot(fn,fftshift(xpsd);grid on;圖15 修正的協(xié)方差法 圖16 修正的協(xié)方差法說明：(1) 修正的協(xié)方差法較協(xié)方差法而言，譜估計(jì)圖大致相同，但前者在信號(hào)源頻率處的譜峰更加突出、尖銳，易于辨別。(2) 通過圖15和16的對(duì)比，我們可發(fā)現(xiàn)，在階次較高的情況下能得到非常滿意的結(jié)果。3.4 Burg法Burg算法是較早提出的建立在數(shù)據(jù)基礎(chǔ)上的AR系數(shù)求解的有效算法。它基于

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

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