Signal feature extraction based on an improved EMD method

Measurement 42 (2009) 796–803

Contents lists available at ScienceDirect

Measurement journal homepage: www.elsevier.com/locate/measurement

Technical note

Signal feature extraction based on an improved EMD method Li Lin *, Ji Hongbing School of Electronic Engineering, Xidian University, Xi’an 710071, PR China

a r t i c l e

i n f o

Article history: Received 25 August 2008 Received in revised form 20 December 2008 Accepted 5 January 2009 Available online 10 January 2009

Keywords: Feature extraction Empirical mode decomposition Inverse EMD filter Sifting stop criterion Optimal envelopes mean

a b s t r a c t Empirical mode decomposition (EMD) is a powerful tool for analyzing composite, nonlinear and non-stationary signals. Intrinsic mode functions (IMFs) obtained by EMD are bandlimited, which can represent the features of signal and reserve the local information. Because of lacking theoretical foundation, EMD has many problems in scale mixing, stop criterion, etc. In this paper, we use an improved EMD method for signal feature extraction. Optimal envelopes mean is obtained by an inverse EMD filter scheme in the improved method. A new sifting stop criterion is proposed to guarantee the orthogonality of the sifting results. And finally, two different methods are utilized for performance evaluation of the decomposition results. Numerical simulation and experimental result demonstrate the validity of the improved method. Ó 2009 Elsevier Ltd. All rights reserved.

1. Introduction The analysis of composite, nonlinear and non-stationary signals has become a hot spot of research in the field of signal measurement and feature extraction recently. A number of new methods have been proposed to analyze these signals. One of the promising methods is the empirical mode decomposition (EMD). EMD was proposed by Huang [1], which is the core of Hilbert–Huang transform (HHT). Since 1998, EMD has been widely researched and proved effective for many actual signals, such as earthquake [1], vibration [2], flow [3], etc. The Hilbert–Huang transformdata processing system (HHT-DPS) has been used and commercialized successfully [4]. However, EMD may generate undesirable IMFs. The undesirable IMFs will cause misinterpretation to the result. The IMFs may cover too wide a frequency range,

Abbreviations: EMD, empirical mode decomposition; IMF, intrinsic mode function; HHT-DPS, Hilbert–Huang transform-data processing system; GA, genetic algorithm; SD, standard deviation; AM–FM, amplitude modulation–frequency modulation; IF, instantaneous frequency; OC, orthogonality criterion; MSE, mean square error; ISM, index of scale mixing; DFT, discrete Fourier transform; EEG, electroencephalogram. * Corresponding author. Tel.: +86 13609183196. E-mail address: [email protected] (L. Lin). 0263-2241/$ - see front matter Ó 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.measurement.2009.01.001

viz. scale mixing. The information of extrema is very important in EMD, and many works focus on the acquirement of efficient extrema. A novel version of EMD is proposed in [5] based on the selection of interpolation points, called doubly-iterative EMD. Kopsinis [6] shows that there are specific extrema, which are related to the analysis signal and are able to lead to much improved decomposition performance. It uses a genetic algorithm (GA) to estimate the specific interpolation points. Furthermore, sifting stop criteria, i.e. the stop conditions are of great importance. First Huang [1] proposes the standard deviation (SD) criterion

SD ¼

T X jhi1 ðtÞ  hi ðtÞj2 2

t¼0

hi1 ðtÞ

:

ð1Þ

where h(t) is the detail information in the sifting process. The three-threshold criterion [7] was proposed as an improvement to SD criterion. It utilizes three thresholds h1, h2 and a, and defines r(t) = jm(t)/a(t)j, m(t) = (ex(t) + em(t))/2, a(t) = (ex(t)  em(t))/2. ex(t) and em(t) are the upper envelope and lower envelope, respectively. EMD sifts until r(t) < h1 for fraction (1  a) and r(t) < h2 for remaining fraction. In this paper, an improved EMD method is proposed for measurement and analysis of composite or multicom-

797

L. Lin, J. Hongbing / Measurement 42 (2009) 796–803

ponent signals. With this new method, composite signals can be decomposed to several IMFs. The IMFs represent the intrinsic and reality feature of the analyzed signal. It enhances the orthogonality of the sifting results and reduces the scale mixing effectively. The improved method is very fit for detection of transient signal and fault signal, extraction of the intrinsic feature and prediction of interested information. It has a broad application potential in fileds of mechanical and vibration, flow, speech, EEG, radar, earthquake, tide, ocean, etc. The remainder of this paper is organized as follows. Section 2 presents the standard EMD algorithm. Section 3 discusses the improved EMD method. Section 4 is the experiments and results. And finally the conclusions are drawn in Section 5.

sðtÞ ¼ mK ðtÞ þ

K X

First, Huang [1] defined the IMF which should satisfy two conditions: (1) in the whole data set, the number of extrema and the number of zero crossings must either equal or differ at most by one; and (2) at any point, the mean value of the upper envelope and lower envelope is zero. Then the sifting process is: (a) For a given signal s(t), construct its upper envelope and lower envelope by connecting all local maxima and local minima with cubic spline functions [8], respectively. (b) Compute the envelopes mean m(t) and extract the detail h(t) = s(t)  m(t). Regard h(t) as new s(t) and repeat the operation above until h(t) satisfies the IMF conditions, then obtain the first-order IMF imf1(t) = h(t). (c) Let the residual r(t) = s(t)  h(t) be a new signal, repeat (a), (b) and (c), obtain the other orders IMFs. Compared with Fourier and Wavelet analysis, EMD generates a set of adaptive and data-driven basis. The basis functions are a series of amplitude modulation–frequency modulation (AM–FM) monocomponent signals. Fourier is the simplest model for stationary signals and linear systems based on sinusoid. Wavelet decomposition is conventionally achieved by repeating application of two filtering operations, high frequency detail superimposed on lowfrequency components. It splits the signal based on predetermined spectral basis by the use of linear time-invariant filters and so precludes the possibility of adapting to local variations in the signals. EMD works like an adaptive high pass filter. It sifts out the fastest changing component of a composite signal first [9]. The cutoff frequency of the high pass filter is adaptive and data-driven. After the first filtering, the signal s(t) can be expressed as

sðtÞ ¼ m1 ðtÞ þ imf1 ðtÞ:

ð2Þ

imf1(t) is the first-order IMF. It is the high frequency component of the signal. Compared with imf1(t), m1(t) is the low frequency component. Repeating the decomposition on m1(t), then we obtain

ð3Þ

The total number of the IMF is K. It should be noted that in the whole sifting process, imfk(t) and mk1(t) have the same frequency feature and extrema number. 3. The improved EMD method 3.1. The optimal envelopes mean First we consider a composite signal s(t) (t = 1, 2, . . . , N) which consists of K AM–FM monocomponent signals

sðtÞ ¼

K X

si ðtÞ ¼

i¼1

2. EMD algorithm

imfk ðtÞ:

k¼1

K X

ai ðtÞ cosðui ðtÞÞ:

ð4Þ

i¼1

The instantaneous frequency (IF) of ith component is fi(t), fi ðtÞ ¼ 21p d/dti ðtÞ and fi(t) > fj(t) for any i < j at all time point t. In the sifting process, the optimal envelopes mean (mk(t))opt can be obtained by the sum of the signals with the lower IF.

ðmk ðtÞÞopt ¼

K X

si ðtÞ:

ð5Þ

i¼kþ1

Unfortunately, because the optimal envelopes mean consists of multicomponent signals, it cannot be obtained exactly by the cubic spline interpolation in the iterative process. To decrease the estimation error of the optimal mean, we use the higher order extrema. Then the optimal mean can be expressed as,

ðmk ðtÞÞopt ¼

n X

ci ðtÞ

ð6Þ

i¼1

where, ci (t),i = 1, 2, . . . , n is the cubic spline interpolation curve to the ith-order extrema of the residual signal PK i¼kþ1 si ðtÞ. This process works like an inverse filter scheme compared with the standard EMD. 3.2. Higher order extrema and inverse EMD filter The first-order maxima u1(tk) and minima l1(tk) of the discrete signal s(t) are the same as the maxima and minima in standard EMD, defined as

(

2

2

u1 ðt k Þ ¼ ftk : ðdsðtÞ=dtÞ ¼ 0; ðd sðtÞ=dt Þ < 0g 2

2

l1 ðtk Þ ¼ ft k : ðdsðtÞ=dtÞ ¼ 0; ðd sðtÞ=dt Þ > 0g

:

ð7Þ

The second-order maxima u2(tk) and minima l2(tk) of s(t) are the maxima of u1(tk) and minima of l1(tk), respectively. Consecutively, the jth-order maxima {uj(tk)} and minima {lj(tk)} are defined as maxima and minima of the (j  1)th’s, respectively. Obviously, the maxima sets and minima sets with different orders satisfy U1  U2      Un and L1  L2      Ln, respectively. n is the maximum order. Besides {uj(tk)} and {lj(tk)}, we also define the jth-order max minima {ulj(tk)} and min maxima {luj(tk)}, which are the maxima of {lj1(tk)} and the minima of {uj1(tk)}, respectively. Usually, n is larger

798

L. Lin, J. Hongbing / Measurement 42 (2009) 796–803

than (K  1). Otherwise, not all the K signal components can be correctly obtained. The process of the inverse EMD filter is: (a) Determine the maximal extrema order n. (b) Find the nth-order extrema {un(tk)}, {ln (tk)}, {uln(tk)} and {lun(tk)}, respectively. (c) Construct the nth-order envelopes un(t), ln(t), uln(t) and lun(t) by cubic spline interpolation, respectively. Make sure that all the extrema locate at the points between un(t) and ln(t). (d) Compute the three different means m1(t) = (un(t) + ln(t))/2, m2(t) = (uln(t) + lun(t))/2 and m3(t) = (un(t) + ln(t) + uln(t) + lun(t))/4. Then judge whether the envelop means satisfy a given condition. If yes, choose the best envelop mean m(t) from the three means, and go to Step (e). Otherwise, jump (f). (e) Let c1(t) = c1(t) + m(t) and s(t) = s(t)  m(t), repeat the above operation until m(t) dissatisfies the given condition. Then obtain the lowest frequency IMF c1(t). (f) Let n = n  1 and repeat Steps (a) to (f), obtain the other orders IMFs until n = 1. The highest frequency IMF is the residual. Notes that when n 6 2, m1(t) = m2(t) = m3(t). Compared with the standard EMD sifting process, the inverse filter scheme is a low pass filter and the lowest frequency component is sifted out first. The cutoff frequency of the low pass filter is also adaptive and data-driven.

3.4. The improved EMD sifting scheme The standard EMD method usually brings redundant low frequency IMFs with very small amplitudes. They are useless for feature extraction and analysis of the original signal. Actually, the components number of signal is directly related to the extrema order, without regarding the influence of sampling interval and noise. And theoretically, the maximum order n is equal to (K  1). The inverse filter scheme of EMD has good decomposition performance for low frequency components. Compared with the IMF to be sifted out, the optimal envelopes mean is also the low frequency component. So, the improved EMD in this paper is: (a) For a given signal s(t), first the residual r(t) is obtained by the inverse EMD filtering. It is the first-order IMF with highest frequency. (b) Let s(t) = s(t)  r(t), repeat Step (a) and Step (b) to obtain other orders IMFs of the signal until residual signal s(t) satisfies the given conditions. The difference between the standard EMD and the improved EMD methods lies in the fitting of the optimal envelopes mean (mk(t))opt in Eq. (5). Furthermore, the orthogonality criterion in Eq. (9) can guarantee the orthogonality of the sifting results that is not considered in the standard EMD sifting process.

3.3. The sifting stop criterion Many stop criteria are based on the local symmetry of the IMF to be sifted out, such as the Cauchy-type criterion [1] and the three-threshold criterion [7]. In theory, an IMF should satisfy N X

imf ðtÞðsðtÞ  imf ðtÞÞ ¼ 0:

ð8Þ

t¼1

which is called orthogonality. For the signal model in Eq. (4), we propose the orthogonality criterion (OC) as the given condition in Step (d), which is defined as

   N X mðtÞ  sðtÞ   OC ¼  :  t¼1 mðtÞ  ðsðtÞ  mðtÞÞ

ð9Þ

If m(t) satisfies m(t) = b  si(t), 0 < b < 1, then OC = b/ (b  b2) = 1/(1  b) > 1. If m(t) is an ‘‘artifact” and irrelative P to s(t), then mðtÞ  sðtÞ ¼ 0, OC < 1. In Step (d), we choose the envelop mean m(t) from {mi(t), i = 1, 2, 3} with the largest OC. Considering the influence of sampling interval and data length, we set OC > 1.05 as the sifting condition in the inverse EMD filter. This means that m(t) can be regarded as an IMF with the amplitude ratio b P 1/21. This condition guarantees the orthogonality of the sifting result.

Fig. 1. (a) Comparisons between the standard EMD and the improved EMD with MSE criterion. (b) Comparisons between the standard EMD and the improved EMD with ISM criterion.

799

L. Lin, J. Hongbing / Measurement 42 (2009) 796–803

MSE ¼ maxfMSEi ;

3.5. Performance evaluation Rilling [10] discusses the conditions of frequency and amplitude when a composite signal can be decomposed completely. The standard mean square error (MSE) criterion is described as

P MSEi ¼

2 t jimfi ðtÞ  si ðtÞj : P 2 t jsi ðtÞj

ð10Þ

If MSEi > 1 then simply let MSEi = 1. And good decomposition is obtained when MSEi approaches to 0. Considering every component estimation error, we define the total standard mean square error

a

The MSE criterion denotes the error between the decomposition results and the original signal components in time domain. The MSE criterion in the time–frequency domain can be represented as,

P FMSEi ¼

t jhfi ðtÞ

P

 fi ðtÞj2

t jfi ðtÞj

2

ð12Þ

:

Where, hfi(t) is the instantaneous frequency of imfi(t) by Hilbert transform. Usually, Hilbert spectrum tends to generate undesired vibration that will influence the evaluation. For single-fre-

s1

0

s2

2

s3

2

s(t) IMF 1 IMF 2 IMF 3 Residual

IMF 1 IMF 2 IMF 3

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0

0.1

0.2

0.3

0.4

0.5 Time (s)

0.6

0.7

0.8

0.9

1

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0

0.1

0.2

0.3

0.4

0.5 Time (s)

0.6

0.7

0.8

0.9

1

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0

0.1

0.2

0.3

0.4

0.5 Time (s)

0.6

0.7

0.8

0.9

1

0 -2

Residual

0

0 -2

c

ð11Þ

2 -2

b

i ¼ 1; 2; . . . ; ng:

2 0 -2

2 0 -2 2 0 -2 2 0 -2 2 0 -2

2 0 -2 2 0 -2 2 0 -2 2 0 -2

Fig. 2. (a) The original signal. (b) Result by the standard EMD method. (c) Result by the improved EMD method.

800

L. Lin, J. Hongbing / Measurement 42 (2009) 796–803

Amplitude

1 0 -1

0

0.02

0.04

0.06

0.08 Time (s)

0.1

0.12

0.14

0.16

Fig. 3. Sound vibration signal.

quency sinusoidal signals, we define the index of scale mixing (ISM) to measure the scale mixing of an IMF.

Sðfi Þ ISMi ¼ 1  P f Sðf Þ

ð13Þ

where, S() is the frequency spectrum of imfi by DFT and fi is the center frequency of si(t). Then the total scale mixing (ISM) is

IMF1

0.5 0 -0.5

IMF2

0.5 0 -0.5

IMF3

a

0.5 0 -0.5

IMF4

ð14Þ

0.5 0 -0.5

ð15Þ

i1

IMF5

i ¼ 1; 2; . . . ; ng:

0.5 0 -0.5

The performance of EMD is only sensitive to the relative parameters of the components, thus leading to the simpler representation,

Residual

ISM ¼ maxfISMi ;

decomposition performance is very poor when u1 > 0.7, whatever the value of u2. The decomposition results show that the MSE of the improved EMD method is a little smaller than that of the standard EMD method, but the ISM of the improved EMD method is distinctly smaller than that of the standard EMD method. For non-stationary signal in Fig. 2a, the decomposition results are shown in Fig. 2b and c. The sampling frequency

0.5 0 -0.5

Obviously, ISM 2 [0, 1]. For discrete signals, the point number of DFT should satisfy that the center frequency locates at the integer points. 4. Experiments and results 4.1. Simulation signals First we consider a simple signal consisting of three sinusoidal components

sðtÞ ¼ s1 ðtÞ þ s2 ðtÞ þ s3 ðtÞ ¼

3 X

ai cosð2pfi t þ /i Þ:

0

0.05

0.1

0.15

0

0.05

0.1

0.15

0

0.05

0.1

0.15

0

0.05

0.1

0.15

0

0.05

0.1

0.15

0

0.05

0.1

0.15

sðtÞ ¼ cosð2pf1 tÞ þ b1 cosð2pu1 f1 t þ h1 Þ þ b2 IMF1 IMF2 IMF3

where, b1 = a2/a1, b2 = a3/a2, u1 = f2/f1, u2 = f3/f2, h1 = /2  /1 and h2 = /3  /2. Due to symmetry, we just consider the situation when 0 < u1 6 1, 0 < u2 6 1 and u1 P u2, hence s1(t) is the highest frequency component and s3(t) the lowest frequency component in the composite signal. To analyze the decomposition performance under different frequency conditions, we set u1 = 0.1, 0.2, . . . , 1 and u2 = 0.1, 0.2, . . . , 1. In this example, the signal has a sampling number of 1000 with f1 = 0.1. As shown in Fig. 1a and b, the performance comparisons between the standard EMD and the improved EMD are evaluated with MSE and ISM criterions, respectively. They are the average of 20 runs for random parameters h1, h2, b1 and b2, where h1, h2 2 [0, 2p] and b1, b2 2 [0.5, 1.5]. The abscissa 100(u1  0.1) + 10u2 is the frequency relations of the three components. For example, abscissa value 61 means u1 = 0.7 and u2 = 0.1, i.e. f1 = 10f2/7 and f2 = 10f3. Only the first three IMFs are considered in this experiment. If the number of IMFs is less than 3, simply set MSE = 1 and ISM = 1. That is why the

b

IMF4

ð16Þ

Residual

 cosð2pu2 f2 t þ h2 Þ

Time (s) 0.5 0 -0.5

0

0.05

0.1

0.15

0

0.05

0.1

0.15

0

0.05

0.1

0.15

0

0.05

0.1

0.15

0

0.05

0.1

0.15

0.5 0 -0.5 0.5 0 -0.5 0.5 0 -0.5

0.5 0 -0.5 Time (s)

Fig. 4. (a) Results by the standard EMD for the vibration signal. (b) Results by the improved EMD for the vibration signal.

801

L. Lin, J. Hongbing / Measurement 42 (2009) 796–803

is 1 kHz. Both the standard EMD method and the improved EMD method have good decomposition performance. The IMFs are approaching the original signal components. However, the improved method reduces the scale maxing, especially for IMF3.

to obtain every harmonics components by EMD without scale mixing. Here, we mainly discuss the decomposition results of the principal component, viz. the first (base) harmonic component which represents the significant features of the signal. It can be found that IMF5 in Fig. 4a and IMF4 in Fig. 4b stand for the base component, respectively. The IMFs obtained by the improved EMD will not cause misinterpretation of the signal. The improved EMD method not only reserves the significant features of the signal, but also extracts these features efficiently.

4.2. Vibration signals HHT method, including EMD is widely used for vibration signals analysis [3]. Fig. 3 shows a segment of sound vibration signal. The sampling frequency is Fs = 44.1 kHz. This signal can be represented as the sum of harmonics components with different coefficients. Because the frequencies of the harmonics components are very adjacent and their intensities differ from each other greatly, it is impossible

a

4.3. EEG signals The electroencephalogram (EEG) is widely used by physicians for interpretation and identification of physiologi-

1

0

-1

0

5

10

0

5

10

15

20 Time (s)

25

30

35

40

b1 0

-1

15

20 Time (s)

25

30

35

40

Fig. 5. (a) EEG signal collected from healthy person. (b) EEG signal collected from epileptic patient.

a

signal IMF1 IMF2 IMF3 IMF4

-20

-10

-30

-40 -50 -60

-40 -50 -60

-70

-70

-80

-80

-90

0

10

20 30 Frequency (Hz)

40

signal IMF1 IMF2 IMF3 IMF4

-20

Power/frequency (dB/Hz)

-30 Power/frequency (dB/Hz)

b

-10

50

-90

0

10

20 30 Frequency (Hz)

40

50

Fig. 6. (a) Power spectrum of the IMFs of EEG signal in Fig. 5a by the standard EMD. (b) Power spectrum of the IMFs of EEG signal in Fig. 5a by the improved EMD.

802

L. Lin, J. Hongbing / Measurement 42 (2009) 796–803

a

signal IMF1 IMF2 IMF3 IMF4

-20

-10

-30

-40 -50 -60

-40 -50 -60

-70

-70

-80

-80

-90

0

10

20 30 Frequency (Hz)

40

signal IMF1 IMF2 IMF3 IMF4

-20

Power/frequency (dB/Hz)

-30 Power/frequency (dB/Hz)

b

-10

50

-90

0

10

20 30 Frequency (Hz)

40

50

Fig. 7. (a) Power spectrum of the IMFs of EEG signal in Fig. 5b by the standard EMD. (b) Power spectrum of the IMFs of EEG signal in Fig. 5b by the improved EMD.

cal and pathological phenomena [11–13]. As shown in Fig. 5, the EEG signal is a nonlinear and non-stationary signal, which contains different kinds of background rhythms: d-wave (1–4 Hz), h-wave (4–8 Hz), a-wave (8– 13 Hz), b-wave (14–30 Hz). The sampling frequency of the EEG system is 100 Hz. The EEG signals are often interfered by noises such as power interference noise and the electromyography (EMG) noise caused by the muscle activity during collecting. These artifacts strongly influence the utility of recorded EEG signals. The IMFs obtained by EMD are band-limited, which can represent the main feature of the signals and reserve the local information. So, it can help physicians to understand EEG signals for identifying and predicting epileptic seizures. In this experiment, we utilize the standard EMD and the improved EMD to decompose the EEG signals, respectively. We only consider the first four IMFs which represent the main features of the signals. Figs. 6 and 7 are the power spectrum of the decomposition results. For healthy person’s EEG signal in Fig. 5a, the principal frequency peaks locate at about 2 Hz and 13 Hz. For epileptic patient’s EEG signal in Fig. 5b, the principal frequency peaks locate at about 4 Hz and 16 Hz. From the results we can find that the frequency band of IMFs obtained by the improved EMD method is narrower than that of the standard EMD method. And the improved EMD has less scale mixing over the standard EMD method. For example, the IMF2 by the standard EMD in Fig. 7a is an artifact frequency component, which will cause misinterpretation to the physiological phenomena. 5. Conclusions In this paper, we propose an improved EMD method for signal feature extraction. The improved method reduces

the scale maxing, thus possessing a broad application potential in real time processing of vibration signal, speech signal, EEG signal and radar signal, etc. Although it may destroy the adaptive and data-driven nature of EMD, the method proposed in this paper is efficient for the decomposition of multicomponent signals. Usually, the observed signals are corrupted by various noises, for example, Figs. 3 and 5 in Section 4. As described in this paper, the EMD method reveals a Wavelet-like decomposition structure for broadband noise. This property allows us to propose a new EMD-based signal noise reduction method. And feature extraction of a signal can be achieved after the noise reduction. Furthermore, the inverse EMD filter has a good decomposition performance for lower frequency components. So, we can combine the inverse EMD and EMD sifting processes to achieve a much better decomposition performance. Finally, there are some issues needed to be further investigated to improve the stability and validity of EMD method, such as theoretical foundation, bound effect in the cubic spline interpolation, influence of sampling rate, interpolation method, etc. Acknowledgments The authors thank the editors and reviewers for their valuable comments and good suggestions. References [1] N.E. Huang, Z. Shen, S.R. Long, M.L. Wu, H.H. Shih, Q. Zheng, N.C. Yen, C.C. Tung, H.H. Liu, The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis, Proceedings of the Royal Society of London Series A 454 (1998) 903–995. [2] H. Ding, Z.Y. Huang, Z.H. Song, Y. Yan, Hilbert–Huang transform based signal analysis for the characterization of gas–liquid twophase flow, Flow Measurement and Instrumentation 18 (2007) 37–46.

L. Lin, J. Hongbing / Measurement 42 (2009) 796–803 [3] Z.K. Peng, Peter W. Tse, F.L. Chu, An improved Hilbert–Huang transform and its application in vibration signal analysis, Journal of Sound and Vibration 286 (2005) 187–205. [4] S. Kizhner, K. Blank, T. Flatley, N.E. Huang, D. Petrick, P. Hestnes, On certain theoretical developments underlying the Hilbert–Huang transform, Proceedings of IEEE Aerospace Conference, Big Sky, MT, USA, March 2006. [5] Y. Kopsinis, S. McLaughlin, Enhanced empirical mode decomposition using a novel sifting-based interpolation points detection, Proceedings of IEEE Statistical Signal Processing, Madison, WI, USA, 2007, pp. 725–729. [6] Y. Kopsinis, S. McLaughlin, Investigation and performance enhancement of the empirical mode decomposition method based on a heuristic search optimization approach, IEEE Transactions on Signal Processing 56 (2008) 1–13. [7] G. Rilling, P. Flandrin, P. Goncalves, On empirical mode decomposition and its algorithms, Proceedings of IEEE EURASIP Workshop on Nonlinear Signal and Image Processing, Grado(I), June 2003.

803

[8] J. Stoer, R. Bulirsch, Introduction to Numerical Analysis, second ed., Springer Press, Berlin, 1992. [9] P. Flandrin, G. Rilling, P. Goncalves, Empirical mode decomposition as a filter bank, IEEE Signal Processing Letters 11 (2004) 112–114. [10] G. Rilling, P. Flandrin, One or two frequencies? The empirical mode decomposition answers, IEEE Transactions on Signal Processing 56 (2008) 85–95. [11] Z.Q. Zuo, S. Puthusserypady, Analysis of schizophrenic EEG synchrony using empirical mode decomposition, Proceedings of International Conference on Digital Signal Processing, Cardiff, 2007, pp. 131–134. [12] T. Wu, G.Z. Yan, B.H. Yang, H. Sun, EEG feature extraction based on wavelet packet decomposition for brain computer interface, Measurement 41 (2008) 618–625. [13] D. Garrett, D.A. Peterson, C.W. Anderson, M.H. Thaut, Comparison of linear, nonlinear and feature selection methods for EEG signal classification, IEEE Transactions on Neural Systems and Rehabilitation Engineering 11 (2003) 141–144.