www.gusucode.com > 局部均值分解源代码 难得的matlab程序代码源码 > lmd/envelope.m

    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;
end

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

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

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

if min([size(x),size(t)]) > 1
	error('x and t must be vectors')
end
s = size(x);
if s(1) > 1
	x = x';
end
s = size(t);
if s(1) > 1
	t = t';
end
if length(t) ~= length(x)
	error('x and t must have the same length')
end

lx = 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 interpolation

envmax = 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)

	
	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];
end
%---------------------------------------------------------------------------------------------------
% 极值点和过零点位置提取
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
% 当连续多个采样值相同时，把最中间的一个值作为极值点，处理方式与连0类似
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  
end

end