EMD Algorithm

  • Published on
    27-Nov-2014

  • View
    76

  • Download
    6

Embed Size (px)

Transcript

The analysis of the EmpiricalMode DecompositionMethodVesselin Vatchev, USCNovember 20, 2002The Empirical Mode Decompo-sition MethodFor a given real signal s(t) we look for a de-composition into simpler signals(modes)s(t) =Mj=1aj(t) cos j(t),where aj(t) is the amplitude and j(t) is thephase of the j-th component. Each of the com-ponents has to have physical and mathemat-ical meaning. Let s(t) be mono componentsignal, i.e. we can nd representation of theforms(t) = a(t) cos (t)that is both physically(

(t) 0) and mathe-matically meaningful. There are innitely manyways to construct such representations but itis often advantageous to write the signal incomplex formS(t) = s(t) +isi(t) = A(t) expi(t)and to take the actual signal to be the realpart of the complex signal. The imaginary partsi(t) of S(t) has to be chosen to achieve asensible physical and mathematical description.If we can x the imaginary parts we can thenunambiguously dene the amplitude and thephase byA(t) =

s2+s2i ; (t) = arctan ssi.How to dene the imaginary part? quadrature method analytic signal method - si(t) is the Hilberttransform of s(t).instantaneous frequency of s(t) is

(t)instantaneous bandwidth is |a

(t)a(t)|.narrow band instead of mono-component.The most popular example iss(t) = Acos (t).Some paradoxes exist. To avoid them, Huang,et al ( N.E. Huang, Z. Shen, and S.R. Long, Anew view of nonlinear water waves: the Hilbertspectrum, Annu. Rev. Fluid Mech. 31 (1999),417457) have developed a method, termedthe Hilbert view, for studying nonstationaryand nonlinear data in nonlinear mechanics.Main tools Empirical Mode Decomposition Method (EMD) Intrinsic Mode Functions (IMFs) Hilbert Transform.Main ideaFirst identify the appropriate time scales thatwill reveal the physical characteristics of thestudied system and then to decompose thefunction into modes intrinsic to the function.The EMD method was motivated from thesimple assumption that any data consists ofdierent simple intrinsic mode oscillations. the time between successive zero-crossings; the time between successive extrema; the time between successive curvature ex-trema.Intrinsic Mode Function (IMF)(a) The number of local extrema of andthe number of its zero-crossings must ei-ther be equal or dier at most by one.(b) At any time t, the mean value of theupper envelope (determined by the localmaxima) and the lower envelope (deter-mined by the local minima) is zero.The instantaneous frequency obtained by ap-plying the Hilbert transform to functions ofthe above type is very well localized in thetime-frequency domain and reveals importantcharacteristics of the signal. Figure 1 is re-produced from Huang et al and compares theHilbert, Wavelet and Fourier spectra for a fre-quency and amplitude modulated signals thatare IMFs. Figure 2 illustrates the EMD methodfor a physical signal of length of day.Comparison of SpectraHilbert, Wavelet, FourierEMF Decompositonfor Length of Day Data(Huang, et al)The EMD Method to decompose f as alinear combination of IMFs nSet n := 0 and f0 := f.Step 1. Set h0 := fn and k := 0.Step 2. Construct the Upper Envelope forhkIdentify all the local extrema, then t allthe local maxima by a cubic spline inter-polant for use as the upper envelope U(t).Step 3. Construct the Lower Envelope forhkProceed in a similar way as Step 2, butreplacing the local maxima as knots, by thelocal minima.Step 4. SiftIdeally, the functions U and L should en-velop the data between them, i.e., L(t) hk(t) U(t), for all t. Their mean is de-noted by mk(t), and the k-th component isdened ashk+1 := hk mk.If hk+1 is not an IMF,then increment k, return to Step 2 andrepeat the procedure (i.e., withhk+1 in place of hk).else dene the IMF n as hk+1 and theresidual fn+1 as fn n. If a con-vergence criteria is not met, in-crement n and return to Step 1.Convergence criteria typically con-sist of testing if the residual is eithersmaller than a predetermined valueor is a monotonic function)In view of the possible deciency of the en-velopes and to make the algorithm more ef-cient Huang et al suggest that the stoppingcriterion on the inner loop be changed fromthe condition that the resulting function tobe an IMF to the single condition that thenumber of extrema equals to zero-crossingsalong with visual assessment of the iterates.f(t) =Nn=1n +rN+1where rN+1 = fN+1 is the residual.f(t) =

Nn=1An(t) exp(i

n(t) dt)

.The residue rN is left on purpose, for it is eithera monotonic function or a constant. Here nis the instantaneous frequency dened as dndtwhere n := arctan(Hn/n)Analysis and commentsThe rst step is to choose a time scale whichis inherent in the function f(t) and has physicalmeaning. the set of inection points; the set of zeros of the function (f(t) cos kt).The second step is to extract some special(with respect to the already chosen time scale)functions, which in the EMD case are calledIMFs. The denition of an IMF, althoughsomewhat vague, has two parts (a) the num-ber of the extrema equals the number of thezeros and (b) the upper and lower envelopesshould have the same absolute value.If we restrict to the larger class of functionswhich just satisfy condition (a), then we areable to provide a characterization in terms ofsolutions of self-adjoint ODEs.Properties of the solutions of a self-adjointODEs.An ODE is called self-adjoint if can be writtenin the formddx

Pdudx

+Qu = 0, (1)for x (a, b) (a and b nite or innite), whereP > 0 and Q is continuous. These equationarise in vibration problems of continuum me-chanics. The solutions are harmonic standingwaves and it is commonly assumed that anywave motion can be resolved into simple har-monic waves (which can be proved mathemat-ically).I.(Max and Min) If Q > 0, then any solution of(1) has exactly one maximum or minimum be-tween successive zeros. (This is the importantcondition for IMFs.)II.(The Prufer substitution) Let Pu

= r(x) cos (x),u(x) = r(x) sin(x), then (1) is equivalent to alinear system of rst order ODEs:ddx= Q(x) sin2 + 1P(x)cos2()drdx= 12

1P(x) Q(x)

r sin2.If Q(x) is positive, then the instantaneous fre-quency dened through the Prufer substitutionis always a positive number. If we want to posesome conditions on the envelopes (on the am-plitude r in the above substitution), then theywill aect only the second equation.Relation between the self-adjoint ODEsand IMFsTheorem. Let f be a given C2[a, b], (a, b R) and all of the zeros of f and f

be simplezeros. Then the following three conditions areequivalent:(i) The number of the zeros of f and the num-ber of the extrema of f on [a, b] dier at mostby 1;(ii) There exist positive, continuously dieren-tiable functions P and Q s.t. f is a solution ofthe self- adjoint ODE(P(t)f

(t))

+Q(t)f(t) = 0;(iii) There exists a C2[a, b] function g s.t.f(t) = 1Q(t)g

(t), g(t) = P(t)f

(t)for some positive, continuously dierentiablefunctions P and Q.Idea of the Proof:(ii) > (i) Follows from the properties of theself-adjoint ODEs.(i) > (ii) # of zeros = # of extrema = M.Then consider the interpolation tableT =

{tj, |f(tj)|, 0}Mj=1, {zj, |aj|, 0}Mj=1

,tj are the extremal points,zj are the zeros (t1 < z1 < t2 < z2 < . . .),and aj are arbitrary positive numbers.Hermite interpolant E(f, t)a for T s.t.(a) The only extrema of E(f, t)a are at thepoints tj, zj for j = 1, 2, . . . , M;(b) The only points in common of the func-tions E(f, t)a and f(t) are tj, j = 1, 2, . . . , M.In general, any piece-wise function E(f, t)a madeby pieces j(t) that satisfy the conditions(y1) = v1, (y2) = v2

(y1) =

(y2) = 0, |

(t)| > 0 for t (y1, y2)will work. (Meyers wavelets.)By choosing aj(j = 1, . . . , M) we can insurethat E(f, t)a satises (b) as well.We haveg(t) := f(t)E(f, t)a= sin(t)for some function (t) s.t. (tj) = 2kj+12 and(zj) = lj for some integers kj and lj.Moreover we can adjust the ajs so that (t) isC1function and

(t) > 0 on the whole interval(a, b).Let r(t) := E(f, t)a be an envelope for f withproperties (a)-(b), and

(t) > 0. Then thefunctionh(t) := r(t) cos (t)has zeros at the points tj, j = 1, 2, . . . , M, andno other zeros on [a, b].The only zeros of h

(t) on [a, b] are the pointszj, j = 1, 2, . . . , M.By observing the changes of the signs of thefunctions f, f

, h, and h

, we can prove thatP(t) := h(t)f

(t), Q(t) := h

(t)f(t)are well dened, strictly positive(f(t1) > 0)C1[a, b] functions.Compa