%我現(xiàn)在百思不得其解的還是篩分過(guò)程中幅值改為amp=2./mean(abs(envmax-envmin),1);的原因%另外，其實(shí)我也懷疑是否寫錯(cuò)了，似乎應(yīng)該是：amp=mean(abs(envmax-envmin),1)/2;%又有新的修改了嗎，sxp=2*(abs(envmoyp))./(abs(envmaxp-envminp))這邊改回來(lái)%EMD  computes Empirical Mode Decomposition%%%   Syntax%%% IMF = EMD(X)% IMF = EMD(X,...,'Option_name',Option_value,...)% IMF = EMD(X,OPTS)% [IMF,ORT,NB_ITERATIONS] = EMD(...)%%%   Description%%% IMF = EMD(X) where X is a real vector computes the Empirical Mode% Decomposition [1] of X, resulting in a matrix IMF containing 1 IMF per row, the% last one being the residue. The default stopping criterion is the one proposed% in [2]:%%   at each point, mean_amplitude < THRESHOLD2*envelope_amplitude%   &%   mean of boolean array {(mean_amplitude)/(envelope_amplitude) > THRESHOLD} < TOLERANCE%   &%   |#zeros-#extrema|<=1%% where mean_amplitude = abs(envelope_max+envelope_min)/2% and envelope_amplitude = abs(envelope_max-envelope_min)/2% % IMF = EMD(X) where X is a complex vector computes Bivariate Empirical Mode% Decomposition [3] of X, resulting in a matrix IMF containing 1 IMF per row, the% last one being the residue. The default stopping criterion is similar to the% one proposed in [2]:%%   at each point, mean_amplitude < THRESHOLD2*envelope_amplitude%   &%   mean of boolean array {(mean_amplitude)/(envelope_amplitude) > THRESHOLD} < TOLERANCE%% where mean_amplitude and envelope_amplitude have definitions similar to the% real case%% IMF = EMD(X,...,'Option_name',Option_value,...) sets options Option_name to% the specified Option_value (see Options)%% IMF = EMD(X,OPTS) is equivalent to the above syntax provided OPTS is a struct % object with field names corresponding to option names and field values being the % associated values %% [IMF,ORT,NB_ITERATIONS] = EMD(...) returns an index of orthogonality%                       ________%         _  |IMF(i,:).*IMF(j,:)|%   ORT = \ _____________________%         /%         ?       || X ||?%        i~=j%% and the number of iterations to extract each mode in NB_ITERATIONS%%%   Options%%%  stopping criterion options:%% STOP: vector of stopping parameters [THRESHOLD,THRESHOLD2,TOLERANCE]% if the input vector's length is less than 3, only the first parameters are% set, the remaining ones taking default values.% default: [0.05,0.5,0.05]%% FIX (int): disable the default stopping criterion and do exactly <FIX> % number of sifting iterations for each mode%% FIX_H (int): disable the default stopping criterion and do <FIX_H> siftiamp=2./mean(abs(envmax-envmin),1);ng % iterations with |#zeros-#extrema|<=1 to stop [4]%%  bivariate/complex EMD options:%% COMPLEX_VERSION: selects the algorithm used for complex EMD ([3])% COMPLEX_VERSION = 1: "algorithm 1"% COMPLEX_VERSION = 2: "algorithm 2" (default)% % NDIRS: number of directions in which envelopes are computed (default 4)% rem: the actual number of directions (according to [3]) is 2*NDIRS% %  other options:%% T: sampling times (line vector) (default: 1:length(x))%% MAXITERATIONS: maximum number of sifting iterations for the computation of each% mode (default: 2000)%% MAXMODES: maximum number of imfs extracted (default: Inf)%% DISPLAY: if equals to 1 shows sifting steps with pause% if equals to 2 shows sifting steps without pause (movie style)% rem: display is disabled when the input is complex%% INTERP: interpolation scheme: 'linear', 'cubic', 'pchip' or 'spline' (default)% see interp1 documentation for details%% MASK: masking signal used to improve the decomposition according to [5]%%%   Examples%%%X = rand(1,512);%%IMF = emd(X);%%IMF = emd(X,'STOP',[0.1,0.5,0.05],'MAXITERATIONS',100);%%T=linspace(0,20,1e3);%X = 2*exp(i*T)+exp(3*i*T)+.5*T;%IMF = emd(X,'T',T);%%OPTIONS.DISLPAY = 1;%OPTIONS.FIX = 10;%OPTIONS.MAXMODES = 3;%[IMF,ORT,NBITS] = emd(X,OPTIONS);%%%   References%%% [1] N. E. Huang et al., "The empirical mode decomposition and the% Hilbert spectrum for non-linear and non stationary time series analysis",% Proc. Royal Soc. London A, Vol. 454, pp. 903-995, 1998%% [2] G. Rilling, P. Flandrin and P. Gon鏰lves% "On Empirical Mode Decomposition and its algorithms",% IEEE-EURASIP Workshop on Nonlinear Signal and Image Processing% NSIP-03, Grado (I), June 2003%% [3] G. Rilling, P. Flandrin, P. Gon鏰lves and J. M. Lilly.,% "Bivariate Empirical Mode Decomposition",% Signal Processing Letters (submitted)%% [4] N. E. Huang et al., "A confidence limit for the Empirical Mode% Decomposition and Hilbert spectral analysis",% Proc. Royal Soc. London A, Vol. 459, pp. 2317-2345, 2003%% [5] R. Deering and J. F. Kaiser, "The use of a masking signal to improve % empirical mode decomposition", ICASSP 2005%%% See also%  emd_visu (visualization),%  emdc, emdc_fix (fast implementations of EMD),%  cemdc, cemdc_fix, cemdc2, cemdc2_fix (fast implementations of bivariate EMD),%  hhspectrum (Hilbert-Huang spectrum)%%% G. Rilling, last modification: 3.2007% gabriel.rilling@ens-lyon.frfunction [imf,ort,nbits] = emd(varargin)[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{:});if display_sifting  fig_h = figure;end%main loop : requires at least 3 extrema to proceedwhile ~stop_EMD(r,MODE_COMPLEX,ndirs) && (k < MAXMODES+1 || MAXMODES == 0) && ~any(mask)  % current mode  m = r;  % mode at previous iteration  mp = m;  %computation of mean and stopping criterion  if FIXE    [stop_sift,moyenne] = stop_sifting_fixe(t,m,INTERP,MODE_COMPLEX,ndirs);  elseif FIXE_H    stop_count = 0;    [stop_sift,moyenne] = stop_sifting_fixe_h(t,m,INTERP,stop_count,FIXE_H,MODE_COMPLEX,ndirs);  else    [stop_sift,moyenne] = stop_sifting(m,t,sd,sd2,tol,INTERP,MODE_COMPLEX,ndirs);  end  % in case the current mode is so small that machine precision can cause  % spurious extrema to appear  if (max(abs(m))) < (1e-10)*(max(abs(x)))    if ~stop_sift      warning('emd:warning','forced stop of EMD : too small amplitude')    else      disp('forced stop of EMD : too small amplitude')    end    break  end  % sifting loop  while ~stop_sift && nbit<MAXITERATIONS    if(~MODE_COMPLEX && nbit>MAXITERATIONS/5 && mod(nbit,floor(MAXITERATIONS/10))==0 && ~FIXE && nbit > 100)      disp(['mode ',int2str(k),', iteration ',int2str(nbit)])      if exist('s','var')        disp(['stop parameter mean value : ',num2str(s)])      end      [im,iM] = extr(m);      disp([int2str(sum(m(im) > 0)),' minima > 0; ',int2str(sum(m(iM) < 0)),' maxima < 0.'])    end    %sifting    m = m - moyenne;    %computation of mean and stopping criterion    if FIXE      [stop_sift,moyenne] = stop_sifting_fixe(t,m,INTERP,MODE_COMPLEX,ndirs);    elseif FIXE_H      [stop_sift,moyenne,stop_count] = stop_sifting_fixe_h(t,m,INTERP,stop_count,FIXE_H,MODE_COMPLEX,ndirs);    else      [stop_sift,moyenne,s] = stop_sifting(m,t,sd,sd2,tol,INTERP,MODE_COMPLEX,ndirs);    end    % display    if display_sifting && ~MODE_COMPLEX      NBSYM = 2;      [indmin,indmax] = extr(mp);      [tmin,tmax,mmin,mmax] = boundary_conditions(indmin,indmax,t,mp,mp,NBSYM);      envminp = interp1(tmin,mmin,t,INTERP);      envmaxp = interp1(tmax,mmax,t,INTERP);      envmoyp = (envminp+envmaxp)/2;      if FIXE || FIXE_H        display_emd_fixe(t,m,mp,r,envminp,envmaxp,envmoyp,nbit,k,display_sifting)      else        sxp=2*(abs(envmoyp))./(abs(envmaxp-envminp));        sp = mean(sxp);        display_emd(t,m,mp,r,envminp,envmaxp,envmoyp,s,sp,sxp,sdt,sd2t,nbit,k,display_sifting,stop_sift)      end    end    mp = m;    nbit=nbit+1;    NbIt=NbIt+1;    if(nbit==(MAXITERATIONS-1) && ~FIXE && nbit > 100)      if exist('s','var')        warning('emd:warning',['forced stop of sifting : too many iterations... mode ',int2str(k),'. stop parameter mean value : ',num2str(s)])      else        warning('emd:warning',['forced stop of sifting : too many iterations... mode ',int2str(k),'.'])      end    end  end % sifting loop  imf(k,:) = m;  if display_sifting    disp(['mode ',int2str(k),' stored'])  end  nbits(k) = nbit;  k = k+1;  r = r - m;  nbit=0;end %main loopif any(r) && ~any(mask)  imf(k,:) = r;endort = io(x,imf);if display_sifting  closeendend%---------------------------------------------------------------------------------------------------% tests if there are enough (3) extrema to continue the decompositionfunction stop = stop_EMD(r,MODE_COMPLEX,ndirs)if MODE_COMPLEX  for k = 1:ndirs    phi = (k-1)*pi/ndirs;    [indmin,indmax] = extr(real(exp(i*phi)*r));    ner(k) = length(indmin) + length(indmax);  end  stop = any(ner < 3);else  [indmin,indmax] = extr(r);  ner = length(indmin) + length(indmax);  stop = ner < 3;endend%---------------------------------------------------------------------------------------------------% computes the mean of the envelopes and the mode amplitude estimatefunction [envmoy,nem,nzm,amp] = mean_and_amplitude(m,t,INTERP,MODE_COMPLEX,ndirs)NBSYM = 2;if MODE_COMPLEX  switch MODE_COMPLEX    case 1      for k = 1:ndirs        phi = (k-1)*pi/ndirs;        y = real(exp(-i*phi)*m);        [indmin,indmax,indzer] = extr(y);        nem(k) = length(indmin)+length(indmax);        nzm(k) = length(indzer);        [tmin,tmax,zmin,zmax] = boundary_conditions(indmin,indmax,t,y,m,NBSYM);        envmin(k,:) = interp1(tmin,zmin,t,INTERP);        envmax(k,:) = interp1(tmax,zmax,t,INTERP);      end      envmoy = mean((envmin+envmax)/2,1);      if nargout > 3        amp = mean(abs(envmax-envmin),1)/2;      end    case 2      for k = 1:ndirs        phi = (k-1)*pi/ndirs;        y = real(exp(-i*phi)*m);        [indmin,indmax,indzer] = extr(y);        nem(k) = length(indmin)+length(indmax);        nzm(k) = length(indzer);        [tmin,tmax,zmin,zmax] = boundary_conditions(indmin,indmax,t,y,y,NBSYM);        envmin(k,:) = exp(i*phi)*interp1(tmin,zmin,t,INTERP);        envmax(k,:) = exp(i*phi)*interp1(tmax,zmax,t,INTERP);      end      envmoy = mean((envmin+envmax),1);      if nargout > 3        amp = mean(abs(envmax-envmin),1)/2;      end  endelse  [indmin,indmax,indzer] = extr(m);  nem = length(indmin)+length(indmax);  nzm = length(indzer);  [tmin,tmax,mmin,mmax] = boundary_conditions(indmin,indmax,t,m,m,NBSYM);  envmin = interp1(tmin,mmin,t,INTERP);  envmax = interp1(tmax,mmax,t,INTERP);  envmoy = (envmin+envmax)/2;  if nargout > 3    amp = mean(abs(envmax-envmin),1)/2;  endendend%-------------------------------------------------------------------------------% default stopping criterionfunction [stop,envmoy,s] = stop_sifting(m,t,sd,sd2,tol,INTERP,MODE_COMPLEX,ndirs)try  [envmoy,nem,nzm,amp] = mean_and_amplitude(m,t,INTERP,MODE_COMPLEX,ndirs);  sx = abs(envmoy)./amp;  s = mean(sx);  stop = ~((mean(sx > sd) > tol | any(sx > sd2)) & (all(nem > 2)));  if ~MODE_COMPLEX    stop = stop && ~(abs(nzm-nem)>1);  endcatch  stop = 1;  envmoy = zeros(1,length(m));  s = NaN;endend%-------------------------------------------------------------------------------% stopping criterion corresponding to option FIXfunction [stop,moyenne]= stop_sifting_fixe(t,m,INTERP,MODE_COMPLEX,ndirs)try  moyenne = mean_and_amplitude(m,t,INTERP,MODE_COMPLEX,ndirs);  stop = 0;catch  moyenne = zeros(1,length(m));  stop = 1;endend%-------------------------------------------------------------------------------% stopping criterion corresponding to option FIX_Hfunction [stop,moyenne,stop_count]= stop_sifting_fixe_h(t,m,INTERP,stop_count,FIXE_H,MODE_COMPLEX,ndirs)try  [moyenne,nem,nzm] = mean_and_amplitude(m,t,INTERP,MODE_COMPLEX,ndirs);  if (all(abs(nzm-nem)>1))    stop = 0;    stop_count = 0;  else    stop_count = stop_count+1;    stop = (stop_count == FIXE_H);  endcatch  moyenne = zeros(1,length(m));  stop = 1;endend%-------------------------------------------------------------------------------% displays the progression of the decomposition with the default stopping criterionfunction display_emd(t,m,mp,r,envmin,envmax,envmoy,s,sb,sx,sdt,sd2t,nbit,k,display_sifting,stop_sift)subplot(4,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(4,1,2)plot(t,sx)hold onplot(t,sdt,'--r')plot(t,sd2t,':k')title('stop parameter')set(gca,'XTick',[])hold offsubplot(4,1,3)plot(t,m)title(['IMF ',int2str(k),';   iteration ',int2str(nbit),' after sifting']);set(gca,'XTick',[])subplot(4,1,4);plot(t,r-m)title('residue');disp(['stop parameter mean value : ',num2str(sb),' before sifting and ',num2str(s),' after'])if stop_sift  disp('last iteration for this mode')endif display_sifting == 2  pause(0.01)else  pauseendend%---------------------------------------------------------------------------------------------------% displays the progression of the decomposition with the FIX and FIX_H stopping criteriafunction display_emd_fixe(t,m,mp,r,envmin,envmax,envmoy,nbit,k,display_sifting)