一種新的估計瞬時頻率的方法-經(jīng)驗包絡(luò)法

鄭近德，程軍圣，楊 宇

（湖南大學(xué) 汽車車身先進(jìn)設(shè)計制造國家重點(diǎn)實(shí)驗室，長沙 410082）

頻率在信號處理、通信、物理學(xué)等領(lǐng)域都是一個很重要的概念，它是刻畫波形的周期性質(zhì)和振蕩模式的一種屬性。頻率在物理上定義為周期的倒數(shù)，據(jù)此定義，如果要定義頻率，必須有一個完整的波形才能有周期。然而對一些平穩(wěn)或非平穩(wěn)信號而言，不存在固定的周期，但它卻有一定的振蕩模式，其頻率隨時間不斷變化，傳統(tǒng)的頻率的定義所具有的物理意義無法明確地描述其頻率瞬變現(xiàn)象。因此，需要一個類似于頻率的物理量來反映和刻畫信號這一性質(zhì)。于是相關(guān)學(xué)者提出了瞬時頻率的概念。Carson等提出了瞬時頻率的概念，并對其定義進(jìn)行了詳細(xì)研究。Gabor［1］給出了解析信號的概念，Ville等［2］提出了現(xiàn)在普遍接受的一般實(shí)信號的瞬時頻率（instantaneous frequency，簡稱IF）的定義。即，實(shí)信號的瞬時頻率定義為該信號所對應(yīng)的解析信號的相位函數(shù)關(guān)于時間的導(dǎo)數(shù)。其中，解析信號是基于希爾伯特變換而定義的。

常用的求信號的瞬時頻率的方法有：相位差分法［3］，鎖相環(huán)法［4-5］，基于希爾伯特變換的解析信號法［6］，Teager能量算子法［7-8］，以及 LMD 方法中的反余弦法［9-10］。這些方法即有自己的優(yōu)勢，但有各自的缺點(diǎn)。相位差分法包括前向差分，后向差分和中心差分，效果比較好的是中心差分法，其優(yōu)點(diǎn)在于對線性調(diào)頻信號是無偏的，且具有零延遲；而且它對應(yīng)為一系列時頻分布的一階矩。但其缺點(diǎn)是對噪聲太敏感，對含噪信號表現(xiàn)很大的方差［3］。Teager能量算子法是數(shù)學(xué)表達(dá)式的近似變換，原理簡單，并且無需進(jìn)行復(fù)數(shù)計算，計算量減小，因此近年來在求取信號的瞬時特征方面得到廣泛應(yīng)用，取得了很好的效果［8］。由于Teager能量算子法基于假設(shè)信號具有線性相位（即瞬時頻率為常數(shù)）且假設(shè)瞬時幅值近似為常數(shù)，因此，對于大部分窄帶信號有很好的效果，但如果瞬時幅值是時間的瞬變函數(shù)，或者波形有波內(nèi)調(diào)制（intrawave modulations）或諧波失真（harmonics distortion），能量算子法將會有很大的誤差，甚至不能使用。

對于單分量實(shí)信號，瞬時頻率最一般的求法是，通過希爾伯特變換定義其解析信號，瞬時頻率定義為其解析信號相位的導(dǎo)數(shù)。這是學(xué)術(shù)界普遍接受的定義瞬時頻率的方法。但是科恩指出，這種定義有幾種自相矛盾的佯謬結(jié)論［11］。首先，瞬時頻率可以不是頻譜中的頻率之一。其次，如果有只有少數(shù)的明顯的頻率組成的一個線性頻譜，那么瞬時頻率可以是連續(xù)的，而且可以在無數(shù)個值范圍內(nèi)變化。第三，雖然解析信號的頻譜對于負(fù)頻率為零，但瞬時頻率卻可以是負(fù)的。第四，對于一個帶限信號，它的瞬時頻率可以在頻帶之外。第五，為了計算某一時刻的解析信號，卻要必須知道全部時間的信號。事實(shí)上，不是所有的單分量信號都可以進(jìn)行希爾伯特變換，必須要滿足Bedrosian定理和 Nuttal定理［12］。

以上的討論一般只適合單分量信號，如果是多分量信號，一般是先通過 EMD［6］，LMD［9-10］等方法將其分解為若個單分量信號和一個殘余項之和。再對得到的單分量信號求瞬時頻率，進(jìn)而求得多分量信號的瞬時頻率。

基于希爾伯特變換的以上條件限制和問題，本文嘗試避開希爾伯特變換求瞬時頻率，提出了一種新的估計瞬時頻率的新方法-經(jīng)驗包絡(luò)法（empirical envelope method，簡稱 EE）。

1 經(jīng)驗包絡(luò)法

經(jīng)驗包絡(luò)法的核心是經(jīng)驗調(diào)幅調(diào)頻分解（簡稱經(jīng)驗AMFM分解），因此，下文先介紹經(jīng)驗AMFM分解。

1.1 經(jīng)驗AMFM分解［12］

一般單分量信號x（t）未必是形如x（t）=a（t）cosφ（t）的形式，而把x（t）寫作調(diào)幅部分和調(diào)頻部分的乘積的形式對于求瞬時頻率是非常重要的，能量算子法和反余弦法都基于此假設(shè)。因此有必要研究把任意的單分量信號唯一地寫作調(diào)幅和調(diào)頻部分的乘積形式的一般方法。諾頓·黃等［12］提出了一種把一般單分量信號特別是內(nèi)稟模態(tài)函數(shù)寫成調(diào)幅調(diào)頻形式的方法，稱之為經(jīng)驗調(diào)幅調(diào)頻分解。具體步驟如下：

（2）三次樣條擬合極值點(diǎn)（τk，xk）（k=1，2，…，M），得到信號的經(jīng)驗包絡(luò)函數(shù)a11（t）。對于一般實(shí)信號極值點(diǎn)是固定的，因此經(jīng)驗包絡(luò)函數(shù)也是唯一定義的。得到經(jīng)驗包絡(luò)函數(shù)后，可以用它來標(biāo)準(zhǔn)化數(shù)據(jù)，x（t）除以經(jīng)驗包絡(luò)函數(shù)a11（t），得：

（3）x1（t）是標(biāo)準(zhǔn)化的數(shù)據(jù)，理論上x1（t）的包絡(luò)估計函數(shù)a12（t），應(yīng)該小于等于1，否則，對x1（t）重復(fù)上述步驟經(jīng)過n次迭代：

直到a1n（t）≤1，即x1n（t）為一純調(diào)頻函數(shù)，迭代結(jié)束。記純調(diào)頻信號x1n（t）為F（t），則存在φ（t），使得：

（4）x（t）的調(diào)幅部分定義為：

至此，x（t）被分解為a（t）cosφ（t）的形式。

上述過程經(jīng)驗地實(shí)現(xiàn)了信號的調(diào)幅調(diào)頻分解。一般地，上述過程收斂過程比較快，迭代次數(shù)不會太多，1～3次即可實(shí)現(xiàn)數(shù)據(jù)的標(biāo)準(zhǔn)化。如經(jīng)驗?zāi)B(tài)分解一樣，上述方法只是經(jīng)驗分解，沒有解析的表達(dá)式和嚴(yán)格的數(shù)學(xué)證明。另外，分解過程可能會引起原始信號波形的失真，但失真的總和是可以忽略不記的，因為上述過程有過零點(diǎn)周期性地嚴(yán)格控制著，而過零點(diǎn)的位置是不變的。

本文在經(jīng)驗AMFM分解的基礎(chǔ)上提出了一種新的估算瞬時頻率的方法-經(jīng)驗包絡(luò)法。

1.2 經(jīng)驗包絡(luò)法

由經(jīng)驗AMFM分解過程易知，任意單分量信號x（t）可以近似寫作x（t）=a（t）cosφ（t）的形式，F(xiàn)（t）=cosφ（t）的瞬時頻率即為原始信號的瞬時頻率，因此，只要求出F（t）的瞬時頻率即可。由此提出了如下的求瞬時頻率的經(jīng)驗包絡(luò)法：

（1）首先由經(jīng)驗AMFM分解，任意單分量信號x（t）可寫作：x（t）=a（t）cosφ（t）。

（2）再令F（t）=cosφ（t），并對其兩邊求導(dǎo)，得：

式中，由于φ'（t）=2πf（t）一般為線性或相對載波部分sinφ（t）變換緩慢的函數(shù)，因此可視φ'（t）=2πf（t）為F'（t）的包絡(luò)部分。

（3）再對F'（t）進(jìn)行經(jīng)驗AMFM分解，得：

這里的b（t）近似為上式中的φ'（t）。因此，原信號的瞬時頻率為：

經(jīng)驗包絡(luò)法是基于信號的經(jīng)驗AMFM分解提出的，計算簡單方便，不需要繁雜的程序和極值點(diǎn)處的特殊處理，只要應(yīng)用兩次經(jīng)驗AMFM分解和一次求導(dǎo)即可。經(jīng)驗包絡(luò)法的核心是經(jīng)驗AMFM分解，分解的效果直接決定了求得的瞬時頻率的準(zhǔn)確性。這里有個矛盾，為了得到的純調(diào)頻信號，則需要增加迭代的次數(shù)，但由于采用三次樣條擬合包絡(luò)，迭代的次數(shù)增加，擬合的誤差會增大，求得的瞬時頻率誤差也會隨之增大。如果迭代的次數(shù)過少，三次樣條擬合的誤差減小，得到的純調(diào)頻信號的個別點(diǎn)仍大于1。實(shí)驗表明，經(jīng)驗包絡(luò)法對信號的純調(diào)頻程度要求不高，即允許有個別點(diǎn)的值大于1，對結(jié)果影響很小。因此，一般迭代次數(shù)越少，經(jīng)驗包絡(luò)法求得的瞬時頻率效果會更好。而如果得到的調(diào)頻信號有大于1的點(diǎn)，反余弦法則會出現(xiàn)突變，甚至不能使用，這也是經(jīng)驗包絡(luò)法優(yōu)于反余弦法之處。經(jīng)驗包絡(luò)法求取瞬時頻率的流程如圖1所示。

圖1 經(jīng)驗包絡(luò)法流程圖Fig.1 Flow-process diagram of emprical envelope method

事實(shí)上，經(jīng)過經(jīng)驗AMFM分解，信號被分解為調(diào)幅部分和純調(diào)頻部分的乘積，對于純調(diào)頻信號F（t）=cosφ（t），可以采取類似于LMD方法中的反余弦法求取瞬時頻率，即對F（t）=cosφ（t）兩邊求反余弦，得：

φ（t）相位以2π展開，再對其求導(dǎo)得到瞬時頻率：

基于反余弦法求瞬時頻率，直接方便，不需要希爾伯特變換，而且由于是基于信號本身的反余弦計算，計算量較小，精確性較高。但反余弦法的缺點(diǎn)是其求得的瞬時頻率在信號的極值點(diǎn)處會有不穩(wěn)定的突變，需要平滑處理。并且對信號的標(biāo)準(zhǔn)化要求較高，如果標(biāo)準(zhǔn)化后的數(shù)據(jù)有大于1的點(diǎn)，則在該處反余弦法會有很大的誤差。

文獻(xiàn)［12］提出另一種基于經(jīng)驗AMFM分解的估計瞬時頻率的方法，即標(biāo)準(zhǔn)希爾伯特變換法（Normalized Hilbert transform，簡稱NHT）。內(nèi)稟模態(tài)函數(shù)經(jīng)過經(jīng)驗AMFM分解，得到純調(diào)頻信號F（t）=cosφ（t），此時F（t）瞬時幅值為1，不再受Bodrosian定理的限制，因而可以對其希爾伯特變換，求取瞬時頻率。由于信號滿足Bodrosian定理，克服了希爾伯特變換會出現(xiàn)負(fù)頻率的缺陷，比直接希爾伯特變換有了很大的提高。但由于仍然采用希爾伯特變換，在端點(diǎn)處會產(chǎn)生能量泄露，因此端點(diǎn)效應(yīng)仍無法避免。

2 仿真信號分析

為了比較希爾伯特變換，標(biāo)準(zhǔn)希爾伯特變換，反余弦法和經(jīng)驗包絡(luò)法四種方法求取瞬時頻率的效果和優(yōu)缺點(diǎn)，先考察信號：

由于信號簡單和篇幅關(guān)系，信號時域波形略。分別用上述四種方法求式（14）所示信號的瞬時頻率，結(jié)果如圖2所示。其中，HT表示希爾伯特變換法，NHT表示標(biāo)準(zhǔn)希爾伯特變換法，ACOS表示反余弦法，EE表示經(jīng)驗包絡(luò)法，Truth表示信號的真實(shí)頻率，無特殊說明，下同。

圖2 不同方法求得的x（t）的瞬時頻率與其真實(shí)頻率Fig.2 IFs get by different methods of x（t）and its true IF

為了說明各種方法求得瞬時頻率的效果和精確性，誤差分析采用相對誤差來作為評價指標(biāo)。相對誤差定義為：

其中，A是真實(shí)值，這里為信號的真實(shí)頻率。Δ是絕對誤差，這里為各種方法求得的瞬時頻率與真實(shí)頻率差值的絕對值。

上述四種方法計算式（14）信號的瞬時頻率的相對誤差如表1。

表1 瞬時頻率的相對誤差Tab.1 The relative error of IFs

由圖2和表1知，希爾伯特變換求取的瞬時頻率有明顯的端點(diǎn)效應(yīng)，而且還出現(xiàn)了負(fù)頻率，相對誤差較大；經(jīng)過標(biāo)準(zhǔn)化后求取的瞬時頻率負(fù)頻率消失，相比希爾伯特變換有了很大的提高，但仍然有端點(diǎn)效應(yīng)；反余弦法求得的瞬時頻率與理論瞬時頻率吻合，但在兩端部分有突變現(xiàn)象；而由本文提出的經(jīng)驗包絡(luò)法求得的瞬時頻率和真實(shí)瞬時頻率最為接近，而且相對誤差比其他方法較小，此例初步說明了經(jīng)驗包絡(luò)法的有效性和優(yōu)越性。

上述考察的是單分量信號，初步表明了本文方法有很好的效果，進(jìn)一步考察多分量信號：

信號x（t）是調(diào)幅調(diào)頻信號x1（t）和調(diào)幅正弦信號x2（t）的疊加。其時域波形及其EMD分解結(jié)果如圖3所示。

圖3 x（t）及其EMD分解結(jié)果Fig.3 x（t）and its EMD decomposition results

分別用上述四種方法求取分量IMF1，IMF2的瞬時頻率，結(jié)果如圖4、圖5所示。需要說明的是，由于分量IMF1和IMF2理論上分別對應(yīng)為調(diào)幅調(diào)頻信號x1（t）和調(diào)幅正弦信號x2（t），故圖4和圖5中的真實(shí)頻率Truth分別表示x1（t）和x2（t）的瞬時頻率。

上述四種方法求得的第一個分量IMF1的瞬時頻率的相對誤差如表2所示。

圖4 不同方法求得的IMF1的瞬時頻率Fig.4 IFs of the IMF1calculated by different methods

表2 不同方法求得IMF1的瞬時頻率的相對誤差Tab.2 Relative errors of the IFs of IMF1 calculated by different methods

由圖5和表2知，希爾伯特變換和標(biāo)準(zhǔn)希爾伯特變換法求得的IMF2瞬時頻率端點(diǎn)效應(yīng)較嚴(yán)重。而反余弦法和經(jīng)驗包絡(luò)法求得的瞬時頻率結(jié)果都非常接近理論瞬時頻率，二者都有非常好的效果。但經(jīng)驗包絡(luò)法求得的瞬時頻率的相對誤差比反余弦法更小。

圖5 不同方法求得的IMF2的瞬時頻率Fig.5 IFs of the IMF2calculated by different methods

上述四種方法求得的第二個分量IMF2的瞬時頻率的相對誤差如表3所示。

表3 不同方法求得IMF2的瞬時頻率的相對誤差Tab.3 Relative errors of the IFs of IMF2 calculated by different methods

由圖5和表3可知，由希爾伯特變換求得的IMF2的瞬時頻率出現(xiàn)了負(fù)頻率，而標(biāo)準(zhǔn)希爾伯特變換克服了這一缺點(diǎn)，有了很大的提高，但端點(diǎn)效應(yīng)很明顯。反余弦法求得的瞬時頻率都有輕微的波動。而由經(jīng)驗包絡(luò)法求得的瞬時頻率相對誤差最小，最接近理論值，與上述方法相比，經(jīng)驗包絡(luò)法有明顯的優(yōu)越性。

需要說明的是，上述圖4，圖5中求得的IMF1，IMF2的瞬時頻率的誤差和以及端點(diǎn)效應(yīng)，不僅僅是因為經(jīng)驗AMFM分解，還有一部分是因為EMD分解，因為EMD分解也會出現(xiàn)端點(diǎn)效應(yīng)。

為更好地說明上述四種方法的分解效力，下面以實(shí)測信號為例來說明。圖6是一外圈有凹槽的6311型球滾動軸承的振動加速度信號及其EMD分解的前四個IMF分量的時域波形。采樣頻率為4 096 Hz，采樣點(diǎn)數(shù)為1 024，轉(zhuǎn)頻為25 Hz。

圖6 軸承振動加速度信號時域波形及其EMD分解的前四個IMF分量波形Fig.6 Bearing vibration acceleration signal waveform and its first four components generated by EMD

由于前兩個分量的瞬時頻率較高，波動較大，不方便比較，因此，分別用希爾伯特變換，標(biāo)準(zhǔn)希爾伯特變換，反余弦法和經(jīng)驗包絡(luò)法求第三個分量IMF3和第四個分量IMF4的瞬時頻率，如下圖7，8所示。

從圖7和圖8中可以看出，希爾伯特變換求得的瞬時頻率有負(fù)頻率，而且有端點(diǎn)效應(yīng)，而反余弦法和經(jīng)驗包絡(luò)法求得的瞬時頻率較為理想。但從圖8可以看出，反余弦求得的瞬時頻率有突變點(diǎn)，經(jīng)驗包絡(luò)法要優(yōu)于反余弦法。

3 結(jié)論

提出了一種新的估計瞬時頻率方法，經(jīng)驗包絡(luò)法，并將其與希爾伯特變換，標(biāo)準(zhǔn)希爾伯特變換和反余弦法進(jìn)行了比較，并通過仿真信號和實(shí)測信號分析，結(jié)果表明經(jīng)驗包絡(luò)法具有一定的優(yōu)越性。經(jīng)驗包絡(luò)法的關(guān)鍵在于求取信號包絡(luò)，經(jīng)驗調(diào)幅調(diào)頻分解中是采用三次樣條擬合信號絕對值的極大值而得，因此會有一定的擬合誤差。如果提高求取包絡(luò)的方法，經(jīng)驗包絡(luò)法的準(zhǔn)確性也將會有很大的提高。實(shí)例已表明有經(jīng)驗包絡(luò)法有很好效果，但對此方法的使用范圍，誤差產(chǎn)生的較小等問題還需要進(jìn)一步研究，以使其更完善和精確，應(yīng)用更廣泛。

［1］ Gabor D.Theory of communication［J］.Journal of the IEEE，1946，93：429-457.

［2］Ville J.Theorie et application de la notion de signal analyticque［J］.Cables et Transmissions，1948，2A（1）：61-74.

［3］張賢達(dá).非平穩(wěn)信號分析與處理［M］.北京：國防工業(yè)出版社，1998.

［4］Boashash B.Estimating and interpreting the instantaneous frequency of a signal-part 1：fundamentals［J］.Proceedings of the IEEE，1992，80（4）：520-538.

［5］Boashash B.Estimating and interpreting the instantaneous frequency of a signal-part 2：algorithms and applications［J］.Proceedings of The IEEE，1992，80（4）：540-568.

［6］Huang N E，Wu Z H，A review on Hilbert-Huang transform：method and its applications to geophysical studies［J］.Reviews of Geophysics，2008，46（2）：1-23.

［7］MaragesP， KaiserJF. Energy separation in signal modulation to speech analysis［J］.IEEE Transactions on Signal Process，1993，41（10）：3024-3051.

［8］ Kaiser J F.Some useful properties of Teager＇s energy operators［J］.Proc.IEEE Int.Conf.Acoust.Speech Signal Processing，1993（3）：140-158.

［9］Smith J S.The localmean decomposition and its app lication to EEG perception data［J］.Journal of the Royal Society Interface，2005，2（5）：443-454.

［10］程軍圣，張 亢，楊 宇，等.局部均值分解與經(jīng)驗?zāi)B(tài)分解的對比研究［J］.振動與沖擊，2009，28（5）：13-16.

［11］科恩著，白居憲，譯.時頻分析：理論與應(yīng)用［M］.西安：西安交通大學(xué)出版社，1998.

［12］ Huang N E，Wu Z，Long S R，et al.On the frequency［J］.Advances in Adaptive Data Analysis.2009，1（2）：177-229.