emd.m

EMD的matlab例程

% EMD.M
%
% G. Rilling, last update: May 2005
%
% 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
%
% default 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)
%		- opts (optional): struct object with (optional) fields:
%			- t: sampling times (line vector) (default: 1:length(x))
%			- stop: threshold, threshold2 and tolerance (optional)
%				for sifting stopping criterion 
%				default: [0.05,0.5,0.05]
%			- display: if equals to 1 shows sifting steps with pause
%				if equals to 2 shows sifting steps without pause (movie style)
%			- maxiterations: maximum number of sifting steps for the computation of each mode
%			- fix (int): disable the stopping criterion and do exactly
%				the value of the field number of sifting steps for each mode
%			- maxmodes: maximum number of imfs extracted
%			- interp: interpolation scheme: 'linear', 'cubic' or 'spline' (default)
%			- fix_h (int): do <fix_h> sifting iterations with |#zeros-#extrema|<=1 to stop 
%				according to N. E. Huang et al., "A confidence limit for the Empirical Mode 
%				Decomposition and Hilbert spectral analysis",
%				Proc. Royal Soc. London A, Vol. 459, pp. 2317-2345, 2003
%			- mask: masking signal used to improve the decomposition
%				according to R. Deering and J. F. Kaiser, "The use of a masking signal to 
%				improve empirical mode decomposition",
%				ICASSP 2005
%
% outputs: 
%		- imf: intrinsic mode functions (last line = residual)
%		- ort: index of orthogonality
%		- nbits: number of iterations for each mode
%
% calls:   
%		- io: computes the index of orthogonality
%
%examples:
%
%>>x = rand(1,512);
%
%>>imf = emd(x);
%
%>>imf = emd(x,struct('stop',[0.1,0.5,0.05],'maxiterations',100));
%Remark: the following syntax is equivalent
%>>imf = emd(x,'stop',[0.1,0.5,0.05],'maxiterations',100);
%
%>>options.dislpay = 1;
%>>options.fix = 10;
%>>options.maxmodes = 3;
%>>[imf,ort,nbits] = emd(x,options);



function [imf,ort,nbits] = emd(varargin);

[x,t,sd,sd2,tol,display_sifting,sdt,sd2t,ner,nzr,lx,r,imf,k,nbit,NbIt,MAXITERATIONS,FIXE,FIXE_H,MAXMODES,INTERP,mask] = init(varargin{:});	

if display_sifting
  figure
end

% maximum number of iterations
% MAXITERATIONS=2000;


%main loop : requires at least 3 extrema to proceed
while ~stop_EMD(r) & (k < MAXMODES+1 | MAXMODES == 0) & ~any(mask)

	% current mode
	m = r;
	

	% mode at previous iteration
	mp = m;

			
	if FIXE
		[stop_sift,moyenne] = stop_sifting_fixe(t,m,INTERP);
	elseif FIXE_H
		stop_count = 0;
		[stop_sift,moyenne,stop_count] = stop_sifting_fixe_h(t,m,INTERP,stop_count,FIXE_H);
		stop_count = 0;
	else
		[stop_sift,moyenne] = stop_sifting(m,t,sd,sd2,tol,INTERP);
	end	

  	if (max(m) - min(m)) < (1e-10)*(max(x) - min(x))
    	if ~stop_sift
		   	warning('forced stop of EMD : too small amplitude')
    	else
      		disp('forced stop of EMD : too small amplitude')
    	end
    	break
  	end

    
  	% sifting loop
  	while ~stop_sift & nbit<MAXITERATIONS
    
    	if(nbit>MAXITERATIONS/5 & mod(nbit,floor(MAXITERATIONS/10))==0 & ~FIXE & nbit > 100)
      		disp(['mode ',int2str(k),', iteration ',int2str(nbit)])
			if exist('s')
	      		disp(['stop parameter mean value : ',num2str(s)])
			end
            [im,iM] = extr(m);
            disp([int2str(sum(m(im) > 0)),' minima > 0; ',int2str(sum(m(iM) < 0)),' maxima < 0.'])
		end
		
   		%sifting
    	m = m - moyenne;
		
		%computation of mean and stopping criterion  
		if FIXE
			[stop_sift,moyenne] = stop_sifting_fixe(t,m,INTERP);
		elseif FIXE_H
			[stop_sift,moyenne,stop_count] = stop_sifting_fixe_h(t,m,INTERP,stop_count,FIXE_H);
		else
			[stop_sift,moyenne,s] = stop_sifting(m,t,sd,sd2,tol,INTERP);  
		end		
	    
    	% display
        
    	if display_sifting
			[envminp,envmaxp,envmoyp] = envelope(t,mp,INTERP);
			if FIXE |FIXE_H
				display_emd_fixe(t,m,mp,r,envminp,envmaxp,envmoyp,nbit,k,display_sifting)
			else
			   	sxp=2*(abs(envmoyp))./(abs(envmaxp-envminp));
			   	sp = mean(sxp);
				display_emd(t,m,mp,r,envminp,envmaxp,envmoyp,s,sp,sxp,sdt,sd2t,nbit,k,display_sifting,stop_sift)
			end
    	end


    	mp = m;
    	nbit=nbit+1;
    	NbIt=NbIt+1;

    	if(nbit==(MAXITERATIONS-1) & ~FIXE & nbit > 100)
            if exist('s')
                warning(['forced stop of sifting : too many iterations... mode ',int2str(k),'. stop parameter mean value : ',num2str(s)])
            else
                warning(['forced stop of sifting : too many iterations... mode ',int2str(k),'.'])
            end
    	end
  
  	end % sifting loop
  	imf(k,:) = m;
  	if display_sifting
  		disp(['mode ',int2str(k),' stored'])
  	end
  	nbits(k) = nbit;
  	k = k+1;

	
  	r = r - m;
  	nbit=0;


end %main loop

if sum(r.^2) & ~any(mask)
	imf(k,:) = r;
end

ort = io(x,imf);

if display_sifting
  close
end

%---------------------------------------------------------------------------------------------------

function stop = stop_EMD(r)
	[indmin,indmax,indzer] = extr(r);
	ner = length(indmin) + length(indmax);
	stop = ner <3;	


%---------------------------------------------------------------------------------------------------


function [stop,envmoy,s]= stop_sifting(m,t,sd,sd2,tol,INTERP)
	try
	   	[envmin,envmax,envmoy,indmin,indmax,indzer] = envelope(t,m,INTERP);
	   	nem = length(indmin) + length(indmax);
	   	nzm = length(indzer);

	   	% evaluation of mean zero
	   	sx=2*(abs(envmoy))./(abs(envmax-envmin));
	   	s = mean(sx);

		stop = ~((mean(sx > sd) > tol | any(sx > sd2) | (abs(nzm-nem)>1)) & (nem > 2));
	catch
		stop = 1;
		envmoy = zeros(1,length(m));
		s = NaN;
	end


%---------------------------------------------------------------------------------------------------
function [stop,moyenne]= stop_sifting_fixe(t,m,INTERP)
	try
		[envmin,envmax,moyenne] = envelope(t,m,INTERP);
		stop = 0;
	catch
		moyenne = zeros(1,length(m));
		stop = 1;
	end


%---------------------------------------------------------------------------------------------------
function [stop,moyenne,stop_count]= stop_sifting_fixe_h(t,m,INTERP,stop_count,FIXE_H)
	try
		[envmin,envmax,moyenne,indmin,indmax,indzer] = envelope(t,m,INTERP);
	   	nem = length(indmin) + length(indmax);
	   	nzm = length(indzer);

		if (abs(nzm-nem)>1)
			stop = 0;
			stop_count = 0;
		else
			stop_count = stop_count+1;
			stop = (stop_count == FIXE_H);
		end
	catch
		moyenne = zeros(1,length(m));
		stop = 1;
	end


%---------------------------------------------------------------------------------------------------

function display_emd(t,m,mp,r,envmin,envmax,envmoy,s,sb,sx,sdt,sd2t,nbit,k,display_sifting,stop_sift)
	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(sb),' before sifting and ',num2str(s),' after'])
	if stop_sift
		disp('last iteration for this mode')
	end
    if display_sifting == 2
      pause(0.01)
    else
      pause
    end

		       
%---------------------------------------------------------------------------------------------------

function display_emd_fixe(t,m,mp,r,envmin,envmax,envmoy,nbit,k,display_sifting)
	subplot(3,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(3,1,2)
    plot(t,m)
    title(['IMF ',int2str(k),';   iteration ',int2str(nbit),' after sifting']);
    set(gca,'XTick',[])
    subplot(3,1,3);
    plot(t,r-m)
    title('residue');
    if display_sifting == 2
      pause(0.01)
    else
      pause
  end
		       
%---------------------------------------------------------------------------------------------------

function [envmin, envmax,envmoy,indmin,indmax,indzer] = envelope(t,x,INTERP)
%computes envelopes and mean with various interpolations
	
NBSYM = 2;		
DEF_INTERP = 'spline';
	
	
if nargin < 2
	x = t;
	t = 1:length(x);
	INTERP = DEF_INTERP;
end

if nargin == 2
	if ischar(x)
		INTERP = x;
		x = t;
		t = 1:length(x);
	end
end

if ~ischar(INTERP)
	error('interp parameter must be ''linear'''', ''cubic'' or ''spline''')
end

if ~any(strcmpi(INTERP,{'linear','cubic','spline'}))
	error('interp parameter must be ''linear'''', ''cubic'' or ''spline''')
end

if min([size(x),size(t)]) > 1
	error('x and t must be vectors')
end
s = size(x);
if s(1) > 1
	x = x';
end
s = size(t);
if s(1) > 1
	t = t';
end
if length(t) ~= length(x)
	error('x and t must have the same length')
end

lx = length(x);
[indmin,indmax,indzer] = extr(x,t);


%boundary conditions for interpolation
		
[tmin,tmax,xmin,xmax] = boundary_conditions(indmin,indmax,t,x,NBSYM);

% definition of envelopes from interpolation

envmax = interp1(tmax,xmax,t,INTERP);	
envmin = interp1(tmin,xmin,t,INTERP);

if nargout > 2
    envmoy = (envmax + envmin)/2;
end


%---------------------------------------------------------------------------------------

function [tmin,tmax,xmin,xmax] = boundary_conditions(indmin,indmax,t,x,nbsym)
% computes the boundary conditions for interpolation (mainly mirror symmetry)