emd_online.m

一部分希尔伯特-黄变换（HHT）的经验模式分解（EMD）程序

function [imf,ort,nbit] = emd_online(x,t,stop,nbpresift,tst,tst2)

% EMD_ONLINE (On Line Empirical Mode Decomposition) computes on-line an EMD
%
% stopping criterion for sifting : 
%   at each point : mean amplitude < threshold*envelope amplitude 
%                   if mean amplitude > max(envelope amplitude)/tolerance
%   &
%   at each point : mean amplitude < threshold2*enveloppe amplitude 
%                   if mean amplitude > max(envelope amplitude)/tolerance2
%
%[imf,ort,nbits] = EMD_ONLINE(x,t,stop,nbpresift,tst,tst2)
% inputs: 
%         - x : analyzed signal 
%         - t (optional) : sampling times (default : 1:length(x))
%         - stop (optional) : threshold, and threshold2 (optional)
%                             tolerance, and tolerance2 (both optional)
%                             for sifting stopping criterion 
%                             default : [0.05,0.5,20,100]
%         - nbpresift (optional) : number of sifting by pieces iterations (default 4)
%         - tst (optional) : if equals to 1 shows sifting steps
%         - tst2 (optional) : if equals to 1 shows sifting by pieces steps                   
%
% 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
% 
% G. Rilling, July 2002


DEFSTOP = [0.05,0.5,20,100];% default parameters for sifting stop

NBPRESIFT = 4;%number of sifting iterations per block

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

if(nargin==2)
  stop = DEFSTOP;
  tst = 0;
  tst2 = 0;
end

if (nargin==3)
  tst=0;
  tst2 = 0;
end

if (nargin==4)
  tst=0;
  tst2 = 0;
end

if (nargin==5)
  tst2 = 0;
end

if nargin > 3
    NBPRESIFT = nbpresift;
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) > 4) | (S(2) > 4) | (length(S) > 2)
  error('stop must have only one row or one column of max four elements')
end

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

if S(2) < 4
    stop(4) = DEFSTOP(4);
end

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

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

if S(2) == 1
  stop=[stop, DEFSTOP(2)];
end


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

if tst
  figure
  figures(1) = gcf;
  figure
  figures(3) = gcf;
end

if tst2
  figure
  figures(2) = gcf;
end

  
MAXITERATIONS=10000;
LARGMIN = 5;

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

LARGTRANS = 10;
LARGTRANSPS = 5;

PAS = 20;

STEP = 5;% maximal number of iterations on a mode




LX = length(x);

% for display
sdt(LX) = 0;
sdt = sdt+sd;
sd2t(LX) = 0;
sd2t = sd2t+sd2;


% number of minima and maxima on the considered zone
lm = 0;
lM = 0;

% number of minima and maxima right of the considered zone, 
% after "stop" or "stopps"
lmr = 0;
lMr = 0;

% same, but left before "start"
lml = 0;
lMl = 0;

% total number of extrema, left and right
nem = 0;
nemr = 0;
neml = 0;

k = 1;
nbit = 0;

% number of modes, and number of modes on which block siftings are completed
nbmodes = 1;
nbmodes_psdone = 0;

start = 1;

% end of the constant part of the window
stop = min(PAS+1,LX);

% start and end of the considered zone
stopr = 1;
startl = 1;

% end of available data on the considered zone
fin = 1;

% start of the considered zone for block sifting
limpsl(1,1:NBPRESIFT) = 1;

% start and end of segment to which sifting is applied
startps = 1;
stopps(1,NBPRESIFT) = 0;
stopps = stopps + 1;

% end of available data for block sifting
finps(1,1:NBPRESIFT) = 1;
finps(1,1) = 10*PAS;

% tests if all data are available for an iteration of block sifting
lafinps(1,NBPRESIFT) = 0;

% allows to interrupt a mode extraction to process to the next mode
% interrupts also if not enough data available
suspps(1,NBPRESIFT) = 0;

% tests for the termination of one iteration of block sifting 
stoptestps(1,NBPRESIFT) = 0;

indmin = [];
indmax = [];

% tests if all data are available for on-line sifting
lafin = 0;

% tests if a mode is entirely extracted
stoptest = 0;

% allows to interrupt a mode extraction to process to the next mode
% interrupts also if not enough data available
susp = 0;

% tells if the considered zone has to be moved forward for having enough extrema
needextr = 1;

% idem for block sifting
needextrps(1,1:NBPRESIFT) = 1;

% tells how many iterations of block sifting have been initiated
nbstartedpresift = 1;

% modes concerned by block sifting
mps = x;

% mode concerned by on-line sifting
m(LX) = 0;

trig = 0;

if tst | tst2
  disp('appuyer sur une touche pour commencer')
  pause
end
  
while sum(stoptest) < nbmodes % global loop

  for k = 1:nbmodes

	  nsteps = 0;
    waittest = 0;

    if k == 1 & trig
        suspps(1,1) = 0;
        trig = 0;
    end
    
    
    while  sum(stoptestps(k,:)) < NBPRESIFT & ~waittest & sum(suspps(k,:)) < nbstartedpresift(k) % boucle de presifting
      for i = 1:nbstartedpresift(k)

          
	if needextrps(k,i) == 1
	  [indmintmp,indmaxtmp] = extr(mps(k,max([(limpsl(k,i)-1),1]):finps(k,i),i));
	  nb = sum(indmintmp > stopps(k,i))+sum(indmaxtmp > stopps(k,i));
	  stoprps(k,i) = finps(k,i);
      
	  if nb < 8*LARGTRANSPS & finps(k,i) < LX
          suspps(k,i) = 1;
          if k == 1 & i == 1
              finps(1,1) = finps(1,1) + 10*PAS;
              [indmintmp,indmaxtmp] = extr(mps(k,max([(limpsl(k,i)-1),1]):finps(k,i),i));
              nb = sum(indmintmp > stopps(k,i))+sum(indmaxtmp > stopps(k,i));
              stoprps(k,i) = finps(k,i);
              trig = 1;
          end
      else

	    lmt = length(indmintmp);
	    lMt = length(indmaxtmp);
	    if lmt > 0
	      indminps(k,1:lmt,i) = indmintmp + max([(limpsl(k,i)-1),1])-1;
	    end
	    if lMt > 0
	      indmaxps(k,1:lMt,i) = indmaxtmp + max([(limpsl(k,i)-1),1])-1;
	    end
	    if lmt < size(indminps,2)
	      indminps(k,length(indmintmp)+1:end,i) = 0;
	    end
	    if lMt < size(indmaxps,2)
	      indmaxps(k,length(indmaxtmp)+1:end,i) = 0;
	    end
	    needextrps(k,i) = 0;
	    
	  end
	  
	  if stoprps(k,i) >= LX
	    lafinps(k,i) = 1;
	    needextrps(k,i) = 0;
	  end
	end
	
	if ~suspps(k,i)
	  curindminps = indminps(k,find(indminps(k,:,i) >= limpsl(k,i)),i);
	  curindmaxps = indmaxps(k,find(indmaxps(k,:,i) >= limpsl(k,i)),i);
	  nemps = length(curindminps) + length(curindmaxps);
	  
	end
      
	% loop of block (pre)sifting
	while (~needextrps(k,i) | lafinps(k,i)) & ~stoptestps(k,i) & ~waittest & ~suspps(k,i)

	  if nemps < 3 & lafinps(k,i)
	    stoptestps(k,:) = 1;
	    stoptest(k) = 1;
	    m(k,:) = mps(k,:,i);
	    if i > 1
	      m(k+1,:) = x - sum(m(1:k,:));
	    end
	    break
	  end
	  
	  if limpsl(k,i) == 1
	    startps(k,i) = 1;
	  else
	    startps(k,i) = stopps(k,i);
	  end
	  if lafinps(k,i)
	    stopps(k,i) = LX;
	    stoptestps(k,i) = 1;
	  else
	    stopps(k,i) = min(curindminps(max([1,end - LARGTRANSPS+1])),curindmaxps(max([1,end - LARGTRANSPS+1]))); % si ~lafinps(k,i)
	  end
	  
      if startps(k,i) == stopps(k,i)
          pause   
          needextrps(k,i) = 1;
          break
      end
      
	  lmr = sum(curindminps > stopps(k,i));
	  lMr = sum(curindmaxps > stopps(k,i));
	  nemrps(k,i) = lmr + lMr;
	  
	  if nemrps(k,i) < 8*LARGTRANSPS
	    needextrps(k,i) = 1;
	  end
	  
	  
	  if limpsl(k,i) == 1
	    margeml = 0;
	    margeMl = 0;
	      
       if curindmaxps(1) < curindminps(1)
		if mps(k,1,i) > mps(k,curindminps(1),i)
	      lmax = fliplr(curindmaxps(2:min(end,NBSYM+1)));
	      lmin = fliplr(curindminps(1:min(end,NBSYM)));
	      lsym = curindmaxps(1);
		else
	      lmax = fliplr(curindmaxps(1:min(end,NBSYM)));
	      lmin = [fliplr(curindminps(1:min(end,NBSYM-1))),1];
	      lsym = 1;
		end
	      else
		if mps(k,1,i) < mps(k,curindmaxps(1),i)
	      lmax = fliplr(curindmaxps(1:min(end,NBSYM)));
	      lmin = fliplr(curindminps(2:min(end,NBSYM+1)));
	      lsym = curindminps(1);
		else
	      lmax = [fliplr(curindmaxps(1:min(end,NBSYM-1))),1];
	      lmin = fliplr(curindminps(1:min(end,NBSYM)));
	      lsym = 1;
		end
	      end
	      
	    tlmin = 2*t(lsym)-t(lmin);
	    tlmax = 2*t(lsym)-t(lmax);
	    
    % in case symmetrized parts do not extend enough
    if tlmin(1) > t(1) | tlmax(1) > t(1)
      if lsym == curindmaxps(1)
	lmax = fliplr(curindmaxps(1:min(end,NBSYM)));
      else
	lmin = fliplr(curindminps(1:min(end,NBSYM)));
      end
      if lsym == 1
	error('bug')
      end
      lsym = 1;
      tlmin = 2*t(lsym)-t(lmin);
      tlmax = 2*t(lsym)-t(lmax);
    end
	    	    
	  else % limpsl ~= 1
	    
	    lmin = curindminps(find(curindminps <= startps(k,i)));
	    lmax = curindmaxps(find(curindmaxps <= startps(k,i)));
	    if length(lmin) < 5 |length(lmax) < 5
	      error('souci')
	    end
	    tlmin = t(lmin);
	    tlmax = t(lmax);
	    margeml = length(lmin);
	    margeMl = length(lmax);
	    
	  end % if limpsl...
	    
	  if lafinps(k,i)
	    
	    margemr = 0;
	    margeMr = 0;
	      
	      if curindmaxps(end) < curindminps(end)
		if mps(k,LX,i) > mps(k,curindmaxps(end),i)
	rmax = [LX,fliplr(curindmaxps(max(end-NBSYM+2,1):end))];
	rmin = fliplr(curindminps(max(end-NBSYM+1,1):end));
	rsym = LX;
		else
	rmax = fliplr(curindmaxps(max(end-NBSYM+1,1):end));
	rmin = fliplr(curindminps(max(end-NBSYM,1):end-1));
	rsym = curindminps(end);
		end
	      else
		if mps(k,LX,i) < mps(k,curindminps(end),i)