pauseendend%---------------------------------------------------------------------------------------------------% 演示分解過程（FIX和FIX_H停止準則）function display_emd_fixe(t,m,mp,r,envmin,envmax,envmoy,nbit,k,display_sifting)subplot(3,1,1)plot(t,mp);hold on;plot(t,envmax,'--k');plot(t,envmin,'--k');plot(t,envmoy,'r');title(['IMF ',int2str(k),';   iteration ',int2str(nbit),' before sifting']);set(gca,'XTick',[])hold  offsubplot(3,1,2)plot(t,m)title(['IMF ',int2str(k),';   iteration ',int2str(nbit),' after sifting']);set(gca,'XTick',[])subplot(3,1,3);plot(t,r-m)title('residue');if display_sifting == 2  pause(0.01)else  pauseendend%---------------------------------------------------------------------------------------% 處理邊界條件（鏡像法）function [tmin,tmax,zmin,zmax] = boundary_conditions(indmin,indmax,t,x,z,nbsym)% 實數情況下，x = zlx = 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;  endelse  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;  endend% 序列末尾情況與序列開頭類似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;  endelse  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;  endend    % 將序列根據對稱中心，鏡像到兩邊tlmin = 2*t(lsym)-t(lmin);tlmax = 2*t(lsym)-t(lmax);trmin = 2*t(rsym)-t(rmin);trmax = 2*t(rsym)-t(rmax);    % 如果對稱的部分沒有足夠的極值點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 % 延拓點上的取值       zlmax = z(lmax); zlmin = z(lmin);zrmax = z(rmax); zrmin = z(rmin);     % 完成延拓tmin = [tlmin t(indmin) trmin];tmax = [tlmax t(indmax) trmax];zmin = [zlmin z(indmin) zrmin];zmax = [zlmax z(indmax) zrmax];end    %---------------------------------------------------------------------------------------------------% 極值點和過零點位置提取function [indmin, indmax, indzer] = extr(x,t)if(nargin==1)  t = 1:length(x);endm = length(x);if nargout > 2  x1 = x(1:m-1);  x2 = x(2:m);  indzer = find(x1.*x2<0);	% 尋找信號符號發生變化的位置  if any(x == 0)	% 考慮信號采樣點恰好為0的位置    iz = find( x==0 );  % 信號采樣點恰好為0的位置    indz = [];    if any(diff(iz)==1) % 出現連0的情況      zer = x == 0;	% x=0處為1，其它地方為0      dz = diff([0 zer 0]);	% 尋找0與非0的過渡點      debz = find(dz == 1);	% 0值起點      finz = find(dz == -1)-1;  % 0值終點      indz = round((debz+finz)/2);	% 選擇中間點作為過零點    else      indz = iz;	% 若沒有連0的情況，該點本身就是過零點    end    indzer = sort([indzer indz]);	% 全體過零點排序  endend% 提取極值點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;	% 最大值% 當連續多個采樣值相同時，把最中間的一個值作為極值點，處理方式與連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]);  endendend%---------------------------------------------------------------------------------------------------function ort = io(x,imf)% ort = IO(x,imf) 計算正交指數%% 輸入 : - x    : 分析信號%        - imf  : IMF信號n = size(imf,1);s = 0;% 根據公式計算for i = 1:n  for j = 1:n    if i ~= j      s = s + abs(sum(imf(i,:).*conj(imf(j,:)))/sum(x.^2));    end  endendort = 0.5*s;end%---------------------------------------------------------------------------------------------------% 函數參數解析function [x,t,sd,sd2,tol,MODE_COMPLEX,ndirs,display_sifting,sdt,sd2t,r,imf,k,nbit,NbIt,MAXITERATIONS,FIXE,FIXE_H,MAXMODES,INTERP,mask] = init(varargin)x = varargin{1};if nargin == 2  if isstruct(varargin{2})    inopts = varargin{2};  else    error('when using 2 arguments the first one is the analyzed signal X and the second one is a struct object describing the options')  endelseif nargin > 2  try    inopts = struct(varargin{2:end});  catch    error('bad argument syntax')  endend% 默認停止條件defstop = [0.05,0.5,0.05];opt_fields = {'t','stop','display','maxiterations','fix','maxmodes','interp','fix_h','mask','ndirs','complex_version'};% 時間序列，停止參數，是否演示，最大迭代次數，每一輪迭代次數，IMF個數，插值方法，每一輪迭代次數（帶條件），mask信號，方向數，是否采用復數模式defopts.stop = defstop;defopts.display = 0;defopts.t = 1:max(size(x));defopts.maxiterations = 2000;defopts.fix = 0;defopts.maxmodes = 0;defopts.interp = 'spline';defopts.fix_h = 0;defopts.mask = 0;defopts.ndirs = 4;defopts.complex_version = 2;opts = defopts;if(nargin==1)  inopts = defopts;elseif nargin == 0  error('not enough arguments')endnames = fieldnames(inopts);for nom = names'  if ~any(strcmpi(char(nom), opt_fields))    error(['bad option field name: ',char(nom)])  end  if ~isempty(eval(['inopts.',char(nom)]))     eval(['opts.',lower(char(nom)),' = inopts.',char(nom),';'])  endendt = opts.t;stop = opts.stop;display_sifting = opts.display;MAXITERATIONS = opts.maxiterations;FIXE = opts.fix;MAXMODES = opts.maxmodes;INTERP = opts.interp;FIXE_H = opts.fix_h;mask = opts.mask;ndirs = opts.ndirs;complex_version = opts.complex_version;if ~isvector(x)  error('X must have only one row or one column')endif size(x,1) > 1  x = x.';endif ~isvector(t)  error('option field T must have only one row or one column')endif ~isreal(t)  error('time instants T must be a real vector')endif size(t,1) > 1  t = t';endif (length(t)~=length(x))  error('X and option field T must have the same length')endif ~isvector(stop) || length(stop) > 3  error('option field STOP must have only one row or one column of max three elements')endif ~all(isfinite(x))  error('data elements must be finite')endif size(stop,1) > 1  stop = stop';endL = length(stop);if L < 3  stop(3) = defstop(3);endif L < 2  stop(2) = defstop(2);endif ~ischar(INTERP) || ~any(strcmpi(INTERP,{'linear','cubic','spline'}))  error('INTERP field must be ''linear'', ''cubic'', ''pchip'' or ''spline''')end% 使用mask信號時的特殊處理if any(mask)  if ~isvector(mask) || length(mask) ~= length(x)    error('masking signal must have the same dimension as the analyzed signal X')  end  if size(mask,1) > 1    mask = mask.';  end  opts.mask = 0;  imf1 = emd(x+mask, opts);  imf2 = emd(x-mask, opts);  if size(imf1,1) ~= size(imf2,1)    warning('emd:warning',['the two sets of IMFs have different sizes: ',int2str(size(imf1,1)),' and ',int2str(size(imf2,1)),' IMFs.'])  end  S1 = size(imf1,1);  S2 = size(imf2,1);  if S1 ~= S2	% 如果兩個信號分解得到的IMF個數不一致，調整順序    if S1 < S2      tmp = imf1;      imf1 = imf2;      imf2 = tmp;    end    imf2(max(S1,S2),1) = 0;	% 將短的那個補零，達到長度一致  end  imf = (imf1+imf2)/2;endsd = stop(1);sd2 = stop(2);tol = stop(3);lx = length(x);sdt = sd*ones(1,lx);sd2t = sd2*ones(1,lx);if FIXE  MAXITERATIONS = FIXE;  if FIXE_H    error('cannot use both ''FIX'' and ''FIX_H'' modes')  endendMODE_COMPLEX = ~isreal(x)*complex_version;if MODE_COMPLEX && complex_version ~= 1 && complex_version ~= 2  error('COMPLEX_VERSION parameter must equal 1 or 2')end% 極值點和過零點的個數ner = lx;nzr = lx;r = x;if ~any(mask) % 如果使用了mask信號，此時imf就已經計算得到了  imf = [];endk = 1;% 提取每個模式時迭代的次數nbit = 0;% 總體迭代次數NbIt = 0;end%---------------------------------------------------------------------------------------------------