matlab中emd未定义,EMD的Matlab程序

emd.m文件

function imf = emd(x)

% Empiricial Mode Decomposition (Hilbert-Huang Transform)

% EMD分解或HHT变换

% 返回值为cell类型，依次为一次IMF、二次IMF、...、最后残差

x   = transpose(x(:));            %将x转置为行矩阵

imf = [];

while ~ismonotonic(x)            %当x(残余函数)不是单调函数时，分解终止条件

x1 = x;

sd = Inf;

while (sd > 0.1) || ~isimf(x1)    %当标准偏差系数sd大于0.1或X1不是固有模态函数时，

分量终止条件(类柯西收敛准则：当sd介于0.2—0.3之

间时，筛选过程停止)

s1 = getspline(x1);        % 极大值点样条曲线

s2 = -getspline(-x1);       % 极小值点样条曲线

x2 = x1-(s1+s2)/2;        %h(t)=x(t)-m(t)其中m(t)为包络线的平均值

sd = sum((x1-x2).^2)/sum(x1.^2);  %筛选停止准则

x1 = x2;  %x2即h(t)

end

imf{end+1} = x1;                  %n=n+1   Cn=h(t),得imf分量

x = x-x1;                         %r=x(t)-Cn

end

imf{end+1} = x;

% 是否单调

function u = ismonotonic(x)            %判断单调性函数

u1 = length(findpeaks(x))*length(findpeaks(-x));

if u1 > 0

u = 0;                           %u = 0表示x不是单调函数

else

u = 1;                           %u =1表示x是单调函数

end

% 是否IMF分量

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   %imf条件：(1)过零点数量与极值点数量必须相等或最多相差不能多于

一个。

(2)极大值、极小值包络均值为零u(t)+m(t)=0

u = 0;                   %u=0表示x不是固有模态函数

else

u = 1;                   %u=1表示x是固有模态函数

end

% 据极大值点构造样条曲线

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

N = length(x);

p = findpeaks(x);

s = spline([0 p N+1],[0 x(p) 0],1:N);  %yy=spline[X,Y,xx]用三次样条的方法得到yy