HHT變換講義EMD

1、HHT變換講義1.1簡介傳統(tǒng)的信號處理方法，如傅立葉分析是一種純頻域的分析方法。它用頻率不同的各復(fù)正弦分量的疊加來擬合原函數(shù)，也即用來分辨。而在有限頻域上的信息不足以確定在任意小范圍內(nèi)的函數(shù)，特別是非平穩(wěn)信號在時間軸上的任何突變，其頻譜將散布在整個頻率軸上。而且，非平穩(wěn)動態(tài)信號的統(tǒng)計特性與時間有關(guān)，對非平穩(wěn)信號的處理需要進行時頻分析，希望得到時域和頻域中非平穩(wěn)信號的全貌和局域化結(jié)果。在傅立葉變換中，人們?nèi)粝氲玫叫盘柕臅r域信息，就得不到頻域信息。反之亦然。后來出現(xiàn)的小波（Wavelet）變換通過一種可伸縮和平移的小波對信號變換，從而達到時頻局域化分析的目的。但這種變換實際上沒有完全擺脫傅立葉變換

2、的局限，它是一種窗口可調(diào)的傅立葉變換，其窗內(nèi)的信號必須是平穩(wěn)的。另外，小波變換是非適應(yīng)性的，小波基一旦選定，在整個信號分析過程中就只能使用這一個小波基了。HHT（Hilbert-Huang Transform）技術(shù)是(1998年由NASA的Norden E Huang 等提出的新的信號處理方法。該方法適用于非線性非平穩(wěn)的信號分析，被認為是近年來對以傅立葉變換為基礎(chǔ)的線性和穩(wěn)態(tài)譜分析的一個重大突破。目前HHT技術(shù)已用于地球物理學(xué)和生物醫(yī)學(xué)等領(lǐng)域的研究，并取得了較好的結(jié)果。存在的問題盡管HHT技術(shù)在處理非線性、非穩(wěn)態(tài)信號方面有很大的優(yōu)勢，但是這個方法本身還是有許多的問題有待進一步研究。正如Huan

3、g 在文章中指出的那樣，對于這種新的信號處理方法，其基的完備性還需要嚴(yán)密的證明。另外，在做Hilbert變換時出現(xiàn)的邊界效應(yīng)也需要更好的方法來解決。但是，HHT技術(shù)中最嚴(yán)重，也是現(xiàn)今研究的最多的是EMD 分解中的包絡(luò)過程。從對EMD分解方法的介紹可以看出，包絡(luò)線的構(gòu)造影響著整個分解的結(jié)果，也決定了后面的Hilbert變換。Huang 采用的三次樣條插值來擬和包絡(luò)線。在實際應(yīng)用中，發(fā)現(xiàn)這樣做會產(chǎn)生嚴(yán)重的邊界效應(yīng)，污染了原始數(shù)據(jù)。特別是對短數(shù)據(jù)而言，這種影響使分析所得的結(jié)果失去了原有的意義。對此，Huang等采用的是根據(jù)信號端點處的振幅和頻率，分別增加兩組特征波的方式進行數(shù)據(jù)延拓。Huang的這種

4、延拓方法已經(jīng)向NASA申請了專利。除此之外，還有人提出了其它方法進行端點延拓。比如國家海洋局的黃大吉等提出的鏡像閉合延拓法和極值延拓法。其中鏡像閉合延拓法是根據(jù)信號的分布特性，把鏡子放在具有對稱性的極值位置，通過鏡像法把鏡內(nèi)信號映射成一個周期性的環(huán)形信號，不存在端點，從根本上避免了端點效應(yīng)。而極值延拓法具有鏡像閉合延拓法相當(dāng)?shù)男Ч辉黾有盘栃蛄斜旧淼拈L度，計算較快，尤其在處理非對稱波形信號時比較優(yōu)越。還有青島海洋大學(xué)的鄧擁軍等，利用神經(jīng)網(wǎng)絡(luò)分析方法來延拓數(shù)據(jù)端。這些方法對邊界效應(yīng)的抑制都有一定的作用，但是也都需要更深一步的研究。1.2 瞬時頻率頻率是個極其重要的物理量，定義為信號周期倒數(shù)，

5、其物理含義顯而易見。對于正弦信號，它的頻率為恒值。但是對于大部分信號，它的頻率是隨時間變化的函數(shù)，故提出瞬時頻率概念。瞬時頻率即表征信號在局部時間點上瞬態(tài)頻率特性，整個持續(xù)期上的瞬時頻率反映了信號頻率的時變規(guī)律。對于隨機時間序列X(t)，對其進行Hilbert變換39，可以得到Y(jié)(t)如下： (1-1)其中，PV為柯西主值（Cauchy principal value)。該式表示Y(t)是X(t)與的卷積。通過這個定義，X(t)和Y(t)組成了一個共軛復(fù)數(shù)對，于是可以得到一個解析信號Z(t)如下 (1-2)其中 (1-3)從理論上講，虛部的定義方法有很多種。但是Hilbert變換為其提供了一個

6、唯一的虛部值，這就使得其結(jié)果成為一個解析函數(shù)。得到了相位，就可以得到瞬時頻率，因為瞬時頻率就是相位導(dǎo)數(shù)。 (1-4)從本質(zhì)上說，公式(1-1)將Hilbert變換定義為X(t)與的卷積，因此它強調(diào)了X(t)的局部特性。公式(1-2)的極坐標(biāo)表達式更進一步澄清了這個表達式的局部特性：它的幅值和相位的最佳局部匹配隨三角函數(shù)變化。1.3 固有模態(tài)函數(shù)IMF（Intrinsic Mode Function）由瞬時頻率的物理意義可知，并不是任意的信號都能用瞬時頻率來討論。只有當(dāng)信號滿足只包括一種振動模式，而沒有復(fù)雜疊加波的情況時才行。實際上，定義一個有意義的瞬時頻率的必要條件就是要求函數(shù)關(guān)于局部零平均值

7、對稱，并且零交叉點和極值點數(shù)量相同?；诖朔N原因，提出了固有模態(tài)函數(shù)的概念。固有模態(tài)函數(shù)滿足以下兩個條件：(1)整個數(shù)據(jù)范圍內(nèi)，極值點和過零點的數(shù)量相等或者相差一個；(2)在任意點處，所有極大值點形成的包絡(luò)線和所有極小值點形成的包絡(luò)線的平均值為零。第一個條件是顯而易見的，它類似于平穩(wěn)過程中傳統(tǒng)的穩(wěn)定且滿足高斯分布正態(tài)分布高斯白噪聲：如果一個噪聲，它的瞬時值服從高斯分布，而它的功率譜密度又是均勻分布的，則稱它為高斯白噪聲。的窄帶信號條件。第二個條件把傳統(tǒng)的全局條件調(diào)整到局部情況。只有滿足了這個條件，得到的瞬時頻率才不會因為不對稱波形的存在而引起不規(guī)則波動。所以這一點是得到正確瞬時頻率的必要條件。

8、而這一點是必須的，因為這樣瞬時頻率就可以不包含由于不對稱波形造成的波動。為了使用瞬時頻率定義，必須要把隨機數(shù)據(jù)歸結(jié)為IMF組件，這樣才可以為每個IMF組件定義瞬時頻率。為了將數(shù)據(jù)歸結(jié)為所需的IMF組件，接下來引入經(jīng)驗?zāi)Ｊ椒纸夥椒ā?.4 經(jīng)驗?zāi)Ｊ椒纸釫MD（Empirical Mode Decomposition）經(jīng)驗?zāi)B(tài)分解方法的大體思路是利用時間序列上下包絡(luò)的平均值確定“瞬時平衡位置”，進而提取固有模態(tài)函數(shù)39。這種方法基于如下假設(shè)：（1）信號至少有兩個極點個極大值和一個極小值；（2）信號特征時間尺度是由極值間的時間間隔來確定的；（3）如果數(shù)據(jù)沒有極值而僅有拐點，可以通過微分、分解、再積分

9、的方法獲得IMF。在此假設(shè)基礎(chǔ)上，Huang等人進一步指出：可以用經(jīng)驗?zāi)B(tài)分解方法將信號的固有模態(tài)篩選出來。經(jīng)驗?zāi)Ｊ椒纸膺^程就是個篩選過程，實現(xiàn)振動模式的提取。該方法的基本思路是用波動上、下包絡(luò)的平均值去確定“瞬時平衡位置”，進而提取出固有模態(tài)函數(shù)。上、下包絡(luò)線是由三次樣條函數(shù)對極大值點和極小值點進行擬合得到的。經(jīng)驗?zāi)Ｊ椒纸膺^程的基本過程可概括如下：(1)尋找信號x(t)所有局部極大值和局部極小值，為更好保留原序列的特性，局部極大值定義為時間序列中的某個時刻的值，它只要滿足既大于前一時刻的值也大于后一時刻的值即可。局部極小值的提取同理，即該時刻的值滿足既小于前一時刻的值也小于后一時刻的值。使用

10、三次樣條函數(shù)進行擬合，獲得上包絡(luò)線xmax(t)和下包絡(luò)線xmin(t)；(2)計算上、下包絡(luò)線的均值m(t)=xmax(t)+xmin(t)/2；(3)用原信號x(t)減去均值m(t)，得到第一個組件h(t)=x(t)-m(t)；由于原始序列的差異，組件h(t)不一定就是一個IMF，如果h(t)不滿足固有模態(tài)函數(shù)兩個條件，就把h(t)當(dāng)成原始信號，重復(fù)(1)-(3)，直到滿足條件為止，這時滿足固有模態(tài)函數(shù)條件的h(t)作為一個IMF，令I(lǐng)1(t)=h(t)，至此第一個IMF已經(jīng)成功的提取了。由于剩余的r(t)=x(t)-I1(t)仍然包含具有更長周期組件的信息，因此可以把它看成新的信號，重復(fù)

11、上述過程，依次得到第二個I2(t)，第三個I3(t)，當(dāng)r(t)滿足單調(diào)序列或常值序列條件時，終止篩選過程，可以認為完成了提取固有模態(tài)函數(shù)的任務(wù)，最后的r(t)稱為余項，它是原始信號的趨勢項。由此可得x(t)的表達式，即原始序列是由n個IMF與一個趨勢項組成。如上所述，整個過程就像篩選過程，根據(jù)時間特性把固有模態(tài)函數(shù)從信號中提取出來。1.5 端點延拓HHT方法的分析質(zhì)量很大程度上取決于EMD分解的質(zhì)量，而在應(yīng)用EMD方法時的一個非常棘手的問題是，由于信號兩端不可能同時處于極大值和極小值，這樣，“篩”過程中構(gòu)成上下包絡(luò)的三次樣條函數(shù)在數(shù)據(jù)序列的兩端就會出現(xiàn)發(fā)散現(xiàn)象（如圖1.1(a)所示，原始信號

12、是中心頻率為1KHz、5個波峰、以100個數(shù)據(jù)描繪一周期的窄帶信號），并且這種發(fā)散的結(jié)果會隨著“篩”過程的不斷進行逐漸向內(nèi)“污染”整個數(shù)據(jù)序列而使所得結(jié)果嚴(yán)重失真（如圖1.1(b)所示）。另外，在進行Hilbert變換時，信號的兩端也會出現(xiàn)嚴(yán)重的端點效應(yīng)。對于一個較長的數(shù)據(jù)序列來講，可以根據(jù)極值點的情況不斷拋棄兩端的數(shù)據(jù)來保證所得到的包絡(luò)的失真度達到最小。但對于一個數(shù)據(jù)點數(shù)少的序列來講，這樣的操作就變得完全不可行。因此，必須對信號或其極值向外進行延拓，以確保包絡(luò)線抵達端點。端點延拓的目的是確保上、下包絡(luò)都與端點相交，以便有與每一個信號點相對應(yīng)的局部平均值。而上、下包絡(luò)是由極大值和極小值連結(jié)而成

13、的，因此只要對極大值和極小值進行延拓，而不必對信號本身進行延拓。極大值和極小值是相間分布的，同時考慮到樣條插值的要求，所以只要在信號左、右兩端分別延拓兩個極大值和兩個極小值即可。由于端點以外沒有信號，任何延拓都是人為的，通過對不同形式的信號的多次試驗，發(fā)現(xiàn)一種端點延拓的方法以端點的一個特征波為依據(jù)進行延拓具有較好的效果，具體做法如下40：（b）EMD分解(a) 擬合結(jié)果原始信號上下包絡(luò)線圖3.1 未作端點延拓的分析效果(a) 擬合結(jié)果（b）EMD分解原始信號上下包絡(luò)線圖3.2 以端點的一個特征波為依據(jù)進行延拓的分析效果設(shè)離散信號： (1-5) (1-6)其采樣步長為,有個極大值和個極小值，對應(yīng)

14、的序列下標(biāo)（,）、時間（ ,）和函數(shù)值（，）記為： (1-7) (1-8) (1-9) (1-10)具體作法如下：1.左端信號左端第一個特征波包含的信號點數(shù)為 (1-11)向外延拓的兩個極值的位置（,）和數(shù)值（，）為： (1-12) (1-13) (1-14) (1-15)2.右端信號右端第一個特征波包含的信號點數(shù)為k2 (1-16)向外延拓的兩個極值的位置（,）和數(shù)值（，）為： (1-17) (1-18) (1-19) (1-20)3.當(dāng)端點的數(shù)值比近端點的第一個極大值大或極小值小時，要進行特殊的處理，以避免信號落到包絡(luò)線之外， (1-21) (1-22) (1-23) (1-24)由圖1.

15、2可見，以端點的一個特征波為依據(jù)進行延拓，分別在信號兩端增加兩個極大值和兩個極小值，從而使原始信號被延長的包絡(luò)線所限制，有效地抑制了端點效應(yīng)，使得EMD分解得到合理的各個IMF模態(tài)。另外，國家海洋局的黃大吉等還提出了鏡像閉合延拓法，即根據(jù)信號的分布特性，把鏡子放在具有對稱性的極值位置，通過鏡像法把境內(nèi)信號映射成一個周期性的環(huán)形信號。這樣做就不存在端點，從根本上避免了端點效應(yīng)40。青島海洋大學(xué)的鄧擁軍等人利用神經(jīng)網(wǎng)絡(luò)分析的方法來進行端點延拓，在數(shù)據(jù)的兩端各得到兩個附加的極大值點和兩個附加的極小值點，由此有效的抑制了端點效應(yīng)41。1.6 EMD結(jié)束準(zhǔn)則根據(jù)理論分析IMF必須滿足兩個條件，但是在實際

16、篩選過程中，嚴(yán)格滿足這兩個條件的信號幾乎不存在，因此如果以這兩個條件為依據(jù)判定IMF，可能根本就得不到結(jié)果或者要以冗長的程序執(zhí)行時間為代價。因此篩選過程要小心進行，為了保證IMF組件的調(diào)幅和調(diào)頻都具有物理意義，同時考慮到程序可行性，必須制定一個結(jié)束篩選的準(zhǔn)則。習(xí)慣上，可以通過標(biāo)準(zhǔn)偏差來完成，它可以由連續(xù)兩個篩選結(jié)果得到18： (1-25)h1(k-1)(t)和h1k(t)分別表示兩個連續(xù)的篩選結(jié)果，一般來說，SD值越小，所得的IMF分量的線性和穩(wěn)定性就越好，典型的SD值為0.2-0.3，這樣對于不同篩選過程就有嚴(yán)格的SD限制了。但是在篩選過程中發(fā)現(xiàn)，SD的這個范圍過于死板，對很多信號而言，以0

17、.2-0.3為限制的值并不科學(xué)，因此在實際篩選過程中，可以以0.2-0.3為參考值，根據(jù)實際情況進行適當(dāng)調(diào)整。另外，如果不是進行信號精確處理，可以考慮控制篩選次數(shù)來決定篩選結(jié)果，這樣不僅保障了程序的有效性，也控制了程序的執(zhí)行時間，在具體問題中取得了很好的效果。1.7 Hilbert譜完成了經(jīng)驗?zāi)Ｊ椒纸膺^程就得到了所有可提取的固有模態(tài)函數(shù)，在Hilbert變換的基礎(chǔ)上，只要根據(jù)公式(1-1)-(1-4)計算瞬時頻率就可以了。對固有模態(tài)函數(shù)進行Hilbert變換后，可以用下面方式表示信號 (1-26)在這兒沒有考慮余項rn，因為它只是單調(diào)函數(shù)或常量。盡管Hilbert變換可以把單調(diào)函數(shù)看成是振動的

18、一部分，但是其余項中的能量很小，一般情況可以不用考慮。由希爾伯特變換得出的振幅和頻率都是時間的函數(shù)，如果用三維圖形表達幅值、頻率和時間之間的關(guān)系，或者把振幅用灰度的形式顯示在頻率時間平面上，就可以得到Hilbert譜。如果把對時間積分，就得到希爾伯特邊際譜： (1-27)邊際譜提供了對每個頻率的總振幅的量測，表達了整個時間長度內(nèi)累積的振幅。另外，作為希爾伯特邊際譜的附加結(jié)果，可以得到如式（1-28）定義的希爾伯特瞬時能量： (1-28)瞬時能量提供了信號能量隨時間的變換情況。事實上，如果振幅的平方對時間積分，可以得到希爾伯特能量譜： (1-29)希爾伯特能量譜提供了對于每個頻率的能量的量測，表

19、達了每個頻率在整個時間長度內(nèi)所累積的能量。1.8 HHT方法的軟件實現(xiàn)HHT變換的整個過程包括EMD與Hilbert變換，整個過程是通過matlab軟件編寫的M文件實現(xiàn)的，其程序流程圖如圖1.3所示。開始輸入信號x(t)r(t)=x(t),n=0x(t)求出x(t)的所有極值點擬和上下包絡(luò)線計算包絡(luò)線的平均值m(t)h(t)=x(t)-m(t)h是否滿足IMF條件n=n+1,Cn (t)=h(t),r (t)=r (t)-Cn (t)Cn(t)是否滿足EMD中止條件對Cn(t)進行Hilbert變換構(gòu)造解析函數(shù)求瞬時頻率邊際譜,希爾波特能量譜瞬時能量譜等結(jié)束是是x(t)=h(t)否否x(t)=

20、h(t)圖1.3 經(jīng)驗?zāi)Ｊ椒纸馀cHilbert變換流程圖1.9 HHT方法的優(yōu)越性1、HHT方法是一種全新的信號分析方法，能夠描繪出信號的時頻譜圖、邊際譜，能量譜等，是一種更具有適應(yīng)性的時頻局域化分析方法28。而傳統(tǒng)的信號處理方法，如傅立葉變換是一種純時域的分析方法，它用頻率從零到無窮大的各復(fù)正弦分量的疊加來擬合原函數(shù)f(t)在每個時刻的值，也即用F()在有限頻域上的信息不足以確定在很小范圍內(nèi)的函數(shù)f(t)，特別是非平穩(wěn)信號在時間軸上的任何突變，其頻譜將散布在整個頻率軸上39。所以，這種分析方法適用于確定性的平穩(wěn)信號，而在非線性、非平穩(wěn)過程的處理上，傅立葉變換將只作為一種數(shù)學(xué)變化手段。而且，非

21、平穩(wěn)信號的統(tǒng)計特性與時間有關(guān)，對非平穩(wěn)信號的處理需要進行時頻分析，希望得到時域和頻域中非平穩(wěn)信號的全貌和局域化結(jié)果。但在傅立葉變換中，若想得到信號的時域信息，就得不到頻域信息，反之亦然。后來出現(xiàn)的小波變換通過一種可伸縮和平移的基小波對信號變換，從而達到時域局域化分析的目的，但這種變化實際上并沒有完全擺脫傅立葉變換的局限，它是一種窗口可調(diào)的傅立葉變換，其窗內(nèi)的信號必須是平穩(wěn)的12。2、HHT局部性能良好而且是自適應(yīng)的，對平穩(wěn)信號和非平穩(wěn)信號都能進行分析。它沒有固定的先驗基底，分解完全基于數(shù)據(jù)本身進行。固有模態(tài)函數(shù)是基于序列數(shù)據(jù)的時間特征尺度得出的不同的數(shù)據(jù)有不同的固有模態(tài)函數(shù)，每個固有模態(tài)函數(shù)可以認為是信號中固有的一個模態(tài)，所以通過Hilbert變換得到的瞬時頻率具有清晰的物理意義，能夠表達信號的局部特征。而傅立葉變換是以余弦函數(shù)為基底進行信