基于遞歸希爾伯特變換的振動信號解調和瞬時頻率計算方法

胡志祥， 任偉新

(合肥工業大學 土木與水利工程學院，合肥　230009)

?

基于遞歸希爾伯特變換的振動信號解調和瞬時頻率計算方法

胡志祥， 任偉新

(合肥工業大學 土木與水利工程學院，合肥230009)

摘要：精確地提取振動信號的瞬時幅值和瞬時頻率對結構的參數識別和健康監測有重要作用。希爾伯特變換是一種常用的信號解調及瞬時頻率計算方法，但在信號不滿足Bedrosian乘積定理的條件時會造成較大誤差。針對這一問題，提出了一種遞歸希爾伯特變換方法，用前一步希爾伯特變換計算出的純調頻信號作為新的信號，遞歸地使用希爾伯特變換以進行信號解調，理論分析表明遞歸希爾伯特變換能夠快速地收斂。最后采用仿真信號對比了遞歸希爾伯特變換與單次希爾伯特變換、經驗調幅調頻分解及Teager能量算子法在信號解調及瞬時頻率計算中的結果，結果表明了遞歸希爾伯特變換方法的實用性及精確性。

關鍵詞：振動信號; 瞬時頻率; 信號解調; 希爾伯特變換; 經驗調幅調頻分解

振動信號蘊涵了動力系統當前的狀態信息，是結構健康監測系統需要監測的重要數據。對線性時不變系統，常用基于傅里葉變換的頻譜分析獲得各階振動頻率；而對非線性或時變結構，傅里葉變換的運用受到局限[1]。為此，眾多學者提出使用時頻分析方法來進行結構模態參數識別。

常用的時頻分析方法包括兩類：第一類是傅里葉變換的發展，如短時傅里葉變換、小波變換、同步擠壓小波變換、S變換等，這些方法能部分地克服傅里葉變換的缺點，在信號處理領域得到了廣泛的應用，但都是用預先設定的基函數來表達原始信號，存在自適應性不足、時頻分辨率低的問題[2-6]；第二類是基于信號分解的方法，如希爾伯特振動分解、廣義解調、迭代希爾伯特變換、希爾伯特-黃變換等，先將復雜信號分解為許多單分量信號，再計算各單成分信號的瞬時幅值和頻率，這幾種分解方法各自有不同的適用范圍，其中希爾伯特-黃變換應用最為廣泛[7-10]。

Huang等創造性地提出了經驗模式分解(Empirical Mode Decomposition,EMD)方法將復合信號分解為一系列本征模態函數(Intrinsic Mode Func-tion,IMF)，再利用希爾伯特變換計算各本征模態函數的幅值函數與瞬時頻率，此即希爾伯特-黃變換[10]。經驗模式分解沒有使用預設的基函數，而是取決于信號本身的特征，因而比較適合非線性非平穩信號分析。但是當IMF不滿足希爾伯特變換Bedrosian乘積定理的成立條件時，利用希爾伯特變換求本征模態函數的幅值和頻率將會存在誤差[11-12]。Huang等[13]為解決這一問題提出了經驗調幅調頻分解方法。張亢等[14]提出用分段波形的方法來計算單分量信號的瞬時頻率。經驗調幅調頻分解受到樣條曲線擬合誤差的影響，并且樣條曲線可能與信號本身相割而使最終結果在部分區間取值>1(或<-1)，將影響后續的頻率計算[15]。而分段波形法適用于幅值波動很小的情況，當相鄰半波內信號的幅值相差較大時，用階梯狀的不連續函數替代信號幅值函數是不合理的。Teager能量算子法也常用于信號解調與瞬率估計，但只適用于對幅值、瞬時頻率變化緩慢的情況[16]。

本文提出一種基于遞歸希爾伯特變換的單分量信號解調和瞬時頻率計算方法，用前一步希爾伯特變換獲得的純調頻函數作為新的信號，遞歸地使用希爾伯特變換計算新的純調頻信號。理論分析發現該遞歸過程收斂迅速，最后得出的純調頻信號有兩個特點：① 零點與原信號位置相同；② 其希爾伯特變換與其正交信號相等。遞歸希爾伯特變換不要求幅值函數與振蕩項頻譜互不重疊，因此在信號解調方面有更大的適用范圍。而所得純調頻信號零點位置與原信號相同，可使求得的瞬時頻率能夠接近真實瞬時頻率。最后通過算例比較了基于遞歸希爾伯特變換的方法與單次希爾伯特變換、經驗調幅調頻分解及Teager能量算子法在信號解調及瞬時頻率計算方面的結果，分析結果表明了遞歸希爾伯特變換的實用性及精確性。

1遞歸希爾伯特變換

對信號x0(t)進行希爾伯特變換即計算x0(t)與1/πt的卷積，其公式為

(1)

(2)

即解析函數的虛部是其實部的希爾伯特變換。該解析函數還可以表示為

z0(t)=A0exp(-jφ0)

(3)

其中

(4)

(5)

因此，可將原信號表示為幅值函數A0和純調頻信號cosφ0的乘積

x0(t)=A0cosφ0

(6)

若以純調頻函數x1(t)=cosφ0作為新的信號，并對其進行希爾伯特變換，可得到新的幅值函數和純調頻信號。不斷重復進行遞歸計算，遞歸公式為

(7)

(8a)

(8b)

在遞歸過程中，每次得到的新信號形狀不斷變化，直到幅值函數An趨于1，該遞歸過程的收斂性在下一節進行分析。xN+1(t)=cosφN為遞歸終止時得到的純調頻信號，它的希爾伯特變換與其正交信號相等，即

H[xN+1(t)]=H[cosφN]=sinφN

(9)

綜合所有的遞歸步驟，可將原信號表示為幅值函數與振蕩項的乘積

(10)

可見，利用遞歸希爾伯特變換可將原信號分解為幅值函數和純調頻函數兩部分，實現信號解調。容易證明，純調頻信號cosφN具有以下性質：

(1) cosφN的零點與原信號x0(t)的零點位置相同；

(2) cosφN的希爾伯特變換與其正交信號sinφN相等，而函數zN=cosφN+i·sinφN為解析函數。

2收斂性分析

(11)

令復信號z1(t)=cosφ0+i·sinφ0，并設其傅里葉變換為W1(ω)。根據Nuttal定理[16]，正交誤差函數的總能量可下式給出

(12)

對x1(t)進行希爾伯特變換，其幅值可由下式計算

(13)

推導中利用了式(11)中正交誤差函數的定義，進一步利用泰勒級數公式，并忽略高次項，可得

A1≈1+e1sinφ0

(14)

因此，進行第二次希爾伯特變后可得到新的純調頻信號x2(t)=cosφ1及其正交信號y2(t)=sinφ1，根據式(11)和式(14)可知

(15)

(16)

最后，再次計算x2(t)的希爾伯特變換，注意到e1為慢變函數，近似地利用Bedrosian乘積定理[11]，可得

(17)

根據式(16)和式(17)，可得出x2(t)對應的正交誤差函數，即

e2=Hcosφ1-sinφ1≈

(18)

可見，每一步遞歸希爾伯特變換使所得的純調頻信號對應的正交誤差函數取值減半。盡管在推導過程中采用了一些近似處理，兩相鄰遞歸步驟得到的純調頻信號分別對應的正交誤差函數的取值不一定正好相差2倍，但收斂的趨勢是顯而易見的。利用式(12)，設各遞歸步驟中en(t)的能量為En，以rn作為正交誤差函數能量遞減的指標，則

rn=En/En+1

(19)

后文算例將通過En和rn的變化來考察遞歸過程的收斂情況。

3仿真信號分析

為驗證遞歸希爾變換在信號解調及瞬時頻率估計方面的性能，將以下信號作為待處理信號：

信號1:x1=[1+0.1cos(2πt)]cos[4πt+

sin(0.5πt)]

信號2:x2=[2+0.2cos(6πt)]cos[4πt+

sin(0.5πt)]

顯然，信號1和信號2中振蕩項相同，瞬時頻率皆為f=2+0.25cos(0.5πt)，而兩信號幅值函數不同。與振蕩項相比，信號1的幅值函數為慢變函數，而信號2的幅值函數包含快變成分。因此，單次希爾伯特變換可分離信號1的幅值函數和振蕩項，但不能成功地對信號2進行解調。圖1給出了幾種信號解調方法分解出的幅值函數，并與實際幅值函數(Truth)進行了對比，圖中采用HT(Hilbert Transform)表示單次希爾伯特變換方法、EAMFMD(Empirical Amplitude Modulation Frequency Modulation Decomposition)表示經驗調幅調頻分解方法、TEAGER表示Teager能量算子法、RHT(Recursive Hilbert Transform)表示遞歸希爾伯特變換法(無特殊說明，下同)。圖2給出了上述解調方法計算出的瞬時頻率，并與實際瞬時頻率進行了對比。

為定量說明不同方法對信號幅值和瞬時頻率的識別性能，定義誤差指標為

(20)

式中：aI為識別出的物理量(幅值或頻率)，a為實際物理量。識別結果與真實結果越接近，誤差指標值越小。表1給出了不同計算方法得到的幅值函數和瞬時頻率識別對應的誤差指標。

表1　幅值和頻率識別的誤差指標

圖1　不同方法計算信號1和信號2的幅值函數對比Fig.1 Comparison of the amplitude functions of signal 1 and signal 2 by different methods

圖2　不同方法計算信號1和信號2的瞬時頻率對比Fig.2 Comparison of the instantaneous frequencies of signal 1 and signal 2 by different methods

對于信號1，RHT和HT方法計算的幅值函數誤差指標分別為0.006和0.010，瞬時頻率誤差指標分別為0.013和0.017，這表明計算值與實際值都比較接近。瞬時頻率在信號端部包含計算誤差，這是希爾伯特變換的端部效應造成的。實際應用中可通過信號鏡像或延拓的方法抑制端部效應造成的誤差。采用EAMFMD方法所得的幅值函數在信號端部包含較大誤差，這是樣條插值在信號端部不穩定導致的；EAMFMD方法估算的幅值與實際值比較接近；由于信號瞬時頻率時采用反余弦法，而解調所得的純調頻函數極值點不一定為1或-1，使計算出的瞬時頻率在極值點處有較大偏差，使瞬時頻率誤差指標高達32.590，實際使用中必須進行濾波處理。TEAGER方法所得幅值函數與實際值接近，但包含高頻誤差，同時瞬時頻率計算結果與實際瞬時頻率起伏趨勢相同，但誤差指標大于HT和RHT方法所得結果。事實上，只有在信號幅值和頻率都為常數的情況下TEAGER方法才能計算出準確結果，雖然幅值和頻率為時變函數但變化較慢時該方法仍可使用，但本文算例表明該方法并不適用于解調幅值函數和瞬時頻率隨時間快速變化的信號。

信號2的幅值函數含有高頻振蕩成分，利用HT方法不能有效解調此類信號。HT方法所得信號幅值和瞬時頻率與實際幅值和瞬時頻率都相差很大，誤差指標分別增至0.051和0.058。EAMFMD和TEAGER方法得出的幅值和瞬時頻率誤差都更大。所以，HT、EAMFMD和TEAGER方法都不能有效分解信號2。由計算結果可知，只有RHT方法能有效地解調信號2，誤差指標與信號1計算結果相同，這體現了RHT方法的穩定性，也表明RHT不受Bedrosian乘積定理成立條件的限制，幅值包含高頻振蕩成分時也能有效進行信號解調。

為驗證遞歸希爾伯特變換的收斂性及第3節提出的瞬時頻率計算方法，考察具有時變瞬時頻率的純調頻信號x=cos[10πt+9sin(πt)]，其瞬時頻率f=5+4.5cos(πt)，fmax=9.5 Hz，fmin=0.5 Hz。對該信號進行遞歸希爾伯特變換，根據式(12)和式(19)，可以計算出每次遞歸過程中的正交誤差函數的能量En及遞減指標rn。圖3和圖4分別給出了En和rn的變化趨勢。正交誤差函數的能量En隨遞歸次數的增大而減小，而在前16步遞歸過程中，遞減指標rn取值約等于2，隨后遞減指標rn快速趨近1，但保持rn≥1?？紤]到正交誤差函數的能量En在16步遞歸計算后只有5.8×10-9，此時的正交誤差函數可忽略不計。由此可知，本例中遞歸希爾伯特變換的收斂速度較快，在實際計算時可根據數據長度及精度需求合理選擇遞歸步數。

圖3　正交誤差函數能量變化趨勢Fig.3 Variation trend of the energy in the quadrature error signal

圖4　遞減指標變化趨勢Fig.4 Variation trend of the declining indicator

4結論

提出了一種基于遞歸希爾伯特變換的信號解調及瞬時頻率計算方法，與單次希爾伯特變換、經驗調幅調頻分解及Teager能量算子法等相比，具有一定的優越性。主要結論有：

(1) RHT遞歸計算得到的信號具有零點與原信號相同、其希爾伯特變換與正交信號相等的特點。

(2) RHT方法可用于信號解調，不受Bedrosian定理成立條件的限制，可提取包含高頻成分的幅值函數，擴展了希爾伯特變換的應用范圍。

(3) 理論分析及算例表明RHT具有較好的收斂性，正交誤差函數的能量快速下降。

總之，遞歸希爾伯特變換是對希爾伯特變換這一應用廣泛的信號處理方法的一個擴展。對任意給定的信號，利用RHT可以計算出零點與它相同、且希爾伯特變換與正交信號相等的純調頻信號。盡管本文對RHT方法進行了一定的理論與實例探討，但還有較多問題值得研究，例如RHT收斂結果唯一性的證明、RHT在解析信號構造方面的應用等，希望這些問題能引起更多學者關注。

參 考 文 獻

[ 1 ] 高維成, 劉偉, 鄒經湘. 基于結構振動參數變化的損傷探測方法綜述[J]. 振動與沖擊, 2004, 23(4):1-7.

GAO Wei-cheng, LIU Wei, ZOU Jing-xiang. Damage detection methods based on changes of vibration parameters: A summary review[J]. Journal of Vibration and Shock, 2004, 23(4): 1-7.

[ 2 ] Qian S, Chen D. Joint time-frequency analysis[J].IEEE Journals and Magazings, 1999, 16(2):52-67.

[ 3 ] Ruzzene M, Fasana A, Garibaldi L, et al. Natural frequencies and dampings identification using wavelet transform: Application to real data[J].Mechanical Systemsand Signal Processing, 1997, 11(2):207-218.

[ 4 ] Wang Chao, Ren Wei-xin, Wang Zou-cai, et al. Instantaneous frequency identification of time-varying structures by continuous wavelet transform[J].Engineering Structures, 2013, 52:17-25.

[ 5 ] Ditommaso R, Mucciarelli M, Ponzo F C. Analysis of non-stationary structural systems by using a band-variable filter[J]. Bulletin of Earthquake Engineering,2012,10(3):895-911.

[ 6 ] Daubechies I, Lu J, Wu H. Synchrosqueezed wavelet transforms: An empirical mode decomposition-like tool[J].Applied and Computational Harmonic Analysis, 2011, 30(2):243:361.

[ 7 ] 朱可恒,宋希庚, 薛冬新. 希爾伯特振動分解在滾動軸承故障診斷中應用[J]. 振動與沖擊,2014,33(14):160-164.

ZHU Ke-heng, SONG Xi-geng, XUE Dong-xin. Roller bearing fault diagnosis using Hilbert vibration decomposition[J]. Journal of Vibration and Shock, 2014, 33(14): 160-164.

[ 8 ] 程軍圣,楊宇,于德介.基于廣義解調時頻分析的多分量信號分解方法[J].振動工程學報,2007,20(6):563-569.

CHENG Jun-sheng, YANG Yu, YU De-jie. A multi-component signal decomposition method based on the general ized demodulation time-frequency analysis[J]. Journal of Vibration Engineering, 2007, 20(6): 563-569.

[ 9 ] Qin Yi, Qin Shu-ren, Mao Yong-fang. Research on iterated Hilbert transform and its application in mechanical fault diagnosis[J].Mechanical Systems and Signal Processing, 2008, 22(8):1967-1980.

[10] Huang N E, Shen Z, Long S R,et al. The empirical mode decomposition and Hilbert spectrum for nonlinear andnonstationary time series analysis[J].Proceedings of The Royal Society A Mathematical Physical and Engineering Sciences, 1998, 454(1971):903-995.

[11] Bedrosian E. A product theorem for Hilbert transforms[J].Proceedings of IEEE, 1963, 51(5): 868-869.

[12] Picinbono B. On instantaneous amplitude and phase of signals[J].IEEE Transactions on Signal Processing, 1997, 45(3): 552-560.

[13] Huang N E, Wu Z. A review on Hilbert-Huang transform: method and its applications to geophysical studies[J].Reviews of Geophysics,2008, 46(2):1-23.

[14] 張亢, 程軍圣, 楊宇,等. 基于分段波形的信號瞬時頻率計算方法[J].湖南大學學報:自然科學版, 2011, 38(11): 54-59.

ZHANG Kang, CHENG Jun-sheng, YANG Yu, et al. A piece-wise based signal instantaneous frequency computing method[J]. Journal of Hunan University:Natural Sciences, 2011, 38(11): 54-59.

[15] 戴豪民, 許愛強. 瞬時頻率計算方法的比較研究和改進[J].四川大學學報:自然科學版, 2014,51(6):1197-1204.

DAI Hao-min, XU Ai-qiang. Comparative research and improvement on the calculation method of instantaneous frequency[J]. Journal of Sichuan University:Natural Sciency Edition, 2014, 51(6):1197-1204.

[16] Alexandros P, Petros M. A comparison of the energy operator and the Hilbert transform approach to signal and speech demodulation[J].Signal Processing,1994,37(1):95-120.

Vibration signal demodulation and instantaneous frequency estimation based on recursive Hilbert transformation

HUZhi-xiang,RENWei-xin(School of Civil Engineering, Hefei University of Technology,Hefei 230009, China)

Abstract:Accurately extracting instantaneous amplitude and instantaneous frequency is important in structure parametic identification and health monitoring. Hilbert transformation is one of the most commonly used methods for signal demodulation and instantaneous frequency computation. However, it may cause larger errors when vibration signals do not satisfy the conditions of Bedrosian prodact theorem. Aiming at this problem, a recursive Hilbert transformation method was proposed. With this method, a pure frequency modulation signal derived in the previous step was taken as a new signal, it was modulated using Hilbent transformation recursively. The theoretical analysis showed that the recursive HirBert transformation can converge rapidly. The proposed method was compared with Hilbert transformation, the empirical AM-FM decomposition, and Teager energy method for simulated signal demodulation and instantaneous frequency computation. The results showed that the recursive Hilbert transformation.

Key words:vibrating signal; instantaneous frequency; signal demodulation; Hilbert transformation; empirical AM-FM decomposition

中圖分類號:TH165.3； TN911.7

文獻標志碼：A

DOI:10.13465/j.cnki.jvs.2016.07.006

通信作者任偉新 男，博士，教授，1960年生

收稿日期：2015-06-25修改稿收到日期：2015-10-23

基金項目：國家自然科學基金(51408177);中國博士后科學基金(2014M551802)

第一作者 胡志祥 男，博士，講師，1985年生