% 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�alv�s
% "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)
