% emprical mode decomposition （只對(duì)信號(hào)進(jìn)行EMD分解，以矩陣形式返回各個(gè)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  <ivan@ieee.org>% function imf = emd(x);tic;if size(x,1)==1    c=x;else    c = x'; % copy of the input signal (as a row vector)endN = length(x);%-------------------------------------------------------------------------% loop to decompose the input signal into successive IMFimf = []; % Matrix which will contain the successive IMF, and the residuewhile (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)等于零，則是極值點(diǎn)的位置。            maxmin = [maxmin, i];         elseif sign(d(i))~=sign(d(i+1))   % we are straddling a zero so注：sign是用來(lái)判斷一個(gè)值與零的關(guān)系，返回值有－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%若成立說(shuō)明極值點(diǎn)只剩一個(gè)或是沒(méi)有極值點(diǎn)了。         break      end            % divide maxmin into maxes and mins注意此處區(qū)別極大值與極小值的方法。很簡(jiǎn)單。      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加上兩端點(diǎn)重新構(gòu)成新的行向量。二者起點(diǎn)和終點(diǎn)相同。                  %-------------------------------------------------------------------------      % spline interpolate to get max and min envelopes; form imf此處用三次樣條曲線擬合得到上下包絡(luò)      maxenv = spline(maxes,h(maxes),1:N);      minenv = spline(mins, h(mins),1:N);%此處能直接用spline命令來(lái)得到上下包絡(luò)嗎？能用，也是三次樣條插值。            m = (maxenv + minenv)/2; % mean of max and min enveloppes簡(jiǎn)單的求上下包絡(luò)的平均值      prevh = h; % copy of the previous value of h before modifying it后面求r(t)殘差時(shí)還用的到原始信號(hào)。      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   endtoc;return