EMD 程序
function imf = emd(x)
% Empiricial Mode Decomposition (Hilbert-Huang Transform)
% imf = emd(x)
% Func : findpeaks
x = transpose(x(:));%转置
imf = [];
while ~ismonotonic(x) %当 x 不是单调函数，分解终止条件
x1 = x;
sd = Inf;%均值
%直到 x1 满足 IMF 条件，得 c1
while (sd > 0.1) | ~isimf(x1) %当标准偏差系数 sd 大于 0.1 或 x1 不是固有模态函数时，分量终止条
件
s1 = getspline(x1);%上包络线
s2 = -getspline(-x1);%下包络线
x2 = x1-(s1+s2)/2;%此处的 x2 为文章中的 h
sd = sum((x1-x2).^2)/sum(x1.^2);
x1 = x2;
end
imf{end+1} = x1;
x = x-x1;
end
imf{end+1} = x;
% FUNCTIONS
function u = ismonotonic(x)
%u=0 表示 x 不是单调函数，u=1 表示 x 为单调的
u1 = length(findpeaks(x))*length(findpeaks(-x));
if u1 > 0, u = 0;
else, u = 1; end
function u = isimf(x)
%u=0 表示 x 不是固有模式函数，u=1 表示 x 是固有模式函数
N = length(x);
u1 = sum(x(1:N-1).*x(2:N) < 0);
u2 = length(findpeaks(x))+length(findpeaks(-x));
if abs(u1-u2) > 1, u = 0;
else, u = 1; end
function s = getspline(x)
%三次样条函数拟合成元数据包络线
N = length(x);
p = findpeaks(x);
评论0