function [imf,ort,nbits] = emd_local(x,t,stop,tst);
% EMD_LOCAL ("Local" Empirical Mode Decomposition) computes a
% local version of EMD
%
% 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
%
%
% [imf,ort,nbits] = EMD_LOCAL(x,t,stop,tst)
% inputs:
% - x : analyzed signal (line vector)
% - t (optional) : sampling times (line vector) (default : 1:length(x))
% - stop (optional) : threshold, threshold2 and tolerance (optional)
% for sifting stopping criterion
% default : [0.05,0.5,0.05]
% - tst (optional) : if equals to 1 shows sifting steps with pause
% if equals to 2 no pause
%
% 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,0.05];
if(nargin==1)
t = 1:length(x);
stop = defstop;
tst = 0;
end
if(nargin==2)
stop = defstop;
tst = 0;
end
if (nargin==3)
tst=0;
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) > 3) | (S(2) > 3) | (length(S) > 2)
error('stop must have only one row or one column of max three elements')
end
if S(1) > 1
stop = stop';
S = size(stop);
end
if S(2) <2
stop(2)=defstop(2);
end
if S(2) < 3
stop(3)=defstop(3);
end
sd = stop(1);
sd2 = stop(2);
tol = stop(3);
if tst
figure
end
% maximum number of iterations
MAXITERATIONS=2000;
% minimum length of the constant part of curve appended to zones badsx==1
LARGMIN = 5;
% maximum number of points symmmetrized for interpolations
NBSYM = 5;
lx = length(x);
sdt(lx) = 0;
sdt = sdt+sd;
sd2t(lx) = 0;
sd2t = sd2t+sd2;
% number of extrema and zero-crossings in residual
ner = lx;
nzr = lx;
r = x;
imf = [];
k = 1;
% iterations counter for the extraction of 1 mode
nbit=0;
% total iterations counter
NbIt=0;
% loop start : modes extracted while at least 3 extrema
while ner > 2
% current mode
m = r;
% mode at previous iteration
mp = m;
sx = sd + 1;
badsx = 1;
% tests if enough extrema to proceed
test = 0;
[indmin,indmax,indzer] = extr(m);
lm=length(indmin);
lM=length(indmax);
nem=lm + lM;
nzm=length(indzer);
j=1;
% sifting loop
while (mean(badsx) > tol | (abs(nzm-nem)>1))&(test == 0)&nbit<MAXITERATIONS
if(nbit>MAXITERATIONS/5 & mod(nbit,floor(MAXITERATIONS/10))==0)
disp(['mode ',int2str(k),' nombre d iterations : ',int2str(nbit)])
disp(['stop parameter mean value : ',num2str(s)])
end
% boundary conditions for interpolations :
if indmax(1) < indmin(1)
if m(1) > m(indmin(1))
lmax = fliplr(indmax(2:min(end,NBSYM+1)));
lmin = fliplr(indmin(1:min(end,NBSYM)));
lsym = indmax(1);
else
lmax = fliplr(indmax(1:min(end,NBSYM)));
lmin = [fliplr(indmin(1:min(end,NBSYM-1))),1];
lsym = 1;
end
else
if m(1) < m(indmax(1))
lmax = fliplr(indmax(1:min(end,NBSYM)));
lmin = fliplr(indmin(2:min(end,NBSYM+1)));
lsym = indmin(1);
else
lmax = [fliplr(indmax(1:min(end,NBSYM-1))),1];
lmin = fliplr(indmin(1:min(end,NBSYM)));
lsym = 1;
end
end
if indmax(end) < indmin(end)
if m(end) < m(indmax(end))
rmax = fliplr(indmax(max(end-NBSYM+1,1):end));
rmin = fliplr(indmin(max(end-NBSYM,1):end-1));
rsym = indmin(end);
else
rmax = [lx,fliplr(indmax(max(end-NBSYM+2,1):end))];
rmin = fliplr(indmin(max(end-NBSYM+1,1):end));
rsym = lx;
end
else
if m(end) > m(indmin(end))
rmax = fliplr(indmax(max(end-NBSYM,1):end-1));
rmin = fliplr(indmin(max(end-NBSYM+1,1):end));
rsym = indmax(end);
else
rmax = fliplr(indmax(max(end-NBSYM+1,1):end));
rmin = [lx,fliplr(indmin(max(end-NBSYM+2,1):end))];
rsym = lx;
end
end
tlmin = 2*t(lsym)-t(lmin);
tlmax = 2*t(lsym)-t(lmax);
trmin = 2*t(rsym)-t(rmin);
trmax = 2*t(rsym)-t(rmax);
% in case symmetrized parts do not extend enough
if tlmin(1) > t(1) | tlmax(1) > t(1)
if lsym == indmax(1)
lmax = fliplr(indmax(1:min(end,NBSYM)));
else
lmin = fliplr(indmin(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
if trmin(end) < t(lx) | trmax(end) < t(lx)
if rsym == indmax(end)
rmax = fliplr(indmax(max(end-NBSYM+1,1):end));
else
rmin = fliplr(indmin(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
mlmax =m(lmax);
mlmin =m(lmin);
mrmax =m(rmax);
mrmin =m(rmin);
% definition of enveloppes from interpolation
envmax = interp1([tlmax t(indmax) trmax],[mlmax m(indmax) mrmax],t,'spline');
envmin = interp1([tlmin t(indmin) trmin],[mlmin m(indmin) mrmin],t,'spline');
envmoy = (envmax + envmin)/2;
% estimation of mode amplitude
amp = abs(envmax-envmin)/2;
if any(amp == 0)
amp = amp+1e-300; % circumvents divide by zero
end
mamp = max(amp);
sx=(abs(envmoy))./amp;
s = mean(sx);
% definition of f equal to 1 where sifting is needed and
% dcays fast to 0 in the neighbourhood
badsx = (sx > sd & amp > mamp/20) | (sx > sd2 & amp > mamp/100);
d = diff([0 badsx 0]);
debs = find(d==1);
fins = find(d==-1)-1;
if length(debs)~=length(fins)
error('pb avec le regroupement en composantes connexes')
end
indextr = sort([indmin indmax]);
d = diff(indextr);
f = [];
connexe = [];
lc = length(debs);
if lc > 0
connexe(lc,lx) = 0;
for i = 1:lc
connexe(i,debs(i):fins(i)) = 1;
% indices of previous and next extrema
indp = min([nem-1,length(find(indextr < debs(i)))]);
inds = max([2,length(find(indextr <= fins(i)))+1]);
% evaluation of extrema densities left and right
llarg = mean(d(min([nem-1,max([1,indp])]):max([1,min([nem-1,indp + 2])])));
rlarg = mean(d(min([nem-1,max([1,inds - 1])]):max([1,min([nem-1,inds + 1])])));
larg = mean([rlarg,llarg]);
% special cases...
if indp == 0
if inds ~= (nem+1)
larg = (indextr(inds)-indextr(inds-1));
else
larg = round(lx/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(lx,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
