emd.m

EMD经验模态分解,用于信号处理

% emprical mode decomposition （只对信号进行EMD分解，以矩阵形式返回各个imf值。）
% EMD:  Emprical mode decomposition
%
% imf = emd(x)
%
% x   - input signal (must be a column vector)
%
% This version will calculate all the imf's (longer)
%
% imf - Matrix of intrinsic mode functions (each as a row)
%       with residual in last row.
%
% See:  Huang et al, Royal Society Proceedings on Math, Physical, 
%       and Engineering Sciences, vol. 454, no. 1971, pp. 903-995, 
%       8 March 1998
%
% Author: Ivan Magrin-Chagnolleau  <[e-mail]>
% 

function imf = emd(x);
tic;
if size(x,1)==1
    c=x;
else
    c = x'; % copy of the input signal (as a row vector)
end
N = length(x);

%-------------------------------------------------------------------------
% loop to decompose the input signal into successive IMF

imf = []; % Matrix which will contain the successive IMF, and the residue

while (1) % the stop criterion is tested at the end of the loop
   
   %-------------------------------------------------------------------------
   % inner loop to find each imf
   
   h = c; % 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如果d(i)等于零，则是极值点的位置。
            maxmin = [maxmin, i];
         elseif sign(d(i))~=sign(d(i+1))   % we are straddling a zero so注：sign是用来判断一个值与零的关系，返回值有－1，0，1三种情况。
            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];%此处是用maxes/mins加上两端点重新构成新的行向量。二者起点和终点相同。
      
      
      %-------------------------------------------------------------------------
      % spline interpolate to get max and min envelopes; form imf此处用三次样条曲线拟合得到上下包络
      maxenv = spline(maxes,h(maxes),1:N);
      minenv = spline(mins, h(mins),1:N);%此处能直接用spline命令来得到上下包络吗？能用，也是三次样条插值。
      
      m = (maxenv + minenv)/2; % mean of max and min enveloppes简单的求上下包络的平均值
      prevh = h; % copy of the previous value of h before modifying it后面求r(t)残差时还用的到原始信号。
      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) );%判断imf的均值是否为零。粗略的判断。
      
   end
   
   imf = [imf; h]; % store the extracted IMF in the matrix imf
   % if size(maxmin,2)<2, then h is the residue
   
   % stop criterion of the algo.
   if size(maxmin,2) < 2
      break
   end
   
   c = c - h; % substract the extracted IMF from the signal
   
end
toc;
return