Time-varying frequency-shifting signal-assisted empirical mode decomposition method for AM–FM signals

ARTICLE IN PRESS Mechanical Systems and Signal Processing 23 (2009) 2458–2469

Contents lists available at ScienceDirect

Mechanical Systems and Signal Processing journal homepage: www.elsevier.com/locate/jnlabr/ymssp

Time-varying frequency-shifting signal-assisted empirical mode decomposition method for AM–FM signals Xu Guanlei a,c,, Wang Xiaotong a,c, Xu Xiaogang b,c a b c

Dalian Naval Academy, Department of Navigation, Dalian of China 116018, China Dalian Naval Academy, Department of Automatization, Dalian of China 116018, China Dalian Naval Academy, Institute of Photoelectric Technology, Dalian of China 116018, China

a r t i c l e i n f o

abstract

Article history: Received 7 October 2008 Received in revised form 12 June 2009 Accepted 19 June 2009 Available online 27 June 2009

A time-varying frequency-shifting signal-assisted empirical mode decomposition (EMD) method is presented in this paper. For AM–FM signals, if the extrema density (average number of extrema per unit length) is higher, the frequency is higher approximately, and vice versa. The assisted signal in this method is directly derived from the local extrema using the least square method, and avoids the difficulty that masking signals are not adaptively examined so that the masking signal may not always enable the EMD to generate true mono-component intrinsic mode functions (IMFs). Therefore, the constructed assisted signal is adaptive and suitable for the analysis of frequencyvarying components. Compared with the traditional EMD, the method mainly solves the following problems of the traditional EMD: the traditional EMD fails to separate the modes whose frequencies are time-varying and lie within an octave, and fails to remove the time-varying frequency mode-mixing. The results of three experiments are given to show our method’s performance and efficiency. & 2009 Elsevier Ltd. All rights reserved.

Keywords: Empirical mode decomposition (EMD) Time-varying frequency-shifting signal (TFS) Assisted signal

1. Introduction and background The empirical mode decomposition (EMD) has been developed as a valuable time–frequency–magnitude resolution tool for non-stationary signals [1] and been widely used [2–4,7,11,16]. However, the original EMD has some problems in engineering work. The original EMD technique fails to separate modes whose frequencies lie within an octave [5,6] and the mode-mixing [1,7] is common in intrinsic mode functions (IMFs). For these problems, some solutions have been presented. Assisted signals (also called as masking signals in some papers [8–10]) are the most presented solutions. The using of a masking signal is firstly proposed by Deering and Kaiser [8] to improve the filtering characteristics of the EMD for the separation of modes whose frequencies lie within an octave, called masking signal-based EMD (MSEMD). However, a procedure for choosing appropriate masking signals is not adequately examined in [8]. There is also no solution proposed to the possibility that the masking signal may not always enable the EMD to generate truly mono-component intrinsic mode functions. Senroy and Suryanarayanan [9] and Senroy et al. [10] improved this algorithm and used it in power system, in which fast Fourier transform (FFT) is employed to achieve the appropriate frequencies of these stationary equivalents so as to construct the masking signals. However, the previous masking signal-based methods are only efficient for the stationary signals whose frequencies are constant. When the signals have time-varying frequencies, these methods will not work.

 Corresponding author at: Dalian Naval Academy, Department of Navigation, Institute of Photoelectric Technology, Dalian of China 116018, China. Tel.: +86 13889511675. E-mail address: [email protected] (X. Guanlei).

0888-3270/$ - see front matter & 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.ymssp.2009.06.006

ARTICLE IN PRESS X. Guanlei et al. / Mechanical Systems and Signal Processing 23 (2009) 2458–2469

2459

The other one major problem of the original EMD is the frequent appearance of mode-mixing [1,3,7], which is defined as a single intrinsic mode function either consisting of signals of widely disparate scales, or a signal of a similar scale residing in different IMF components. Mode-mixing is a consequence of signal intermittency. As discussed by Huang et al. [1,7], the intermittence could not only cause serious aliasing in the time–frequency distribution, but also make the individual IMF devoid of physical meaning. To alleviate this from occurring, Huang et al. [7] proposed the intermittence test, which can indeed ameliorate some of the difficulties. However, the approach itself has its own problems. Firstly, the intermittence test is based on a subjectively selected scale. With this subjective intervention, the EMD ceases to be totally adaptive. Secondly, the subjective selection of scales only works if there are clearly separable and definable time scales in the data. Similarly, Xu et al. [2–4] proposed the neighborhood limited empirical mode (NLEMD) for this problem through controlling the density of extrema by selected scale. NLEMD still has the problem that the decomposition ceases to be totally adaptive. Deering et al. also addressed this problem in [8] using masking signals. However, a procedure for choosing appropriate masking signals is still not adequately examined [8]. In addition to the masking signal-based signal-assisted methods, the noise-assisted methods are also proposed. Gledhill [11] used noise to test the robustness of the EMD algorithm. Flandrin et al. [12] used added noise to overcome the difficulty of the original EMD method that the method ceases to work if the data lack the necessary extrema. The noise-assisted method (EEMD) to overcome the mode-mixing is proposed by Wu and Huang [13]. EEMD utilizes the full advantage of the statistical characteristics of white noise [5] and the filter-like characteristics [14] of EMD to perturb the signal in its true solution neighborhood and to cancel itself out after serving its purpose. However, EEMD still fails to separate modes whose frequencies lie within an octave. Meanwhile, unfortunately all the mode-mixing removal methods pointed out earlier only hold for the stationary signals. Still, there are some papers, such as [17–19], involving the separation of AM–FM signals. However, they focused on the different problems with ours. In addition, some theoretical problems and the derivations are not involved in these papers [17–19]. The method in [18] needs to detect the spectral peak in the Fourier transformed domain, thus it only adapts to the stationary signals or sub-stationary signals. If all the signals are time-varying and nonlinear, there is no peak in the spectral domain. That is, in such case that method ceases to work. Similarly, the method in [19] misfits the general AM–FM signals, and the usage situation can be found in [19, Section 4.4, p. 22] i.e., the few assumptions from line -13. To say, the two improved methods are not suitable for the general AM–FM signals used by here. In this paper, a time-varying frequency-shifting signal (TFS) is derived from the extrema, which are employed for the assisted signal to separate the modes whose frequencies lie within an octave and overcome the problem of mode-mixing. The following of this paper is organized as follows. The second section of this paper briefly reviews the MSEMD and EEMD. The third section is the construction of assisted signal. In the fourth section, our improved decomposition method is derived. The mode-mixing removal is introduced in Section 5. Our experiments are given in Section 6. Finally the last section concludes this paper.

2. MSEMD and EEMD Now we will review MSEMD and EEMD briefly, respectively. MSEMD [9,10] (for the sake of similarity, here only one component’s decomposition is shown) is defined as follows: (1) Perform FFT on the signals to be decomposed to estimate frequency components fk(k ¼ 1,y,n). (2) Construct the masking signals: maskk(t) ¼ 5.5Am sin(2p(fkfk1)), where Am is the magnitude of fk obtained from FFT spectrum. (3) Obtain two signals s(t)+maskk(t) and s(t)maskk(t). Perform EMD on both signals to obtain their first IMFs only, imf+ and imf. Then ck(t) ¼ (imf++imf)/2. P (4) Obtain the residue, r k ðtÞ ¼ sðtÞ  ni¼k ci ðtÞ. (5) Perform steps (1)–(5) iteratively using all the masking signals until the final residue only contains the remaining component f1. Finally get sðtÞ ¼

n X

ci ðtÞ þ rðtÞ

(1)

i¼1

Since MSEMD uses FFT to achieve the efficient frequency components of masking signals, the non-stationary or the timevarying frequency components will not work. That is to say, in FFT spectrum, if all the frequency components overlap together, efficient frequency components are inaccessible. Another drawback of MSEMD is that it fails to separate the modes perfectly, including the stationary components, whose frequencies lie within an octave. For example, as shown in Fig. 1, there is no suitable masking signal found for the perfect separation, and it shows that MSEMD fails to separate the two components wholly. Note that there are still about 20% error (difference) at least between the real mode and the decomposed IMF.

ARTICLE IN PRESS X. Guanlei et al. / Mechanical Systems and Signal Processing 23 (2009) 2458–2469

1.2

1.2

1

1

0.8

0.8 Error

Error

2460

0.6

0.6

0.4

0.4

0.2

0.2

0 140

168

196 224 252 Freqency /Hz

280

0 170

204

238 272 306 Freqency /Hz

340

Fig. 1. The results using MSEMD: The time support is 0.5 s. In (a), the two components: c1(t) ¼ cos 280pt and c2(t) ¼ cos 200pt. The masking signal changes from 5.5 cos 280pt to 5.5 cos 560pt. In (b), the two components: c1(t) ¼ cos 340pt and c2(t) ¼ cos 200pt. The masking signal changes from P P 5.5 cos 340pt to 5.5 cos 680pt. Using MSEMD [9,10] to get the first component imf1(t). The error is defined as error l ¼ t jcl ðtÞ  imf l ðtÞj= t jcl ðtÞj. The errors can be found from the vertical axles.

EEMD [13] is defined as follows:

(1) (2) (3) (4)

Add a white noise series to the targeted data s(t). Decompose the data s(t) with added white noise into IMFs using EMD. Repeat step (1) and step (2) again and again, but with different white noise series each time. Obtain the (ensemble) means of corresponding IMFs of the decompositions as the final result.

The effects of the decomposition using the EEMD are that the added white noise series cancel each other, and the mean IMFs stays within the natural dyadic filter [14] windows, significantly reducing the chance of mode-mixing and preserving the dyadic property. Unfortunately, since EEMD holds the dyadic property of EMD, it still fails to separate modes whose frequencies lie within an octave.

3. TFS construction That EMD fails to separate the modes whose frequencies lie within an octave has been explored by Wu and Huang [5], Xu et al. [6], and Rilling and Flandrin [14,15], respectively. In the following, we will discuss the detailed approach, which avoid the drawback of previous methods.

3.1. Construction of assisted signal For AM–FM signals, the extrema are the appearance of the frequencies of the signals. If the extrema density (average number of extrema per unit length) is higher, the frequency is higher, and vice versa. EMD is adaptive in that it fully P P explores the local extrema. For the multi-component sðtÞðsðtÞ ¼ 2l¼1 sl ðtÞ ¼ 2l¼1 al cos ðfl ðtÞÞÞ, if the number or the extrema density of extrema of s(t) is exactly the same as that of the component s1(t), the approximate phase of s1(t) might be estimated using the extrema of s(t). The following is the phase estimation using the least square method from these extrema. First, find all the time locations of extrema, including the maxima and minima: t ex ¼ ft ex;1 ; t ex;2 ; . . . ; t ex;m g Pn

(2) l

Let x(t) ¼ cos f(t), fðtÞ ¼ l¼0 bl t . For tex ¼ {tex,1, tex,2,y,tex,m}, set x(tex,n) ¼ 1 and x(tex,n+1) ¼ 1, where s(tex,n) is a maximum and s(tex,n+1) is a minimum. Then we have

fðtex;n Þ ¼ ðn þ 1Þp

(3)

 2 Therefore through minimizing Ab  F2 to get the optimal coefficients b ¼ ðAT AÞ1 AT F

(4)

ARTICLE IN PRESS X. Guanlei et al. / Mechanical Systems and Signal Processing 23 (2009) 2458–2469

2461

where 3 2 bl tl 6b 7 6 ex;1 6 l1 7 7 6 tl 6 7 6 ex;2 6. b ¼ 6 .. 7; A ¼ 6 7 6 6 4 6b 7 4 1 5 t lex;m b0 2

t l1 ex;1

t ex;1



t l1 ex;2



t ex;2







t l1 ex;m



t ex;m

1

2

3

7 1 7 7 7 7 5 1

and

2p 6 3p 6 6. F¼6 6 .. 6 6 mp 4

ðm þ 1Þp

3 7 7 7 7 7 7 7 5

Thus one signal y(t) ¼ cos f(t) is achieved, which is employed as the assisted signal for the separation of components. If the signal is too long, the order of phase can be very high and thus the accuracy will decrease. To avoid this problem, the simple and better way is to estimate the phase piecewise. Since it works in the same manner, the piecewise estimation can be performed easily and is neglected here. 3.2. The estimation of amplitudes PN P For one multi-component signal sðtÞ ¼ N l¼1 sl ðtÞ ¼ l¼1 al ðtÞ cos jl ðtÞ (al(t)AR, 1rlrN, f1(t) ¼ (1/2p) djl(t)/dt4(1/ 2p) djl+1(t)/d ¼ fl+1(t)). Now the phase, jl(t), of the component sl(t) ¼ al(t) cos jl(t)(1rlrN) is known. Perform Hilbert P PN jjl ðtÞ ð1  l  NÞ (where the symbol ‘‘cpx’’ denotes the complex transform on s(t), and scpx ðtÞ ¼ N l¼1 scpx;l ðtÞ ¼ l¼1 al ðtÞe signals) are obtained. If in s(t) there is a complex component xcpx;k ðtÞ ¼ ejjk ðtÞ ð1  k  NÞ, then multiply them scpx ðtÞxcpx;l ðtÞ ¼

k1 X

al ðtÞej½jl ðtÞjk ðtÞ þ ak ðtÞ

l¼1

þ

N X

al ðtÞej½jl ðtÞjk ðtÞ

(5)

l¼kþ1

Perform Fourier transform on this Eq. (5), the spectrum can be separated into three distinct parts: ( ) k1 X spðf Þ ¼ FTfscpx ðtÞ  xcpx;l ðtÞg ¼ FT al ðtÞej½jl ðtÞjk ðtÞ ( þ FTfak ðtÞg þ FT

l¼1 N X

al ðtÞe

) j½jl ðtÞjk ðtÞ

l¼kþ1

(6) ¼ spþ ðf Þ þ sp0 ðf Þ þ sp ðf Þ Pk1 PN j½jl ðtÞjk ðtÞ j½jl ðtÞjk ðtÞ where spþ ðf Þ ¼ FTf l¼1 al ðtÞe g, sp0(f) ¼ FT{ak(t)}, sp ðf Þ ¼ FTf l¼kþ1 al ðtÞe g. FT is Fourier transform operator and IFT the reverse Fourier transform operator. In the Fourier spectrum jsp(f)j, generally there are two distinct spectral valleys between the three parts, then the three parts of signals can be separated efficiently. We can easily obtain the sp0(f) ¼ FT{ak(t)} and thus the component sk ðtÞ ¼ ak ðtÞxcpx;k ðtÞ ¼ ak ðtÞejjk ðtÞ (where ‘‘*’’ is the conjugated operator). For example, there is a complex multi-component that contains 3 mono-components (chirps): P P xðtÞ ¼ 3l¼1 xl ðtÞ ¼ 3l¼1 al ðtÞejfl ðtÞ . Note that the spectrums of the three components have overlapping in Fourier transform domain (see Fig. 2(c)). However, if one of the phases of them is known, through Fourier spectrum the according amplitude can be separated (see Fig. 2(d)): a1(t) ¼ IFT{sp0(f)}, where 12 Hzrfr500 Hz. How to get the minimum frequency, 12 Hz, in Fig. 2(d)? It is not practical that one points out the frequency and tries it every time after Fourier transformation. Note one fact that these AM–FM components have different frequencies and therefore they are orthogonal. Thus we can obtain the following relations: X ½ðxðtÞ  aðtÞ cos jl ðtÞÞðaðtÞ cos jl ðtÞÞ t

¼

X ½ðal ðtÞ cos jl ðtÞ  aðtÞ cos jl ðtÞÞðaðtÞ cos jl ðtÞÞ ¼ 0

(7)

t

In practice, through minimizing the following equation one gets the optimal estimation of the amplitude: X _ _ ½ðsðtÞ  aðtÞ cos jl ðtÞÞðaðtÞ cos jl ðtÞÞ

(8)

t

where _

al ðtÞ ¼ IFTfsp0 ðf Þgðf 0;min  f  f 0;max Þ

(9)

Through searching the two appropriate frequencies, f0,min, f0,max, in the frequency domain, one achieves the minimum of _ Eq. (9) and thus the optimal estimation of the amplitude al ðtÞ.

ARTICLE IN PRESS 2462

X. Guanlei et al. / Mechanical Systems and Signal Processing 23 (2009) 2458–2469

400

5

300

-5 5

200

-5 5

Frequency /Hz

0

0

0 -5

100

10 0 -10

0 0

0.1

0.2

0.3

0.4

0

0.5

500

1800

450

1600

400

1400

350 300 250 200 150

0.1

0.2

0.3

0.4

0.5

Time / s

Spectrum energy

Spectrum energy

time /s

1200 1000 800 600 400

100

200

50 0 500 400 300 200 100

0 0

-100 -200 -300 -400 -500

frequency / Hz

500 400 300 200 100

0

-100 -200 -300 -400 -500

frequency f /Hz

Fig. 2. (a) The actual time–frequency distributing of the three complex components: xl ðtÞ ¼ al ðtÞejfl ðtÞ ðl ¼ 1; 2; 3Þ. (b) The real parts of the three components al(t) cos fl(t)(l ¼ 1,2,3) and their sum. (c) The Fourier spectrum of the three complex components. (d) The Fourier spectrum FTfxðtÞejfl ðtÞ g.

Since the estimated component has the max frequencies of all the components, i.e., f0,max can be chosen as max{FRQ[sp(f)]}. Thus only getting the frequency f0,min(min{FRQ[sp(f)]}rf0,minr0) through minimizing Eq. (8) is sufficient. Here FRQ is the bandwidth achieving operator. Finally, in Fourier transformed domain, Eq. (8) turns to the searching of the frequency f0,min: ( ) X f 0;min ¼ arg min ½ðsðtÞ  IFTðsp0 ðf ÞÞ cos jl ðtÞÞðIFTðsp0 ðf ÞÞ cos jl ðtÞÞ (10) |fflfflfflfflffl{zfflfflfflfflffl} t f 2½0;f min \FRQ ½sp0 ðf Þ

where fmin ¼ min{FRQ[sp(f)]}. Once the frequency f0,min is determined, the optimal amplitude is also determined through Eq. (9). 4. Signal-assisted empirical mode decomposition (SAEMD) For one multi-component, more components (Z3) will greatly enlarge the difference between the extrema’s time locations of s(t) and the extrema’s time locations of s1(t). The best way is to use the IMF derived from EMD to decrease the number of components. Therefore, the classical EMD has to be employed. The following is our signal-assisted empirical mode decomposition: S1. Extract the extrema, including the local maxima and minima, of the multi-component signal s(t). S2. Use EMD to obtain imf1(t).

ARTICLE IN PRESS X. Guanlei et al. / Mechanical Systems and Signal Processing 23 (2009) 2458–2469

2463

S3. Based on imf1(t), use the LSM-based method defined in part II. A to construct TFS: cos j(t). _ S4. Perform amplitude estimation, aðtÞ, using the formulations pointed out earlier in part II. B based on the assisted signal TFS: ( ) X arg min ½ðimf 1 ðtÞ  IFTðsp0 ðf ÞÞ cos jl ðtÞÞðIFTðsp0 ðf ÞÞ cos jl ðtÞÞ (11) f 0;min ¼ |fflfflfflfflffl{zfflfflfflfflffl} t f 2½0;f min \FRQ ½sp0 ðf Þ

_

where fmin ¼ min{FRQ[sp(f)]}, al ðtÞ ¼ IFTfsp0 ðf Þgðf 0;min  f  f 0;max Þ. _ S5. The first component: c1 ðtÞ ¼ aðtÞ cos jðtÞ. S6. s(t) ¼ s(t)c1(t). S7. Repeat steps S1S6 until the residue r(t) is lower than a threshold value of error tolerance.

Note that SAEMD has to obtain the new component from the IMFs derived from traditional EMD. If the IMF is one actual AM–FM component without mixing other component, SAEMD still performs an estimation without affecting the final results.

10

350

0

300

-10 5 Frequency f/hz

250

0 -5 5 0

200 150 100

-5 20

50

0 -20

0 0

0.2

0.4 0.6 Time t/s

0.8

1.0

0

0.2

0.4

0.6

0.8

1.0

Time t/s

10

350

0

300

-10

250

5.5 5

Frequency f/hz

4.5 5 0 -5

4 3.5

200

3 2.5

150

2

5

100

0

50

-5

0

1.5 1 0.5

0

0.2

0.4

0.6

Time t/s

0.8

1.0

0

0.2

0.4

0.6

0.8

1.0

Time t/s

Fig. 3. The first simulation results for the decomposition of three components. (a) The original three components and their sum. (b) The real timefrequency distribution of the three components. (c) The result of our method. (d) The time-frequency distribution of our result. (e) The result of traditional EMD. (f) The time-frequency distribution of traditional EMD. (g) The result of EEMD. (h) The time-frequency distribution of EEMD.

ARTICLE IN PRESS 2464

X. Guanlei et al. / Mechanical Systems and Signal Processing 23 (2009) 2458–2469

10

350

0

300

10

-10

9 8

Frequency f/hz

250

5 0 -5

7 200

6 5

150

4

5

100

0

50

-5

0

3 2 1

0

0.2

0.4 0.6 Time t/s

0.8

1.0

0

10

350

0

300

-10

250

0.2

0.4 0.6 Time t/s

0.8

1.0

10 9

Frequency f/hz

8

5 0 -5

7 6

200

5 150

5

100

0

50

4 3 2

-5

1 0

0 0

0.2

0.4

0.6 Time t/s

0.8

1.0

0

0.2

0.4 0.6 Time t/s

0.8

1.0

Fig. 3. (Continued)

5. Mode-mixing removal In order to remove the mode-mixing, the construction of assisted signal should be changed accordingly. Since the distances between the adjacent extrema determine the scale of the modes, naturally they are the better selections for the control of different scale mode. For the sake of mode-mixing removal, we must obey the following rules: (1). The difference of any two distances between three adjacent extrema can not change abruptly. (2). The mode whose extrema’s distances are the smallest is the high-frequency mode to be decomposed. To obtain appropriate (2), every time before finding the extrema, we give one reference value N, the difference of any two distances between three adjacent extrema is limited by maxfjðt ex;l  t ex;lþ1 Þ=ðt ex;l1  t ex;l Þj; jðtex;l1  t ex;l Þ=ðt ex;l  t ex;lþ1 Þjg  N

(12)

If (12) is not satisfied, then the 1D case method defined in [15] is employed to decrease the difference of the adjacent two distances. Note that our method is adaptive because the ratio of the adjacent two distances is adopted instead of the constant distances adopted in [15]. In most cases, that the N is 2–5 is sufficient through our empirical simulations. In our paper, we set N ¼ 3. Once the tex ¼ {tex,1,tex,2,y,tex,m} (2) is obtained, the assisted signal can be constructed using Eqs. (3) and (4). Then the separation of the high-frequency component is the same as part III.

ARTICLE IN PRESS X. Guanlei et al. / Mechanical Systems and Signal Processing 23 (2009) 2458–2469

2465

6. Experiments and discussion 6.1. General AM–FM signal-based simulations In this section, three simulations are shown to demonstrate our method. Our results are compared with the traditional EMD by Huang et al. [1] and EEMD [13]. The sifting number of traditional EMD is 20 times. For the sake of simplicity of comparison, here we use three components and the results only include the first three components. Our sampling frequency is 1000 Hz, and the time support is 1 s. The EEMD uses noise 30 times and the variance of noise is 1. The other parameters of EEMD used are the same as our method. For the fist simulation, the parameters of the three signals (x1(t), x2(t) x3(t)) are as follows. The amplitude of x1(t) is 4.5+cos(10pt), and the frequency is decreasing from 300 to 100 Hz. The amplitude of x2(t) is 4+0.8 cos(15pt), and the frequency is decreasing from 240 to 80 Hz. The amplitude of x3(t) is 3.5+cos(4pt), and the frequency is decreasing from 60 to 20 Hz. For the second simulation, the parameters of the three signals (x1(t), x2(t), x3(t)) are as follows. The amplitude of x1(t) is 4.5+cos(10pt), and the frequency is decreasing from 200 to 50 Hz (at t ¼ 0.5 s) and increasing from 50 to 200 Hz. The amplitude of x2(t) is 4+0.8 cos(15pt), and the frequency is decreasing from 160 to 40 Hz (at t ¼ 0.6 s) and increasing from 40

10

250

0 -10

200 Frequency f/hz

5 0 -5 5 0

150

100

-5 50

20 0

0

-20 0

0.2

0.4

0.6

0.8

1.0

0

0.2

0.4

Time t/s

0.6

0.8

1.0

Time t/s

10

250

5.5

0

5 4.5

200

-10

4 Frequency f/hz

5 0 -5

3.5

150

3 2.5

100

2

5

1.5 50

1

0

0.5 -5

0 0

0.2

0.4 0.6 Time t/s

0.8

1.0

0

0.2

0.4

0.6

0.8

1.0

Time t/s

Fig. 4. The second simulation results for the decomposition of three components. (a) The original three components and their sum. (b) The real timefrequency distribution of the three components. (c) The result of our method. (d) The time-frequency distribution of our result. (e) The result of traditional EMD. (f) The time-frequency distribution of traditional EMD. (g) The result of EEMD. (h) The time-frequency distribution of EEMD.

ARTICLE IN PRESS 2466

X. Guanlei et al. / Mechanical Systems and Signal Processing 23 (2009) 2458–2469

10

250

10 9

0 200

8

Frequency f/hz

-10 5 0 -5

7 150

6 5

100

4 3

5 50

2

0

1

-5

0 0

0.2

0.4

0.6

0.8

1.0

0

0.2

Time t/s

0.4

0.6

0.8

1.0

Time t/s

10

250

10 9

0 200

8

Frequency f/hz

-10 10 0 -10

7 150

6 5

100

4 3

5 50

2

0

1

-5

0

0 0

0.2

0.4

0.6

0.8

1.0

0

Time t/s

0.2

0.4

0.6

0.8

1.0

Time t/s Fig. 4. (Continued)

to 100 Hz. The amplitude of x3(t) is 3.5+cos(4pt), and the frequency is increasing from 10 to 15 Hz (at t ¼ 0.5 s) and decreasing from 15 to 10 Hz. Note the three frequency functions are quadratic functions of time t. The first simulation (Fig. 3) includes three AM–FM components whose frequencies are linear. Since there are two (the first two) components whose frequencies lie within an octave, the traditional EMD fails to separate them. The results can also be seen from the errors that the three components by the traditional EMD and EEMD, respectively, have error1 ¼ 0.8731, error2 ¼ 0.9784, error3 ¼ 0.1454, error01 ¼ 0.8610, error02 ¼ 0.9689, and error03 ¼ 0.2995, which are much higher than ours: error1 ¼ 0.0481, error2 ¼ 0.0727, and error3 ¼ 0.0526. The visual results are shown clearly in Fig. 3(c)–(h). The second simulation (Fig. 4) includes three AM–FM components whose frequencies are quadric. Since there are two components whose frequencies lie within an octave, the traditional EMD fails to separate them. The errors that the three components by the traditional EMD and EEMD are, respectively: error1 ¼ 0.6002, error2 ¼ 0.7529, error3 ¼ 0.5595, error01 ¼ 0.7986, error02 ¼ 0.8497, and error03 ¼ 1.0123, which are much higher than ours: error1 ¼ 0.1034, error2 ¼ 0.1240, and error3 ¼ 0.0738. The visual results are shown clearly in Fig. 4. Note that the components whose frequencies lie within an octave are often cooperated into a new component or their parts are cooperated into a new component. The traditional EMD breaks them into new components. This not only leads to the incomplete decomposition, but also leads to some mode-mixing in the IMFs. For example, in the first simulation (Fig. 3), the first two components are mostly included in the first IMF (the first IMF of Fig. 3(e)) by the traditional EMD. However, in the second simulation (Fig. 4), the second IMF (the middle IMF of Fig. 4(e)) by the traditional EMD not only includes a part of the second component, but also includes a part of the third component, which are different modes at all. Instead, EEMD only performs the action of one filter bank [12,14].

ARTICLE IN PRESS X. Guanlei et al. / Mechanical Systems and Signal Processing 23 (2009) 2458–2469

2467

For the third simulation, the parameters of the three signals (x1(t), x1(t), x1(t)) are as follows. The amplitude of x1(t) is 3e(t0.5)2/0.95, and the frequency is decreasing from 200 to 100 Hz. The amplitude of x2(t) is 4+0.8 cos(15pt), and the frequency is decreasing from 50 to 25 Hz. The amplitude of x3(t) is 3.5+cos(4pt), and the frequency is decreasing from 10 to 5 Hz. The last simulation (Fig. 5) includes three AM–FM components whose frequencies are linear. Differently, there is an intermittent mode. Using our method, the ideal result (Fig. 5(c)) is obtained rather than the mode-mixing result of traditional EMD (Fig. 5(e)). The according results are also seen from the time–frequency distributions (Fig. 5(d) and (f)). The errors that the three components by the traditional EMD and EEMD are, respectively: error1 ¼ 0.3607, error2 ¼ 0.1716, error3 ¼ 0.1268, error01 ¼ 0.5343, error02 ¼ 0.9410, and error03 ¼ 1.2735, which are much higher than ours: error1 ¼ 0.0377, error2 ¼ 0.0221, and error3 ¼ 0.0209. 7. Conclusion In this paper, a time-varying frequency-shifting signal-assisted EMD is developed, which overcomes the shortcomings that other signal-assisted EMD methods have. The assisted signal is derived from the local extrema, which is adaptive and suitable for the analysis of frequency-varying components. Since other signal-assisted EMD methods will not work in the case that the frequencies are time-varying, the comparison with other signal-assisted EMD methods is not given except for the traditional EMD and EEMD. Our method mainly solves the problem that traditional EMD and EEMD fail to separate the modes whose frequencies are time-varying and lie within an octave. The time-varying frequency mode-mixing, previous

5

250

0 -5 5 Frequency f/hz

200

0 -5 5 0 -5 20

150

100

50

0 -20

0 0

0.2

0.4 0.6 Time t/s

0.8

1.0

0

0.2

0.4

0.6

0.8

1.0

Time t/s

250

5

4.5 0 4

200 -5 Frequency f/hz

3.5

5 0 -5

3

150

2.5 2

100

1.5

5 50

1

0

0.5

-5

0 0

0.2

0.4

0.6 Time t/s

0.8

1.0

0

0.2

0.4

0.6

0.8

1.0

Time t/s

Fig. 5. The simulation for mode-mixing removal. (a) The original three components and their sum. (b) The real time-frequency distribution of the three components. (c) The result of our method. (d) The time-frequency distribution of our result. (e) The result of traditional EMD. (f) The time-frequency distribution of traditional EMD. (g) The result of EEMD. (h) The time-frequency distribution of EEMD.

ARTICLE IN PRESS 2468

X. Guanlei et al. / Mechanical Systems and Signal Processing 23 (2009) 2458–2469

250

10

6

0 200

5 Frequency f/hz

-10 5 0 -5

150

4 3

100

2

5 50

1

0 0

-5 0

0.2

0.4

0.6

0.8

1.0

0

0.2

0.4

0.6

0.8

1.0

Time t/s

Time t/s

250

5

4.5

0

4 200 3.5 Frequency f/hz

-5 5 0 -5

3

150

2.5 2

100

1.5

5

1

50 0

0.5 0

-5 0

0.2

0.4

0.6

0.8

1.0

0 0

Time t/s

0.2

0.4

0.6

0.8

1.0

Time t/s Fig. 5. (Continued)

methods fail to remove, may be removed using our method as well. Finally, three simulations are given to show our method’s efficiency.

Acknowledgement First we will thank the reviewers for their suggestions and the editors for their work very much. This work was partly supported by the NSF of Chinese Liaoning province and the Kaifang Foundation of Zhejiang University of China. References [1] N.E. Huang, Z. Shen, S.R. Long, M.C. Wu, 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 454 (1998) 903–995. [2] G.L. Xu, X.T. Wang, X.G. Xu, et al., Image enhancement algorithm based on neighborhood limited empirical mode decomposition, Acta Electronica Sinica 34 (3) (2006) 99–103. [3] G.L. Xu, X.T. Wang, X.G. Xu, Neighborhood limited empirical mode decomposition and application in image processing, IEEE Proceedings of the ICIG (2007) 149–154. [4] G.L. Xu, X.T. Wang, X.G. Xu, T. Zhu, Multi-band image fusion algorithm based on neighborhood limited empirical mode decomposition, Journal of Infrared and Millimeter Waves 25 (3) (2006) 225–228. [5] Z. Wu, N.E. Huang, A study of the characteristics of white noise using the empirical mode decomposition method, Proceedings of the Royal Society of London A 460 (2004) 1597–1611. [6] G.L. Xu, X.T. Wang, X.G. Xu, et al., The conditions and principle for the direct decomposition from multi-component to mono-component using EMD, Progress in Natural Science of China 16 (10) (2006) 1356–1360 (in Chinese).

ARTICLE IN PRESS X. Guanlei et al. / Mechanical Systems and Signal Processing 23 (2009) 2458–2469

2469

[7] N.E. Huang, Z. Shen, S.R. Long, A new view of nonlinear water waves: the Hilbert spectrum, Annual Review of Fluid Mechanics 31 (1999) 417–457. [8] R. Deering, J.F. Kaiser, The use of a masking signal to improve empirical mode decomposition, ICASSP 2005, vol. IV, IEEE, pp. 485–488. [9] N. Senroy, S. Suryanarayanan, Two techniques to enhance empirical mode decomposition for power quality applications, in: Power Engineering Society General Meeting, IEEE, Tampa, FL, USA, 2007, pp. 1–6. [10] N. Senroy, S. Suryanarayanan, Paulo F. Ribeiro, An improved Hilbert–Huang method for analysis of time-varying waveforms in power quality, IEEE Transactions on Power Systems 22 (4) (2007) 1843–1850. [11] R.J. Gledhill, Methods for investigating conformational change in biomolecular simulations, A dissertation for the degree of Doctor of Philosophy, Department of Chemistry, The University of Southampton, 2003, 201pp. [12] P. Flandrin, P. Gonc- alve`s, G. Rilling, EMD equivalent filter banks, from interpretation to applications, in: N.E. Huang, S.S.P. Shen (Eds.), Hilbert–Huang Transform: Introduction and Applications, World Scientific, Singapore, 2005, pp. 67–87, 360pp. [13] Z. Wu, Norden E. Huang, Ensemble empirical mode decomposition: a noise assisted data analysis method, COLA Technical Report, 2005. Available from: /ftp://grads.iges.org/pub/ctr/ctr_193.pdfS. [14] P. Flandrin, G. Rilling, P. Gonc-alve`s, Empirical mode decomposition as a filter bank, IEEE Signal Processing Letters 11 (2004) 112–114. [15] G. Rilling, P. Flandrin, One or two frequencies? The empirical mode decomposition answers, IEEE Transactions on Signal Processing 56 (1) (2008) 85–95. [16] J. Cheng, D. Yu, Y. Yang, A fault diagnosis approach for roller bearings based on EMD method and AR model, Mechanical Systems and Signal Processing 20 (2006) 350–362. [17] Y. Qin, S. Qin, Y. Mao, Research on iterated Hilbert transform and its application in mechanical fault diagnosis, Mechanical Systems and Signal Processing 22 (8) (2008) 1967–1980. [18] W. Yang, Interpretation of mechanical signals using an improved Hilbert–Huang transform, Mechanical Systems and Signal Processing 22 (5) (2008) 1061–1071. [19] M. Feldman, Time-varying vibration decomposition and analysis based on the Hilbert transform, Journal of Sound and Vibration 295 (3–5) (2006) 518–530.