經(jīng)典功率譜估計及其仿真
2008-04-12 00:00:00宋寧關(guān)華
現(xiàn)代電子技術(shù) 2008年11期
摘 要：介紹了各種經(jīng)典功率譜估計方法，不僅從理論上對各種方法的譜估計質(zhì)量進行了分析比較，而且通過Matlab實驗仿真驗證了理論分析的正確性。著重對使用比較廣泛的Welch法進行了深入的研究，給出了窗函數(shù)選擇的一般要求，通過仿真分析了不同的窗函數(shù)對Welch法譜估計質(zhì)量的影響，比較了他們的優(yōu)缺點。最后分析了采樣點數(shù)較少即短數(shù)據(jù)對Welch法譜估計質(zhì)量的影響。
關(guān)鍵詞：經(jīng)典譜估計；估計質(zhì)量；Welch法；窗函數(shù)；短數(shù)據(jù)
中圖分類號：TN911.7 文獻標識碼：B
文章編號：1004-373X(2008)11-159-03
Classical Power Spectrum Density Estimation and Its Simulation
SONG Ning，GUAN Hua
(School of Information Electric Engineering，Shandong Jianzhu University，Jinan，250101，China)
Abstract：Various classical Power Spectrum Density (PSD) estimation methods are introduced，estimation quality of each method is analyzed and compared in both theory and simulation using the software Matlab.Then further study is made in Welch method which is used most widely.General selecting criterion of window function is presented and estimation quality of Welch method using different window function is compared.Finally，the impact of fewer data on estimation quality of Welch method is analyzed.
Keywords：classical PSD estimation;estimation quality;Welch method;window function;fewer data
1 引 言
信號的頻譜分析是研究信號特性的重要手段之一，對于確定性信號，可以用Fourier變換來考察其頻譜性質(zhì)，而對于廣義平穩(wěn)隨機信號，由于它一般既不是周期的，又不滿足平方可積，嚴格來說不能進行Fourier變換，通常是求其功率譜來進行頻譜分析。
功率譜反映了隨機信號各頻率成份功率能量的分布情況，可以揭示信號中隱含的周期性及靠得很近的譜峰等有用信息，應(yīng)用極其廣泛，例如，在語音信號識別、雷達雜波分析、地震勘探信號處理、水聲信號處理、系統(tǒng)辨識中非線性系統(tǒng)識別、物理光學(xué)中透鏡干涉、流體力學(xué)的內(nèi)波分析、太陽黑子活動周期研究等許多領(lǐng)域，發(fā)揮了重要作用。
然而，實際應(yīng)用中的平穩(wěn)隨機信號通常是有限長的，只能根據(jù)有限長信號估計原信號的真實功率譜，這就是功率譜估計問題。功率譜估計分為經(jīng)典譜估計和現(xiàn)代譜估計，本文主要研究經(jīng)典譜估計。經(jīng)典譜估計又稱非參數(shù)模型譜估計，主要方法有：周期圖法，相關(guān)圖法及改進的周期圖估計法。
2 周期圖法
Schuster于1899年首先提出周期圖法，也稱直接法，取平穩(wěn)隨機信號X(n)的有限個觀察值x(0)，x(1)，…，x(N－1)對功率譜S(ω)進行估計：
(ω)=1N|XN(ejω)|2=1N|∑N－1n=0x(n)e－jωn|2
(1)
主要性能指標有：
(1) 估計的均值:
E［[AKS^](ω)］=12πS(ω)*W(ω)
W(ω)是窗函數(shù)ω(n)的Fourier變換。當(dāng)N→∞時，E［[AKS^](ω)］→S(ω)，是無偏估計；N是有限值時，是有偏估計，偏差為：
bia([AKS^](ω))=12πS(ω)W(ω－λ)dλ－S(ω)
(2) 估計的方差:
var([AKS^](ω))=|12πN∫π－πS(λ)D0(ω－λ)D0(ω+λ)dλ|2+
E2([AKS^](ω))
其中D0(ω)是矩形窗d0(n)=1，|n|≤|N－1|的Fourier變換，可見[AKS^](ω)不是S(ω)的一致估計；
隨著N的增大，譜估計起伏增大，N→∞時，var([AKS^](ω))S2(ω)。
(3) 估計的分辨率：數(shù)據(jù)窗為長度為N的矩形窗時，Re s{S(ω)}=0.892πN，N增大時，分辨率提高，但會使[AKS^]x(ω)的起伏加劇，可見方差與分辨率是一對矛盾。
周期圖法應(yīng)用比較廣泛，主要是由于它與序列的頻譜有對應(yīng)關(guān)系，可以采用FFT快速算法來計算。但是，這種方法需要對無限長的平穩(wěn)序列進行截斷，相當(dāng)于對其加矩形窗，使之成為有限長數(shù)據(jù)。同時，這也意味著對自相關(guān)函數(shù)加三角窗，使功率譜與窗函數(shù)卷積，從而產(chǎn)生頻譜泄漏，容易使弱信號的主瓣被強信號的旁瓣所淹沒，造成頻譜的模糊和失真，使得譜分辨率較低。
3 相關(guān)圖法
相關(guān)圖法是1958年由Blackman與Tukey首先提出的，理論基礎(chǔ)是維納-辛欽定理，基本思想是通過改善對相關(guān)函數(shù)的估計方法，對周期圖進行平滑處理以改善周期圖譜估計的方差性能，亦稱BT法、間接法。首先，由信號的有限個觀察值x(0)，x(1)，…，x(N－1)估計出自相關(guān)函數(shù)[AKR^](m)=1N∑N－1－|m|n=0x(n)x(n+m)，然后在(－M，M)內(nèi)對[AKR^](m)作Fourier變換，得到功率譜[AKS^](ω)=∑Mm=－M[AKR^](m)ω(m)e－jωn，ω(n)為窗函數(shù)。
當(dāng)M=N－1時，BT法與周期圖法估計出的功率譜是一樣的；當(dāng)M BT法的缺點在于當(dāng)M→N時，[AKR^](m)的方差很大，使譜估計質(zhì)量下降；由[AKR^](m)得到的[AKS^](ω)不一定為正值，從而可能失去功率譜的物理意義。 4 經(jīng)典譜估計方法的改進 4.1 經(jīng)典譜估計性能分析 平穩(wěn)隨機信號X(n)的自相關(guān)函數(shù)為： R(n)=E［X(n)X(n+m)］ (2) 功率譜為： S(ω)=∑∞m=－∞R(m)e－jωm (3) 然而，實際應(yīng)用中求不出總集意義上的自相關(guān)函數(shù)和功率譜，因而假定X(n)是各態(tài)歷經(jīng)的，x(n)是他的一個樣本，則： R(m)=limN→∞12N+1∑Nn=－Nx(n)x(n+m) (4) S(ω)=limN→∞E{12N+1|∑Nn=－Nx(n)e－jωn|2} (5) 由以上兩式可以看到，盡管自相關(guān)函數(shù)可以用時間平均來代替總集平均，但功率譜必須用總集平均來計算。這是因為隨機過程X(n)的每一次實現(xiàn)x(n)，其Fourier變換仍然是一個隨機過程，在每個頻率ω處都是一個隨機變量，因此，求數(shù)學(xué)期望是必要的，也就是說，對R(m)作Fourier變換后，S(ω)不具有各態(tài)歷經(jīng)性，所以，真實譜S(ω)應(yīng)在總集意義上求。另外，如果沒有求數(shù)學(xué)期望運算，式(4)的求極限運算在任何統(tǒng)計意義上都不會收斂。而式(1)與式(5)相比，既無求均值計算，也無求極限運算，只能看作是對真實譜S(ω)作均值運算時的一個樣本，缺少了統(tǒng)計平均，當(dāng)然也就產(chǎn)生了方差。 為了改進周期圖法的估計性能，常用的方法有兩種：一是平滑，就是用適當(dāng)?shù)拇昂瘮?shù)對譜進行平滑；二是平均，就是對同一信號做多次周期圖估計后再平均，在一定程度上彌補上述所缺的求均值運算。平滑的方法在BT譜估計中已作了介紹，下面介紹兩種平均方法。 4.2 Bartlett法 Bartlett譜估計法的主要目標是改善方差，基本原理是：將長度為N的數(shù)據(jù)分為L段，每段長度為M，先對每段數(shù)據(jù)用周期圖法進行譜估計，然后對L段求平均得到長度為N的數(shù)據(jù)的功率譜。由上述原理可得功率譜為： [AKS^](ω)=1L∑Li=1i^(ω)=1ML∑Li=1|∑M－1n=0xim(n)e－jωn|2 (6) (1) 估計的均值：E［[AKS^](ω)］=12πS(ω)*W(ω)，和周期圖法一樣。 (2) 估計的方差：當(dāng)N→∞時，var([AKS^](ω))1LS2(ω)，是周期圖法的1/L。 (3) 估計的分辨率： Re s{S(ω)}=0.892πM=0.89L2πN 由于MN，故其分辨率比周期圖法低，可見，Bartlett法方差的改善是以犧牲分辨率為代價的。 4.3 Welch法 Welch譜估計法是對Bartlett法的改進，旨在保持Bartlett法方差性能的同時，改善其分辨率，又稱加權(quán)交疊平均法，其基本原理是：對數(shù)據(jù)分段時，使每一段有部分重疊，然后對每一段數(shù)據(jù)用一個合適的窗函數(shù)進行平滑處理，最后對各段譜求平均。由上述原理可得功率譜為： [AKS^](ω)=1MUL∑Li=1|∑M－1n=0xim(n)e－jωn|2 (7) 其中U=∑nω(n)，ω(n)是窗函數(shù)。 (1) 估計的均值: E［[AKS^](ω)］=12πMUS(ω)*W(ω)，W(ω)是ω(n)的Fourier變換。 (2) 估計的方差： 當(dāng)N→∞時，var([AKS^](ω))9L16NS2(ω)。 (3) 估計的分辨率：當(dāng)窗函數(shù)為矩形窗、重疊50%時，Re s{S(ω)}=1.282πM。 因為Welch法允許各段數(shù)據(jù)交疊，所以數(shù)據(jù)段數(shù)L會增加，使方差得到更大的改善，但是數(shù)據(jù)的交疊有減小了每一段數(shù)據(jù)的不相關(guān)性，使方差的減小不會達到理論程度。另外，采用合適的窗函數(shù)可以減小信號的頻譜泄漏，同時也可以增加譜峰的寬度，從而提高譜分辨率。 5 Matlab仿真 取平穩(wěn)隨機信號為： x(n)=2sin(2π#8226;200#8226;n+ω1)+ 3cos(2π#8226;300#8226;n+ω2)+e(n) 其中ω1，ω2為均勻分布的隨機相位，e(n)為零均值的白噪聲，分別用周期圖法、BT法、Bartlett法、Welch法對其進行功率譜估計。 5.1 仿真結(jié)果 在各種譜估計方法中，采樣頻率1 000 Hz，采樣點數(shù)1 000，F(xiàn)FT點數(shù)512，Bartlett法中矩形窗長度為100，Welch法中Hanning窗長度為255，重疊數(shù)據(jù)為128。仿真結(jié)果如圖1所示。 圖1 仿真結(jié)果(一) 5.2 仿真結(jié)果分析 由圖1可以直觀地看出： (1) 周期圖法和BT法的特點是離散性大，曲線粗糙，方差較大，但是分辨率較高； (2) Bartlett法和Welch法的收斂性較好，曲線平滑，方差較小，但是功率譜主瓣較寬，分辨率低，這是由于對隨機序列加窗截斷所引起的Gibbs效應(yīng)造成的； (3) 與Bartlett法相比，Welch法的估計曲線比較粗糙，但是分辨率較好，原因是Welch法中對數(shù)據(jù)進行截斷時加的是Hanning窗，而在Bartlett法中使用的是矩形窗，相對于矩形窗，Hanning窗的主瓣包含更多的能量，因而使功率譜的主瓣較窄，分辨率較高。 由上述理論分析及仿真實驗可知，Welch法采用加窗交疊求功率譜，可以有效減小方差和偏差，一般情況下能接近一致估計的要求，因而得到廣泛應(yīng)用。同時還可以發(fā)現(xiàn)，對信號加不同的窗函數(shù)，譜估計的質(zhì)量是不同的，下面就通過進一步的仿真來分析研究Welch法中使用不同的窗函數(shù)時譜估計的質(zhì)量。 5.3 不同窗函數(shù)的Welch譜估計 在選擇窗函數(shù)時，一般有如下要求： (1) 窗口寬度M要遠小于樣本序列長度N，以排除不可靠的自相關(guān)值； (2) 當(dāng)平穩(wěn)信號為實過程時，為保證平滑周期圖和真是功率譜也是實偶函數(shù)，平滑窗函數(shù)必須是實偶對稱的； (3) 平滑窗函數(shù)應(yīng)當(dāng)在m=0出游峰值，并且m隨絕對值增加而單調(diào)下降，使可靠的自相關(guān)值有較大的權(quán)值； (4) 功率譜是頻率的非負函數(shù)且周期圖是非負的，因而要求窗函數(shù)的Fourier變換是非負的。 仍以上述平穩(wěn)隨機信號為例，采樣頻率、采樣點數(shù)、FFT點數(shù)、窗長度及重疊數(shù)據(jù)不變，窗函數(shù)采用矩形窗、Blackman窗、Hamming窗，仿真結(jié)果如圖2所示。 圖2 仿真結(jié)果(二) 由圖2可以看出，使用不同的窗函數(shù)譜估計的質(zhì)量是不一樣的，矩形窗的主瓣較窄，分辨率較好，但方差較大，噪聲水平較高；而Blackman窗和Hamming窗的主瓣較寬，分辨率較低，但方差較小，噪聲水平較低。因此，在進行譜分析時選擇何種窗函數(shù)，要視具體情況而定。如果強調(diào)高分辨率，能精確讀出主瓣頻率，而不關(guān)心幅度的精度，例如測量震動物體的自震頻率，可以選用主瓣寬度比較窄的矩形窗；對受到強干擾的窄帶信號，若干擾靠近信號，則可選用旁瓣幅度較小的窗函數(shù)，若離開通帶較遠，則可選用漸近線衰減速度比較快的窗函數(shù)?？傊槍Σ煌男盘柡筒煌奶幚砟康膩磉x擇合適的窗函數(shù)，這樣才能得到良好的效果。 5.4 短數(shù)據(jù)的Welch譜估計 以上研究了不同窗函數(shù)的Welch譜估計質(zhì)量，下面研究采樣數(shù)據(jù)較少時Welch譜估計的質(zhì)量。仍以上述平穩(wěn)隨機信號為例，采樣頻率1 000，采樣點數(shù)100，F(xiàn)FT點數(shù)256，采用矩形窗、Blackman窗、Hamming窗，窗長度50，重疊數(shù)據(jù)25，仿真結(jié)果如圖3所示。 圖3 仿真結(jié)果(三) 比較圖2和圖3可知，在短數(shù)據(jù)情況下，得到的信號功率譜的譜峰不夠尖銳，若信號的功率更小，很容易被淹沒在噪聲中而無法識別。 6 結(jié) 語 在經(jīng)典譜估計中，無論是周期圖法還是其改進的方法，都存在著頻率分辨率低、方差性能不好的問題，原因是譜估計時需要對數(shù)據(jù)加窗截斷，用有限個數(shù)據(jù)或其自相關(guān)函數(shù)來估計無限個數(shù)據(jù)的功率譜，這其實是假定了窗以外的數(shù)據(jù)或自相關(guān)函數(shù)全為零，這種假定是不符合實際的，正是由于這些不符合實際的假設(shè)造成了經(jīng)典譜估計分辨率較差。另外，經(jīng)典譜估計的功率譜定義中既無求均值運算又無求極限運算，因而使得譜估計的方差性能較差，當(dāng)數(shù)據(jù)很短時，這個問題更為突出。如何選取最佳窗函數(shù)、提高頻譜分辨率，如何在短數(shù)據(jù)情況下提高信號譜估計質(zhì)量，還需要進一步研究。 參 考 文 獻 ［1］金連文，韋崗.現(xiàn)代數(shù)字信號處理簡明教程［M］.北京：清華大學(xué)出版社，2004. ［2］飛思科技產(chǎn)品研發(fā)中心.Matlab7輔助信號處理技術(shù)與應(yīng)用［M］.北京：電子工業(yè)出版社，2005. ［3］王曉峰，王炳和.周期圖及其改進方法中譜分析率的Matlab分析［J］.武警工程學(xué)院學(xué)報，2003(6):75-77. ［4］姚武川，姚天任.經(jīng)典譜估計方法的Matlab分析［J］.華中理工大學(xué)學(xué)報，2000，28(4)：45-47. 作者簡介 宋 寧 男，1983年出生，碩士研究生。研究方向為自動化裝置的集成化與智能化。 關(guān) 華 男，1963年出生，碩士生導(dǎo)師，博士，副教授。主要從事信號處理、圖像處理、DSP等方面的研究。 注：本文中所涉及到的圖表、注解、公式等內(nèi)容請以PDF格式閱讀原文。
