emd_improve.m

用于信号处理的emd改进算法,供大家讨论学习

function [imf,r,ort,nbits] = EMD_CE(M,S,x,t)
% function for extrema-based sifting(Components from Extrema,CE)
% with M as the maximum number of sifting steps allowed,and S as the 
% stopping criterion;

% "A confidence limit for the  empirical mode decomposition and 
% Hilbert spectral analysis",The Royal Society,2003

if(nargin==3)
  t = 1:length(x);
  fs=1;      %采样频率:Hz
end
fs=1/(t(2)-t(1));%采样频率:Hz
%%%%%%%检验 x %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
T = size(x);
if ((T(1) > 1) & (T(2) > 1)) | (length(T) > 2)%如果x不止一行不止一列,或者三维三维以上
  error('x must have only one row or one column')
end
if T(1) > 1
  x = x';
end
%%%%%%%检验 t %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
T = size(t);
if ((T(1) > 1) & (T(2) > 1)) | (length(T) > 2)
  error('t must have only one row or one column')
end
if T(1) > 1
  t = t';
end
clear T;
%%%%%%% x与t %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
lx = length(x);
if (length(t)~=lx)
  error('x and t must have the same length')
end

% maximum number of iterations
MAXITERATIONS=M;
%%%%%%%%%%%%%%%%初始化%%%%%%%%%%%%%%%%
% maximum number of extrema and zero-crossings in residual
nem = lx;
nzm = lx;

r = x;
imf = [];
imftest=0;
k = 1;%%记录第k个IMF

% iterations counter for extraction of 1 mode
nbit=0;

% total iterations counter
NbIt=0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%抽取过程%%%%%%%%%%%%
m=r;
[mins,maxs,zers] = extrema(r);
lm=length(mins);
lM=length(maxs);
nem=lm + lM;
nzm=length(zers);

if imftest~=0;return;end

while nem > 2          
% current mode
m = r;

% mode at previous iteration
mp = m;

% tests if enough extrema to proceed
test = 0;
imftest=0;

% sifting loop
while (imftest<S) & (test == 0) 
    % boundary conditions for interpolations :
            if mins(1)<maxs(1);
                tc0=maxs(1)-mins(1);tlmax=mins(1)-tc0;mlmax=mean(m(maxs));
                if tlmax>=1;tlmax=1;end;
                if m(1)>m(maxs(1));mlmax=m(1);tlmax=1;end;
                tc1=mins(1)-tlmax;  tlmin=tlmax-tc1;  mlmin=mean(m(mins));
                if tlmin>=1;tlmin=1;end;
                if m(1)<m(mins(1));mlmin=m(1);tlmin=1;end;
            else
                tc0=mins(1)-maxs(1);tlmin=maxs(1)-tc0;mlmin=mean(m(mins)); 
                if tlmin>=1;tlmin=1;end;
                if m(1)<m(mins(1));mlmin=m(1);tlmin=1;end;
                tc1=maxs(1)-tlmin;  tlmax=tlmin-tc1;  mlmax=mean(m(maxs));
                if tlmax>=1;tlmax=1;end;
                if m(1)>m(maxs(1));mlmax=m(1);tlmax=1;end;
            end
            
            if mins(end)>maxs(end);
                tc0=mins(end)-maxs(end);trmin=mins(end)+tc0;mrmin=mean(m(mins));
                if trmin<=lx;trmin=lx;end;
                if m(end)<m(mins(end));mrmin=m(end);trmin=lx;end;
                tc1=trmin-maxs(end);    trmax=trmin+tc1;    mrmax=mean(m(maxs));
                if trmax<=lx;trmax=lx;end;
                if m(end)>m(maxs(end));mrmax=m(end);trmax=lx;end;
            else
                tc0=maxs(end)-mins(end);trmax=maxs(end)+tc0;mrmax=mean(m(maxs));
                if trmax<=lx;trmax=lx;end;
                if m(end)>m(maxs(end));mrmax=m(end);trmax=lx;end;
                tc1=trmax-mins(end);    trmin=trmax+tc1;    mrmin=mean(m(mins));
                if trmin<=lx;trmin=lx;end;
                if m(end)<m(mins(end));mrmin=m(end);trmin=lx;end;
            end
            
            % definition of envelopes from interpolation
            envmax = spline(([tlmax maxs trmax])/fs,[mlmax m(maxs) mrmax],t);
            envmin = spline(([tlmin mins trmin])/fs,[mlmin m(mins) mrmin],t);
            envmoy = (envmax + envmin)/2;
            m = m - envmoy;
            [mins,maxs,zers] = extrema(m);
            lm=length(mins);
            lM=length(maxs);
            nem = lm + lM;
            nzm = length(zers);
            
            %判断m是否为IMF
            if abs(nem-nzm)<=1
                imftest=1+imftest;
            else
                imftest=0;
            end
            
            % end loop : stops if not enough extrema
            if nem < 3
                test = 1;
            end
            nbit=nbit+1; 
            NbIt=NbIt+1;
            
            if(nbit==(MAXITERATIONS-1))
                warning(['forced stop of sifting : too many iterations... mode ',int2str(k)]);%,'. stop parameter mean value : ',num2str(s)])
                test=1;
            end
        end
        
        imf(k,:) = m;
        
        nbits(k) = nbit;
        k = k+1;
        r = r - m;
        [mins,maxs,zers] = extrema(r);
        nem = length(mins) + length(maxs);
        nzm = length(zers);
        nbit=1;
        
        if (max(r) - min(r)) < (1e-10)*(max(x) - min(x))
            if nem > 2
                warning('forced stop of EMD : too small amplitude')
            else
                disp('forced stop of EMD : too small amplitude')
            end
            break;
        end
    end
    %end;
%%%%%%%%%%%%抽取过程结束%%%%%%%%%%%%
ort = io(x,imf);
return;