頻率估計的時頻分析方法研究

郭 振，郭建濤

（1.信陽師范學院 計算機與信息技術學院，河南 信陽 464000；2.信陽師范學院 物理電子工程學院，河南 信陽 464000）

信號處理與分析中最重要也是最基本的兩個變量是時間和頻率。隨著現(xiàn)代電子技術的發(fā)展使其信號形式越來越復雜，信號的非線性、非平穩(wěn)時變特性已經(jīng)成為信號的主要特征，單純的時域或者頻域處理已經(jīng)不能滿足需要。采用時頻分析可以在時頻二維平面給出信號的時變頻譜，是解決上述問題的一種有效方法。

現(xiàn)代通信，特別是跳頻擴頻通信系統(tǒng)中，通過頻率和相位調(diào)制來展寬信號的頻譜。在雷達、聲納以及電子對抗中的信號探測系統(tǒng)中，系統(tǒng)和目標之間的相對移動，傳播媒質(zhì)的擾動都會使頻率發(fā)生改變。信號頻率的時變已經(jīng)成為其非平穩(wěn)性的主要表現(xiàn)之一，因此在很多實際信號的處理中，估計一個非平穩(wěn)過程的瞬時頻率是一項重要的工作。利用時頻分析研究瞬時頻率估計及各種方法的性能對國民經(jīng)濟的發(fā)展和國防建設有重要的意義。筆者闡述了瞬時頻率時頻分析理論研究和應用發(fā)展概況，給出了幾種典型的估計算法以及相應的性能評價。

1 瞬時頻率

瞬時頻率這個概念起源于通信中的調(diào)頻研究。Armstrong[1]發(fā)現(xiàn)進行頻率調(diào)制可以有效地壓制噪聲，從而引起人們對調(diào)頻信號中“頻率”這個物理量的興趣，促使人們?nèi)ド钊胙芯款l率調(diào)制和瞬時頻率的概念及其數(shù)學描述。

1）解析信號的相位差分

Gabor[2]提出了從實信號產(chǎn)生相應復信號的方法，這對瞬時頻率概念的進一步發(fā)展起了重要作用。實信號s（t）=a（t）cos（φ（t））的瞬時頻率定義式為：

其中，z（t）為實信號的解析信號，φ（t）為解析信號相位。

2）基于時頻分布矩（特定時間的頻率平均）：許多TFD例如WVD的一階矩給出瞬時頻率值，另外一些TFD例如STFT的一階矩給出瞬時頻率的近似值。定義式為：

應用時頻分布WVD的時間一階條件矩計算得到的瞬時頻率和應用解析信號有限差分法得到的瞬時頻率結(jié)果完全相同。估計過程存在噪聲時，估計性能下降，統(tǒng)計上次優(yōu)化，而且計算耗時，所以一階矩方法應用較少。

對于實際應用需要將定義式離散化，通常采用差分形式，且以中心有限差分最為常見。對于多分量信號，瞬時頻率的定義就失去了它的意義。這種情況下，一般根據(jù)應用的需要，分離各個分量，再分別求出每個分量的瞬時頻率。

瞬時頻率提供一種作為時間函數(shù)的頻域能量集中的表現(xiàn)形式，在信號識別、跟蹤、估計和建模方面有重要意義。針對瞬時頻率的時頻分析方法主要從估計精度、計算復雜性、信號包含噪聲的特點以及信號本身的多分量特性等方面進行研究。

2 估計方法

對于一個合理的時頻分布，其能量峰脊應出現(xiàn)在瞬時頻率的跡線附近，利用時頻分布可以從兩個角度估計瞬時頻率：基于時間一階條件矩的瞬時頻率估計和基于譜峰檢測的瞬時頻率估計。基于能量檢測方法進行頻率估計實際上是依據(jù)式（1）進行的。已經(jīng)知道這兩種定義對于WVD分布是等價的，但基于一階矩方法由于有較大的計算量，而沒有額外的性能優(yōu)勢，所以文獻中多采用能量檢測的瞬時頻率估計方法。

2.1 線性時頻方法

常見的線性時頻表示主要有短時傅里葉變換、Gabor變換以及小波變換。從考察信號的頻率成分隨時間的演化特性角度來說，小波變換的結(jié)果令人費解，因為它本質(zhì)上是一種時間—尺度的多分辨分析。而且一旦母小波選擇不當，應用效果大受影響。

短時傅里葉變換方法又稱加窗傅里葉變換算法，它的基本思想是假定非平穩(wěn)信號在分析窗函數(shù)的一個短的時間間隔內(nèi)是準平穩(wěn)，并移動分析窗函數(shù)，計算出各個不同時刻的功率譜。對頻譜進行峰值檢測即完成了對信號瞬時頻率的估計。

該時頻方法所加的時窗平滑了噪聲的影響，使STFT有較好的抗噪聲干擾的能力。同時該方法是線性，不會產(chǎn)生交叉項干擾，而且方法簡單，計算量小。主要缺陷是：對于一定的時刻，只是對附近窗口內(nèi)的信號作分析，存在所謂的“窗效應”，其時間和頻率分辨率受窗長的約束，不能同時優(yōu)化。通常將其作為粗略估計值，作為精確估計的基礎。

2.2 雙線性時頻表示

克服時頻分辨率受限的一種方法是引入雙線性時頻表示。這類分布多是魏格納分布的變形，可以用統(tǒng)一的形式表示，習慣稱之為Cohen類時頻分布。

2.2.1 Wigner-Ville分布（WVD）

WVD最早由Wigner[3]于1932年提出，Ville[4]于1948年把它應用到信號處理領域。WVD定義為：

它是雙線性乘積核或瞬時自相關函數(shù)關于時延的傅里葉變換。

對于常數(shù)振幅二次相位的線性調(diào)頻信號，WVD方法具有理想的時頻聚集性能，可以使瞬時頻率曲線和時頻能量峰脊重合，得到無偏的瞬時頻率估計。文獻[5]指出存在高斯白噪聲的情況下，線性調(diào)頻信號需要SNR在3 dB以上，估計性能才是優(yōu)化的，并隨著SNR的下降而緩慢下降；在3 dB以下，性能顯著變差。

對于其他頻率特征的信號，通過選擇適當?shù)拇昂瘮?shù)，使窗內(nèi)的頻率特征近似線性，也可得到無偏估計的效果，但是信號必須限制在較高的信噪比水平。當信號信噪比較低時，這種方法表現(xiàn)出較高的方差。由此可以看出，WVD適應于線性時頻信號，或者高信噪比的高階調(diào)頻信號；相應地研究WVD分布，需要克服強噪聲、非線性調(diào)頻以及算法固有的交叉項問題。這些問題的改善分散在下述不同的時頻分析中。

WVD計算量大，沒有合適的快速算法，離實時處理還有一定的距離。已有的WVD快速算法，一般采用快速傅里葉變換的方法，包括Boashash[6]等提出的重排核函數(shù)法和Martin[7]等人提出共軛對消法。

2.2.2 多項式魏格納分布

Boashash[8]等提出了多項式魏格納分布（PolynomialWigner-Ville Distribution，PWVD），目的在于解決多項式相位調(diào)頻情況下的頻率估計。PWVD定義：

其中 Kz（t，τ）為核函數(shù)：

在c=0.5，b=1，q=2的情況下，PWVD轉(zhuǎn)化為WVD。

非線性調(diào)頻問題轉(zhuǎn)化為確定上述核函數(shù)中系數(shù)使多項式核可以將二次或多次調(diào)頻信號轉(zhuǎn)換成正弦波，從而使PWVD對多項式調(diào)頻信號具有理想的時頻聚集性。但是對于多分量信號，PWVD的交叉項將變得更加復雜，將噪聲的影響在時頻平面彌散化，加劇了噪聲的干擾。因此，在信噪比較低時，基于PWVD的瞬時頻率估計的方差較高，即使采用優(yōu)化的核設計，與線性調(diào)頻情況下的方差也要高出4 dB以上。另外，PWVD在實現(xiàn)時，必須對原信號插值，多項式核相乘使計算變得復雜。

2.2.3 交叉魏格納分布

與WVD類似，可以定義兩個信號或者一個信號的多分量的交叉魏格納分布（Cross Wigner-Ville Distribution XWVD），即瞬時互相關函數(shù)關于時延的傅里葉變換，表達式為：

應用XWVD[9]估計瞬時頻率的基本原理是在每次迭代過程這中，應用基于XWVD的譜峰檢測估計出的瞬時頻率來產(chǎn)生新的參考信號，新的參考信號和待估計的信號生成一個新的XWVD，再從新生成的XWVD估計出瞬時頻率。通過迭代使信號能量聚集程度增加，從而使在噪聲中估計出瞬時頻率的概率增加，而且算法具有收斂性。選擇合適的幅度變化率，達到C-R界的估計方差，信噪比最小可以約為-2 dB。XWVD的定義不僅克服了WVD中交叉項問題，而且提高了頻率估計在低信噪比情況下的性能。缺點是不能處理多分量信號，而且是迭代算法。

2.2.4 偽 WVD（Ps.WVD）

WVD本質(zhì)上是一種二次型時頻分布，對多分量信號將出現(xiàn)交叉項，在時頻平面引入了模糊，這是WVD的一個重大缺陷，主要通過核函數(shù)的設計來加以抑制。其中偽Wigner-Ville分布通過與窗函數(shù)的卷積實現(xiàn)平滑操作，減少交叉項影響，變換結(jié)果出現(xiàn)負值。雖然采用平滑偽魏格納分布，可以得到正值結(jié)果，但是由于平滑作用使得時頻面出現(xiàn)模糊，影響時頻分辨率。另外，窗函數(shù)設計應用于STFT或者WVD，從本質(zhì)上存在一個前提，即信號在窗口內(nèi)滿足平穩(wěn)性的假設，對于快速時變信號，這個條件很難在任意時刻得到滿足。因此，該類分布在時頻估計方面應用較少，并出現(xiàn)了大量的基于自適應時頻表示的頻率估計方法。

2.3 自適應時頻表示

2.3.1 自適應窗口長度選擇

短時傅里葉變換和WVD變形中采用固定長度窗函數(shù)，很難在所有時刻滿足平穩(wěn)性的假設，從而引起頻率估計偏差或者方差的增大，相應于大的或小的窗口長度。為了解決這個問題，可以采用自適應于估計頻率的窗口長度。這種方法可以應用于SFTF或WVD的多種變形形式。基本方法一般是首先推導基于時頻頻譜分布進行頻率估計的偏差和方差的表達式，通過均方誤差的表達式在偏差和估計方差之間進行折中，從而得到自適應的窗口長度。

為解決STFT的“窗效應”問題，文獻[10]中提出了一種自適應窗長的STFT方法，該窗函數(shù)為含兩個控制函數(shù)的高斯窗函數(shù)，該方法對于瞬態(tài)信號和長時間信號分量有一定的自適應能力。文獻[11]中提出短時線性調(diào)頻窗函數(shù)的STFT對信號的瞬時頻率進行估計。文獻[12]利用自適應方法估計各種調(diào)頻信號，通過蒙特卡羅仿真，指出除了三角形式的調(diào)頻信號以外，包括階梯狀、三角函數(shù)形式的調(diào)頻，其估計的累積偏差和方差都明顯低于偽魏格納分布的結(jié)果。

2.3.2 自適應核函數(shù)方法

為了更好地在抑制交叉項的同時保持很好的信號時頻聚集性，對于具有不同時變頻率特性的信號，核函數(shù)也要根據(jù)每個信號不同的特性來選取，每一步的最佳參數(shù)通過基函數(shù)與信號殘余內(nèi)積的模的最大化確定。信號自適應擴展得到基函數(shù)后，再利用不同的時頻分布，利用能量峰值探測求較好的瞬時頻率估計。

S.Mallat和Z.Zhang[13]提出在更大的范圍內(nèi)尋找匹配原始信號的基函數(shù)，利用Gauss函數(shù)的尺度變換、時頻移位組成一個Gabor函數(shù)集。然后根據(jù)最大匹配投影原理尋找最佳基函數(shù)的線性組合，以達到自適應分解之目的。但是該迭代算法對時頻平面的劃分是一種格型分割（2j），當待分析的信號是chirp信號時，這種匹配相當于零階逼近，勢必會造成分解過程存在許多截斷和分量之間的混合畸變。為克服這一缺陷，L.Angrisani等[14]提出了“小調(diào)頻波變換”，簡稱 CT。該算法在核函數(shù)上乘以一個二次項的調(diào)頻信號，由兩個參數(shù)分布表示切率（chirprate）和時頻平面單元的曲率，并且采用高斯函數(shù)的開方形式使其具有歸一化能量。對于多項式相位的頻率估計誤差在1～2%以內(nèi)；多分量情況下也給出了類似的結(jié)果，但是STFT和偽WVD卻出現(xiàn)了嚴重的干擾，只有部分瞬時頻率可以估計出來。當待分析信號的瞬時頻率具有一定的周期性時，文獻[15]給出與CT類似地稱之為WT的時頻變換。在高斯母函數(shù)的基礎上，乘上一個用初始幅度、頻率、相位中心頻率的線性調(diào)頻信號，與仿真驗證該方法在處理瞬時頻率具有時間周期性信號時頻率估計的均方誤差根的大小比CT變換要低3 dB以上。同樣的信號利用STFT處理，多分量之間的相互干擾使得即使經(jīng)過重排處理，其可讀性也較差。

文獻[16]采用高斯函數(shù)的尺度變換、旋轉(zhuǎn)、時移和頻移構(gòu)成一簇基函數(shù)。其中旋轉(zhuǎn)通過分數(shù)階傅里葉變換來實現(xiàn)。通過旋轉(zhuǎn)角度的調(diào)節(jié)，可使基函數(shù)在頻帶中心的鄰域內(nèi)更好地匹配信號。該方法適應于局部波形具有高斯形狀的特征，具有較好的時頻聚集性和抗噪聲干擾的能力。但由于基于高斯函數(shù)的分數(shù)階傅里葉變換的自適應時頻分布不滿足時頻邊緣特性，因此不能用時間一階條件矩的方法估計瞬時頻率。文獻[17]在2 dB噪聲的環(huán)境中，準確地估計出瞬時頻率，而應用WVD估計出的瞬時頻率由于受噪聲的影響而出現(xiàn)了鋸齒，特別是在兩端，估計值出現(xiàn)振蕩，且遠小于實際值。

N.E.Huang提出了著名的Hilbert-Huang變換方法（簡稱為 HHT），其主要內(nèi)容是經(jīng)驗模型分解 （Empirical Mode Decomposition，簡稱為EMD）和Hilbert譜。在這一理論中，Huang等人引入了一種單分量的數(shù)學模型，他稱之為內(nèi)模函數(shù)（Intrinsic Mode Function，簡稱為 IMF），并給出了將一個任意信號自適應地分解成IMFs之和的分解算法，即EMD。這種方法主要是利用基于經(jīng)驗的模式分解方法，把復雜的數(shù)據(jù)序列分解為簡單的、有限個基本模式分量，得到的基本模式分量具有很好的Hilbert變換特性，使得瞬時頻率具有實際的物理意義。同時對復雜數(shù)據(jù)序列的瞬時頻率的引入，消除了傳統(tǒng)的信號分析方法如傅里葉分析所產(chǎn)生的表示非線性，非平穩(wěn)信號的偽諧波。由于它是基于信號局部特征的，因此它適用于非線性、非平穩(wěn)信號的處理。

3 研究與發(fā)展趨勢

從考察信號的頻率成分隨時間的演化特性角度來說，時頻分析技術力求在計算復雜性和時頻分辨率之間找到最佳平衡點。客觀地說，各種時頻分析技術難分軒，關鍵是其適合何種類型的信號。從時頻分析技術應用于瞬時頻率估計的研究過程中，不難發(fā)現(xiàn)以下的發(fā)展趨勢：

1）瞬時頻率與時頻分布關系密切，但時頻分布的計算量非常大，尋找WVD的快速計算方法，提高WVD的實現(xiàn)速度是一件十分有意義的跟蹤，必將進一步推動時頻分布在實際中的應用。

2）高階時頻分布由于較好的時頻聚集性能，成為估計非線性瞬時頻率的有力工具，然而對于多分量信號，高階時頻分布的交叉項將變得非常復雜。如何在盡可能保留高階時頻分布好的時頻聚集性能前提下減少交叉項的干擾是一個重要的發(fā)展方向。

3）自適應時頻分布的優(yōu)越性能使其必將成為瞬時頻率估計的趨勢。如何實現(xiàn)自適應信號分解，以及在信號長度增加的情況下，如何實現(xiàn)參數(shù)優(yōu)化是很有前途的研究方向。

[1]ArmstrongE.Amethodofreducingdistributesinradiosignaling by a method of frequency modulation[J].Proc.IEE， 1963（24）:669-740.

[2]Gabor D.Theory of communication[J].Proc.IEE，1946，93（3）:429-457.

[3]Qian S，CHEN Dang-pang.Joint time-frequency analysis[J].Prentice Hall，1996:41-52，181-187.

[4]Ville J.Theorie et application de al notion de signal analytique[J].Cable et Transmission，1948，2A（1）:61-74.

[5]Rao P，Taylor F J.Estimation of the instantaneous frequency using the discrete wigner distribution[J].Electronics Letters，1990， 26（4）:246-248.

[6]Boashash B.An efficient real-time implementation of the Wigner-Ville Distribution[J].IEEE Transaction on Acoustics Speech and Signal processing， 1987（35）:1611-1618.

[7]Martin W，F(xiàn)landrin P.Wigner spectral analysis of nonstationary process[J].IEEE Transaction on Acoustics Speech and Signal processing， 1987（33）:1461-1470.

[8]Barkat B，Boashash B.Design of higher order polynomial wigner-ville distributions[J].IEEE transaction on signal processing， 1999，47（9）:2608-2611.

[9]Boashash B.Time-frequency signal analysis[M].In Advances in Spectrum Analysis and Array processing，1990.

[10]Czerwinski R N.Adaptive short-time Fourier analysis[J].IEEE Signal Processing Letters， 1997， 4（2）:42-45.

[11]Kwok H K，Jones D L.Improved Instantaneous frequency estimation using an adaptive short-time Fourier transform[J].IEEE Trans.Signal Processing， 2000，48（10）：2964-2972.

[12]Changdra Sekhar S，Sreenias T V.Adaptive spectrogram vs.adaptive pseudo-Wigner-Ville distribution for instantaneous frequency estimation[J].Signal Processing， 2003（83）:1529-1543.

[13]Mallat S，Zhang Z.Adaptive time-frequency transform[J].IEEE transaction on Signal processing，1993（III），241-244.

[14]Angrisani L， D’Arco M.A measurement method based on an improved version ofthe Chirplettransform for instantaneous frequency estimation[C]//IEEE Instrumentation and Measurement Technology Conference，2001:1123-1130.

[15]Angrisani L， D’Arco M， Moriello Schiano L，et al.On the use of the Warblet transform for instantaneous frequency estimation[J].IEEE，IMTC 2004：935-941.

[16]馬世偉.非平穩(wěn)信號的參數(shù)自適應時頻表示及其應用的研究[D].上海：上海大學，2006.

[17]CHEN Guang-hua，CAO Jia-lin，WANG Jian.Use adaptive time-frequency distribution ofinstantaneous frequency estimation[J].Systems Engineering and Electronics，2002，24（1）:31-35.