A NOVEL FORWARD-PDE APPROACH AS AN ALTERNATIVE TO EMPIRICAL MODE DECOMPOSITION

A Novel Forward-Pde Approach as an Alternative to Empirical Mode Decomposition

Abstract

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 involved in the EMD method have been obstacles for mathematical modelling. Here, we propose a novel forward heat equation approach to represent the mean envelope and sifting process. This new model can provide a better mathematical analysis of classical EMD as well as identifying its limitations. Our approach achieves a better performance for a “mode-mixing” signal as compared to the classical EMD approach and is more robust to noise. Furthermore, we discuss the ability of EMD to separate signals and possible improvements by adjusting parameters.

A NOVEL FORWARD-PDE APPROACH AS AN ALTERNATIVE TO EMPIRICAL MODE DECOMPOSITION

 Heming Wang†, Richard Mann⋆ and Edward R. Vrscay† †Department of Applied Mathematics, ⋆David R. Cheriton School of Computer Science University of Waterloo, Ontario, Canada N2L3G1

Index Terms—  Empirical Mode Decomposition, Spectral Analysis, Partial Differential Equation, Heat Equation

1 Introduction

Huang et al. [1] introduced the empirical mode decomposition (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 process to separate different modes of oscillation which are referred to as Intrinsic Mode Functions (IMF). EMD is essentially a decomposition algorithm that extracts the highest local frequency components 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 frequencies. The amplitudes and instantaneous frequencies may then be combined to produce a local time-frequency 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 earthquakes [1], medical experiments [3] and rotating machinery [4]. There are some kinds of signals, however, for which the sifting process fails to separate into different oscillatory modes. Because most of the work on EMD has focused on algorithms as opposed to mathematical analysis, there has been very little work on developing a rigorous theoretical basis for EMD as well as an understanding of why it fails for certain kinds of signals. The need for a mathematical model which explains the principle of EMD and provides a description of the region where it can work effectively has been the motivation for this work.

One major obstacle for mathematical modelling of EMD is the interpolation process employed by the algorithm. In perhaps the first effort to model the EMD interpolation procedure [5], a rather large number of variables were encountered. This, plus the fact that iteration is involved, makes it difficult to arrive at an accurate expression in the model. In this paper, we propose a forward heat PDE approach to solve these problems. Instead of taking the average of two envelopes, a forward heat equation is introduced to construct a mean envelope. This mean envelope can be viewed as the result of passing a signal through a smooth filter - in this case, a Gaussian filter. After obtaining the mean envelope, 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, but provides a solid mathematical basis for the method.

The remainder of the paper is organized as follows. In Section 2, we briefly review the details of classical EMD method and some related work. In Section 3, we introduce the new forward-heat PDE approach. In Section 4, mathematical interpretation of EMD is provided, and the limitations are analyzed. In Section 5, numerical implementation and experimental results are presented.

2 Related 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 .
2. Interpolate between maximal points to obtain upper envelope function and between minimal points to obtain lower envelope function (x).
3. Computer the local mean: .
4. Extract mean from signal: .
5. If is not an IMF, iterate 3) and 4) until it is.
6. After finding IMF, subtract it from and repeat Step 2 to obtain the residual.
There are, however, several drawbacks [6]:
Vague Definition of IMF which presents obstacles in implementation. 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 necessary to choose appropriate stopping criteria for the sifting process.
Boundary Effects: Proper boundary conditions are necessary in order to minimize errors at the boundaries. Otherwise, there can be “tweaking” at the endpoints.
Mode Mixing: Whenever the signal contains riding waves, some frequency components will vanish 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 , the upper and lower envelopes of a function are defined as follows,

 Uδ(x)=sup|y|<δh(x+y),   Lδ(x)=inf|y|<δh(x+y). (1)

After Taylor expansions are applied to the envelopes, the mean envelope is defined as

 mδ(x)=12(Uδ(x)+Lδ(x))≈h(x)+δ22h′′(x). (2)

The sifting process – the process to obtain the Intrinsic Mode Function (IMF) – is then defined as follows,

 hn+1(x)=hn(x)−mδ(x),   h0(x)=S(x). (3)

Using the following Taylor expansion in ,

 hn+1=h(x,t+Δt)=hn+Δt∂h∂t+O(Δt2), (4)

the authors arrive at the following PDE,

 ∂h∂t+1δ2h+12∂2h∂x2=0,     h(x,0)=S(x), (5)

which is known as a backward heat equation since the diffusivity constant is negative. (Note that the initial condition, , to this PDE is the original signal .) Unfortunately, there are several drawbacks to this approach:
1. The parameter , which is chosen empirically, has a significant influence on the result. For a generalized signal , the solution is

 h(x,T)=∑ke(ω2k2−1δ2)TAkcos(ωkx+ϕk). (6)

As increases, the amplitudes of components with lower frequencies will be decreased at each step and therefore vanish at the end of the algorithm. Only the higher frequencies survive. Therefore, choosing requires an additional knowledge of the signal.
2. Even if we extract correct the frequency component from the signal, we cannot guarantee that the amplitude of the component is correct. In order to distinguish two frequency components, one sometimes has to decrease their amplitudes to very small values.

3. As mentioned earlier, Eq.  (5) is a backward heat/diffusion equation. Because the diffusivity parameter is negative, the evolution of a signal will be opposite to that of a signal under the standard (forward) diffusion PDE – signals become less smooth and local amplitudes grow exponentially. As expected, numerical methods also suffer from instability.

3 Forward Pde Algorithm

3.1 Introduction

We now outline our PDE-based method to perform a new type of interpolation in the EMD algorithm. The idea is very simple: Instead of taking the average of two envelope functions of a signal to produce a mean (Step 3 in Section 2.1), we proceed as follows. For prescribed values of the diffusivity constant and time (which can be adjusted, as will be discussed below), solve the following initial value problem for the heat/diffusion equation,

 ∂h∂t=a∂2h∂x2,    h(x,0)=S(x) (7)

and define the mean curve of to be . (Of course, is equivalent to the convolution of with the Gaussian function with standard deviation .) One of the primary motivations for this definition is that the time rate of change of is zero at spatial inflection points of . An example is shown in Figure 1. This is the basis of the following modified EMD algorithm applied to a signal :
1. Initialize: Let and set .
2. Find mean of : Solve the PDE in (7) for for . Then define .
3. Extract mean: Define .
4. If is not an IMF, let , and go to Step 2.

3.2 Parameter selection

3.2.1 How to choose parameter a

Parameter is crucial to determine the mean envelope so it must be carefully chosen. To ensure that the mean envelope is always within the range of the signal amplitude, it is necessary that . As mentioned in [10], if the sampling rate is not sufficiently large, sampling effects will cause a loss of accuracy. We must assume that , the maximum frequency to be extracted, satisfies . This implies that . It is then safe to set . Ideally the parameter should be set to . There are two practical approach to estimate : (1) autocorrelation, (2) zero-crossing rate.

3.2.2 The pair of parameters T and N

The parameters and determine the shape of the mean curve, and represents the number of iterations. In order to be able to separate the high-frequency component from the other components, we impose the following condition,

 [1−e−f20T1−e−T]N=δ, (8)

where , an adjustable parameter, is close to zero. Here, is the cutoff frequency ratio: If we let denote the ratio between two frequency components (lower/higher), then the algorithm may fail to separate the components , i.e., if the two components are too close to each other. When , the ratio of the norms of the lower- and higher-frequency components will satisfy . Therefore, finding the optimal parameters reduces to the following two problems,

 f0=⎡⎢⎣log(1−(δα)1N(1−ϵ))logϵ⎤⎥⎦1/2 ,  (1−e−T)N=1−δ. (9)

Under these restrictions, we seek to maximize . Assume that , and let . Also assume that the two signals are sufficiently separated, i.e., .

The results of one numerical experiment, with and , are shown in Figure 2. From these results, we conclude that the theoretical cutoff frequency should be .

4 Mathematical explanation of EMD and its limitations

4.1 Forward PDE interpretation of EMD

Here we consider the following simple model which is sufficiently general to represent many realistic signals,

 S(x)=K∑k=1Akcos(ωkx+ϕk)+C=K∑k=1sk(x)+C. (10)

Solving Eq. (7) for first mean envelope PDE, we obtain

 ma(x,T)=K∑k=1e−aω2kTsk+C (11)

After iterations, our modified EMD algorithm yields the following result for the th cosine component ,

 hk,N=(1−e−aω2kT)hk,N−1=(1−e−aω2kT)Nsk. (12)

Now suppose, without loss of generality, that . It is easy to show that for sufficiently large,

 hN=N∑k=1(1−e−aω2kT)Nsk≃(1−e−aω2KT)NsK≃sK, (13)

where the final approximation is valid for sufficiently large. By choosing the appropriate set of parameters, the IMF extracted after iterations will be (at least approximately) the highest-frequency component .

4.2 Border effect

Border effects arise mostly at the mean curve procedure , which involves the solution of the heat PDE (Eq. (5)). Unlike the classical approach, we can control the boundary effect by imposing appropriate boundary condition on the PDE. For example, by assuming the signal is periodic, or fixed at both ends.

4.3 View of Filter

As stated by Fladrin [11], the EMD algorithm is equivalent 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 equivalent to passing it through a high-pass filter.

5 Numerical Results

5.1 Two-mode mixing

This experiment addresses the mode mixing separation problem. The signal is built by concatenating two sinusoids with different frequencies, as shown in Figures 5, 4 and 3. Unlike classical 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 frequency details. Classical EMD fails to separate these different modes.

5.2 Performance Measure of Separation Ability

As shown 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 two-component signal, , is the amplitude ratio and and is the frequency ratio, with . Denote the frequency components as and . Now use the following performance measure for separation capability:

 PM=||IMF1−sin(2πx)||L2||sin(2πx)||L2 (14)

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 classical EMD algorithm. The mean curve of a signal is obtained by evolving the signal with the heat/diffusion equation and therefore avoids 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 according to the properties of the signal being analyzed.

References

• [1] N. E. Huang, S. Zheng, S. R. Long, M. C. Wu, H. H. Shih, Q. Zheng, N. C. Yen, C. C. Tung, and H. H. Liu, “The empirical mode decomposition and the hilbert spectrum for nonlinear and non-stationary time series analysis,” in Proceedings 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. Yan, and Z. Wu, “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 decomposition,” IEEE Transactions on Information Technology in Biomedicine, vol. 16, no. 6, pp. 1135–1142, 2012.
• [4] Y. Lei, J. Lin, Z. He, and M. J. Zuo, “A review on empirical mode decomposition in fault diagnosis of rotating machinery,” Mechanical Systems and Signal Processing, vol. 35, no. 1, pp. 108–126, 2013.
• [5] H. El, S. Diop, R. Alexandre, and A. O. Boudraa, “Analysis of intrinsic mode functions: a pde approach,” IEEE Signal Processing 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,” Interdisciplinary Mathematics, vol. 5, pp. 1–26, 2005.
• [7] Z. Wu 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 framework 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 Processing, 2009. ICASSP 2009. IEEE International Conference on. IEEE, 2009, pp. 3429–3432.
• [10] G. Rilling and P. Flandrin, “On the influence of sampling on the empirical mode decomposition,” in Acoustics, Speech and Signal Processing, 2006. ICASSP 2006 Proceedings. 2006 IEEE International Conference on. IEEE, 2006, vol. 3, pp. III–III.
• [11] P. Flandrin, G. Rilling, and P. Goncalves, “Empirical mode decomposition as a filter bank,” IEEE signal processing 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 transactions on signal processing, vol. 56, no. 1, pp. 85–95, 2008.
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters