emd.asv

经改进的希尔伯特——黄分解

function [imf,res] = emd(x,t,stop);

% default for stopping
defstop = [0.05,0.5,0.05];

if(nargin==1)
  t = 1:length(x);
  stop = defstop;

end

if(nargin==2)
  stop = defstop;

end



S = size(x);
if ((S(1) > 1) & (S(2) > 1)) | (length(S) > 2)
  error('x must have only one row or one column')
end

if S(1) > 1
  x = x';
end

S = size(t);
if ((S(1) > 1) & (S(2) > 1)) | (length(S) > 2)
  error('t must have only one row or one column')
end

if S(1) > 1
  t = t';
end

if (length(t)~=length(x))
  error('x and t must have the same length')
end

S = size(stop);
if ((S(1) > 1) & (S(2) > 1)) | (S(1) > 3) | (S(2) > 3) | (length(S) > 2)
  error('stop must have only one row or one column of max three elements')
end

if S(1) > 1
  stop = stop';
  S = size(stop);
end

if S(2) < 3
  stop(3)=defstop(3);
end

if S(2) < 2
  stop(2)=defstop(2);
end

sd = stop(1);
sd2 = stop(2);
tol = stop(3);


% maximum number of iterations
MAXITERATIONS=2000;

% maximum number of symmetrized points for interpolations
NBSYM = 2;

lx = length(x);

sdt(lx) = 0;
sdt = sdt+sd;
sd2t(lx) = 0;
sd2t = sd2t+sd2;

% maximum number of extrema and zero-crossings in residual
ner = lx;
nzr = lx;

r = x;
imf = [];
k = 1;

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

% total iterations counter
NbIt=0;

while ner > 2
        
  % current mode
  m = r;
  
  % mode at previous iteration
  mp = m;
  
  sx = sd+1;
  
  % tests if enough extrema to proceed
  test = 0;
  
  [indmin,indmax,indzer] = extr(m);
  lm=length(indmin);
  lM=length(indmax);
  nem=lm + lM;
  nzm=length(indzer);
  
  j=1;
  
  % sifting loop
  while ( mean(sx > sd) > tol | any(sx > sd2) | (abs(nzm-nem)>1)) & (test == 0) & nbit<MAXITERATIONS
    
    if(nbit>MAXITERATIONS/5 & mod(nbit,floor(MAXITERATIONS/10))==0)
      disp(['mode ',int2str(k),' nombre d iterations : ',int2str(nbit)])
      disp(['stop parameter mean value : ',num2str(s)])
    end
   
    % boundary conditions for interpolations :
 
if  indmax(1)<indmin(1)           %left 
    k1=indmax(2)-indmax(1)
elseif indmax(1)>indmin(1)
    k1=indmin(2)-indmin(1)
elseif length(indmax)==length(indmin)==1
    k1=2*abs(indmax(1)-indmin(1))
end
    
    tt1=t(indmax(1))-2*add*k1       %left extend points
    uu1=x(indmam(1))
    tt2=t(indmax(1))-add*k1
    uu2=m(indmax(1))
    ttt1=t(indmin(1))-2*k1*add
    vv1=m(indmin(1))
    ttt2=t(indmin(1))-k1*add
    vv2=m(indmin(1))
     
    
if indmin(end)<indmax(end)          %right
    k2=indmax(end)-indmax(end-1)
elseif indmax(end)<indmin(end)
    k2=indmin(end)-indmin(end-1)
elseif length(indmax)==length(indmin)==1
    k2=2*abs(indmax(end)-indmin(end))
end
    
    tt3=t(indmax(end))+add*k2     %right extend points
    uu3=m(indmax(end))
    tt4=t(indmax(end))+2*add*k2
    uu4=m(indmax(end))
    ttt3=t(indmin(end))+add*k2
    vv3=m(indmin(end))
    ttt4=t(indmin(end))+2*add*k2
    vv4=m(indmin(end))
    
    %对数据进处理
 
    if m(1)>m(indmax(1))
        tt2=t(1)
        uu2=m(1)
    end
    if m(1)<m(indmin(1))
        ttt2=t(1)
        vv2=m(1)
    end
    if m(end)>m(indmax(end))
        tt3=t(end)
        uu3=m(end)
    end
    if m(end)<m(indmin(end))
        ttt3=t(end)
        vv3=m(end)
    end  
        ttmin=[ttt1 ttt2 t(indmin) ttt3 ttt4]
        ttmax=[tt1 tt2 t(indmax) tt3 tt4]
        ttxmin=[vv1 vv2 m(indmin) vv3 vv4]
        ttxmax=[uu1 uu2 m(indmax) uu3 uu4]
    
    maxbl= interp1(ttmax,ttxmax,t,'spline');	
    minbl = interp1(ttmin,ttxmin,t,'spline');
    midbl = (maxbl + minbl)/2;

    m = m - midbl;
   

    % evaluation of mean zero
     sx=2*(abs(midbl))./(abs(maxbl-minbl));
    s = mean(sx);
   % display_sift(t,r,m,mp,k,sx,s,nbit,sdt,sd2t,minbl,maxbl, midbl)
    if nem < 3
      test = 1;
    end

    mp = m;
    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)])
    end
  
  end
  imf(k,:) = m;
  
  nbits(k) = nbit;
  k = k+1;
  r = r - m;
  [indmin,indmax,indzer] = extr(r);
  ner = length(indmin) + length(indmax);
  nzr = length(indzer);
  nbit=1;

  if (max(r) - min(r)) < (1e-10)*(max(x) - min(x))
    if ner > 2
      warning('forced stop of EMD : too small amplitude')
    else

      disp('forced stop of EMD : too small amplitude')
    end
    break
  end
  
end

imf(k,:) = r;