希爾伯特-黃變換
https://blog.csdn.net/lvsehaiyang1993/article/details/80875679
轉載
在我們正式開始講解Hilbert-Huang Transform之前，不妨先來了解一下這一偉大算法的兩位發明人和這一算法的應用領域
Section I 人物簡介
希爾伯特：公認的數學界“無冕之王”，1943年去世于瑞士蘇黎世。除此之外，自不必過多介紹。
黃鍔：1937年出生于湖北省；1975年進入NASA（美國國家宇航局）；美國國家工程院院士。
Section II Hilbert-Huang的應用領域
醫學領域：探測心率不齊、登革熱的擴散、血壓的變化
交通領域：探測公路橋梁安全
安全領域：辨識發言者的身份
地理領域：地震工程
航天領域：衛星資料分析
在了解了這一偉大發明的背景后，下面我們要正式的開始入手希爾伯特-黃變換了，我將嘗試以盡可能簡要的語言向大家介紹這一發明，并盡可能的避免不必要的數學推導。
Section III Hilbert-Huang的算法詳細介紹
如下圖所示，在希爾伯特-黃的運算步驟中，原始腦電信號/其他時序信號被作為Huang的算法的輸入，在經過huang的算法處理過后被當做Hilbert的輸入進行處理。這便是Hilbert-Huang最簡單明了的運算步驟。在這里為了繼續往下的講解更加方便，我們先來介紹兩個概念。上文中，我們提到了“huang的算法”，在正式的書面語言中，我們并不這么稱呼它，而是將“huang的算法”稱為EMD(Empirical mode decomposition,經驗模式分解)。而另外一個概念IMF在這里直接講解或許會使大家暈頭轉向（或許有人注意到圖中的IMF后面有一個‘s’，而這里卻沒有加‘s’，對英語只有基礎了解的人也應猜到IMF不止一個）。
當將到這里的時候，大部分人或許會萌生出一個念頭——“難道Huang僅僅是對Hilbert的錦上添花嗎”。好吧，至少本人當初就是這樣想的，畢竟Hilbert比huang更早出名，而且Hilbert是數學史上公認的“大牛”，哦，不對，是“大王”。用當前時興的話來說就是“huang有可能抱了Hilbert的大腿”。但當我真正了解了這一偉大的發明之后，我才徹徹底底打消了這個十分愚蠢的念頭。
我個人并不喜歡吊人胃口，這里把結論說在前面“Huang的算法幾乎是Hilbert使用的前提條件，Hilbert Transform則是Hilbert-Huang算法的精要所在”（注意句中出現了“幾乎”一詞）。下面我就給大家講一下這句話的由來。比如我們造了一款叫做“榔頭”的手機，“榔頭”手機對用戶的使用提出了下列要求：1.晚上不能使用。2.下雨天不能打。3.室內不能打。4.室外的偏遠郊區也不能打。實際上，Hilbert正是這樣一款“榔頭”手機，它對用戶的使用提出了近乎苛刻的要求。Hilbert變換算法要求輸入信號只能是線性穩態的。請注意這里是兩個詞“線性”“穩態”。無論是在自然界還是在人類社會中，絕大部分的信號要么是“線性非穩態”，要么是“非線性穩態”，要么干脆是“非線性非穩態”。我們關心的重點——EEG信號正是這樣一類“非線性非穩態”的信號。這也就導致了絕大部分信號不能夠愉快的進入Hilbert的“碗里”來。此時，Huang的EMD算法起到了這樣的作用，它能夠將所有的時域信號轉化為“線性穩態”，解了Hilbert算法的軟肋。
首先，我們先說一說Huang的EMD算法。為了講解清晰起見，我將對照下圖予以講解：
上圖中，深藍色的線條是EEG信號（截取自瑞士聯邦理工學院DEAP數據庫 s01 trail1 channel1的前200個數據點）。圖中，紅線上的紅點是該EEG signal的極大值點，綠線上的綠點是該signal的極小值點。我們分別為極大值點和極小值點做三次包絡線做好的包絡線分別是紅色包絡線和綠色包絡線兩條線。為這兩條線做出均值線即為圖中圍繞y=0軸（注意y=0軸的位置，并非是圖中的坐標系的x軸，x軸所代表的線是y=-10）震蕩的淺藍色均值線。之后，我們將原始EEG信號減去均值線，得到疑似IMF線（圖中未標出）（這里終于出現IMF了，可是我還是無法讓你直觀理解，大家先暫且忍忍，強行記憶一下）。之后，我們對疑似IMF進行判斷（需要同時滿足兩個條件，下面講）。如果滿足條件，則疑似IMF升級為正式IMF。然后將原始信號減去正式IMF的結果賦值給原始信號，說白了就是讓這一IMF從原始信號里“滾蛋”。就好比蒙面歌王的某一期的獲勝者一樣，都贏了，不滾干嘛，難道還要和加時賽選手（減完后剩下的原始信號，也即新原始信號）在一起嗎？另一個方面，如果疑似IMF未能通過檢驗，則將當前IMF作為原始信號，并回到做極值點的包絡線那一步重新開始。現在講一下兩個重要的條件：
條件1：均值線（總得有很多數構成吧）的平均值趨近于0（一般和0做差<0.1）
條件2：原始信號的極值點個數（包括極大值點個數+極小值點個數）和原始信號同y=0的交點個數之差不能大于1（小于等于1）
那么這樣一個程序什么時候可以循環結束呢，答案是，當某一次IMF被發現是單調函數或者是缺少極大/小值點即可讓程序結束。下圖是程序流程圖。
空口無憑，我們處理一段真正的腦電試試看（程序會在之后給出）
圖中共有4*2個圖，位于（1,1）這個位置的是腦電的原始信號。之后從（1,2）->(4,2)均為IMF。其中，除了（1,1）自身，每一副圖都是（1,1）的一個IMF（現在知道什么是IMF了吧）。通過觀察不難發現。一個典型的IMF分量的上下包絡線肯定是對稱的。最后一個IMF（4,2）被稱為余項用r表示。觀察即可知該IMF分量沒有極小值點（端點除外），所以程序才會結束。通常來講，別的書上會這樣用數學公式告訴你：
其中ξ(t)就是原始信號，IMFi就是K個固有模態函數。rK就是原始信號減完IMF后剩下的余項。
下面是求解EMD算法的Matlab源程序，特此聲明，本程序為我本人在網上找到的，除了注釋外，其他版權皆歸屬原作者，由于不清楚原作者是誰，未能標出，如果侵犯權利，請聯系我刪除源碼。
%非主函數，被調用
function n = findpeaks(x)%用于尋找極值點，該函數只會求極大值
% Find peaks.
% n = findpeaks(x)
n = find(diff(diff(x)>0)<0);%一階導數大于0二階導數小于0的點
u = find(x(n+1)>x(n));
n(u) = n(u) + 1;
end
%非主函數，被調用<br>%判斷x是否單調，返回0代表不是單調，返回1代表是單調
function u = ismonotonic(x)
u1 = length(findpeaks(x))*length(findpeaks(-x));%如果最大/最小值有一個為0即可判斷程序滿足退出條件了
if u1 > 0
u = 0;
else
u = 1;
end
end
%非主函數，被調用。判斷當前x是不是真IMF
function u = isimf(x)
N = length(x);
u1 = sum(x(1:N-1).*x(2:N) < 0);%求x與y=0軸交點的個數
u2 = length(findpeaks(x))+length(findpeaks(-x));%求極值點個數
if abs(u1-u2) > 1
u = 0;
else
u = 1;
end
end
%非主函數，被調用，作用是獲得x的包絡線
function s = getspline(x)
N = length(x);
p = findpeaks(x);
s = spline([0 p N+1],[0 x(p) 0],1:N);
end
%主函數
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
imf(end+1,:) = x1;
x = x-x1;
end
imf(end+1,:) = x;
end
Section IV Hilbert算法的介紹
在上一章中，我們介紹了EMD算法，在這一部分中，我會介紹Hilbert算法，這一節有些許數學趣味，對數學趣味不感興趣的直接跳到應用部分。
由最后一步可以知道，當頻率大于0時，相位向左移90度；反之，向右移90度。這便是希爾伯特變換。
一般來講，對于原始信號x(t)的希爾伯特變換H[x(t)]，通常被寫為
z(t)=x(t)+j H[x(t)]
其中，x(t)被稱為復信號z(t)的實部，H[x(t)]被稱為復信號z(t)的虛部， z(t)被稱為x(t)的解析信號
一般情況下，matlab會將z(t)給出，而不直接給出原始信號的希爾伯特變換，所以需要使用imag函數求解z(t)的虛部，這才是真正的希爾伯特變換。
|
QQ好友和群
QQ空間
騰訊微博
騰訊朋友
收藏
淘帖
頂
踩
