|
我发现运行程序后得到的IMF1分量的右端点发散,怎么回事啊???
clc;
clear;
t=0:0.0004:1;
x=sin(2*pi*50*t)+0.3*sin(5.5*pi*50*t);
% x=sin(2*pi*50*t)+0.3*sin(6.5*pi*50*t);
tst=1;
stop = [0.01,0.1,0.01];
NBSYM = 2;
lx = length(x);
sd = stop(1);
sd2 = stop(2);
tol = stop(3);
MAXITERATIONS=2000;
sdt(lx) = 0;
sdt = sdt+sd;
sd2t(lx) = 0;
sd2t = sd2t+sd2;
% maximum number of extrema and zero-crossings in residual
ner = lx;
nzr = lx;
r = x;
imf = [];
k = 1;
% iterations counter for extraction of 1 mode
nbit=0;
% total iterations counter
NbIt=0;
while ner > 2
% current mode
m = r;
% mode at previous iteration
mp = m;
sx = sd+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(sx > sd) > tol | any(sx > sd2) | (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 :
lm=length(indmax);
ln=length(indmin);
if indmax(1) < indmin(1)
if m(1) > m(indmin(1))
lmax = fliplr(indmax(2:min(lm,NBSYM+1)));
lmin = fliplr(indmin(1:min(ln,NBSYM)));
lsym = indmax(1);
else
lmax = fliplr(indmax(1:min(lm,NBSYM)));
lmin = [fliplr(indmin(1:min(ln,NBSYM-1))),1];
lsym = 1;
end
else
if m(1) < m(indmax(1))
lmax = fliplr(indmax(1:min(lm,NBSYM)));
lmin = fliplr(indmin(2:min(ln,NBSYM+1)));
lsym = indmin(1);
else
lmax = [fliplr(indmax(1:min(lm,NBSYM-1))),1];
lmin = fliplr(indmin(1:min(ln,NBSYM)));
lsym = 1;
end
end
if indmax(end) < indmin(end)
if m(end) < m(indmax(end))
rmax = fliplr(indmax(max(lm-NBSYM+1,1):end));
rmin = fliplr(indmin(max(ln-NBSYM,1):end-1));
rsym = indmin(end);
else
rmax = [lx,fliplr(indmax(max(lm-NBSYM+2,1):end))];
rmin = fliplr(indmin(max(ln-NBSYM+1,1):end));
rsym = lx;
end
else
if m(end) > m(indmin(end))
rmax = fliplr(indmax(max(lm-NBSYM,1):end-1));
rmin = fliplr(indmin(max(ln-NBSYM+1,1):end));
rsym = indmax(end);
else
rmax = fliplr(indmax(max(lm-NBSYM+1,1):end));
rmin = [lx,fliplr(indmin(max(ln-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(lm,NBSYM)));
else
lmin = fliplr(indmin(1:min(ln,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(lm-NBSYM+1,1):end));
else
rmin = fliplr(indmin(max(ln-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 envelopes 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;
m = m - envmoy;
[indmin,indmax,indzer] = extr(m);
lm=length(indmin);
lM=length(indmax);
nem = lm + lM;
nzm = length(indzer);
% evaluation of mean zero
sx=2*(abs(envmoy))./(abs(envmax-envmin));
s = mean(sx);
mp = m;
nbit=nbit+1;
NbIt=NbIt+1;
if(nbit==(MAXITERATIONS-1))
warning(['forced stop of sifting : too many iterations... mode ',int2str(k),'. stop parameter mean value : ',num2str(s)])
end
end
imf(k,:) = m;
if tst
disp(['mode ',int2str(k),' enregistre'])
end
nbits(k) = nbit;
k = k+1;
r = r - m;
[indmin,indmax,indzer] = extr(r);
ner = length(indmin) + length(indmax);
nzr = length(indzer);
nbit=1;
if (max(r) - min(r)) < (1e-10)*(max(x) - min(x))
if ner > 2
warning('forced stop of EMD : too small amplitude')
else
disp('forced stop of EMD : too small amplitude')
end
break
end
end
imf(k,:) = r;
if tst
close
end |
|