emd2.m

经验模态分解的扩展运用

% EMD:  Emprical mode decomposition(结合EMD分解和硬阈值法的去噪程序,返回值为去噪后的波形数据，并同时绘出了该去噪后的结果。)
%
% imf = emd(x)
%
% x   - input signal (must be a column vector or a row 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]>
% 基于硬阈值去噪的emd分解。

function imf = emd2(x);
tic;
num=size(x);
if num(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%注：默认步长是1
         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) );%此处好像是用方差之类的东西来得到判断条件。不是！
      
   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
d=[];
for t=1:4
    mad=median(abs(imf(t,:)-median(imf(t,:))));
    p=mad/0.6745;
    t1=p*sqrt(2*log(size(imf,2)));
    t1=t1.*ones(1,size(imf,2));
    c1=sign(imf(t,:)).*abs(imf(t,:)-t1).*(abs(imf(t,:))>=t1);
    d=[d;c1];
end
imf1=add(d,1,4)+add(imf,5,size(imf,1));
imf=imf1;
figure;
plot(imf1);
toc;
return