Copyright © 2008 Yannis Kopsinis and Steve McLaughlin. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract

Empirical mode decomposition (EMD) is a signal analysis method which has received much attention lately due to its application in a number of fields. The main disadvantage of EMD is that it lacks a theoretical analysis and, therefore, our understanding of EMD comes from an intuitive and experimental validation of the method. Recent research on EMD revealed improved criteria for the interpolation points selection. More specifically, it was shown that the performance of EMD can be significantly enhanced if, as interpolation points, instead of the signal extrema, the extrema of the subsignal having the higher instantaneous frequency are used. Even if the extrema of the subsignal with the higher instantaneous frequency are not known in advance, this new interpolation points criterion can be effectively exploited in doubly-iterative sifting schemes leading to improved decomposition performance. In this paper, the possibilities and limitations of the developments above are explored and the new methods are compared with the conventional EMD.

1. Introduction

The empirical mode decomposition (EMD) method [1] is an algorithm for the analysis of multicomponent signals [2] that works by breaking them down into a number of amplitude and frequency modulated (AM/FM) zero mean signals, termed intrinsic mode functions (IMFs). In contrast to conventional decomposition methods, which perform the analysis by projecting the signal under consideration into a number of predefined basis vectors, EMD expresses the signal as an expansion of basis functions which are signal-dependent, and are estimated via an iterative procedure called sifting. This attribute of EMD potentially leads to a number of merits. To name a few: it can be applied regardless of the nonstationary and/or nonlinear characteristics of the signal under consideration. The results are not prejudiced by the predetermined basis, a fact often leading to IMFs which preserve the physical meanings of the intrinsic processes underlying the signal. Moreover, the resulting IMFs are zero-mean narrowband functions well suited for meaningful instantaneous frequency (IF) estimates via the Hilbert transform or other alternative techniques [3]. As a result, EMD in conjunction with IF estimation offers an alternative path towards time-frequency signal representation.

The main drawback of EMD is the lack of a strong theoretical analysis capable of evaluating and predicting EMD behaviour in generalized signal conditions. However, recently, Rilling and Flandrin [4] have made some initial steps in this direction by theoretically analyzing the EMD outcomes in a two-tone signal case. Although the signal can be considered simplistic, the analysis resulted in important conclusions. It was shown that there is quite a wide range of two-tone combinations, for which EMD results do not, at least directly, agree with intuition and physical interpretation. Revealing such a limitation should definitely be a target for future EMD variants.

The lack of theoretical developments on EMD has restricted the potential for improvements of the method itself. Literally, the current popular variation of EMD is roughly the same as the one proposed in the original EMD paper [1], with the attempts for development of novel variants being limited to less than a handful [5–8]. A recent study on EMD [9], even though it is not analytical, revealed specific aspects that offer insight to its performance. More specifically, information about improved criteria for interpolation points selection was extracted based on a genetic algorithm (GA) optimization approach. A recently presented novel EMD variant [8] called doubly-iterative EMD (DI-EMD) (DI-EMD Matlab functions can be downloaded from http://www.see.ed.ac.uk/~ykopsini/emd/emd.html) succeeds in estimating interpolation points in agreement with the optimized criteria derived in [9] leading to enhanced overall decomposition performance. In this paper, the performance and behaviour of DI-EMD is further investigated mainly through the Rilling two-tone signal model [4].

2. Emd Method

EMD [1] adaptively decompose a multicomponent signal [2] into a number of zero-mean, narrowband IMFs , ,(1) Each one of the IMFs, say the th one , is estimated with the aid of an iterative process, called sifting, applied to the residual multicomponent signal:(2)

During the th iteration of the sifting process, an estimate of the th IMF is computed as follows:(3)where, is the temporal estimate of the th IMF at the th iteration and is an estimate of the local mean of .

Usually, the local mean is estimated as the average of two envelopes, an upper envelope and a lower one, which enfold the corresponding IMF estimate . In general, the envelopes are constructed in agreement with the following algorithm:

First, some time instances , called nodes, which correspond to the upper and the lower envelope, respectively, are specified according to some criteria. These time instances indicate the positions , at which the upper and lower envelopes “touch” the temporal IMF estimate as can be seen in Figure 1. In order to succeed in this , serve as interpolation points in a piecewise polynomial interpolation scheme, usually cubic spline interpolation.

Figure 1: Quantities related to the EMD method.

Finally, the current estimate of the local mean is given by(4)

After sifting iterations, which is a number chosen either statically or dynamically according to specific criteria [10, 11], the sifting process is concluded and the th IMF is set equal to . Alternatively, the intermediate estimates of the local mean can be summed up together forming a total mean envelope:(5)which is actually an enhanced estimate of the local mean of the residual signal under consideration . According to such a reformulation, the th IMF can be obtained directly from the expression:(6)Either (3) with (4) or (6) with (5) and (4) are equivalent expressions to the sifting process [9].

In other words, from (6) it can be inferred that EMD considers signals as fast oscillations superimposed on slow oscillations [12] and the sifting process aims to iteratively estimate the slow oscillating signals using (5). As a consequence, the th IMF is an estimate of the fast oscillating component of the signal . Lets say, for example, that consists of AM/FM signals:(7)which have corresponding instantaneous frequencies (IF) . It turns out that the sifting process tries to extract in each time instant this signal among those consisting which has the higher instantaneous frequency. At the same time, the extracted frequency modulated (FM) signal although tends to be narrowband is not necessarily monocomponent permitting in such a way the presence of amplitude modulation (AM). As a result, the fast oscillating signal that the sifting process tries to estimate with IMF is given by (8)where corresponds to . Note that does not necessarily coincide with one of the signals comprising . may consist of parts of these signals depending on which one of them, in specific time instances, has the higher instantaneous frequency. Apparently, the “ideal” local mean of is given by .

In turn, the slow oscillating part is further processed through a number of sifting iterations for its separation to a fast oscillation part (which is the next IMF and a slow oscillating part which serves as input to the next sifting process.

2.1. Example of EMD Application

One of the potential applications of EMD is the analysis of short-duration transient signals. The reason is that essentially, the signal analysis performance of EMD does not depend on the length of the signal and/or the available samples. In principle, the only requirement is the availability of a large enough number of maxima and minima, which depend on the order of the spline interpolation. As long as the above requirement is fulfilled, the full performance of EMD is guaranteed, otherwise, EMD cannot be applied.

The examined transient signal is the one shown in Figure 2(a), which consists of three nonlinear chirp signal components depicted, in the time-frequency plane, in Figure 2(b). The three signal components in the time domain can be seen with solid lines in Figure 3.

Figure 2: (a) A transient signal having 256 samples. (b) The transient signal in the frequency domain.

Figure 3: Three chirp signals and the corresponding IMF estimates.

When the standard EMD with cubic spline interpolation is used for the analysis of the above signal, the result is a number of IMFs with only three of them having significant energy. It turns out that the higher energy estimated IMFs, shown in Figure 3 with dashed lines, roughly coincide with the actual signal components.

Based on the instantaneous frequency estimates of the estimated IMFs, a quite accurate spectrogram can be drawn as it is seen in Figure 4(a). It is important to note that the conventional short-time Fourier transform (STFT) regardless of the adopted window length is not capable of producing such a sharp and detailed spectrogram as it can be seen in Figures 4(a) and 4(b). This is happening due to the transient nature of the specific signal. In other words, if the STFT window is considered long enough to include an adequate number of samples, in order to achieve reliable frequency resolution, the signal itself changes considerable within the time span of the specific window leading to poor time-frequency representation.

Figure 4: (a) Spectrogram using IF estimates of the IMFs. (b) Spectrogram using STFT with kaiser window of length equal to 64 samples. (c) Spectrogram using STFT with kaiser window of length equal to 128 samples.

3. Doubly Iterative EMD

Several variants of EMD can be formed by altering the way that the local mean is obtained. The local mean estimates are determined from the upper and the lower envelope construction, which as discussed above, basically comes down to the adopted criteria for the nodes selection as well as the interpolation scheme used. With respect to the standard version of EMD, usually employed in practice, the maxima and minima, referred to as local extrema, of signal (or in the first iteration) are used as interpolation points and natural cubic splines are used for interpolation. More specifically, and , where the operator denotes the th derivative of function . In [7, 9], it was shown that the local extrema of the IMF estimates, , in each sifting iteration are far from being the optimum choice of interpolation points. It turned out that the decomposition performance was significantly improved if the interpolation points, where set fixed in all the sifting iterations and equal to the extrema of the fast oscillating signal . Hereafter, the extrema above will be called desired and the corresponding nodes are given by and .

At first glance, the observation that the desired extrema perform much better than the standard local extrema may be considered useless in the sense that actually is the signal that the sifting process aims to extract. In that sense, its extrema cannot be known in advance. However, in [9] we saw that it is possible to obtain interpolation points estimates which are closer than the local extrema to the desired ones. This can be realized by adopting, for example, the local extema of a high-pass filtered version of the signal under consideration . The filtering results in a signal with attenuated slow oscillating components leading to improved estimates of the extrema of the desired fast oscillating counterpart. In other words, are preprocessed in each iteration in order to resemble to and then estimate the desired extrema from this. In practice, the above technique exhibits certain difficulties mainly related to the filter cutoff choice which should also be time varying in the case of nonstationary signals. Apart from that, it can be argued that the filtering preprocess compromises the spontaneous, data-driven nature of the EMD operation.

According to DI-EMD, the estimation of the desired extrema is approached from a different viewpoint which is based on the fact that for the estimation of the desired extrema it is not the actual fast oscillating signal that is needed to be known, but its first derivative. As we will see the first derivative of can be effectively estimated directly in each iteration. More importantly, this can be naturally realized within the data-driven framework of EMD.

Figure 5(a) shows an example of a residual signal after IMF extractions (which is denoted as for generalization purposes due to the fact that the procedure described next can be applied not only to the residual signal but to any tentative IMF estimate . The dotted curve depicts the corresponding local mean and the filled circles indicate the local extrema of which lay in the positions where the first derivative of , shown with solid line in Figure 5(b), equals to zero. Moreover, the desired extrema, which are the extrema of the fast oscillating signal component, , are shown with asterisks in Figure 5(b). It can be readily seen that the local extrema are deviated not only from the desired positions but also have been smeared out completely in many cases. can be written as(9)where is the first derivative of the local mean depicted in dotted line. The estimation of the desired extrema can be achieved by estimating and then computing the positions in which it vanishes to zero. Due to the fact that the differentiation operator does not effect the frequency content of the signal, we still expect and to be the fast and the slow oscillating component of respectively. Following the discussion in Section 2, the fast oscillating part and the local mean of a signal can be efficiently estimated through a sifting process. As a result, the application of a predefined number of sifting iterations on produces an estimate of which in turn can be subtracted from leading to an estimate of the first derivative of the fast oscillating part of shown in Figure 5(c). The filled circles indicates the zero-crossing points of which are adopted as estimates of the desired extrema. From now on, the sifting iterations used for the desired extrema estimates will be referred to as internal and the normal EMD sifting iterations used for the current IMF estimate will be called external.

Figure 5: (a) The signal under consideration (solid line) together with its slow oscillating counterpart (dotted line). The corresponding local and desired extrema are depicted with filled circles and asterisks, respectively. (b) Shows the first derivative of the signals in (a). (c) Shows the first derivative of the fast oscillating part, that is,

A summary of the DI-EMD is given next by (1) set and ;(2) apply a number of sifting operations on in order to obtain an estimate of the first derivative of the total local mean ;(3) find the zero-crossings of the first derivative of the fast oscillating part ;(4) in order to estimate , perform a sifting iteration on using as interpolation nodes. the positions of the zero-crossing points estimated at Step (3). The characterization of the extrema as maxima or minima simply result from ;(5) set and return to Step (2) until the stoping criterion is fulfilled. As we are going to see in the simulations section, the estimation of the desired extrema described at Steps (2) and (3) is not necessary to be performed for each external sifting iteration.

4. Method Evaluation

In this section, DI-EMD will be tested in two different simulation scenarios. The first one is based on a two-pure tone signal model [4] and the second one contains a tone of linearly increased amplitude leading to a decomposition problem with gradually increased, with respect to time, difficulty.

4.1. Two-Tone Signal Example

The signal under consideration is given by(10)When takes values in , the term is the higher frequency component and is the lower frequency one. The aforementioned study tried to analytically answer the question: for which ranges of the parameters , the conventional EMD(i) is capable of separating the two signal components, that is, the first IMF equals to ;(ii) it considers the signal as a single component, that is, ;(iii) or is something else.

The latter behaviour can be considered as undesirable since the produced results are not directly related to the analyzed signal and its intrinsic properties. Interestingly, conventional EMD interprets the first extracted component neither as nor as for a quite wide range of and when [4].

In Figure 6, the gray area on the plane corresponds to parameters ranges where the EMD behaviour is undesirable since it does not produces IMFs directly related to the signal components under consideration. The metrics used for this graph are similar to the ones used in [4]. In the specific case, the of the area corresponds to undesirable decisions.

Figure 6: EMD decision areas. The white areas correspond to either equal to or and the gray area corresponds to a different decision.

Figure 7(a) exhibits a comparison between the ill-behaviour area of the conventional EMD and ill-behaviour area that corresponds to the 3rd-order DI-EMD when 30 internal iterations per external iteration are used. In both cases, the number of siftings is set equal to 10. The light gray area corresponds to the ill-behaviour of the conventional EMD, the black area corresponds to ill-behaviour of the 3rd-order DI-EMD and the midtone gray area corresponds to common ill-behaviour area between the two methods. Figure 7(b) compares the conventional EMD with the DI-EMD having 13th-order splines for the external iteration and 3rd-order splines for the internal. From Figure 7 is observed that the ill-behaviour area is reduced in some extent, especially when the high-order DI-EMD method is used. More specifically, using the DI-EMD configuration of Figure 7(a), the ill-behaviour area is reduced to the of the area and using the higher-order DI-EMD (Figure 7(b)) leads to undesirable area. In the two-tone example, the use of high-order splines in the internal iterations offered no further improvement. However, as will be seen in the next example, the latter observation cannot be taken as a general rule.

Figure 7: (a) The light gray area corresponds to ill-behaviour of the conventional EMD, the black area corresponds to ill-behaviour of 3rd-order DI-EMD, and the midtone gray area corresponds to common ill-behaviour area between the two methods. (b) The same as (a), but the 13 th-order splines have been used for the external sifting iterations.

4.2. Increased AM Signal Example

For the second performance evaluation, we adopted the signal shown in Figure 8. It consists of a constant frequency and constant amplitude sinusoid with frequency and amplitude and a sinusoid having a linearly increased amplitude and constant frequency . Two cases are explored. In the first one and in the second one . This simulation example explores the ability of several variants of EMD to extract the faster oscillating signal, that is, the signal , in a progressively “hostile” environment, namely, in the presence of a slow oscillating signal having a gradually increased amplitude. In such an example, the performance can be quantified by the value of the ratio up to which EMD succeeds in resolving the .

Figure 8: Multicomponent signal consisting of a constant amplitude signal and a signal with a linearly increased amplitude.

Starting with , Figure 9 shows in the logarithmic scale the difference between the fast oscillating signal and the IMF estimated after 2000 external sifting iterations using several different variants of EMD. In fact, the error curves undergo rapid oscillations which have been smoothed out here for visualization purposes. We observe that in general the EMD methods exhibits a gradually increased error as a function of as long as the problem of extracting becomes more and more difficult. Moreover, for each method there is a point at which the error is rising abruptly due to the fact that the fast signal can no longer be extracted properly.

Figure 9: Error between the fast oscillating signal and the IMF estimated after 2000 external sifting iterations.

We can consider a specific error value, say 0, that indicates the value up to which a specific EMD method succeeds in resolving . According to this performance measure, we can explore the ability of extracting the fast oscillating signal as a function of the external sifting iterations as it is shown in Figure 10. With respect to the figure legend, st-EMD() corresponds to the standard EMD using local extrema and spline interpolation of th order. Furthermore, the notation DI-EMD() corresponds to doubly-iterative EMD using spline interpolation of orders th and th for the external and the internal sifting processes with the internal sifting process being realized by a preset number of iterations. Moreover, there is the option for the internal sifting process to be performed not after every external sifting iteration, but every ones. The latter option corresponds to reduced complexity DI-EMD.

Figure 10: Performance results of different EMD variants for the case of .

Clearly, in such an example, the proposed doubly-iterative EMD outperforms the standard EMD. More specifically, the standard EMD needs more than 2000 sifting iterations in order to extract up to the point where the slow oscillating signal has amplitude 5 times larger than the amplitude of the fast oscillating signal. On the other hand, all the variants of DI-EMD exceed the value within 200 external sifting iterations. At the same time, the complexity remains low. For example, in the case of DI-EMD() the aforementioned performance is achieved with only a 40% increase in the total number of siftings. Moreover, the performance can be further improved with the use of higher-order splines in the internal iterations. However, it turns out that the higher the order of splines is, the larger the number of external siftings has to be in order to achieve the potential performance. Remember that in the case of the two tone signal the high-order splines in the internal iterations did not lead to performance improvements. The other way around is happening in the case of the increased AM example, namely, the performance deteriorates when high-order splines in the external iterations are used.When the frequency relation is reduced to , the advantages of using DI-EMD instead of the conventional EMD are reduced as it can be seen in Figure 11.

Figure 11: Performance results of different EMD variants for the case of .

In general, based on examples examined here, it can be argued that the use of DI-EMD is not “harmful,” however, the corresponding performance improvements depend on the signal to be analyzed. Moreover, high-order splines although, especially in the internal iterations can help, they do not guarantee enhanced performance compared to the cubic spline case and should be carefully used.

5. Conclusions

In this paper, the doubly-iterative EMD which incorporates an enhanced technique for interpolation points estimates for the EMD method was examined in two different simulation examples. An improvement was demonstrated in the overall decomposition performance which can in some cases enhanced when the doubly-iterative method is combined with envelope estimation using high-order spline interpolation.

Acknowledgments

This work has been presented in part by Y. Kopsinis and S. McLaughlin, “Enhanced Empirical Mode Decomposition Using a Novel Sifting-Based Interpolation Points Detection” at the IEEE Statistical Signal Processing Workshop, SSP 2007. This work was performed as part of the BIAS consortium under a grant funded by the EPSRC under their Basic Technology Programme.

References

1. N. E. Huang, Z. Shen, S. R. Long, et al., “The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis,” Proceedings of the Royal Society of London A, vol. 454, no. 1971, 903 pages, 1998.
2. L. Cohen, Time-Frequency Analysis, Prentice-Hall, Upper Saddle River, NJ, USA, 1995.
3. N. E. Huang and Z. Wu, “An adaptive data analysis method for nonlinear and nonstationary time series: the empirical mode decomposition and Hilbert spectral analysis,” in Proceedings of the 4th International Conference on Wavelet Analysis and Its Applications (WAA '05), Macao, China, November-December 2005.
4. G. Rilling and P. Flandrin, “One or two frequencies? The empirical mode decomposition answers,” IEEE Transactions on Signal Processing, vol. 56, no. 1, 85 pages, 2008.
5. R. Deering and J. F. Kaiser, “The use of a masking signal to improve empirical mode decomposition,” in Proceedings of the IEEE International Conference on Acoustics, Speech, and Signal Processing (ICASSP '05), vol. 4, p. 485, Philadelphia, Pa, USA, March 2005.
6. Y. Washizawa, T. Tanaka, D. P. Mandic, and A. Cichocki, “A flexible method for envelope estimation in empirical mode decomposition,” in Proceedings of the 10th International Conference on Knowledge-Based Intelligent Information and Engineering Systems (KES '06), vol. 4253, p. 1248, Bournemouth, UK, October 2006.
7. Y. Kopsinis and S. McLaughlin, “Investigation of the empirical mode decomposition based on genetic algorithm optimization schemes,” in Proceedings of the IEEE International Conference on Acoustics, Speech, and Signal Processing (ICASSP '07), vol. 3, p. 1397, Honolulu, Hawaii, USA, April 2007.
8. Y. Kopsinis and S. McLauglin, “Enhanced empirical mode decomposition using a novel sifting-based interpolation points detection,” in Proceedings of the14th IEEE/SP Workshop on Statistical Signal Processing (SSP '07), p. 725, Madison, Wis, USA, August 2007.
9. Y. Kopsinis and S. McLaughlin, “Investigation and performance enhancement of the empirical mode decomposition method based on a heuristic search optimization approach,” IEEE Transactions on Signal Processing, vol. 56, no. 1, 1 pages, 2008.
10. G. Rilling, P. Flandrin, and P. Gonçalvès, “On empirical mode decomposition and its algorithms,” in Proceedings of the 6th IEEE/EURASIP Workshop on Nonlinear Signal and Image Processing (NSIP '03), Grado, Italy, June 2003.
11. N. E. Huang, M.-L. C. Wu, S. R. Long, et al., “A confidence limit for the empirical mode decomposition and Hilbert spectral analysis,” Proceedings of the Royal Society A, vol. 459, no. 2037, 2317 pages, 2003.
12. G. Rilling and P. Flandrin, “On the influence of sampling on the empirical mode decomposition,” in Proceedings of the IEEE International Conference on Acoustics, Speech, and Signal Processing (ICASSP '06), vol. 3, p. 444, Toulouse, France, May 2006.