希爾伯特-黃變換說明及程序.docx

Hilbert-Huang 變換 希爾伯特-黃轉(zhuǎn)換希爾伯特-黃轉(zhuǎn)換（Hilbert-Huang Transform），由臺灣中央研究院院士黃鍔（Norden E. Huang）等人提出，將欲分析資料分解為本質(zhì)模態(tài)函數(shù)（intrinsic mode functions, IMF），這樣的分解流程稱為經(jīng)驗(yàn)?zāi)B(tài)分解(Empirical Mode Decomposition, EMD)的方法。然后將IMF作希爾伯特轉(zhuǎn)換(Hilbert Transform)，正確地獲得資料的瞬時頻率。此方法處理對象乃針對非穩(wěn)態(tài)與非線性訊號。與其他數(shù)學(xué)轉(zhuǎn)換運(yùn)算（如傅立葉變換）不同，希爾伯特-黃轉(zhuǎn)換算是一種應(yīng)用在數(shù)據(jù)資料上的算法，而非理論工具。質(zhì)模態(tài)函數(shù)(IMF)任何一個資料，滿足下列兩個條件即可稱作本質(zhì)模態(tài)函數(shù)。 局部極大值(local maxima)以及局部極小值(local minima)的數(shù)目之和必須與零交越點(diǎn)(zero crossing)的數(shù)目相等或是最多只能差1，也就是說一個極值后面必需馬上接一個零交越點(diǎn)。 在任何時間點(diǎn)，局部最大值所定義的上包絡(luò)線(upper envelope)與局部極小值所定義的下包絡(luò)線，取平均要接近為零。因此，一個函數(shù)若屬于IMF，代表其波形局部對稱于零平均值。此類函數(shù)類似于弦波（sinusoid-like），但是這些類似于弦波的部分其周期與振幅可以不是固定。因?yàn)?，可以直接使用希爾伯特轉(zhuǎn)換，求得有意義的瞬時頻率。經(jīng)驗(yàn)?zāi)B(tài)分解(EMD)EMD算法流程圖建立IMF是為了滿足希爾伯特轉(zhuǎn)換對于瞬時頻率的限制條件之前置處理，也是一種轉(zhuǎn)換的過程。我們將IMF來做希爾伯特轉(zhuǎn)換可以得到較良好的特性，不幸的是大部分的資料并不是IMF，而是由許多弦波所合成的一個組合。如此一來，希爾伯特轉(zhuǎn)換并不能得到正確的瞬時頻率，我們便無法準(zhǔn)確的分析資料。為了解決非線性（non-linear）與非穩(wěn)態(tài)（non-stationary）資料在分解成IMF時所遇到的困難，便發(fā)展出EMD。經(jīng)驗(yàn)?zāi)B(tài)分解是將訊號分解成IMF的組合。經(jīng)驗(yàn)?zāi)B(tài)分解是借著不斷重復(fù)的篩選程序來逐步找出IMF。以訊號為例，篩選程序的流程概述如下:步驟 1: 找出中的所有局部極大值以及局部極小值，接著利用三次樣條(cubic spline)，分別將局部極大值串連成上包絡(luò)線與局部極小值串連成下包絡(luò)線。步驟 2: 求出上下包絡(luò)線之平均，得到均值包絡(luò)線。步驟 3: 原始信號與均值包絡(luò)線相減，得到第一個分量。步驟 4: 檢查是否符合IMF的條件。如果不符合，則回到步驟1并且將當(dāng)作原始訊號，進(jìn)行第二次的篩選。亦即重復(fù)篩選次直到符合IMF的條件，即得到第一個IMF分量，亦即步驟 5: 原始訊號減去可得到剩余量，表示如下式步驟 6: 將當(dāng)作新的資料，重新執(zhí)行步驟1至步驟5，得到新的剩余量。如此重復(fù)次.當(dāng)?shù)趥€剩余量已成為單調(diào)函數(shù)(monotonic function)，將無法再分解IMF時，整個EMD的分解過程完成。原始訊號可以表示成個IMF分量與一個平均趨勢(mean trend)分量的組合，亦即如此一來，原始資料便分解成n個IMF和一個趨勢函數(shù)，我們便可將IMF做希爾伯特轉(zhuǎn)換來進(jìn)行瞬時頻率的分析。結(jié)論傅立葉變換是將一個訊號分解成無限多個弦波來分析資料，但是希爾伯特-黃轉(zhuǎn)換則是將一個訊號分解成數(shù)個近似于弦波的訊號（周期、振幅不固定）和一個趨勢函數(shù)來做分析。兩者各有其優(yōu)缺點(diǎn)，整理如下優(yōu)點(diǎn)：1.避免復(fù)雜的數(shù)學(xué)運(yùn)算2.可分析頻率會隨時間變化的訊號3.較適于分析氣候、經(jīng)濟(jì)等具有趨勢的資料4.可以找出一個函數(shù)的趨勢缺點(diǎn)：1.缺乏嚴(yán)謹(jǐn)?shù)奈锢硪饬x2.需要復(fù)雜的遞回，運(yùn)算時間反而比短時距傅立葉變換要長3.希爾伯特轉(zhuǎn)換未必能正確計(jì)算出本質(zhì)模態(tài)函數(shù)之瞬時頻率4.無法使用快速傅立葉變換5.只有在特例（組合較簡單的資料）時使用希爾伯特-黃轉(zhuǎn)換較快2.Matlab代碼/matlabcentral/fileexchange/196812.1 EMD分解代碼（emd.m）function imf = emd(x)% Empiricial Mode Decomposition (Hilbert-Huang Transform)% imf = emd(x)% Func : findpeaks x = transpose(x(:);imf = ;while ismonotonic(x) x1 = x; sd = Inf; while (sd 0.1) | isimf(x1) s1 = getspline(x1); s2 = -getspline(-x1); x2 = x1-(s1+s2)/2; sd = sum(x1-x2).2)/sum(x1.2); x1 = x2; end imfend+1 = x1; x = x-x1;endimfend+1 = x; % FUNCTIONS function u = ismonotonic(x) u1 = length(findpeaks(x)*length(findpeaks(-x);if u1 0, u = 0;else, u = 1; end function u = isimf(x) N = length(x);u1 = sum(x(1:N-1).*x(2:N) 1, u = 0;else, u = 1; end function s = getspline(x) N = length(x);p = findpeaks(x);s = spline(0 p N+1,0 x(p) 0,1:N);2.2 找極值代碼（findpeaks.m）function n = findpeaks(x)% Find peaks.% n = findpeaks(x) n = find(diff(diff(x) 0) x(n);n(u) = n(u)+1;2.3 繪制時-頻曲線以及尺度分解代碼（plot_hht.m）function plot_hht(x,Ts)% Plot the HHT.% plot_hht(x,Ts)% Ts ：time interval （sec） % : Syntax% The array x is the input signal and Ts is the sampling period.% Func : emd % Get HHT.imf = emd(x);for k = 1:length(imf) b(k) = sum(imfk.*imfk); th = angle(hilbert(imfk); dk = diff(th)/Ts/(2*pi);endu,v = sort(-b);b = 1-b/max(b); % Set time-frequency plots.N = length(x);c = linspace(0,(N-2)*Ts,N-1);for k = v(1:2) figure, plot(c,dk,k.,Color,b(k k k),MarkerSize,3); set(gca,FontSize,8,XLim,0 c(end),YLim,0 1/2/Ts); xlabel(Time), ylabel(Frequency);end % Set IMF plots.M = length(imf);N = length(x);c =