摘 要:從經(jīng)典功率譜估計周期圖法原理入手，從理論上分析了其存在的局限性，借助Welch算法對其進行修正。依靠Matlab強大的數(shù)值分析和信號處理能力，進行實驗仿真，比較不同的窗函數(shù)，不同的數(shù)據(jù)長度對Welch法譜估計質(zhì)量的影響，并分析了造成這些影響的原因。

關鍵詞:功率譜估計;周期圖法;Welch算法;Matlab

中圖分類號:TP911 文獻標識碼:A

文章編號:1004-373X(2010)03-007-03

Matlab Simulation Analysis of Power Spectrum Estimation Based on Welch Method

YI Xin，QU Aihua

(Naval Command College，Nanjing，211800，China)

Abstract:The paper mainly introduces the principles of Periodogram method of classical PSD estimation，analyzes the deficiency of Periodogram method in theory，and makes use of Welch to amend Perodogram method.By the use of simulation in Matlab，the impacts of different window function and different lenghth of data on estimation quality of Welch are discussed and the reasons of the impacts are analyzed.

Keywords:power spectrum estimation;periodogram method;Welch method;Matlab

0 引 言

隨機信號在時間上是無限的，在樣本上是無窮多，因此隨機信號的能量是無限的，它應是功率信號。功率信號不滿足傅里葉變換的絕對可積條件，因此嚴格意義上其傅里葉變換是不存在的。因此，對隨機信號的頻域分析，不再是簡單的頻譜，而是功率譜。

功率譜估計(PSD)是利用有限長的數(shù)據(jù)估計信號的功率譜，它涉及信號系統(tǒng)、隨機信號分析、概率統(tǒng)計、隨機過程、矩陣代數(shù)等一系列的基礎學科，廣泛應用于雷達、聲納、通信、地質(zhì)勘探、天文、生物醫(yī)學工程等眾多領域[1]。功率譜估計可以分為經(jīng)典譜估計和現(xiàn)代譜估計[2]，本文主要研究經(jīng)典譜估計中的Welch算法。目前經(jīng)典譜估計中主要采用兩種方法進行估計，分別是周期圖法和自相關法[3]，而Welch算法就是在周期圖法的基礎上進行改進得到的。

1 周期圖法

Schuster于1899年首先提出了周期圖這一概念，因為它是直接由傅里葉變換得到的，又稱之為直接法。該方法是把隨機信號x(n)的N點觀察數(shù)據(jù)xN(n)視為一能量有限信號，直接取xN(n)的傅里葉變換，得xN(ω)，然后再取其幅值的平方，并除以N，作為對x(n)真實的功率譜P(ejω)的估計，即:

(ω)=1NXN(ω)2=1N∑N-1n=0x(n)e-jωn2

(1)

由此可見，式(1)在n→∞時，既不存在均值，也不存在極限，它只能看作是對真實譜做均值運算時的一個樣本。缺少了統(tǒng)計平均，在記錄的信號序列長度一定的條件下，要保證足夠高的譜分辨率，譜估計的方差會很大，譜的正確性會很差。因此周期圖的方差性能不好，而且，當數(shù)據(jù)長度N太大時，譜曲線呈現(xiàn)較大的起伏;當數(shù)據(jù)長度N太小時，譜的分辨率又不好。據(jù)此，直接法進行譜估計不滿足一致性估計條件。因此，必須對周期圖進行改進。

2 直接法譜估計的改進

為了改進周期圖法的估計性能，常用的方法有兩種:一是平均，就是對同一信號做多次周期圖估計后再平均，在一定程度上彌補上述所缺的求均值運算。該方法事實上就是經(jīng)典譜估計的間接法，本文不做討論。二是平滑，就是用適當?shù)拇昂瘮?shù)對譜進行平滑。其思想是把一長度為N的數(shù)據(jù)xN(n)分成L段，每段的長度為M，分別求每一段的功率譜，然后加以平均，以達到所希望的目的。

若對分段的數(shù)據(jù)加矩形窗[1，4]，則第i段的數(shù)據(jù)變?yōu)?

xiN(n)=xN[n+(i-1)M]d1[n+(i-1)M]，

0≤n≤M-1，1≤i≤L

(2)

由此得到修正后的周期圖，即平均周期圖:

P(ω)=1ML∑Li=1∑M-1n=0xiNe-jωn2

(3)

此方法就是Bartlett法，它很好地改善了直接法的方差特性，但是它是以犧牲偏差和分辨率為代價的。

Welch法是對Bartlett法的改進。主要改進在兩個方面:一是在對xN(n)分段時，允許每段數(shù)據(jù)存在部分的交疊;二是每一段的數(shù)據(jù)窗口可以不是矩形窗口。這樣可以改善由于矩形窗所造成的分辨率較差的影響。然后按照Bartlett法求每一段的功率譜，并對結(jié)果進行歸一化，從而得到進一步修正的周期圖，即:

(ω)=1MU∑Li=1∑M-1n=0xiNd2(n)e-jωn2

(4)

式中:U為歸一化因子;d2(n)是數(shù)據(jù)窗口。

因為Welch算法各段允許交疊，從而增大了段數(shù)L，這樣可以更好地改善方差特性。但是，數(shù)據(jù)的交疊又減小了每一段的不相關性，使方差的減小不會達到理論計算的程度。另外，選擇合適的窗函數(shù)可以減小頻譜的泄漏，改善分辨率。

3 Matlab仿真

為了分析Welch算法的性能，利用Matlab中提供的pwelch[5]函數(shù)實現(xiàn)Welch算法的功率譜估計。采樣頻率為1 000，采樣點數(shù)為1 000，F(xiàn)FT點數(shù)為256，交疊數(shù)為20，窗函數(shù)采用矩形窗、海明窗和blackman窗[6-8]，仿真程序如下:

Fs=1000;

n=0:1/Fs:1;

xn=2sin(2*pi*100*n)+4*cos(2*pi*150*n)+randn(size(n));

nfft=256;

[Pxx，f]=pwelch(xn，rectwin(100)，20，nfft，F(xiàn)s，′half′);%矩形窗

[Pxx1，f]=pwelch(xn，hamming(100)，20，nfft，F(xiàn)s，′half′);%海明窗

[Pxx2，f]=pwelch(xn，blackman(100)，20，nfft，F(xiàn)s，′half′);%blackman窗

xpsd=10*log10(Pxx);

xpsd1=10*log10(Pxx1);

xpsd2=10*log10(Pxx2);

figure(1)

plot(f，xpsd);

figure(2)

plot(f，xpsd1);

figure(3)

plot(f，xpsd2);

如圖1所示，由矩形窗處理的譜估計的主瓣寬度最窄，分辨率最好，但是其旁瓣比其他窗函數(shù)的旁瓣要高，因此其正弦譜線附近的旁瓣泄漏比較嚴重，而且其起伏性較大，所以其方差特性最差[9]。由blackman窗和海明窗處理譜估計的主瓣寬度最寬，因此其分辨率相對較差，但其旁瓣較小，大大改善了由矩形窗處理的譜估計旁瓣較大所產(chǎn)生的譜失真。究其原因，選擇不同的窗函數(shù)其主瓣寬度不一樣，造成譜估計的分辨率也不相同;另外，選擇不同的窗函數(shù)旁瓣的衰減速度也不相同，因而譜估計旁瓣的泄漏程度也不一樣。

圖1 長數(shù)據(jù)條件下不同窗函數(shù)的Welch算法譜估計

改變采樣點數(shù)為100，交疊數(shù)為10，窗長度為50，其他參數(shù)不變。由于保持采樣頻率，降低采樣點數(shù)，所以要對信號采取抽取，這里使用的是FIR濾波器，采用的是凱撒窗，其中對于矩形窗使用的是24階，海明窗和blackman窗采用的是50階。采樣點數(shù)的大大降低，造成采樣后的信號已經(jīng)和原信號存在較大的誤差，這時譜變得較平滑，但是各個窗的主瓣寬度變得更寬，造成分辨率降低。換句話說，減小數(shù)據(jù)長度，造成窗函數(shù)主瓣的寬度將變寬，從而減小譜的起伏。同時，由于信號不夠尖銳，在功率較小的情況下就會淹沒在噪聲中。較大的邊瓣掩蓋了功率譜中的較弱的成分，或者產(chǎn)生假的峰值，這些情況在圖2中已經(jīng)出現(xiàn)。

圖2 短數(shù)據(jù)長度條件下不同窗函數(shù)的Welch算法估計

4 結(jié) 語

綜上所述，盡管Welch算法在短數(shù)據(jù)情況下存在局限性，但是由于其原理簡單、實現(xiàn)容易，目前得到了廣

泛的應用[10]。在實際的應用中，Welch算法對周期圖

的平滑和平均是和窗函數(shù)的使用緊密相關聯(lián)的。平滑和平均主要是用來改善周期圖的方差性能，但往往又減小了分辨率和增大了偏差。不存在任何窗函數(shù)能夠綜合地改善估計譜的方差、偏差和分辨率等性能。因此，應該根據(jù)不同的信號、不同的處理目的合理地選擇窗函數(shù)。

參考文獻

[1]宋寧，關華.經(jīng)典功率譜估計及其仿真[J].現(xiàn)代電子技術(shù)，2008，31(11):159-161.

[2]馮磊.經(jīng)典功率譜估計與現(xiàn)代功率譜估計的對比[J].商業(yè)文化，2009(5):239-243.

[3]胡廣書.數(shù)字信號處理[M].北京:清華大學出版社，2007.

[4]姚武川，姚天任.經(jīng)典譜估計方法的Matlab分析[J].華中理工大學學報，2000，28(4):45-47.

[5]寧長春，陳天祿，索郎桑姆，等.數(shù)字信號處理中常用的Matlab工具箱函數(shù)簡介[J].西藏科技，2007(12):75-78.

[6]魏鑫，張平.周期圖法功率譜估計中的窗函數(shù)分析[J].現(xiàn)代電子技術(shù)，2005，28(3):14-15.

[7]邵玉斌.Matlab/Simulink通信系統(tǒng)建模與仿真實例分析[M].北京:清華大學出版社，2008.

[8]范瑜，鄔正義.功率譜估計的Welch方法中的窗函數(shù)研究[J].常熟高專學報，2000，14(7):36-39.

[9]瞿海雁，李鸝，錢小凌.如何在Matlab中優(yōu)化基本周期圖法對隨機信號進行的功率譜估計[J].首都師范大學學報:自然科學版，2006，27(5):33-36.

[10]羅敏，劉嵩.基于Welch算法的功率譜估計的實現(xiàn)[J].北京工商大學學報:自然科學版，2007，25(3):58-59.