www.pudn.com > HHT_power-system_power-quality_disturbances-detect > emdyy.m, change:2010-04-26,size:1133b

```function imf = emd(x)    %EMD

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;

sd = sum((x1-x2).^2)/sum(x1.^2);
x1 = x2;
end

imf{end+1} = x1;
x          = x-x1;
end
imf{end+1} = x;

function u = ismonotonic(x)

u1 = length(findpeaks(x))*length(findpeaks(-x));
if u1 > 0, u = 0;
else      u = 1; end            %u=0表示x不是单调函数，u=1表示x为单调的

function u = isimf(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    %u=0表示x不是固有模式函数，u=1表示x是固有模式函数

function s = getspline(x)       %三次样条函数拟合成元数据包络线

N = length(x);
p = findpeaks(x);
s = spline([0 p N+1],[0 x(p) 0],1:N);
```