一種新的估計瞬時頻率的方法-經驗包絡法

鄭近德，程軍圣，楊 宇

（湖南大學 汽車車身先進設計制造國家重點實驗室，長沙 410082）

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

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

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

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

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

1 經驗包絡法

經驗包絡法的核心是經驗調幅調頻分解（簡稱經驗AMFM分解），因此，下文先介紹經驗AMFM分解。

1.1 經驗AMFM分解［12］

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

（2）三次樣條擬合極值點（τk，xk）（k=1，2，…，M），得到信號的經驗包絡函數a11（t）。對于一般實信號極值點是固定的，因此經驗包絡函數也是唯一定義的。得到經驗包絡函數后，可以用它來標準化數據，x（t）除以經驗包絡函數a11（t），得：

（3）x1（t）是標準化的數據，理論上x1（t）的包絡估計函數a12（t），應該小于等于1，否則，對x1（t）重復上述步驟經過n次迭代：

直到a1n（t）≤1，即x1n（t）為一純調頻函數，迭代結束。記純調頻信號x1n（t）為F（t），則存在φ（t），使得：

（4）x（t）的調幅部分定義為：

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

上述過程經驗地實現了信號的調幅調頻分解。一般地，上述過程收斂過程比較快，迭代次數不會太多，1～3次即可實現數據的標準化。如經驗模態分解一樣，上述方法只是經驗分解，沒有解析的表達式和嚴格的數學證明。另外，分解過程可能會引起原始信號波形的失真，但失真的總和是可以忽略不記的，因為上述過程有過零點周期性地嚴格控制著，而過零點的位置是不變的。

本文在經驗AMFM分解的基礎上提出了一種新的估算瞬時頻率的方法-經驗包絡法。

1.2 經驗包絡法

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

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

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

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

（3）再對F'（t）進行經驗AMFM分解，得：

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

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

圖1 經驗包絡法流程圖Fig.1 Flow-process diagram of emprical envelope method

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

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

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

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

2 仿真信號分析

為了比較希爾伯特變換，標準希爾伯特變換，反余弦法和經驗包絡法四種方法求取瞬時頻率的效果和優缺點，先考察信號：

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

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

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

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

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

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

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

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

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

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

分別用上述四種方法求取分量IMF1，IMF2的瞬時頻率，結果如圖4、圖5所示。需要說明的是，由于分量IMF1和IMF2理論上分別對應為調幅調頻信號x1（t）和調幅正弦信號x2（t），故圖4和圖5中的真實頻率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知，希爾伯特變換和標準希爾伯特變換法求得的IMF2瞬時頻率端點效應較嚴重。而反余弦法和經驗包絡法求得的瞬時頻率結果都非常接近理論瞬時頻率，二者都有非常好的效果。但經驗包絡法求得的瞬時頻率的相對誤差比反余弦法更小。

圖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的瞬時頻率出現了負頻率，而標準希爾伯特變換克服了這一缺點，有了很大的提高，但端點效應很明顯。反余弦法求得的瞬時頻率都有輕微的波動。而由經驗包絡法求得的瞬時頻率相對誤差最小，最接近理論值，與上述方法相比，經驗包絡法有明顯的優越性。

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

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

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

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

從圖7和圖8中可以看出，希爾伯特變換求得的瞬時頻率有負頻率，而且有端點效應，而反余弦法和經驗包絡法求得的瞬時頻率較為理想。但從圖8可以看出，反余弦求得的瞬時頻率有突變點，經驗包絡法要優于反余弦法。

3 結論

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

［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］張賢達.非平穩信號分析與處理［M］.北京：國防工業出版社，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］.振動與沖擊，2009，28（5）：13-16.

［11］科恩著，白居憲，譯.時頻分析：理論與應用［M］.西安：西安交通大學出版社，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.