振动信号EMD分析matlab算法

%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> sifting % 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�alves % "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�alves 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.fr function [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 proceed while ~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-16)*(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;