emd_online.m

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

 	rmax = fliplr(curindmaxps(max(end-NBSYM+1,1):end));
	rmin = [LX,fliplr(curindminps(max(end-NBSYM+2,1):end))];
	rsym = LX;
		else         
	rmax = fliplr(curindmaxps(max(end-NBSYM,1):end-1));
	rmin = fliplr(curindminps(max(end-NBSYM+1,1):end));
	rsym = curindmaxps(end);
		end
	      end
	      
	    trmin = 2*t(rsym)-t(rmin);
	    trmax = 2*t(rsym)-t(rmax);

    if trmin(end) < t(LX) | trmax(end) < t(LX)
      if rsym == indmax(end)
	rmax = fliplr(curindmaxps(max(end-NBSYM+1,1):end));
      else
	rmin = fliplr(curindminps(max(end-NBSYM+1,1):end));
      end
      if rsym == LX
 	error('bug')
      end
      rsym = LX;
      trmin = 2*t(rsym)-t(rmin);
      trmax = 2*t(rsym)-t(rmax);
    end

        
	  else % lafin(k,i) ~= 1
	    
	    rmin = curindminps(end - LARGTRANSPS + 1:end);
	    rmax = curindmaxps(end - LARGTRANSPS + 1:end);
	    trmin = t(rmin);
	    trmax = t(rmax);
	    margemr = length(rmin);
	    margeMr = length(rmax);
	    
    end
	  
	  mpslmax = mps(k,lmax,i);
	  mpslmin = mps(k,lmin,i);
	  mpsrmax = mps(k,rmax,i);
	  mpsrmin = mps(k,rmin,i);
	  
	  
	  envmax = interp1([tlmax t(curindmaxps((margeMl+1):(end-margeMr))) trmax],[mpslmax mps(k,curindmaxps((margeMl+1):(end-margeMr)),i) mpsrmax],t(startps(k,i):stopps(k,i)),'spline');
	  envmin = interp1([tlmin t(curindminps((margeml+1):(end-margemr))) trmin],[mpslmin mps(k,curindminps((margeml+1):(end-margemr)),i) mpsrmin],t(startps(k,i):stopps(k,i)),'spline');
	  
	  envmoy = (envmax + envmin)/2;
	  

	  
	  if i == NBPRESIFT
	    m(k,startps(k,i):stopps(k,i)) = mps(k,startps(k,i):stopps(k,i),i) - envmoy;
	    fin(k) = stopps(k,i);
	    susp(k) = 0;
	    if startps(k,i) == 1
	      nbmodes_psdone = k;
	    end
	  else
	    if nbstartedpresift(k) < i+1
	      nbstartedpresift(k) = i+1;
	    end
	    mps(k,startps(k,i):stopps(k,i),i+1) = mps(k,startps(k,i):stopps(k,i),i) - envmoy;
	    finps(k,i+1) = stopps(k,i);
	    suspps(k,i+1) = 0;
	  end
	  
	  limpsl(k,i) = min([curindminps(max(end - 15*LARGTRANSPS,1)),curindmaxps(max(end - 15*LARGTRANSPS,1))]);
	  
	  % display
	  if tst
	    figure(figures(3))
	    subplot(2*NBPRESIFT,1,2*i-1)
	    plot(t(stopps(k,i):finps(k,i)),mps(k,stopps(k,i):finps(k,i),i))
	    hold on
	    plot(t(1:startps(k,i)),mps(k,1:startps(k,i),i),'k')
	    plot(t(startps(k,i):stopps(k,i)),mps(k,startps(k,i):stopps(k,i),i),'r')
	    plot(t(startps(k,i):stopps(k,i)),envmax,'k--')
	    plot(t(startps(k,i):stopps(k,i)),envmin,'k--')
	    axis([1,LX,min(mps(k,1:finps(k,i),i)),max(mps(k,1:finps(k,i),i))])
	    hold off
	    title(['mode ',int2str(k),'. Presifting # ',int2str(i)])
	    subplot(2*NBPRESIFT,1,2*i)
	    if i == NBPRESIFT
	      plot(t(startps(k,i):finps(k,i)),m(k,startps(k,i):finps(k,i)))
	      hold on
	      plot(t(1:startps(k,i)),m(k,1:startps(k,i)),'k')
	      hold off
	      axis([1,LX,min(m(k,1:finps(k,i))),max(m(k,1:finps(k,i)))])
	    else
	      plot(t(startps(k,i):stopps(k,i)),mps(k,startps(k,i):stopps(k,i),i+1))
	      hold on
	      plot(t(1:startps(k,i)),mps(k,1:startps(k,i),i+1),'k')
	      plot(t(stopps(k,i):finps(k,i+1)),mps(k,stopps(k,i):finps(k,i+1),i+1))
	      hold off
	      axis([1,LX,min(mps(k,1:finps(k,i+1),i+1)),max(mps(k,1:finps(k,i+1),i+1))])
	    end
	    pause(0.01)
	  end
	    
	    
	  
	  
	  
	end % petit while
      end % for
    end % grand while 
    
    
    % interrupts loop on a mode for processing to the next one    
    waittest = 0;
    
    
    
    while ~stoptest(k) & ~waittest & ~susp(k)	

	
      if needextr(k)
	if fin(k) > 0
	  stopr(k) = fin(k);
	  [indm,indM] = extr(m(k,max([(startl(k)-1),1]):stopr(k)));
	else
	  indm = [];
	  indM  = [];
	  susp(k) = 1;
	end
    
	
	if sum([(indm > stop(k)),(indM > stop(k))]) < 4*LARGTRANS & fin(k) < LX %stopr(k) >= fin(k) & fin(k) < LX
	  susp(k) = 1;
	else
	  startr(k) = stopr(k) + 1;
	  [indmintmp,indmaxtmp] = extr(m(k,max([(startl(k)-1),1]):stopr(k)));
	  lmt = length(indmintmp);
	  lMt = length(indmaxtmp);
	  indmin(k,1:lmt) = indmintmp + max([(startl(k)-1),1])-1;
	  indmax(k,1:lMt) = indmaxtmp + max([(startl(k)-1),1])-1;
	  if lmt < size(indmin,2)
	    indmin(k,length(indmintmp)+1:end) = 0;
	  end
	  if lMt < size(indmax,2)
	    indmax(k,length(indmaxtmp)+1:end) = 0;
	  end
	  needextr(k) = 0;
	  if k==2
	  end
	end
	
	if stopr(k) >= LX
	  lafin(k) = 1;
	  stop(k) = LX;
	  needextr(k) = 0;
	end
      end
      
      if ~susp(k)
	curindmin = indmin(k,find(indmin(k,:) >= startl(k)));
	curindmax = indmax(k,find(indmax(k,:) >= startl(k)));
	nem = length(curindmin) + length(curindmax);
	
	lml(k) = sum(curindmin < start(k));
	lMl(k) = sum(curindmax < start(k));
      end

      % loop of local on-line sifting
      while (~needextr(k) | lafin(k)) & ~stoptest(k) & ~waittest & ~susp(k)
		
	if nem < 3 & lafin(k) == 1
	  stoptest(k) = 1;
	  if nbit(k) > 1 
	    m(k+1,:) = x - sum(m(1:k,:));
	  end
	  break
	end
	
	dm = diff(curindmin);
	dM = diff(curindmax);
	
	
	if lml(k)+lMl(k) < 4 | (start(k) < max([curindmin(min(LARGTRANS,end)),curindmax(min(LARGTRANS,end))]) + mean([dm(1:min([end,5])),dM(1:min([end,5]))]))
	  margeml = 0;
	  margeMl = 0;
	    
	    if curindmax(1) < curindmin(1)
	      if m(k,1) > m(k,curindmin(1))
	lmax = fliplr(curindmax(2:min(end,NBSYM+1)));
	lmin = fliplr(curindmin(1:min(end,NBSYM)));
	lsym = curindmax(1);
	      else
	lmax = fliplr(curindmax(1:min(end,NBSYM)));
	lmin = [fliplr(curindmin(1:min(end,NBSYM-1))),1];
	lsym = 1;
	      end
	    else
	      if m(k,1) < m(k,curindmax(1))
	lmax = fliplr(curindmax(1:min(end,NBSYM)));
	lmin = fliplr(curindmin(2:min(end,NBSYM+1)));
	lsym = curindmin(1);
	      else
	lmax = [fliplr(curindmax(1:min(end,NBSYM-1))),1];
	lmin = fliplr(curindmin(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 == curindmax(1)
	lmax = fliplr(curindmax(1:min(end,NBSYM)));
      else
	lmin = fliplr(curindmin(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

      debfen = startl(k);
	  
	else
	  lmin = curindmin(1:LARGTRANS);
	  margeml = LARGTRANS;
	  lmax = curindmax(1:LARGTRANS);
	  margeMl = LARGTRANS;
	  debfen = max([lmin,lmax]);
	  tlmin = t(lmin);
	  tlmax = t(lmax);
	end
	
	if lafin(k) == 1
	  margemr = 0;
	  margeMr = 0;

        if curindmax(end) < curindmin(end)
	      if m(k,LX) > m(k,curindmax(end))
	rmax = [LX,fliplr(curindmax(max(end-NBSYM+2,1):end))];
	rmin = fliplr(curindmin(max(end-NBSYM+1,1):end));
	rsym = LX;
	      else
	rmax = fliplr(curindmax(max(end-NBSYM+1,1):end));
	rmin = fliplr(curindmin(max(end-NBSYM,1):end-1));
	rsym = curindmin(end);
	      end
	    else
	      if m(k,LX) < m(k,curindmin(end))
	rmax = fliplr(curindmax(max(end-NBSYM+1,1):end));
	rmin = [LX,fliplr(curindmin(max(end-NBSYM+2,1):end))];
	rsym = LX;
	      else         
	rmax = fliplr(curindmax(max(end-NBSYM,1):end-1));
	rmin = fliplr(curindmin(max(end-NBSYM+1,1):end));
	rsym = curindmax(end);
	      end
	    end
	  trmin = 2*t(rsym)-t(rmin);
	  trmax = 2*t(rsym)-t(rmax);

    if trmin(end) < t(LX) | trmax(end) < t(LX)
      if rsym == curindmax(end)
	rmax = fliplr(curindmax(max(end-NBSYM+1,1):end));
      else
	rmin = fliplr(curindmin(max(end-NBSYM+1,1):end));
      end
      if rsym == LX
 	error('bug')
      end
      rsym = LX;
      trmin = 2*t(rsym)-t(rmin);
      trmax = 2*t(rsym)-t(rmax);
    end

    
      finfen(k) = LX;
	  
	else
	  rmin = curindmin((end - LARGTRANS + 1):end);
	  rmax = curindmax((end - LARGTRANS + 1):end);
	  margemr = LARGTRANS;
	  margeMr = LARGTRANS;
	  finfen(k) = min([rmin,rmax]);
	  trmin = t(rmin);
	  trmax = t(rmax);
	end
	    
	mlmin = m(k,lmin);
	mlmax = m(k,lmax);
	mrmin = m(k,rmin);
	mrmax = m(k,rmax);
    
	lcur = finfen(k) - debfen + 1;
	
	
	envmax = interp1([tlmax t(curindmax(margeMl+1:(end-margeMr))) trmax],[mlmax m(k,curindmax(margeMl+1:(end-margeMr))) mrmax],t(debfen:finfen(k)),'spline');
	envmin = interp1([tlmin t(curindmin(margeml+1:(end-margemr))) trmin],[mlmin m(k,curindmin(margeml+1:(end-margemr))) mrmin],t(debfen:finfen(k)),'spline');
		      
	envmoy = (envmax + envmin)/2;

	
	amp = abs(envmax-envmin)/2;
	mamp = max(amp);
	sx = abs(envmoy)./amp;
	s = mean(sx);
	
	
	
    % definition of f equal to 1 where sifting is needed and decays 
    % fast to 0 in the neighbourhood

	badsx = (sx > sd & amp > mamp/tol) | (sx > sd2 & amp > mamp/tol2);
	d = diff([0 badsx 0]);
	debs = find(d==1);
	fins = find(d==-1)-1;
	if length(debs)~=length(fins)
	  error('pb avec les composantes connexes')
	end
	indextr = sort([curindmin curindmax]);
	d = diff(indextr);
	connexe = [];
	lc = length(debs);
	f = [];
	if lc > 0
	  connexe(lc,lcur) = 0;
	  for i = 1:lc
	    connexe(i,debs(i):fins(i)) = 1;
	    indp = min([lc-1,length(find(indextr < debs(i)))]);
	    inds = max([2,length(find(indextr <= fins(i)))+1]);
	    llarg = mean(d(min([max([nem-1,1]),max([1,indp])]):max([1,min([nem-1,indp + 2])])));
	    rlarg = mean(d(min([max([nem-1,1]),max([1,inds - 1])]):max([1,min([nem-1,inds + 1])])));
	    larg = mean([rlarg,llarg]);
	    if indp == 0
	      if inds ~= (nem+1)
		larg = (indextr(inds)-indextr(inds-1));
	      else
		larg = round(lcur/2);
	      end
	    else
	      if  inds ~= (nem+1)
		larg = max([(indextr(inds)-indextr(inds-1)),(indextr(indp+1)-indextr(indp))]);
	      else
		larg =(indextr(indp+1)-indextr(indp));
	      end
	    end
	    larg = 2 * max([round((fins(i)-debs(i))/4),larg,LARGMIN]);
	    w(1:round(larg/2)) = [1:round(larg/2)]/round(larg/2);
	    w((2*larg+2-round(larg/2)):(2*larg+1)) = fliplr(w(1:round(larg/2)));
	    w(round(larg/2):(2*larg+2-round(larg/2))) = 1;
	    indd = max(1,debs(i)-larg);
	    indf = min(lcur,fins(i)+larg);
	    connexe(i,indd:debs(i)) = w((larg+1-debs(i)+indd):(larg+1));
	    connexe(i,fins(i):indf) = w((larg+1):(larg+1+indf-fins(i)));
	  end
	  
	  f = max(connexe,[],1);
	  
	else
	  f(lcur) = 0;
	end
	    
    %                                   ___
    % definition of f2, window of form /   \ used as a multiplier on f
	f2 = [];
	f2((start(k)-debfen+1):(stop(k)-debfen+1)) = 1;
	if debfen == 1
	  dl = round(max([mean([dm(1:min([end,3])),dM(1:min([end,3]))]),start(k)/2]));
	  f2(1:start(k)) = max([dl-start(k)+1:dl]/dl,0);
	else
	  f2(1:(start(k)-debfen+1)) = [0:(start(k)-debfen)]/(start(k)-debfen);
	end
	if lafin(k) == 1
	  f2(stop(k)-debfen+1:LX-debfen+1) = 1;
	else
	  f2((stop(k)-debfen+1):(finfen(k)-debfen+1)) = fliplr([0:(finfen(k)-stop(k))])/(finfen(k)-stop(k));
	end
	
	f1 = f;
	f = f.*f2;
	    
	    
    mp = m(k,:);
	m(k,debfen:finfen(k)) = m(k,debfen:finfen(k)) - f.*envmoy;
	    
	% display
	    
	if tst == 1
	  figure(figures(1))
	  subplot(4,1,1)
	  plot(t(1:fin(k)),mp(1:fin(k)));hold on;
	  axis([t(1),t(end),-max(abs(mp)),max(abs(mp))])
	  title(['IMF ',int2str(k),';   iteration ', int2str(nbit(k)),' before sifting']);
	  plot(t(debfen:finfen(k)),envmax,'--k');plot(t(debfen:finfen(k)),envmin,'--k');plot(t(debfen:finfen(k)),envmoy,'r');
	  hold  off
      set(gca,'XTick',[])
	  subplot(4,1,2)
	  plot(t(debfen:finfen(k)),sx);
      hold on
      plot(t,sdt,'--r')
      plot(t,sd2t,':k')
      title('stop parameter')
      set(gca,'XTick',[])
	  axis([t(1),t(end),0,max(sx)])
	  hold off
	  subplot(4,1,3);
	  plot(t(debfen:finfen(k)),f);
	  hold on;
	  plot(t(debfen:finfen(k)),f2,'--k');
	  hold off
	  axis([t(1),t(end),0,1.1])
	  hold off
      title('window')
      subplot(4,1,4);
	  plot(t,m(k,:));
	  axis([t(1),t(end),-max(abs(m(k,:))),max(abs(m(k,:)))])
	  title(['IMF ',int2str(k),';   iteration ', int2str(nbit(k)),' after sifting']);
      set(gca,'XTick',[])
      pause(0.01)
	  % clf
	end
	
	
	
	% if start change
	if sum(badsx(round((start(k)-debfen+1)/2):min([(start(k)+5*PAS-debfen+1),(stop(k) - debfen +1)])))==0
	  if start(k) >= LX
	    stoptest(k) = 1;
	    startl(k) = LX;
	    suspps(k+1,1) = 0;
	    disp(['mode ',int2str(k),' termin