EMD as dyadic filter bank

Soon after EMD was proposed researchers started analysing it and testing it for different sorts of applications. And obviously for papers, but it follows the other. One of the first discovered behaviours was acting as a dyadic filter bank on random data. This was shown in 2004 both on fractal Gaussian noise [1] and on white noise [2]. Those are very interesting papers and definitely worth reading.

Quick example of such behaviour is presented in Figures 1 and 2. The first one shows input data – 256 data samples of normal random noise with 0 mean and variance 1 – and its EMD decomposition. Each row relates to different IMF. Figure 2 presents Fourier spectrum for signals in Figure 1, i.e. first row is spectrum for whole signal, second is for first IMF, third for second IMF and so on. Vertical lines relates middle (red) and end (black) of dyadic (powers of two) bands. Values of those lines are given as 256/n^2 and 256/(n+1)^2 respectively for black and red. Although they are not ideal they fit quite nicely and closely.

What is the reason for such behaviour? Actually it’s not understood well. It has been show empirically that it exists, but as many EMD related features it lacks of mathematical explanation (again, [1] and [2] discuss this very nicely). Definitely it has to do with the type of input data (EMD on signal with 2 modes will replicate those modes, see previous post) and… spline types. Figures 3 and 4 shows exactly he same features as Figures 1 and 2, but where produced using Akima spline technique instead of Cubic natural splines. It can be seen that the vertical lines are not as close matching in Fig. 4 as they are in Fig. 2.

1) EMD decomposition of normal random data. First row is the input signal and the following relate to respective IMF. Plots are not in the same y scale.


1) EMD decomposition of normal random data. First row is the input signal and the following relate to respective IMF. Plots are not in the same y scale.

2) Fourier spectrum made for original data and each IMF. Red and black vertical lines indicate respectively mean and end of filter bands. The exact values are powers of 2 - 256/(n+1)^2  for red and 256/n^2 for black, where n is the number of IMF. Plots are not in the same y scale.

2) Fourier spectrum made for original data and each IMF. Red and black vertical lines indicate respectively mean and end of filter bands. The exact values are powers of 2 – 256/(n+1)**2 for red and 256/n**2 for black, where n is the number of IMF. Plots are not in the same y scale.

3) EMD, using Akima's spline instead of natural cubic, performed on normal random data. First row is the input signal and the following relate to respective IMF. Plots are not in the same y scale.


3) EMD, using Akima’s spline instead of natural cubic, performed on normal random data. First row is the input signal and the following relate to respective IMF. Plots are not in the same y scale.

4) Fourier spectrum made for original data and each IMF from Fig. 3 (using Akima's splines). Red and black vertical lines indicate respectively mean and end of filter bands. The exact values are powers of 2 - 256/(n+1)**2  for red and 256/n**2 for black, where n is the number of IMF. Plots are not in the same y scale.


4) Fourier spectrum made for original data and each IMF from Fig. 3 (using Akima’s splines). Red and black vertical lines indicate respectively mean and end of filter bands. The exact values are powers of 2 – 256/(n+1)^2 for red and 256/n^2 for black, where n is the number of IMF. Plots are not in the same y scale.

[1] P. Flandrin, G. Rilling, and P. Goncalves, “Empirical Mode Decomposition as a Filter Bank,” IEEE Signal Process. Lett., vol. 11, no. 2, pp. 112–114, Feb. 2004.
[2] Z. Wu and N. E. Huang, “A study of the characteristics of white noise using the empirical mode decomposition method,” Proc. R. Soc. A Math. Phys. Eng. Sci., vol. 460, no. 2046, pp. 1597–1611, Jun. 2004.

Advertisements

Mode mixing in EMD

One of the problems researchers claim that empirical mode decomposition (EMD) has is the difficulty of unmixing mixed modes [1]. Usually it means that EMD cannot (always) untangle two or more sine mono frequency components s(t) = \sum_{n=0}^{N} A_n\sin(\omega_n t). Some say that this is big issue and propose tweaks to EMD [2], additional methods [3] or provide restrictions on when it is possible to extract [4]. The most widely analysed case is for N=2, i.e. s(t) = \sin(t) + b \sin(\omega t), where b and \omega are relative amplitude and frequency respectively. This has been studied thoroughly in [4].

Why is this not necessarily an issue? EMD is supposed to extract intrinsic mode functions (IMFs), i.e. functions with a single mode/frequency. Although, when hearing oscillations many probably think about sinusoid, but it actually can be in any form. Wavelets are an example of weirdly looking oscillations (small waves that have beginning and end). As for the definition of oscillatory behaviour there is not a single best one. Huang et al. [5] defined IMF as a function that:

  1. has number of local extrema and zero-crossings equal or different by 1,
  2. mean of its top and bottom envelopes is zero.

Authors were aiming to meta-define a function that is suitable for Hilbert transform, by which they have defined \omega(t) instantaneous frequency. Components of decomposition, i.e. IMFs, are thus described in formula c_j (t) = a(t) cos( \phi(t) ), where a(t)>0 is an amplitude and \phi(t) = \int \omega(t) dt is a phase. Representation of a function by multiplication of two other functions is ambiguous, same as number 5 can be presented as 5 = 2 \cdot 2.5 or 5 = 0.2 \cdot 25.
Here is where the Hilbert transformation plays a role. It provide a recipe for representing signal in amplitude times phase functions, but it is still one of the methods. Nevertheless, once we are allowing signal to have modulation in amplitude mode mixing makes little sense. Taking two sinusoids as an example, they can be either in written in form of s_{I} (t) = \sin( \omega_1 t) + \sin( \omega_2 t) or, using trigonometric properties, as s_{II} (t) = 2 \cos(\frac{\omega_1 - \omega_2}{2} t ) \sin(\frac{\omega_1 + \omega_2}{2} t) . In the first form we have two sines with constant amplitudes, whereas in the second one there is a single sine function, of frequency 0.5(\omega_1+\omega_2), which amplitude is modulated by 2\cos(\frac{\omega_1 - \omega_2}{2} t ). Both forms give exactly the same result – the difference is in written representation. It is up to user to decide how he sees it.

A numerical example. Given two sines EMD is behaving choosy on parameters such as a relative amplitude and a relative frequency. Figure 1 presents EMD decomposition on signal S_1 = \cos(15\pi t ) + \cos(39\pi t). It can be seen that two extracted IMFs have exactly the frequencies (number of extrema divided by 2) of which the signal has been synthesised. This is due to the big difference between used frequencies. If one uses close values, i.e. allowing for beating, then the most information is contained in the first IMF (Figure 2). To produce this Figure sines with frequencies 17.5 Hz and 19.5 Hz (S_2 = \cos(35\pi t) + \cos(39\pi t)) were used. Second row of Figure 2 shows the difference between input signal and first IMF – barely visible on the same amplitude scale.

Fig. 1. Input signal for EMD consists of two sines with frequency 7.5 and 19.5 Hz. Those modes have been extracted in IMFs - first contains 19.5 Hz and second 7.5 Hz.

Fig. 1. Input signal for EMD consists of two sines with frequency 7.5 and 19.5 Hz. Those modes have been extracted in IMFs – first contains 19.5 Hz and second 7.5 Hz.

Fig. 2. EMD performed on signal made of two sines with 17.5 and 19.5 Hz frequencies. In first row original signal (blue) and first IMF (green) are overlapping. Second row shows difference between them.

Fig. 2. EMD performed on signal made of two sines with 17.5 and 19.5 Hz frequencies. In first row original signal (blue) and first IMF (green) are overlapping. Second row shows difference between them.

References:
[1] N. E. Huang and Z. Wu, “A review on Hilbert-Huang transform: Method and its applications to geophysical studies,” Rev. Geophys., vol. 46, no. 2, p. RG2006, Jun. 2008.
[2] Z. Wu and N. E. Huang, “Ensemble Empirical Mode Decomposition: A Noise-assisted Data Analysis Method’,” Adv. Adapt. Data Anal., vol. 01, no. 01, pp. 1–41, Jan. 2009.
[3] R. Deering and J. J. F. Kaiser, “The Use of a Masking Signal to Improve Empirical Mode Decomposition,” in Proceedings. (ICASSP ’05). IEEE International Conference on Acoustics, Speech, and Signal Processing, 2005., 2005, vol. 4, pp. 485–488.
[4] G. Rilling and P. Flandrin, “One or Two Frequencies? The Empirical Mode Decomposition Answers,” IEEE Trans. Signal Process., vol. 56, no. 1, pp. 85–95, Jan. 2008.
[5] N. E. Huang, Z. Shen, S. R. Long, M. C. Wu, H. H. Shih, Q. Zheng, C. C. Tung, H. H. Liu, and N.-C. Yen, “The Empirical Mode Decomposition and the Hilbert Spectrum for Nonlinear and Non-stationary Time Series Analysis,” Proc. R. Soc. A Math. Phys. Eng. Sci., vol. 454, no. 1971, pp. 903–995, Mar. 1998.

Empirical mode decomposition – Introduction

Aside

Empirical mode decomposition (EMD) is a data-driven decomposition method and was originally proposed by Huang et. al in 1998 [1]. Since that time the method has gained a lot of attention in the science community. EMD has been applied in a wide range of different fields, including geophysics, biomedicine, neuroscience, finance and many more.

The method is defined by an algorithm as follows:

  1. Identify all local extrema (both minima and maxima) in input signal s(t).
  2. If the number of extrema is less or equal 2 then s(t) is considered as a trend (r(t) := s(t) ) — finish with one component.
  3. Estimate top (env_max) and bottom(env_min) envelopes by interpolating respectively local maxima and local minima with natural cubic splines.
  4. Calculate local mean (mean of both envelopes) — m(t) = 0.5(env_{max} + env_{min}).
  5. Subtract the mean from the input signal h(t) = s(t) - m(t) .
  6. If h_j fulfills the stopping criteria, then it is considered an intrinsic mode function (IMF) (a component c_j(t) ) and algorithm starts again from step 1 for a signal s(t) := m(t). Otherwise, it starts with s(t) := h(t).

Empirical evidence is that the algorithm converges to finite number of IMFs, with which input signal can be reconstructed :

s(t) = \sum_n^N c_i (t) + r(t)

where N is the number of components.

Image

Top (blue) and bottomo (red) envelopes of signal (green). Their average, described as carrier wave, is presented in dashed black.

Image

Empirical mode decomposition on signal composed of three sine modes.

[1] N. E. Huang, Z. Shen, 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,” Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, vol. 454, pp. 903-995, 1998.

This blog.

It didn’t take long to forget about this blog. Now, hopefully, with revoke of motivation to post more periodically I would like to reactivate it. Content should be related to work I’m doing currently and will try to post at least once a week, probably on Tuesdays.
Yes, this is more a statement for myself. A reminder for future me!
Come on lad!