Hilbert-Huang Transform

In the recent times HHT is gaining huge popularity for spectral analysis of nonlinear and nonstationary time series data (Huang et al. 2016). It is a purely data adaptive method, which can produce physically meaningful representations of time- series data. This method does not require a priori selection of functions, but instead it decomposes the signal into intrinsic oscillation modes derived from the succession of extrema. The HHT involves two major steps (i) the use of EMD method or its variants such as EEMD (Wu and Huang 2005) or CEEMDAN (Torres et al. 2011) to decompose a time series into a collection of orthogonal time series, namely IMFs; (ii) the use of HSA to obtain instantaneous frequency, which may be helpful to identify embedded structures of time series data. This decomposition helps to better understand the internal structure of the signal and the components involved (Klionski et al.

2008). The original data can be reconstructed by the linear addition of each of the IMFs, which are independent to each other.

Each IMF must meet the following conditions:

(a) the number of extreme values in the overall data must match the number of zero crossings, or differ by only one

where N = total number of maxima; N . = total number of minima; and

max 7 mm 7

N = total number of zero crossings; and

(b) at any point of time, the mean value of the envelope of local maxima and envelope of local minima must be zero

where £raax(f) is the envelope of the local maxima; and £roin(0 is the envelope of the local minima, which were obtained by spline interpolation procedures.

Empirical Mode Decomposition

The process of extracting the IMFs from a time series X(t) (which is known as the sifting process) consists of the following steps:

  • 1. Identify all extrema (maxima and minima) of the signal X(t)
  • 2. Connect these maxima with any interpolation function (e.g., cubic spline) to construct an upper envelope (£max(0); use the same procedure for minima to construct a lower envelope (£m|n (t))
  • 3. Compute the mean of the upper and lower envelope, m(t)
  • 4. Calculate the difference time series d(t) = X(t)-m(t)
  • 5. Let d(t) be the new signal and repeat steps (1) to (4) until d(t) becomes a zero- mean series with no riding waves (i.e., there are no negative local maxima and positive local minima) with smoothened amplitudes.

To satisfy step (5), an appropriate criterion is to be applied to stop the sifting iterations in order to guarantee that the IMF retains enough physical sense of both amplitude and frequency modulations (Huang and Wu 2008). A number of stopping criteria have been reported in the literature (Huang et al. 1998; Huang and Wu 2008). One popular criteria is the modified Cauchy type stopping criterion (Huang and Wu 2008) computed from two consecutive sifting results in:

where [1] [2]k’ is the index for IMF; T is the index for iteration of the sifting operation (to get k"‘ IMF); T is the data length; dkj(t) is the deviation of the original time series from the mean in the iteration to evolve the k"' IMF; is the tolerance value specified by the user (normally, 0.2-0.3 as per Huang and Wu 2008).

  • 6. On satisfying the zero-mean condition, d(t) can be designated as the first intrinsic mode function IMF1.
  • 7. Compute the residue R{(t) by subtracting IMF1 from original signal (i.e., Rt(t)=X(t)-IMF((r)), is used as new signal. The sifting process is repeated upon this signal to get IMF2.
  • 8. The higher oscillatory modes are obtained by treating the residue (Rk(t)) as the signal (X (0), iteratively.


The k"‘ residue is defined as Rk (t) = X(t)~^ IMFj(t). The process will be continued


till the resulting residue is a monotonic function or a function having only one extrema.


Then the original signal can be reconstructed as X(t) = ^[lMFk(t)] + RK(t), where

k= I

К is the number of decomposed IMFs.

During the decomposition of a signal into IMFs, sometimes it may fail to assign the signals with similar frequencies into separate IMFs. As a result, each IMF contains different modes of oscillations, which make IMFs to lose physical meaning or falsely represent the physical processes associated with mode. To solve this problem, noise assisted data analysis methods were proposed by different researchers in the past. The Ensemble EMD (EEMD) (Wu and Huang 2005) and the Complete Ensemble EMD with Adaptive Noise (CEEMDAN) (Torres et al. 2011) falls in this category.

Ensemble Empirical Mode Decomposition

The procedure of EEMD is presented below (Wu and Huang 2005):

  • (iii) If m < M, go to step (i) with m = m + 1. Repeat Steps (i) and (ii) multiple times, with different white noise series
  • 3. Obtain an ensemble mean lMFk (?)of the corresponding IMFs; report the mean MFk(t) as the final IMF.

The effect of the added white noise should decrease following statistical rules (Wu and Huang 2005):

where t// is the amplitude of the added noise and i//y is the final standard deviation of error, which is defined as the difference between the input signal and the corresponding IMFs. Equations (2.27) and (2.28) indicate that as the ensemble number increases, the effect of added white noise gets nullified. In short, k"' true IMF (MFk(t)) is the mean of the corresponding IMFs obtained by EMD over ensemble of artificial signals, generated by adding different realizations of white noise (wu(/)) to the original signal X(t).

Complete Ensemble Empirical Mode Decomposition with Adaptive Noise (CEEMDAN)

In CEEMDAN method, a noise is added at each stage of the decomposition to result in a unique residue in each mode from the residue of previous mode (or the true signal, for the first mode) and currently generated IMF.

1. Perform the EMD for Mrealizations Xm(t) = X(t) + e0wm(t) and compute the first mode of CEEMDAN IMF{(t), by averaging the realizations

where m=I,2,...,M the index for realizations; e0 is the noise parameter for the initial step.

This step shows that the first mode obtained by CEEMDAN will be same as that obtained by EEMD.

  • 2. At the first stage = 1) calculate the first residue as R{(t) = X(t)~IMFx(t)
  • 3. Decompose the realizations Ru„(t) = R{(t) + elEl(wm(t)), m =

until their first EMD mode gets evolved. Then compute the second mode


where the operator Ek(.) is an operator represents the evolution of the k'h mode by EMD; £,is the noise parameter for the first stage (k = 1)

4. Calculate the k"' residue as

where IMFk (t) is the IMFs obtained by CEEMDAN

5. Compute the first mode of Rk(t)+ ekEk(wm(t)), and define the (k + l)lh mode by CEEMDAN as

6. Go to step (4) for next к

Steps 4 to 6 are performed until the obtained residue is no longer feasible to be decomposed (i.e., the residue is monotonic or having only one extrema).

к "

The final residue becomes RK(t) = X(t)-£^IMFk(t). Therefore, the above


decomposition process is complete and provides an exact reconstruction of the original data, i.e.,

Most of the geophysical processes are influenced by multiple potential causal variables and it is often necessary to investigate the association of processes with such input (causal) variables. Over the years, researchers developed different variants of EMD (or its noise-assisted versions), which can handle multiple time series simultaneously. For example, Rilling et al., (2007) proposed bidimensional EMD, Rehman and Mandic (2010a) proposed tridimensional EMD, etc. The Multivariate EMD (MEMD) proposed by Rehman and Mandic (2010b) is the most generalized version of EMD to decompose multiple time series signals simultaneously.

Multivariate Empirical Mode Decomposition

Multivariate EMD proposed by Rehman and Mandic (2010b) is an extension of the traditional EMD, which decomposes multiple time series simultaneously after identifying the common scales inherent in different time series of concern. In this method, multiple envelops are produced by taking projections of multiple inputs along different directions in an m-dimensional space.

Assuming V(t) = {v, (?), v2(?).... vm (?)}being the m vectors as a function of time ?,

and X% ={Xlk,x2k.....,xj} denoting the direction vector along different directions

given by angles (pk = | in a direction set X (k = 1,2,3, ... К, К is the

total number of directions). It can be noted that the rotational modes appear as the counterparts of the oscillatory modes in EMD or its variants. The IMFs of m temporal datasets can be obtained by the following steps:

  • 1. Generate a suitable set of direction vectors by sampling on a (m - 1) unit hypersphere
  • 2. Calculate the projection p4** (?) of the datasets V(t) along the direction vector

for all к

  • 3. Find temporal instants ?,% corresponding to the maxima of projection for all к
  • 4. Interpolate [?, V(?,%)] to obtain multivariate envelop curves eft(0 for all к
  • 1 чг
  • 5. The mean of envelope curves (M(t)) is calculated by M(t) = — V (?)

^ k=I

6. Extract the detail D(t) using D(t) = V(t) - M(t). If D(?) fulfills the stoppage criterion for a multivariate IMF, apply the above procedure from step (1) onward upon the residue series (i.e., V(t) - D(t)). Otherwise, repeat the steps from (2) onward upon the series £)(?).

Hammersley sampling sequence can be used for the generation of direction vectors (Huang et al. 2016) and the stoppage criteria proposed for EMD (i.e., the sum squared difference between the deviations of the mode from the mean signal in two consecutive iterations should be less than a specified tolerance) can be used in the implementation of MEMD.

Multivariate dataset could be a collection of totally different signals and MEMD is a technique, which can capture the common scales present in those signals. The extension of EMD to a multivariate dataset involves the following steps: (i) performs multiple real-valued projections of the signals along multiple directions on hyperspheres (n-spheres); (ii) multidimensional envelopes of the signal are created by interpolating (component wise) the extrema of projected signals; (iii) averaging of the envelops to obtain the local mean. In the implementation of the above steps, MEMD uses the mathematical concept of quaternion. The quaternions, first described by Irish mathematician Sir William Rowan Hamilton in 1843, are a number system that extends the complex numbers to higher dimensions.

Any number of the form q = a : + bi + cj + dk, where a, b, c, and d are real numbers, i2 =/ = k2 = -1, and ij = k= -ji, jk = i and ki = /; is a quaternion.

Under addition and multiplication, quaternions have all the properties of a field, except multiplication is not commutative. _

The norm of quaternions is defined as ||r/|| = +/; +c +d and the unit vector

is e"e = cos 0+и sin в.

The conceptual idea of projection of direction vectors and sampling are provided in Figures 2.9 and 2.10.

To generate a suitable set of direction vectors, unit hyperspheres can be sampled based on uniform angular sampling and quasi-Monte Carlo-based low-discrepancy

Concept of projection of a point along different directions in a 3D space

FIGURE 2.9 Concept of projection of a point along different directions in a 3D space: (a) the direction vector OA in 3D space, can also be represented by a point on the surface of a unit sphere: (b) multiple direction vectors represented by points on a particular longitude line.

Concept of projection of a signal along multiple directions by rotation

FIGURE 2.10 Concept of projection of a signal along multiple directions by rotation (a) For projections along longitudinal lines on a sphere, multiple axes represented by a set of vectors are chosen in thex-y plane, with angle taken with respect to x-axis; (b) projections of the input signal can be taken by rotating the input signal along rotation axes represented by a set of unit quaternions.

sequences. A simple and practically convenient choice for a set of direction vectors is to employ uniform angular sampling of a unit sphere in an «-dimensional hyperspherical coordinate system. The resulting set of direction vectors covers the whole (n - 1) sphere. A coordinate system in an «-dimensional Euclidean space can then be defined to serve as a point set (and the corresponding set of direction vectors) on an (n - /) sphere. Let {£?,, 6f,..., 0(nl)} be the (« - 1) angular coordinates, then an «-dimensional coordinate system having {x, }"=| as the « coordinates on a unit (« - 1) sphere is given by:

The point set corresponding to the «-dimensional coordinate system is now very convenient to generate. There are many ways to generate such sequence. In certain applications the quasi-Monte Carlo method and their lower discrepancy (close to uniform sampling) version is used. Rehman and Mandic (2010b) used Hammersley- Halton sequence (involving prime numbers, quasi-Monte Carlo method) to generate the direction vectors. Multiple direction vectors in «-dimensional space is depicted in Figure 2.11.

Direction vectors for taking projections of tri-variate signals on a two-sphere generated by using (a) spherical coordinate system and (b) a low-discrepancy Hammersley sequence

FIGURE 2.11 Direction vectors for taking projections of tri-variate signals on a two-sphere generated by using (a) spherical coordinate system and (b) a low-discrepancy Hammersley sequence.

  • [1] Initialize the number of ensemble realizations M, the amplitude of the addedwhite noise for the mlh realization (wm), and set in = 1
  • [2] Perform m,h trial on the signal (i) Add white noise series with the given amplitude to the investigatedsignal where XJi) indicates the m'h artificial signal; X(t) is the true signal; andwm(t) represents the zero mean unit variance white noise signal for the m'htrial. (ii) Decompose the noise-added signal XJf) into К number of IMFs IMFkm(k= 1,2K) using the EMD method described in the previous section,where IMF. = k'h IMF of the mlh trial KJ)1
< Prev   CONTENTS   Source   Next >