emd.m

EMD的matlab例程

	lx = length(x);
	
	if (length(indmin) + length(indmax) < 3)
		error('not enough extrema')
	end

	if indmax(1) < indmin(1)
    	if x(1) > x(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 x(1) < x(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 x(end) < x(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 x(end) > x(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 

	xlmax =x(lmax); 
	xlmin =x(lmin);
	xrmax =x(rmax); 
	xrmin =x(rmin);

	tmin = [tlmin t(indmin) trmin];
	tmax = [tlmax t(indmax) trmax];
	xmin = [xlmin x(indmin) xrmin];
	xmax = [xlmax x(indmax) xrmax];
  
%---------------------------------------------------------------------------------------------------

function [indmin, indmax, indzer] = extr(x,t);
%extracts the indices corresponding to extrema

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

m = length(x);

if nargout > 2
	x1=x(1:m-1);
	x2=x(2:m);
	indzer = find(x1.*x2<0);
	
	if any(x == 0)
	  iz = find( x==0 );
	  indz = [];
	  if any(diff(iz)==1)
	    zer = x == 0;
	    dz = diff([0 zer 0]);
	    debz = find(dz == 1);
	    finz = find(dz == -1)-1;
	    indz = round((debz+finz)/2);
	  else
	    indz = iz;
	  end
	  indzer = sort([indzer indz]);
	end
end
  
d = diff(x);

n = length(d);
d1 = d(1:n-1);
d2 = d(2:n);
indmin = find(d1.*d2<0 & d1<0)+1;
indmax = find(d1.*d2<0 & d1>0)+1;


% when two or more consecutive points have the same value we consider only one extremum in the middle of the constant area

if any(d==0)
  
  imax = [];
  imin = [];
  
  bad = (d==0);
  dd = diff([0 bad 0]);
  debs = find(dd == 1);
  fins = find(dd == -1);
  if debs(1) == 1
    if length(debs) > 1
      debs = debs(2:end);
      fins = fins(2:end);
    else
      debs = [];
      fins = [];
    end
  end
  if length(debs) > 0
    if fins(end) == m
      if length(debs) > 1
        debs = debs(1:(end-1));
        fins = fins(1:(end-1));

      else
        debs = [];
        fins = [];
      end      
    end
  end
  lc = length(debs);
  if lc > 0
    for k = 1:lc
      if d(debs(k)-1) > 0
        if d(fins(k)) < 0
          imax = [imax round((fins(k)+debs(k))/2)];
        end
      else
        if d(fins(k)) > 0
          imin = [imin round((fins(k)+debs(k))/2)];
        end
      end
    end
  end
  
  if length(imax) > 0
    indmax = sort([indmax imax]);
  end

  if length(imin) > 0
    indmin = sort([indmin imin]);
  end
  
end  
  
%---------------------------------------------------------------------------------------------------

function [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)	

x = varargin{1};
if nargin == 2
	if strcmp(class(varargin{2}),'struct') 
		inopts = varargin{2};
	else
		error('when using 2 arguments the first one is the analysed signal x and the second one is a struct object describing the options')
	end
elseif nargin > 2
	try
		inopts = struct(varargin{2:end});
	catch
		error('bad argument syntax')
	end
end

% default for stopping
defstop = [0.05,0.5,0.05];

opt_fields = {'t','stop','display','maxiterations','fix','maxmodes','interp','fix_h','mask'};

defopts.stop = defstop;
defopts.display = 0;
defopts.t = 1:max(size(x));
defopts.maxiterations = 2000;
defopts.fix = 0;
defopts.maxmodes = 0;
defopts.interp = 'spline';
defopts.fix_h = 0;
defopts.mask = 0;

opts = defopts;



if(nargin==1)
	inopts = defopts;
elseif nargin == 0
	error('not enough arguments')
end


names = fieldnames(inopts);
for nom = names'
	if length(strmatch(char(nom), opt_fields)) == 0
		error(['bad option field name: ',char(nom)])
	end
	eval(['opts.',char(nom),' = inopts.',char(nom),';'])
end

t = opts.t;
stop = opts.stop;
display_sifting = opts.display;
MAXITERATIONS = opts.maxiterations;
FIXE = opts.fix;
MAXMODES = opts.maxmodes;
INTERP = opts.interp;
FIXE_H = opts.fix_h;
mask = opts.mask;

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('option field 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 option field 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('option field stop must have only one row or one column of max three elements')
end

if ~all(isfinite(x))
	error('data elements must be finite')
end

if S(1) > 1
  stop = stop';
  S = size(stop);
end

if S(2) < 3
  stop(3)=defstop(3);
end

if S(2) < 2
  stop(2)=defstop(2);
end


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

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

%special procedure when a masking signal is specified
if any(mask)
	S = size(mask);
	if min(S) > 1
		error('masking signal must have the same dimension as the analyzed signal x')
	end
	if S(1) > 1
		mask = mask';
	end
	if max(S) ~= max(size(x))
		error('masking signal must have the same dimension as the analyzed signal x')
	end
	opts.mask = 0;
	imf1 = emd(x+mask,opts);
	imf2 = emd(x-mask,opts);
	if size(imf1,1) ~= size(imf2,1)
		warning(['the two sets of IMFs have different sizes: ',int2str(size(imf1,1)),' and ',int2str(size(imf2,1)),' IMFs.'])
	end
	S1 = size(imf1,1);
	S2 = size(imf2,1);
	if S1 ~= S2
		if S1 < S2
			tmp = imf1;
			imf1 = imf2;
			imf2 = tmp;
		end
		imf2(max(S1,S2),1) = 0;
	end
	imf = (imf1+imf2)/2;

end


sd = stop(1);
sd2 = stop(2);
tol = stop(3);

lx = length(x);

sdt = sd*ones(1,lx);
sd2t = sd2*ones(1,lx);

if FIXE
	MAXITERATIONS = FIXE;
	if FIXE_H
		error('cannot use both ''fix'' and ''fix_h'' modes') 
	end
end

% number of extrema and zero-crossings in residual
ner = lx;
nzr = lx;

r = x;

if ~any(mask) % if a masking signal is specified "imf" already exists at this stage
	imf = [];
end
k = 1;

% iterations counter for extraction of 1 mode
nbit=0;

% total iterations counter
NbIt=0;

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