www.pudn.com > hilbert.zip > emd_lily02.m, change:2008-05-14,size:6902b


% EMD.M 
%
% G. Rilling, July 2002
%
% computes EMD (Empirical Mode Decomposition) according to:
%
% N. E. Huang et al., "The empirical mode decomposition and the 
% Hilbert spectrum for non-linear and non stationary time series analysis,"  
% Proc. Royal Soc. London A, Vol. 454, pp. 903-995, 1998
%
% with variations reported in:
%
% G. Rilling, P. Flandrin and P. Gon鏰lv鑣
% "On Empirical Mode Decomposition and its algorithms"
% IEEE-EURASIP Workshop on Nonlinear Signal and Image Processing
% NSIP-03, Grado (I), June 2003
%
% stopping criterion for sifting : 
%   at each point : mean amplitude < threshold2*envelope amplitude
%   &
%   mean of boolean array ((mean amplitude)/(envelope amplitude) > threshold) < tolerance
%   &
%   |#zeros-#extrema|<=1
%
% inputs:  - x : analysed signal (line vector)
%          - t (optional) : sampling times (line vector) (default : 1:length(x))
%          - stop (optional) : threshold, threshold2 and tolerance (optional)
%                              for sifting stopping criterion 
%                              default : [0.05,0.5,0.05]
%          - tst (optional) : if equals to 1 shows sifting steps with pause
%                             if equals to 2 no pause
%
% outputs: - imf : intrinsic mode functions (last line = residual)
%          - ort : index of orthogonality
%          - nbits : number of iterations for each mode
%
% calls:   - extr (finds extrema and zero-crossings)
%          - io : computes the index of orthogonality

function [imf,ort,nbits] = emd_lily0(x,t,stop,tst);
fs=0.001;
% default for stopping
defstop = [0.05,0.5,0.05];

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

if(nargin==2)
  stop = defstop;
  tst = 0;
end

if (nargin==3)
  tst=0;
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);

if tst
  figure
end

% 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 (lm==1)&(lM==1)
            k1=2*abs(indmax(1)-indmin(1));
        end
         if indmax(1) < indmin(1)
         k1=indmax(2)-indmax(1);
         else if indmax(1) > indmin(1)
         k1=indmin(2)-indmin(1);
             end
         end
         tlmax2=t(indmax(1))-2*k1*(1/fs);
         mlmax2=m(indmax(1));
         tlmax1=t(indmax(1))-k1*(1/fs);
         mlmax1=m(indmax(1));
         tlmin2=t(indmin(1))-2*k1*(1/fs);
         mlmin2=m(indmin(1));
         tlmin1=t(indmin(1))-k1*(1/fs);
         mlmin1=m(indmin(1));
         
         if (lm==1)&(lM==1)
            k2=2*abs(indmax(1)-indmin(1));
        end
         if indmax(lM) > indmin(lm)
         k2=indmax(lM)-indmax(lM-1);
         else if indmax(lM) < indmin(lm)
         k2=indmin(lm)-indmin(lm-1);
             end
         end
         trmax2=t(indmax(lM))+2*k2*(1/fs);
         mrmax2=m(indmax(lM));
         trmax1=t(indmax(lM))+k2*(1/fs);
         mrmax1=m(indmax(lM));
         trmin2=t(indmin(lm))+2*k2*(1/fs);
         mrmin2=m(indmin(lm));
         trmin1=t(indmin(lm))+k2*(1/fs);
         mrmin1=m(indmin(lm));
    
   
    % definition of envelopes from interpolation
        
    envmax = interp1([tlmax2 tlmax1 t(indmax) trmax1 trmax2],[mlmax2 mlmax1 m(indmax) mrmax1 mrmax2],t,'spline');
    envmin = interp1([tlmin2 tlmin1 t(indmin) trmin1 trmin2],[mlmin2 mlmin1 m(indmin) mrmin1 mrmin2],t,'spline');

    envmoy = (envmax + envmin)/2;

    m = m - envmoy; %极值对称延拓法
   
    [indmin,indmax,indzer] = extr(m);
    lm=length(indmin);
    lM=length(indmax);
    nem = lm + lM;
    nzm = length(indzer);

    % evaluation of mean zero
    sx=2*(abs(envmoy))./(abs(envmax-envmin));
    s = mean(sx);
    
    % display
        
     if tst 
      subplot(4,1,1)
      plot(t,mp);hold on;
      plot(t,envmax,'--k');plot(t,envmin,'--k');plot(t,envmoy,'r');

      title(['IMF ',int2str(k),';   iteration ',int2str(nbit),' before sifting']);
      set(gca,'XTick',[])
      hold  off

      subplot(4,1,2)
      plot(t,sx)
      hold on
      plot(t,sdt,'--r')
      plot(t,sd2t,':k')
      title('stop parameter')
      set(gca,'XTick',[])
      hold off

      subplot(4,1,3)
      plot(t,m)
      title(['IMF ',int2str(k),';   iteration ',int2str(nbit),' after sifting']);
      set(gca,'XTick',[])

      subplot(4,1,4);
      plot(t,r-m)
      title('residue');
      disp(['stop parameter mean value : ',num2str(s)])
      if tst == 2
        pause(0.01)
      else
        pause
      end
      
    end
 
    % end loop : stops if not enough extrema
    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;
  if tst
    disp(['mode ',int2str(k),' enregistre'])
  end
  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;
n=size(imf);
n=n(1);
ort = io(x,imf);

if tst
  close
end