% 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  figureend% maximum number of iterations% MAXITERATIONS=2000;%main loop : requires at least 3 extrema to proceedwhile ~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 loopif sum(r.^2) & ~any(mask)	imf(k,:) = r;endort = io(x,imf);if display_sifting  closeend%---------------------------------------------------------------------------------------------------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;endif nargin == 2	if ischar(x)		INTERP = x;		x = t;		t = 1:length(x);	endendif ~ischar(INTERP)	error('interp parameter must be ''linear'''', ''cubic'' or ''spline''')endif ~any(strcmpi(INTERP,{'linear','cubic','spline'}))	error('interp parameter must be ''linear'''', ''cubic'' or ''spline''')endif min([size(x),size(t)]) > 1	error('x and t must be vectors')ends = size(x);if s(1) > 1	x = x';ends = size(t);if s(1) > 1	t = t';endif length(t) ~= length(x)	error('x and t must have the same length')endlx = 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 interpolationenvmax = 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)