A Novel Foward-PDE Approach as an Alternative to Empirical Mode Decomposition

A NO VEL FOR W ARD-PDE APPR OA CH AS AN AL TERNA TIVE TO EMPIRICAL MODE DECOMPOSITION Heming W ang † , Richar d Mann ? and Edwar d R. Vrscay † † Department of Applied Mathematics, ? David R. Cheriton School of Computer Science Uni versity of W aterloo, Ontario, Canada N2L3G1 ABSTRA CT In this paper we present a mathematical model of the Empirical Mode Decomposition (EMD). Although EMD is a powerful tool for signal processing, the algorithm itself lacks an appropriate theoretical basis. The interpolation and iteration processes inv olved in the EMD method have been obstacles for mathematical modelling. Here, we propose a nov el forward heat equation approach to represent the mean en velope and sifting process. This new model can provide a better mathematical analysis of classical EMD as well as identifying its limitations. Our approach achiev es a better performance for a “mode-mixing” signal as compared to the classical EMD approach and is more robust to noise. Further - more, we discuss the ability of EMD to separate signals and possible improv ements by adjusting parameters. Index T erms — Empirical Mode Decomposition, Spectral Analysis, Partial Dif ferential Equation, Heat Equation 1. INTR ODUCTION Huang et al. [1] introduced the empirical mode decompo- sition (EMD) in 1998 as a tool to analyze linear and non- stationary signals. EMD has been applied quite successful in science and engineering. It treats a signal as a mixture of mono-components and applies a sifting pr ocess to separate different modes of oscillation which are referred to as Intrin- sic Mode Functions (IMF). EMD is essentially a decomposi- tion algorithm that e xtracts the highest local frequenc y com- ponents from the signal for each IMF . A repeated application produces a decomposition of a signal into of components with decreasing frequency . The Hilbert transform is then applied to each component in order to determine instantaneous fre- quencies. The amplitudes and instantaneous frequencies may then be combined to produce a local time-frequenc y analysis of the signal. The ability of EMD to capture intrinsic physical features for non-stationary signals [2] has been demonstrated for real- world signals with limited frequency components, such as signals obtained from earthquak es [1], medical experiments [3] and rotating machinery [4]. There are some kinds of sig- nals, ho wev er, for which the sifting process fails to separate into different oscillatory modes. Because most of the w ork on EMD has focused on algorithms as opposed to mathemat- ical analysis, there has been v ery little work on de veloping a rigorous theoretical basis for EMD as well as an understand- ing of why it fails for certain kinds of signals. The need for a mathematical model which explains the principle of EMD and pro vides a description of the region where it can w ork effecti vely has been the moti vation for this w ork. One major obstacle for mathematical modelling of EMD is the interpolation process employed by the algorithm. In perhaps the first ef fort to model the EMD interpolation proce- dure [5], a rather lar ge number of variables were encountered. This, plus the fact that iteration is in volved, makes it difficult to arri ve at an accurate expression in the model. In this pa- per , we propose a forward heat PDE approach to solve these problems. Instead of taking the av erage of two env elopes, a forward heat equation is introduced to construct a mean en ve- lope. This mean en velope can be viewed as the result of pass- ing a signal through a smooth filter - in this case, a Gaussian filter . After obtaining the mean en velope, we repeat the same steps as those of the original EMD sifting process in order to extract the IMFs. Our approach generates results which are similar to classical EMD, b ut provides a solid mathematical basis for the method. The remainder of the paper is organized as follows. In Section 2, we briefly revie w the details of classical EMD method and some related work. In Section 3, we introduce the new forward-heat PDE approach. In Section 4, mathemat- ical interpretation of EMD is pro vided, and the limitations are analyzed. In Section 5, numerical implementation and exper - imental results are presented. 2. RELA TED WORKS 2.1. CLASSICAL EMD ALGORITHM The classical EMD algorithm may be summarized as follows: 1. Find all local maximal and minimal points of the signal S ( x ) . 2. Interpolate between maximal points to obtain upper en ve- lope function E u p per ( x ) and between minimal points to obtain lower en velope function E l ower (x). 3. Computer the local mean: m ( x ) = 1 2 ( E u p per ( x ) + E l ower ( x )) . 4. Extract mean from signal: c ( x ) = S ( x ) − m ( x ) . 5. If c ( x ) is not an IMF , iterate 3) and 4) until it is. 6. After finding IMF , subtract it from S ( x ) and repeat Step 2 to obtain the residual. There are, howe ver , several dra wbacks [6]: V ague Definition of IMF which presents obstacles in im- plementation. For an IMF , the number of extrema and zero- crossings must differ at most by one. In addition, the “local mean” of the IMF should be close to zero. It is therefore nec- essary to choose appropriate stopping criteria for the sifting process. Boundary Effects: Proper boundary conditions are neces- sary in order to minimize errors at the boundaries . Otherwise, there can be “tweaking” at the endpoints. Mode Mixing: Whenev er the signal contains riding waves, some frequency components will v anish after performing EMD. In an effort to solve this problem, Huang introduced a new method called ensemble empirical mode decomposition [7], in which Gaussian noise is first added, and the signal then denoised. 2.2. Backward Heat Equation As mentioned earlier , interpolation represents an obstacle in the mathematical modelling of EMD [8]. A PDE approach was proposed in [5, 9] to overcome this obstacle. Here, for a prescribed δ > 0, the upper and lo wer en velopes of a function h ( x ) are defined as follows, U δ ( x ) = sup | y | < δ h ( x + y ) , L δ ( x ) = inf | y | < δ h ( x + y ) . (1) After T aylor expansions are applied to the en velopes, the mean en velope is defined as m δ ( x ) = 1 2 ( U δ ( x ) + L δ ( x )) ≈ h ( x ) + δ 2 2 h 00 ( x ) . (2) The sifting process – the process to obtain the Intrinsic Mode Function (IMF) – is then defined as follows, h n + 1 ( x ) = h n ( x ) − m δ ( x ) , h 0 ( x ) = S ( x ) . (3) Using the following T aylor expansion in t , h n + 1 = h ( x , t + ∆ t ) = h n + ∆ t ∂ h ∂ t + O ( ∆ t 2 ) , (4) the authors arriv e at the following PDE, ∂ h ∂ t + 1 δ 2 h + 1 2 ∂ 2 h ∂ x 2 = 0 , h ( x , 0 ) = S ( x ) , (5) which is kno wn as a backwar d heat equation since the dif- fusivity constant is negati ve. (Note that the initial condition, h ( x , 0 ) , to this PDE is the original signal S ( x ) .) Unfortunately , there are sev eral drawbacks to this approach: 1. The parameter δ , which is chosen empirically , has a sig- nificant influence on the result. For a generalized signal s = ∑ k A k cos ( ω k x + φ k ) , the solution is h ( x , T ) = ∑ k e ( ω 2 k 2 − 1 δ 2 ) T A k cos ( ω k x + φ k ) . (6) As T increases, the amplitudes of components with lo wer fre- quencies ω < √ 2 δ will be decreased at each step and therefore vanish at the end of the algorithm. Only the higher frequen- cies ω ≥ √ 2 δ surviv e. Therefore, choosing δ requires an addi- tional knowledge of the signal. 2. Even if we extract correct the frequency component from the signal, we cannot guarantee that the amplitude of the com- ponent is correct. In order to distinguish two frequency com- ponents, one sometimes has to decrease their amplitudes to very small v alues. 3. As mentioned earlier , Eq. (5) is a backward heat/dif fusion equation. Because the diffusi vity parameter is negativ e, the ev olution of a signal will be opposite to that of a signal un- der the standard (forward) diffusion PDE – signals become less smooth and local amplitudes grow exponentially . As ex- pected, numerical methods also suffer from instability . 3. FOR W ARD PDE ALGORITHM 3.1. Introduction W e now outline our PDE-based method to perform a new type of interpolation in the EMD algorithm. The idea is very sim- ple: Instead of taking the a verage of two en velope functions of a signal S ( x ) to produce a mean (Step 3 in Section 2.1), we proceed as follows. For prescribed v alues of the diffusi vity constant a > 0 and time T > 0 (which can be adjusted, as will be discussed belo w), solve the follo wing initial value problem for the heat/diffusion equation, ∂ h ∂ t = a ∂ 2 h ∂ x 2 , h ( x , 0 ) = S ( x ) (7) and define the mean curve of S ( x ) to be m ( x ) = h ( x , T ) . (Of course, m ( x ) is equi valent to the con volution of S ( x ) with the Gaussian function with standard deviation a .) One of the primary motiv ations for this definition is that the time rate of change of h ( x , t ) is zero at spatial inflection points of h . An example is shown in Figure 1. This is the basis of the following modified EMD algorithm applied to a signal S ( x ) : 1. Initialize: Let n = 0 and set h 0 ( x , 0 ) = S ( x ) . 2. Find mean of h n ( x , 0 ) : Solve the PDE in (7) for h n ( x , t ) for 0 ≤ t ≤ T . Then define m n ( x ) = h n ( x , T ) . 3. Extract mean: Define c n ( x ) = h n ( x , 0 ) − h n ( x , T ) . 4. If c n is not an IMF , let h n + 1 ( x , 0 ) = c n ( x ) , n → n + 1 and go to Step 2. Fig. 1 . Mean Env elope Obtained by Forw ard Heat Equation 3.2. Parameter selection 3.2.1. How to choose parameter a Parameter a is crucial to determine the mean en velope so it must be carefully chosen. T o ensure that the mean en velope is al ways within the range of the signal amplitude, it is neces- sary that a ≤ 1 ω 2 . As mentioned in [10], if the sampling rate f s is not sufficiently large, sampling effects will cause a loss of accuracy . W e must assume that f max , the maximum fre- quency to be extracted, satisfies f max < f s . This implies that 1 4 π 2 f 2 s < 1 ω 2 max . It is then safe to set a = 1 4 π 2 f 2 s . Ideally the pa- rameter a should be set to a = 1 ω 2 max . There are two practical approach to estimate a : (1) autocorrelation, (2) zero-crossing rate. 3.2.2. The pair of parameters T and N The parameters P and T determine the shape of the mean curve, and N represents the number of iterations. In order to be able to separate the high-frequenc y component from the other components, we impose the following condition, " 1 − e − f 2 0 T 1 − e − T # N = δ , (8) where δ > 0, an adjustable parameter, is close to zero. Here, f 0 is the cutoff frequency ratio : If we let 0 < f < 1 denote the ratio between two frequenc y components (lower/higher), then the algorithm may fail to separate the components f > f 0 , i.e., if the two components are too close to each other . When f < f 0 , the ratio of the norms of the lower - and higher-frequenc y components will satisfy || S l ower S higher || < δ . Therefore, finding the optimal parameters reduces to the following tw o problems, f 0 = " log ( 1 − ( δ α ) 1 N ( 1 − ε )) log ε # 1 / 2 , ( 1 − e − T ) N = 1 − δ . (9) Under these restrictions, we seek to maximize f 0 . Assume that a = 1 ω max , and let e − T = ε . Also assume that the two signals are sufficiently separated, i.e., || S l ower S higher || < δ . The results of one numerical experiment, with N = 100 . 0 and T = 10 . 0, are shown in Figure 2. From these results, we conclude that the theoretical cutof f frequency should be f 0 ' 0 . 7. Fig. 2 . Cutoff frequency v aried with parameter T and N 4. MA THEMA TICAL EXPLAN A TION OF EMD AND ITS LIMIT A TIONS 4.1. Forward PDE inter pretation of EMD Here we consider the following simple model which is suffi- ciently general to represent many realistic signals, S ( x ) = K ∑ k = 1 A k cos ( ω k x + φ k ) + C = K ∑ k = 1 s k ( x ) + C . (10) Solving Eq. (7) for first mean en velope PDE, we obtain m a ( x , T ) = K ∑ k = 1 e − a ω 2 k T s k + C (11) After N iterations, our modified EMD algorithm yields the following result for the k th cosine component s k ( x ) , h k , N = ( 1 − e − a ω 2 k T ) h k , N − 1 = ( 1 − e − a ω 2 k T ) N s k . (12) Now suppose, without loss of generality , that ω 1 < ω 2 < · · · < ω K = ω max . It is easy to show that for N sufficiently large, h N = N ∑ k = 1 ( 1 − e − a ω 2 k T ) N s k ' ( 1 − e − a ω 2 K T ) N s K ' s K , (13) where the final approximation is valid for T sufficiently large. By choosing the appropriate set of parameters, the IMF ex- tracted after N iterations will be (at least approximately) the highest-frequency component s K . Fig. 3 . Experiment on mode-mixing signal by classical EMD (left) and forward-PDE approach (right) 4.2. Border effect Border effects arise mostly at the mean curve procedure , which in volv es the solution of the heat PDE (Eq. (5)). Un- like the classical approach, we can control the boundary ef- fect by imposing appropriate boundary condition on the PDE. For example, by assuming the signal is periodic, or fixed at both ends. 4.3. V iew of Filter As stated by Fladrin [11], the EMD algorithm is equi valent to a set of filter banks, which is justified in our PDE method. In each iteration of our PDE approach, the mean of the signal is obtained by passing it through a low-pass filter . Subsequent subtraction of the mean from the signal is therefore equi v alent to passing it through a high-pass filter . 5. NUMERICAL RESUL TS 5.1. T wo-mode mixing This experiment addresses the mode mixing separation prob- lem. The signal is built by concatenating two sinusoids with different frequencies, as sho wn in Figures 3 to 5. Unlike clas- sical EMD, the forward-PDE approach can distinguish the two modes and produce a reasonable separation. As such, it can extract features from mode-mixed signals and obtain better instantaneous frequenc y details. Classical EMD f ails to separate these different modes. Fig. 4 . Hilbert-Huang Spectrum for mode-mixing signal us- ing classical EMD approach. Fig. 5 . Hilbert-Huang Spectrum for mode-mixing signal us- ing forward-PDE approach. 5.2. PERFORMANCE MEASURE OF SEP ARA TION ABILITY As sho wn in [12], by applying EMD to mixtures of two cosine signals, we can examine the separation capability for different frequency component rations. Consider the following tw o- component signal, S ( x ) = cos ( 2 π x ) + α cos ( 2 π f x ) , α is the amplitude ratio and α ∈ ( 10 − 2 , 10 2 ) and f is the frequency ratio, with f ∈ ( 0 , 1 ) . Denote the frequency components as S 1 ( x ) = cos ( 2 π x ) and S 2 ( x ) = α cos ( 2 π f x ) . Now use the fol- lowing performance measure for separation capability: PM = || I M F 1 − sin ( 2 π x ) || L 2 || sin ( 2 π x ) || L 2 (14) Fig. 6 . Left: Performance Measure Re garding α and f for classical EMD. Right: Performance Measure Regarding α and f for forward-PDE approach. The results are presented in Figure 6. It should be noted that the performance of our PDE approach is similar to that of the classical EMD method. 6. CONCLUSIONS This paper presents a forward-PDE modification of the clas- sical EMD algorithm. The mean curve of a signal is obtained by e volving the signal with the heat/dif fusion equation and therefore av oids any complicated methods of extracting local maxima or minima. Our approach provides a mathematical interpretation of the EMD algorithm as well as its limitations. It also performs better on “mode-mixed” signals. Our method also allows the parameters in the PDE to be adjusted accord- ing to the properties of the signal being analyzed. 7. REFERENCES [1] N. E. Huang, S. Zheng, S. R. Long, M. C. W u, H. H. Shih, Q. Zheng, N. C. Y en, C. C. T ung, and H. H. Liu, “The empirical mode decomposition and the hilbert spectrum for nonlinear and non-stationary time series analysis, ” in Pr oceedings of the Royal Society of London A: mathematical, physical and engineering sciences . The Royal Society , 1998, vol. 454, pp. 903–995. [2] T . Y . Hou, M. P . Y an, and Z. W u, “ A variant of the emd method for multi-scale data, ” Advances in Adaptive Data Analysis , vol. 1, no. 04, pp. 483–516, 2009. [3] V . Bajaj and R. B. Pachori, “Classification of seizure and nonseizure eeg signals using empirical mode decompo- sition, ” IEEE T ransactions on Information T echnology in Biomedicine , vol. 16, no. 6, pp. 1135–1142, 2012. [4] Y . Lei, J. Lin, Z. He, and M. J. Zuo, “ A revie w on empir- ical mode decomposition in fault diagnosis of rotating machinery , ” Mechanical Systems and Signal Pr ocess- ing , vol. 35, no. 1, pp. 108–126, 2013. [5] H. El, S. Diop, R. Alexandre, and A. O. Boudraa, “ Anal- ysis of intrinsic mode functions: a pde approach, ” IEEE Signal Pr ocessing Letters , vol. 17, no. 4, pp. 398–401, 2010. [6] N. E. Huang et al., “Introduction to the hilbert-huang transform and its related mathematical problems, ” In- ter disciplinary Mathematics , vol. 5, pp. 1–26, 2005. [7] Z. W u and N. E. Huang, “Ensemble empirical mode decomposition: a noise-assisted data analysis method, ” Advances in adaptive data analysis , vol. 1, no. 01, pp. 1–41, 2009. [8] H. El, S. Diop, R. Alexandre, and V . Perrier, “ A pde based and interpolation-free frame work for modeling the sifting process in a continuous domain, ” Advances in Computational Mathematics , vol. 38, no. 4, pp. 801– 835, 2013. [9] H. El, S. Diop, R. Alexandre, and A. O. Boudraa, “ A pde characterization of the intrinsic mode functions, ” in Acoustics, Speech and Signal Pr ocessing, 2009. ICASSP 2009. IEEE International Conference on . IEEE, 2009, pp. 3429–3432. [10] G. Rilling and P . Flandrin, “On the influence of sam- pling on the empirical mode decomposition, ” in Acous- tics, Speech and Signal Pr ocessing, 2006. ICASSP 2006 Pr oceedings. 2006 IEEE International Confer ence on . IEEE, 2006, vol. 3, pp. III–III. [11] P . Flandrin, G. Rilling, and P . Goncalv es, “Empirical mode decomposition as a filter bank, ” IEEE signal pro- cessing letters , vol. 11, no. 2, pp. 112–114, 2004. [12] G. Rilling and P . Flandrin, “One or two frequencies? the empirical mode decomposition answers, ” IEEE trans- actions on signal pr ocessing , v ol. 56, no. 1, pp. 85–95, 2008.