www.pudn.com > duichengyantuo.rar > duichenyantuo.m, change:2009-04-22,size:3098b

t=linspace(0,1,1024); 
x=cos(2*pi*3*t)+0.4*cos(2*pi*10*t)+0.9*sin(2*pi*25*t); 
%figure(1) 
plot(x); 
[indmin,indmax,indzer] = extr(x); 
for i=1:986 
    x1(1,i)=x(10+i-1); 
end 
figure(2) 
plot(x1) 
x2=[x1 x1 x1] 
figure(3) 
plot(x2) 
figure(4) 
t1=1:1:2958; 
plot(t1(1:987),x2(1:987),'-',t1(987:1973),x2(987:1973),'-',t1(1973:2958),x2(1973:2958),'-'); 
 
 
n=4; 
N = length(x2); 
 
%------------------------------------------------------------------------- 
% loop to decompose the input signal into n successive IMFs 
 
imf = []; % Matrix which will contain the successive IMF, and the residue 
 
for k=1:n % loop on successive IMFs 
    
   %------------------------------------------------------------------------- 
   % inner loop to find each imf 
    
   h = x; % at the beginning of the sifting process, h is the signal 
   SD = 1; % Standard deviation which will be used to stop the sifting process 
    
   while SD > 0.3 % while the standard deviation is higher than 0.3 (typical value) 
       
      % find local max/min points 
      d = diff(h); % approximate derivative 
      maxmin = []; % to store the optima (min and max without distinction so far) 
      for i=1:N-2 
         if d(i)==0                        % we are on a zero 
            if sign(d(i-1))~=sign(d(i+1))  % it is a maximum 
               maxmin = [maxmin, i]; 
            end 
         elseif sign(d(i))~=sign(d(i+1))   % we are straddling a zero so 
            maxmin = [maxmin, i+1];        % define zero as at i+1 (not i) 
         end 
      end 
       
      if size(maxmin,2) < 2 % then it is the residue 
         break 
      end 
       
      % divide maxmin into maxes and mins 
      if maxmin(1)>maxmin(2)              % first one is a max not a min 
         maxes = maxmin(1:2:length(maxmin)); 
         mins  = maxmin(2:2:length(maxmin)); 
      else                                % is the other way around 
         maxes = maxmin(2:2:length(maxmin)); 
         mins  = maxmin(1:2:length(maxmin)); 
      end 
       
      % make endpoints both maxes and mins 
      maxes = [1 maxes N]; 
      mins  = [1 mins  N]; 
       
       
      %------------------------------------------------------------------------- 
      % spline interpolate to get max and min envelopes; form imf 
      maxenv = spline(maxes,h(maxes),1:N); 
      minenv = spline(mins, h(mins),1:N); 
       
      m = (maxenv + minenv)/2; % mean of max and min enveloppes 
      prevh = h; % copy of the previous value of h before modifying it 
      h = h - m; % substract mean to h 
       
      % calculate standard deviation 
      eps = 0.0000001; % to avoid zero values 
      SD = sum ( ((prevh - h).^2) ./ (prevh.^2 + eps) ); 
       
   end 
    
   imf = [imf; h]; % store the extracted IMF in the matrix imf 
   % if size(maxmin,2)<2, then h is the residue 
        plot(imf); 
   % stop criterion of the algo. if we reach the end before n 
   if size(maxmin,2) < 2 
      break 
   end 
    
   x= x - h; % substract the extracted IMF from the signal 
 
end 
c=imf'; 
figure(1) 
nt= plot4c(t,c,1);