其中ξ(t)就是原始信號(hào)，IMFi就是K個(gè)固有模態(tài)函數(shù)。rK就是原始信號(hào)減完IMF后剩下的余項(xiàng)。
下面是求解EMD算法的Matlab源程序，特此聲明，本程序?yàn)槲冶救嗽诰W(wǎng)上找到的，除了注釋外，其他版權(quán)皆歸屬原作者，由于不清楚原作者是誰(shuí)，未能標(biāo)出，如果侵犯權(quán)利，請(qǐng)聯(lián)系我刪除源碼。
%非主函數(shù)，被調(diào)用
function n = findpeaks(x)%用于尋找極值點(diǎn)，該函數(shù)只會(huì)求極大值
% Find peaks.
% n = findpeaks(x)
n = find(diff(diff(x)>0)<0);%一階導(dǎo)數(shù)大于0二階導(dǎo)數(shù)小于0的點(diǎn)
u = find(x(n+1)>x(n));
n(u) = n(u) + 1;
end
%非主函數(shù)，被調(diào)用<br>%判斷x是否單調(diào)，返回0代表不是單調(diào)，返回1代表是單調(diào)
function u = ismonotonic(x)
u1 = length(findpeaks(x))*length(findpeaks(-x));%如果最大/最小值有一個(gè)為0即可判斷程序滿足退出條件了
if u1 > 0
u = 0;
else
u = 1;
end
end
%非主函數(shù)，被調(diào)用。判斷當(dāng)前x是不是真IMF
function u = isimf(x)
N = length(x);
u1 = sum(x(1:N-1).*x(2:N) < 0);%求x與y=0軸交點(diǎn)的個(gè)數(shù)
u2 = length(findpeaks(x))+length(findpeaks(-x));%求極值點(diǎn)個(gè)數(shù)
if abs(u1-u2) > 1
u = 0;
else
u = 1;
end
end
%非主函數(shù)，被調(diào)用，作用是獲得x的包絡(luò)線
function s = getspline(x)
N = length(x);
p = findpeaks(x);
s = spline([0 p N+1],[0 x(p) 0],1:N);
end
%主函數(shù)
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算法，在這一部分中，我會(huì)介紹Hilbert算法，這一節(jié)有些許數(shù)學(xué)趣味，對(duì)數(shù)學(xué)趣味不感興趣的直接跳到應(yīng)用部分。