A new method for multicomponent signal decomposition based on self-adaptive filtering

Measurement 44 (2011) 1312–1327

Contents lists available at ScienceDirect

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

A new method for multicomponent signal decomposition based on self-adaptive filtering Yi Qin ⇑, Baoping Tang, Jiaxu Wang, Xiao Ke State Key Laboratory of Mechanical Transmission, Chongqing University, Chongqing 400044, People’s Republic of China

a r t i c l e

i n f o

Article history: Received 26 June 2009 Received in revised form 8 July 2010 Accepted 23 February 2011 Available online 1 March 2011 Keywords: Orthogonal empirical mode decomposition (OEMD) Multicomponent signals Analytic band-pass filtering Algorithm Mechanical vibration signal Fault diagnosis

a b s t r a c t Since the empirical mode decomposition (EMD) lacks strict orthogonality, a new method for multicomponent signal decomposition, orthogonal empirical mode decomposition (OEMD), is proposed by this paper. The essential principle of this method is to obtain the intrinsic mode functions (IMFs) and the residue by self-adaptive band-pass filtering. Firstly, the feasibility of OEMD is theoretically analyzed, then its strict orthogonality and completeness is proved, and the orthogonal basis used in OEMD is generated. Secondly, the method of analytical band-pass filtering which preserves perfect band-pass feature in the frequency domain is presented, then two fast algorithms to implement OEMD are proposed, i.e. IMF sequential searching (ISS) algorithm and IMF binary searching (IBS) algorithm. The speed of IBS is faster than that of ISS, whereas IBS algorithm may obtain much more false IMFs than ISS when signals are of complex spectral constitutions. Finally, OEMD is applied to both synthetic signals and mechanical vibration signals, the results show that compared with EMD, OEMD better solves mode aliasing, avoids the occurrence of false mode, is free of end extension, and can be effectively applied to mechanical fault diagnosis. Ó 2011 Elsevier Ltd. All rights reserved.

1. Introduction Multicomponent signal analysis is still a big problem in the signal processing domain. Hilbert Huang Transform (HHT), is a new method proposed by Huang et al., which is very suitable for the analysis of multicomponent signals [1,2]. Nowadays, it has obtained good results when applied to the field of ocean engineering, biomedical engineering, earthquake signal and structure analysis, bridge and building monitoring, electric net detecting, fault diagnosis, etc. [3–9]. Therefore, it gradually becomes a hot theme point in the signal processing field. The major innovation of HHT is the putting forward of intrinsic mode function (IMF) and the invention of empirical mode decomposition (EMD) [10]. And the essence of EMD is to act self-adaptive filtering on the signal, and obtain a finite number of IMFs

⇑ Corresponding author. Tel.: +86 023 65106974; fax: +86 023 65106641. E-mail address: [email protected] (Y. Qin). 0263-2241/$ - see front matter Ó 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.measurement.2011.02.014

and a residue. All the obtained IMFs are monocomponent signals, which can make the results of Hilbert transform have reasonable physical meaning. However, there are three main problems: (a) End effect problem. When applying EMD method to sift the signal, two ends of time series will suffer from the disperse phenomenon, and this disperse will diffuse towards the inside of the signal and ‘‘empoison’’ the whole time series [11]. Besides, when applying Hilbert transform to each IMF, there are end effects occurred at the two ends. There is not a perfect solution to this problem till now. (b) The problem of mode aliasing. Since EMD can pick out the high frequency component at any time instant, when there is abnormity in the signal, it cannot sift it out [12]. Furthermore, when the frequencies of some components in the signal are close [13], EMD cannot separate them so as to form mode aliasing. This phenomenon indicates in some degree that EMD is not an orthogonal decomposition method. (c) The problem of false mode. As EMD is not a strictly orthogonal decomposition method, there are some meaningless components in the results of EMD for some

1313

Y. Qin et al. / Measurement 44 (2011) 1312–1327

signals. This limitation may influence the accuracy of analysis for the signal. Aiming at these problems, an orthogonal decomposition algorithm needs to be researched. EMD acts essentially as a dyadic filter bank resembling those involved in wavelet decompositions [14,15], thus we can obtain the IMFs by means of a proper adaptive filtering algorithm, and orthogonal empirical mode decomposition (OEMD) is innovationally proposed. In this paper, firstly, the feasibility of obtaining IMFs by band-pass filtering is proved, then the process of OEMD is presented, and the orthogonality and completeness of OEMD are proved theoretically, and the orthogonal basis functions are generated. Then, the method of analytical band-pass filtering is presented, and its perfect band-pass characteristic in the frequency domain is demonstrated. Based on analytical band-pass filtering, two fast algorithms to implement OEMD are proposed, i.e. IMF sequential searching (ISS) algorithm and IMF binary searching (IBS) algorithm, and the two algorithms are compared. ISS can obtain the IMF that have maximum bandwidth, therefore the number of IMFs obtained by the ISS is less than that obtained by IBS, especially for complex signals, but its calculation speed is lower than IBS. The decomposition error and convergence of the two algorithms are also discussed. Finally, the proposed method is applied to synthetic signals and mechanical vibration signals and is compared with EMD. The results show that OEMD can be applied to multicomponent signal decomposition and the analysis of mechanical vibration signals. Furthermore, this proposed method can effectively decompose the signal into a sum of IMFs without end extension, and has better performance in conquering mode aliasing and avoiding the occurrence of false modes than EMD.

m1 ðtÞ ¼ xðtÞ  c1 ðtÞ

ð4Þ

Regarding m1(t) as a new x(t), we can filter out other IMFs of the original signal one by one. Finally, we can write x(t) as

xðtÞ ¼

n X

ci ðtÞ þ rðtÞ

ð5Þ

i¼1

2. Orthogonal empirical mode decomposition

where ci(t) is the ith IMF, r(t) is the residue, i.e. the trend item, n is the number of IMFs in the signal. In order to comprehend the decomposition rule of EMD, let us see the local decomposition of a signal x(t), which is shown in Fig. 1. Let us consider the time instant tC. We can see from Fig. 1 that there are two fluctuant periods in the local range of tA to tB; however, the envelope esup(t) and einf(t) only has less than one fluctuant period. We assume that the signal c(t) obtained by subtracting the mean signal m(t) = [esup(t) + einf(t)]/2 from x(t) is an IMF. And the next IMF c0 (t) is obtained by decomposing m(t), therefore, its instantaneous frequency at the time instant tC is lower than that of c(t), which indicates that EMD has the characteristic of selfadaptive multi-resolution and can self-adaptively separate the component of high frequency from the signal. It then follows that EMD is a process of self-adaptive filtering, and the frequency bands of c1(t), . . . , cn(t) are from high to low. Based on this characteristic, we can also make use of self-adaptive band-pass filtering to obtain IMFs. In the following, a simple and fast method to decompose a signal into the sum of IMFs and the residue will be elaborated. Suppose that the non-negative frequency range of the real signal x(t) is [0, fc), then we can find a frequency f1 via a certain searching algorithm which makes the filtered result after x(t) passing through the band-pass filter

2.1. Principle of OEMD

H1 ðf Þ ¼

In Ref. [1] Huang gave the definition of an IMF as that an IMF is a function that satisfies two conditions: (a) 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 (b) at any time instant, the mean value of the envelope defined by the local maxima and the envelope defined by the local minima is zero. The process to decompose a signal into the sum of IMFs is: first find out all the extreme points of the signal (supposed to be x(t)), then fit upper envelope esup(t) and lower envelope einf(t) respectively with all the maximum points and minimum points, which satisfy

einf ðtÞ 6 xðtÞ 6 esup ðtÞ



1; f1 6 jf j < fc

satisfy the requirements of IMF’s definition, hence it becomes the first IMF, and we note it as c1(t). In order to achieve multicomponent AM-FM demodulation, H1(f) should be the widest band-pass filter which enables us to obtain reasonable IMFs after filtering. Then the second IMF c2(t) is searched in the frequency range of [0, f1). Suppose that the low-limit frequency of current widest bandpass filter H2(f) is f2, then we continue to search for the third IMF in the frequency range of [0, f2). This process is repeated until fn is zero or the signal r(t) obtained by filtering x(t) with the low-pass filter

ð1Þ

B

A C

Let the mean of upper and lower envelopes be m(t)

mðtÞ ¼ ðesup ðtÞ þ einf ðtÞ=2Þ

ð6Þ

0; elsewhere

x(t )

ð2Þ

Subtracting m(t) from x(t), we obtain

cðtÞ ¼ xðtÞ  mðtÞ

esup (t )

einf (t )

ð3Þ

Regarding c(t) as a new x(t), and repeating the process above, we can filter out the first-order IMF c1(t) of the original signal (determined by stop criterion). Note that

o

tA

tC

tB

Fig. 1. Local decomposition of signal with EMD.

t

1314

Y. Qin et al. / Measurement 44 (2011) 1312–1327

 Hd ðf Þ ¼

1; 0 6 jf j < fn 0; elsewhere

n X

is a monotonic function (there are no extrema in the function) or preserves little enough energy. Then the decomposition is terminated and we get n IMFs c1(t), . . . , cn(t) and a residue r(t) (when fn is zero, let r(t) = 0). We know from the above process of decomposition that the band-pass filter used to obtain the kth IMF ck(t), is

 Hk ðf Þ ¼

1; fk 6 jf j < fk1 0; elsewhere

where X(f) = FT{x(t)} (FT stands for Fourier transform), then the kth IMF is

Proof. According to the Parseval theorem, we can obtain

1

¼

Z

1

1 Z 1

C j ðf ÞC k ðf Þ df Hj ðf ÞHk ðf ÞXðf ÞX  ðf Þ df

1

From the decomposition rules of the proposed method we have fj P fk1, therefore according to Eq. (8) we get

Hj ðf ÞHk ðf Þ ¼ 0

1

1

cj ðtÞck ðtÞ dt ¼ 0

ð9Þ

The proposition is proved. h Then let us prove that r(t) and any IMF ck(t) are orthogonal. Proof. According to the Parseval theorem, we can obtain

Z

1

1

! Hi ðf Þ þ Hd ðf Þ Xðf Þ ¼ Xðf Þ

i¼1

it follows that n X

ci ðtÞ þ rðtÞ

ð11Þ

i¼1

Since the components obtained in above decomposition process are orthogonal to each other, the process is called orthogonal empirical mode decomposition (OEMD). Eq. (11) shows that the OEMD have completeness as well. Now let us consider the impulse response of the filter’s transfer function. Using the inverse Fourier transform, the impulse response hk(t) of Hk(f) can be given by

Z

1

Hk ðf Þej2pft df

¼ ½sinðj2pfk1 tÞ  sinðj2pfk tÞ=ðptÞ

ð12Þ

Then by the Parseval theorem, we can obtain

Z

1

1



hj ðtÞhk ðtÞ dt ¼

Z

1

Hj ðf ÞHk ðf Þ df ¼ 0

ð13Þ

1

where j < k. Hence, hk(t) can be regarded as the orthogonal basis function and {hi(t), i = 1, 2, . . .} forms an orthogonal basis in OEMD. In addition, it can be seen from Eq. (12) that hk(t) is not compactly supported, so it is not suitable for filtering in the time domain, and a special filtering approach in the frequency domain needs to be applied for implementing OEMD, as will be shown in Section 3. 2.2. Remarks on OEMD

Hence

Z

¼

n X

1

where IFT stands for inverse Fourier transform. Now let us prove that two arbitrary IMFs cj(t) and ck(t) (j < k) are orthogonal.

cj ðtÞck ðtÞ dt ¼

Hi ðf ÞXðf Þ þ Hd ðf ÞXðf Þ

i¼1

hk ðtÞ ¼

ck ðtÞ ¼ IFTfC k ðf Þg

n X

i¼1

ð8Þ

C k ðf Þ ¼ Hk ðf ÞXðf Þ

1

C i ðf Þ þ Rðf Þ ¼

xðtÞ ¼

When k = 1, there is f0 = fc. Let

Z

And since

ð7Þ

rðtÞck ðtÞ dt ¼ ¼

Z

1

1 Z 1

Rðf ÞC k ðf Þ df Hd ðf ÞHk ðf ÞXðf ÞX  ðf Þ df

1

Since fk P fn, it follows from Eqs. (7) and (8) that

Hj ðf ÞHk ðf Þ ¼ 0

Both EMD and OEMD can decompose a signal into the sum of IMFs, whereas the two decomposition methods are essentially different. The essence of EMD is to separate the finest local mode from the signal based on the characteristic time scale, which is defined by the time interval between the extrema of data; while OEMD obtains the IMFs by self-adaptive band-pass filtering. It is clear that selfadaptive band-pass filtering is based on Fourier analysis, which is a global method. Consequently, the local characteristic of EMD is better than that of OEMD. However, OEMD has strict orthogonality, and overcomes some drawbacks of EMD mentioned in Section 1. In Section 4, some examples will be applied to demonstrate the advantages of OEMD. 3. Fast implementation of OEMD

Hence

Z

1

1

rðtÞck ðtÞ dt

¼0

Then the proposition is proved. h

ð10Þ

It is easy to note that the decomposition is orthogonal and complete if the ideal filtering performance can be assured in the frequency domain from above proof process of orthogonality and completeness. This deduction makes the numeric implementation of OEMD possible. Firstly,

1315

Y. Qin et al. / Measurement 44 (2011) 1312–1327

suppose that we sample the real signal x(t) with the sampling frequency of fs, and we get a discrete time series of length N (N is usually an even number), x(n) (n = 0, 1, . . . , N  1). In the same way, without taking into account the negative frequencies, the discrete form of the transfer function of the filter Hk(f) described in Eq. (8) can be given by

 Hk ðmÞ ¼

1 mk 6 m < mk1 6 N=2 0

0 6 m < mk or mk1 6 m < N

ð14Þ

3.1. Analytical band-pass filtering of discrete signal Using the Hilbert transform [16], the analytical signal of a real signal x(n) can be defined as

ð15Þ

And its discrete spectrum is

8 Xð0Þ > > > < 2XðmÞ Z k ðmÞ ¼ > XðN=2Þ > > : 0

m¼0 0
ð16Þ

N=2 < m < N

where ^ xðnÞ is the discrete Hilbert transform of x(n). The discrete spectrum of the filtered result after z(n) passing through the band-pass filter described in Eq. (14) is

Y k ðmÞ ¼ Z k ðmÞHk ðmÞ

ð17Þ

By calculating the discrete IFT of Yk(m), we get the result of band-pass filtering for the discrete analytical signal

yk ðnÞ ¼ pk ðnÞ þ jqk ðnÞ

ð18Þ

It is clear that yk(n) is a discrete analytical signal, thus qk(n) should be the discrete Hilbert transform of pk(n). Now let us deduce the relation between qk(n) and x(n). By performing Fourier transform on both sides of Eq. (18), we can write

Y k ðmÞ ¼ Pk ðmÞ þ jQ k ðmÞ

ð19Þ

where Pk(m), Qk(m) are respectively the discrete Fourier transform of pk(n) and qk(n). It follows from the definition of Hilbert transform that

8 Pk ð0Þ > > > < 2P ðmÞ k Pk ðmÞ þ jQ k ðmÞ ¼ > P ðN=2Þ k > > : 0

m¼0 0 < m < N=2 m ¼ N=2 N=2 < m < N

Y k ðmÞ ¼ Z k ðmÞHk ðmÞ ¼

2XðmÞ mk 6 m < mk1 6 N=2 0 0 6 m < mk or mk1 6 m < N

ð21Þ

Then according to the characteristics of the discrete spectrum of the real signal, we have

X k ðmÞ ¼ X k ðN  mÞ;

m ¼ 0; 1; . . . ; N  1

ð22Þ

Therefore, the discrete spectrum of pk(n) at the whole interval can be written as

8 0 > > > > > > < XðmÞ P k ðmÞ ¼ 0 > > > XðmÞ > > > : 0

0 6 m < mk mk 6 m < mk1 mk1 6 m 6 N  mk1

ð23Þ

N  mk1 < m 6 N  mk N  mk < m < N

Similarly, when mk = 0, we get

8 > < XðmÞ 0 6 m < mk1 P k ðmÞ ¼ 0 mk1 6 m 6 N  mk1 > : XðmÞ N  mk1 < m < N

ð24Þ

From Eqs. (23) and (24) we can see that pk(n) is the result of acting band-pass filtering on x(n), and its frequency is strictly restricted in the range of [fk, fk1). The above filtering process is named as ‘‘analytical band-pass filtering’’ by the authors. Since pk(n) is a real signal with definite physical meaning, and its frequency is strictly in the range of [fk, fk1), this filtering approach can be applied to implement OEMD. The analytical band-pass filtering not only is a method preserving perfect frequency domain characteristics, but also obtains a useful byproduct, i.e. the discrete Hilbert transform of pk(n), which provides convenience for further Hilbert spectrum analysis. However, it is worth noting that after filtering the signal with above method, though the frequency of pk(n) is strictly restricted in [fk, fk1), it is not always able to completely recover the corresponding component in the original signal. The reason is as follows. For a multicomponent signal, the spectra of some components may have a wide frequency range, therefore, only making use of the frequency information of one pass band to reconstruct the corresponding component will bring error when the bandwidth of the component is not equal to the bandwidth of the filter or the spectra of some components lap over. Specially, for harmonic signals, spectrum leakage will bring error to the filtered result, whereas without spectrum leakage the filtered result is consistent with the ideal result (see Section 4.1). The performance of this filtering method can be tested by the following examples. Consider the two-component AM-FM synthetic signal

xðtÞ ¼ x1 ðtÞ þ x2 ðtÞ ð20Þ

When mk – 0, from Eqs. (14), (16), and (17) we get 

8 0 6 m < mk > <0 P k ðmÞ ¼ XðmÞ mk 6 m < mk1 > : 0 mk1 6 m 6 N=2

Pk ðmÞ ¼ Pk ðN  mÞ;

where mk = N  fk/fs, mk1 = N  fk1/fs. By carrying out IFT to Hk(m), we can obtain its impulse response hk(n). It is easy to note that hk(n) is analytical. Therefore, we may transform the discrete time signal x(n) into corresponding discrete analytical signal z(n) before filtering x(n) by means of the filter described in Eq. (14). In the following, we will introduce the analytical band-pass filtering of the discrete signal.

zðnÞ ¼ xðnÞ þ j^xðnÞ

Since the two sides of Eq. (19) are equal, we can obtain

ð25Þ

with

x1 ðtÞ ¼ cos½2p ð60t þ20t 2 Þ

ð26Þ

x2 ðtÞ ¼ ½1þ0:2sinð2p 10tÞcos½2pf0 t þ0:5sinð2p 10tÞÞ ð27Þ where f0 = 160 Hz. The sampling frequency fs = 500 Hz and the sampling points N = 512. The time domain waveform

1316

Y. Qin et al. / Measurement 44 (2011) 1312–1327

and spectrum of x(t) are shown in Fig. 2. Then, in order to separate the chirp signal x1(t), the proposed filtering approach is applied to this signal (the pass band is [20, 140] Hz). The filtered signal and its spectrum are respectively shown in Fig. 3a, the error between the filtered result and the ideal result is shown in Fig. 3b. From Fig. 3 it can be seen that the filtered result are almost the same as the ideal result (the chirp signal x1(t)), which means this proposed method based on analytical bandpass filtering can effectively separate the chirp signal. In order to further illustrate the performance of this proposed filtering approach, another experiment is carried out. We use the same kind of synthetic signal, but change the carrier frequency f0 of x2(t) from 160 Hz to 145 Hz. Under the same sampling conditions, its time domain waveform and spectrum are shown in Fig. 4. Similarly, analytical band-pass filtering is applied to the signal (the pass band is [20, 140] Hz). The filtered signal and its spectrum is shown in Fig. 5a, the error between the filtered result and the ideal result is shown in Fig. 5b. It can be seen from Fig. 5 that the spectra of the two components lap over in one frequency range which is included in the filtering pass band, thus the error between the filtered result and the ideal result is larger than that of the previous experiment. This is a limitation of analytical band-pass filtering. 3.2. The fast algorithms for OEMD 3.2.1. Basic notations and IMF criterion The notations used in the following two algorithms are listed below.        

Hd(m) – the discrete low-pass filter; H(m) – the discrete band-pass filter; x(n) – original signal; y(n) – filtered signal; N – the length of original signal; ci(n) – the ith IMF; r(n) – residue; M – the index of exploratory point, its corresponding real physical frequency is Mfs/N;

 bmin – lower frequency bin limit of the band range to be searched;  bmax – upper frequency bin limit of the band range to be searched;  ni – lower frequency bin limit of the band range to be searched when searching for the IMF with the widest pass band, i.e. the frequency bin corresponding to the situation that when the signal passes through a band-pass filter with the pass band [ni, bmax), we get a non-IMF;  yi – upper frequency bin limit of the band range to be searched when searching for the IMF with the widest pass band, i.e. the frequency bin corresponding to the situation that when the signal passes through a bandpass filter with the pass band [yi, bmax), we get an IMF. The second condition required in the definition of IMF is hard to be satisfied in actual applications. Hence, we use the following criterion to judge whether a signal y(t) is an IMF. First, we judge whether y(t) satisfies the first condition; then the local mean m(t) is calculated by averaging upper and lower envelopes of y(t). Define the standard deviation (SD) by

PLy 2 m ðtÞ SD ¼ Pt¼0 Ly 2 t¼0 y ðtÞ

ð28Þ

where Ly is the length of y(t). A typical threshold for SD can be set between 0.08 and 0.15 if the SD value is lower than the threshold, we regard y(t) as an IMF. 3.2.2. ISS algorithm Based on the analytical band-pass filtering method, we can discuss how to obtain IMFs and the residue fast. It is known from Section 2 that each IMF represents a monocomponent signal with a certain frequency band. Therefore, it is only possible to obtain each IMF and the residue by exploratory searching. The first method coming to your mind might be sequential searching, and the procedure of sequential searching algorithm can be represented as follows. Initialization Let the frequency bin range of x(n) be [0, N/2); the upper frequency bin limit of the band range

Fig. 2. Synthetic signal x(t) (f0 = 160 Hz) and its spectrum.

Y. Qin et al. / Measurement 44 (2011) 1312–1327

1317

Fig. 3. Performance of the proposed filtering method when f0 = 160 Hz (the pass band is [20, 140] Hz): (a) filtered result and its spectrum, (b) error signal between the filtered result and the ideal result.

Fig. 4. Synthetic signal x(t) (f0 = 145 Hz) and its spectrum.

to be searched is denoted by Mmax, and let Mmax = N/2, M = 0 and i = 0.  Step 1. The discrete low-pass filter used in the decomposition is given by

( Hd ðmÞ ¼

1 0 6 m < M max ð29Þ 0 M max 6 m < N

After applying analytical band-pass filtering to the original signal x(n) with the filter described in Eq. (29), the filtered result y(n) is obtained. Then we calculate the energy ratio between y(n) and x(n)

P jyðnÞj2 e ¼ Pn 2 n jxðnÞj

ð30Þ

if e is less than a certain threshold T, which can be set as 0.05, i.e. the energy of y(n) is little enough compared to

1318

Y. Qin et al. / Measurement 44 (2011) 1312–1327

Fig. 5. Performance of the proposed filtering method when f0 = 145 Hz (the pass band is [20, 140] Hz): (a) filtered result and its spectrum, (b) error signal between the filtered result and the ideal result.

x(n), or if y(n) becomes a monotonic function, we terminate the decomposition and note r(n) = y(n). Otherwise, we further judge whether y(n) is an IMF: if it is an IMF, let i i + 1, ci(n) = y(n) and r(n) = 0, then terminate the decomposition; if not, let M M + 1 and go to the next step.  Step 2. The discrete band-pass filter used in the decomposition is defined as

 HðmÞ ¼

1 M 6 m < M max 0

0 6 m < M or M max 6 m < N

ð31Þ

After acting analytical band-pass filtering upon x(n) with the filter H(m), the filtered result y(n) is obtained. If y(n) is an IMF, let i i + 1, ci(n) = y(n), Mmax = M and M = 0, then return to step 1; Otherwise, let M M + 1, then go back to step 2 and repeat. In this paper, the above algorithm is named as IMF sequential searching (ISS). To demonstrate the decomposition processes clearly, Fig. 6 shows the flowchart of ISS algorithm. From its rule of decomposition, it is easy to note that ISS can assure the obtained IMF preserve maximum bandwidth, so it is very suitable for the decomposition of the multicomponent AM–FM signal, but the efficiency is low.

3.2.3. IBS algorithm Considering that the frequency band reflected by each IMF is arranged from high frequency to low frequency, we propose another approach to obtain IMFs based on binary searching. The procedure of binary searching algorithm is represented as follows. Initialization bmin = 0. and bmax = N/2, and let i = 0.  Step 1. If x(n) is an IMF, then we terminate the decomposition and note c1(n) = x(n), r(n) = 0. Otherwise, take the exploratory point as M ¼ bðbmin þ bmax Þ=2c, where bxc represents the biggest integer without exceeding x, and let ni = 0, then go to the next step.  Step 2. Define the discrete band-pass filter as

HðmÞ ¼



1 M 6 m < bmax 0

0 6 m < M or bmax 6 m < N

ð32Þ

After analytical band-pass filtering, the filtered signal y(n) is obtained, and then the energy ratio e between y(n) and x(n) with Eq. (30). If e is less than the certain threshold value T1 (it depends on the component which have minimum energy in the signal, and it is set as 0.05 in this paper), the lower limit of the searching frequency band should be changed. Let M bðM þ bmin Þ=2c, then come back to step 2, and repeat the steps above. Otherwise, let bmin = M and go to the next step.

Y. Qin et al. / Measurement 44 (2011) 1312–1327

1319

Start M max = N / 2 , M = 0 , i = 0

Construct Hd(m) with Eq. (29)

y(n) is obtained by filtering x(n) with Hd(m)

Calculate e with Eq. (30)

yes

Is e
Is y(n) an IMF? yes

M ←M+1 r ( n) = y ( n)

i ← i +1

Construct H(m) with Eq. (31)

y(n) is obtained by filtering x(n) with H(m)

c i ( n) = y ( n)

no r ( n) = 0

Is y(n) an IMF? yes i ← i + 1 , ci ( n) = y ( n) , M max = M , M = 0

End Fig. 6. The flowchart of ISS algorithm.

 Step 3. If y(n) is an IMF, let yi = M, and go to the next step. If y(n) is not an IMF, we take M ¼ bðbmin þ bmax Þ=2c. Then, if M = bmin, the sequential searching is done in the range of (bmin, bmax) with the bmin + 1 as the starting point, until the filtered result is 0 0 an IMF at the exploratory point M , then let yi = M , ni = 0, and go to the next step; otherwise, take ni = bmin, return to step 2 and repeat.  Step 4. If yi  ni = 1, it indicates that the IMF with the maximum bandwidth is obtained. Then go to the next step. Otherwise, we take M ¼ bðyi þ niÞ=2c, and filter x(n) with the filter given in Eq. (32). It follows that y(n) is obtained. If y(n) is an IMF, let yi = M, otherwise take ni = M, then return to step 4 and repeat.  Step 5. Let i i + 1, note that ci(n) = y(n). Then let bmin = 0 and bmax = yi, and judge the signal y(n) obtained after applying analytical filtering to the signal x(n) with the low-pass filter as follows.

 Hd ðmÞ ¼

1 0 6 m < bmax 0

bmax 6 m < N

ð33Þ

If the energy of y(n) is little enough compared to x(n) (the energy ratio between y(n) and x(n) is calculated by using Eq. (30), and the threshold T2 is set as 0.05) or y(n) becomes a monotonic function, then we terminate the decomposition and note r(n) = y(n). Otherwise, we further judge whether y(n) is an IMF: if it is not, let M ¼ bðbmin þ bmax Þ=2c; ni ¼ 0 and return to step 2, otherwise, let i i + 1, note ci(n) = y(n), r(n) = 0, and terminate the decomposition. Since the exploratory point makes binary movement in the searching range, the algorithm is named as IMF binary searching (IBS). The flowchart of IBS algorithm is shown in Fig. 7. 3.2.4. Discussions It can be proven that both the algorithms are orthogonal decomposition methods. A detailed proof will be given in the Appendix. Furthermore, it is easy to note that the filtered result obtained by analytical band-pass filtering with Eq. (23) must be an IMF if mk = mk1  1 and X(mk) – 0, therefore N/2 IMFs can be obtained by ISS (IBS) at most

1320

Y. Qin et al. / Measurement 44 (2011) 1312–1327

Start bmin = 0 , bmax = N / 2 , i = 0

yes Is x(n) an IMF? no M = ⎣(bmin + bmax ) / 2⎦ , ni = 0

Construct H(m) with Eq. (32)

y(n) is obtained by filtering x(n) with H(m)

Calculate e with Eq. (30) ni = bmin

yes

M ← ⎣( M + bmin ) / 2⎦

e
no M =

no Is y(n) an IMF?

M = bmin ?

⎣(bmin + bmax ) / 2⎦

yes

yes

yi = M

Sequential searching yi = M ′

no yi-ni=1? M = ⎣( yi + ni ) / 2⎦

yes i ← i + 1 , ci ( n) = y ( n)

y(n) is obtained by filtering x(n) with H(m) given by Eq. (32) If y(n) is an IMF, let yi = M, otherwise let ni = M

bmin = 0 , bmax = yi

Construct Hd(m) with Eq. (33)

y(n) is obtained by filtering x(n) with Hd(m)

Calculate e with Eq. (30)

Is e
no Is y(n) an IMF?

no

yes

yes r ( n) = y ( n)

i ← i + 1 , c i ( n ) = y ( n ) , r ( n) = 0

c i ( n ) = y ( n ) , r ( n) = 0

End Fig. 7. The flowchart of IBS algorithm.

when the length of the discrete signal is N. It follows that ISS and IBS must be convergent as long as the signal is finite in the time domain. ISS and IBS are only approximate algorithms for OEMD because ideal filtering cannot be implemented. Thus, it is

necessary to discuss the error between ISS (IBS) and OEMD. The error is mainly generated by analytical band-pass filtering. Assume that c(n) is an arbitrary component of the signal s(n). In this case, C(m) and S(m) respectively represent the discrete spectra of c(n) and s(n), and the frequency

Y. Qin et al. / Measurement 44 (2011) 1312–1327

bin range of c(n) is defined as [fmin, fmax]. After applying ISS (IBS) to s(n), let the pass band of the filter which is used to get c(n) be [hmin, hmax]. Therefore, the estimated error of c(n) can be written as

^cðnÞ ¼

h max X

2p

½SðmÞ þ SðN  mÞej N mn

m¼hmin



hmax X

2p

½CðmÞ þ CðN  mÞej N mn

ð34Þ

m¼hmin

When the pass band of the filter almost accords with the frequency range of c(n) and spectral aliasing hardly exists, i.e. hmin  fmin, hmax  fmax, C(m)  S(m)(m e [fmin, fmax]), the estimated error of c(n) will be quite small. Unfortunately, if the spectrum of c(n) and other components lap over in a wide frequency range, both ISS and IBS can hardly separate c(n) from s(n), while other methods may also have similar pitfall. The exploratory point makes binary movement in the IBS algorithm, while the exploratory point makes sequential movement in the ISS algorithm, hence IBS can more quickly separate components from the signal than ISS. For a harmonic signal, since the spectrum of one component locates in a narrow frequency band, the error between the obtained component and the real component will be small if the used filter’s bandwidth includes the main frequency band of the real component. Therefore, IBS is more suitable for analyzing harmonic signals in order to increase the computing speed. However, since the IMF obtained by ISS have wider frequency band than that obtained by IBS, IBS may obtain more false IMFs than ISS, especially when the signal (such as the noisy signal with small signal-tonoise ratio (SNR)) has complex spectral constitution, as will be shown next. Consider a four-component signal

sðtÞ ¼ c1 ðtÞ þ c2 ðtÞ þ c3 ðtÞ þ nðtÞ

ð35Þ

with

c1 ðtÞ ¼ cos½2p  ð30t þ 10t 2 Þ c2 ðtÞ ¼ 0:7  ½1 þ 0:3sinð2p  5tÞ  cos½2p  110t þ 0:4sinð2p  10tÞ c3 ðtÞ ¼ 0:4  cos½2p  220t þ 0:2sinð2p  10tÞ

ð36Þ ð37Þ ð38Þ

where t e [0, 2], the sampling frequency fs = 500 Hz, and n(t) is a white noise. When the variance of n(t) is 0.1 (i.e. SNR is 19.0 dB), ISS and IBS are applied to this signal, respectively. Both ISS and IBS obtain 3 IMFs, while EMD (Rilling’s algorithm [17]) obtain 9 IMFs. The computation times of the two algorithms are measured, and the rootmean-square (RMS) error between the corresponding IMF and the ideal result is calculated for each component. The results are listed in Table 1. Another experiment is performed by modifying the variance of n(t). Let the variance of n(t) be 0.3 (i.e. SNR is 9.9 dB), ISS and IBS are applied to x(t) respectively (Both the two algorithms obtain 3 IMFs, while EMD obtain 10 IMFs), and the results of this experiment are listed in Table 2. If we increase the variance of n(t) further, the difference between ISS and IBS will be demonstrated. If we take the variance of n(t) as 0.7 (i.e.

1321

SNR is 2.1 dB), the numbers of IMFs obtained by ISS and IBS are 8 and 9, respectively; If we take the variance of n(t) as 1 (i.e. SNR is 0.5 dB), the numbers of IMFs obtained by ISS and IBS are 10 and 12, respectively. It can be seen from these experiments that the decomposition performances of the two algorithms are approximate when the SNR is small, while IBS obtains more false IMFs than ISS when the SNR is large, i.e. when the spectrum of the signal is complex. In order to better select the decomposition algorithm, a high pass FIR filter (the pass band is [fs/4, fs/ 2]) can be applied to the original signal, and the filtered signal is used for estimating the variance of the noise. Then, the estimation of the signal’s SNR can be calculated. If the SNR is small, the ISS algorithm should be selected, otherwise, the IBS algorithm should be selected. In addition, it is worth noting that both the two algorithms cannot correctly separate the components from the signal if the frequencies of some components are close to each other, like other common signal analysis methods (such as wavelet transform and EMD). This is one of the limitations in OEMD. 4. Comparisons and applications 4.1. Comparisons To demonstrate the advantages of OEMD on solving mode aliasing, a transient disturbing signal is considered. The signal is defined as

xðtÞ ¼ A1 sinð2pftÞ þ pðtÞA2 sinð2pmftÞ

ð39Þ

where A1 = 1, A2 = 0.4, m = 3, f = 50; p(t) = 1 when 0.425 6 t 6 0.575 and p(t) = 0 when 0.575 < t 6 1 or 0 < t < 0.425; the sampling frequency is 2000 Hz, and the number of sampling points is 2000. OEMD (IBS algorithm) is applied to decompose the signal x(t) and the result is shown in Fig. 8. The IMF component c1 which is shown in Fig. 9 is the hidden transient disturbing signal in x(t). Then according to the regularity that the amplitude must jump at the starting point and the ending point of the disturbing signal, we can accurately detect that the starting and the ending of the disturbance are respectively at the time instant of 0.425 s and 0.575 s. Furthermore, the frequency, amplitude and phase of the disturbing signal can be obtained based on c1. We also carry out EMD to x(t), and the result is shown in Fig. 10. As can be seen from Fig. 10, the phenomenon of mode aliasing occurs, and every IMF cannot demonstrate the real physical process. Thus it can be seen that EMD cannot be applied to decompose this type of signal, while OEMD can obtain correct analysis result. Consequently, OEMD can be effectively applied to singular detection. 4.2. Application into impulsive feature extraction To illustrate potentials of the proposed decomposition method in the engineering field, a vibration displacement signal sampled from a rotor experiment bench (single disc rotor system) is used for analysis. The rotary frequency of the rotor is 20 Hz, the sampling frequency is 500 Hz, and

1322

Y. Qin et al. / Measurement 44 (2011) 1312–1327

Table 1 Comparison of two algorithms ISS and IBS when the SNR is 19.0 dB. Algorithm

RMS error of c1

RMS error of c2

RMS error of c3

Computation time

ISS IBS

0.0734 0.0734

0.1311 0.1311

0.1344 0.1344

0.5996 0.1720

Table 2 Comparison of two algorithms ISS and IBS when the SNR is 9.9 dB. Algorithm

RMS error of c1

RMS error of c2

RMS error of c3

Computation time

ISS IBS

0.2064 0.1874

0.2151 0.2317

0.2175 0.2173

0.3895 0.1142

Fig. 8. Detecting the transient disturbing signal with OEMD, where x3(t) is the original signal.

Fig. 9. The transient disturbing signal obtained by OEMD.

the number of sampling points is 500. The time domain waveform of the signal is shown in Fig. 11. We make use of a rubbing screw to impact the rotor. A control circuit is used for making the rubbing screw move up and down, in order that the rotor suffers from one impulse per revolution. A displacement sensor is placed on the casing of the roller bearing to pick up vibration signals. Fig. 11 shows the sampled vibration displacement signal. It is hard to see the impulsive signal from Fig. 11. Therefore, OEMD (IBS algorithm) is applied to decompose the vibration displacement signal and two IMFs, c1 and c2, and a residue r are obtained, which are shown in Fig. 12. From the

figure it can be seen that the component c1 contains typical periodic impulsive signal whose peak value is quite clear, while the component c2 corresponds to the periodic vibration signal of the rotor. According to c1 as shown in Fig. 12, we can easily know the time instants each impulse occurs (0.04 s, 0.09 s, 0.14 s, . . .), i.e., the time interval between two impulses is 0.05 s. Obviously, the impulsive frequency (20 Hz) accords with the fact. In order to demonstrate the advantage of the proposed method, we also decompose the original signal with EMD and get five components and a residue, which is shown in Fig. 13. It is easy to note that there are some periodical

Y. Qin et al. / Measurement 44 (2011) 1312–1327

1323

Fig. 10. Detecting the transient disturbing signal with EMD, where x3(t) is the original signal.

Fig. 11. Vibration displacement signal of a rotor which is suffered from a series of impulses.

Fig. 12. OEMD result of the signal shown in Fig. 11.

impulses in the first component. Compared Fig. 12 with Fig. 13, we can see that the impulsive periodicity of the characteristic impulsive signal obtained by OEMD is clearer than that of the characteristic impulsive signal ob-

tained by EMD. In addition, there are some meaningless components in the results of EMD. Hence it indicates that OEMD is superior to EMD for analyzing the impulsive vibration signal.

1324

Y. Qin et al. / Measurement 44 (2011) 1312–1327

Fig. 13. EMD result of the signal shown in Fig. 11.

4.3. Application into gearbox fault diagnosis In order to further demonstrate the efficiency of this proposed approach for the multicomponent AM-FM signal, OEMD is applied to a vibration signal of a faulty secondary gearbox as shown in Fig. 14. An accelerometer which is placed on the shell of the gearbox is used for sampling the vibration signal. The number of the gear’s teeth is 30 at the gearbox input shaft, the numbers of the two gears’ teeth are 69 and 18 respectively at the gearbox intermediate shaft, and the number of the gear’s teeth is 81 at the gearbox output shaft. Hence the transmission ratio is 10.35. The gear at the gearbox input shaft has two broken teeth, the rotary frequency of the input shaft (fr) is 15 Hz and the sampling frequency is 10 kHz. The vibration acceleration signal is decomposed by OEMD (ISS algorithm) and 11 IMFs are obtained in all. For a comparison of IBS and ISS, IBS is applied to the same signal and 19 IMFs are obtained. It follows that ISS is more beneficial to analyzing complex signals than IBS. Fig. 15 shows the first IMF obtained by ISS. From the figure it can be seen that the first IMF is an AM-FM signal, therefore

envelope analysis should be applied to it in order to extract the fault feature. It is worth noting that the analytical signal of the first IMF has been obtained after OEMD, and the spectrum analysis is applied to its amplitude envelope so as to obtain the envelope spectrum which is shown in Fig. 16. From Fig. 16 it can be seen that there are obvious spectral lines at the characteristic frequency fr (15 Hz) of the fault gearbox and its double frequency (30 Hz), which shows that the gear at the input shaft has local fault. This is consistent with the fact. In order to better demonstrate the validity of the proposed method, the conventional envelope analysis method based on Hilbert transform is applied to the same signal. Fig. 17 shows the envelope spectrum of the fault vibration signal as shown in Fig. 15 by the conventional Hilbert envelope analytic method with band filtering (the pass band is [1000, 2000] Hz). From Fig. 17 it can be seen that there is no spectrum line at the characteristic frequency fr (15 Hz) of the fault gearbox. If the pass band is [2000, 3000] Hz, the envelope spectrum of the same signal is shown in Fig. 18, from which it can be seen that two obvious spectral

Fig. 14. Vibration acceleration signal of a secondary gearbox with a broken teeth fault.

Y. Qin et al. / Measurement 44 (2011) 1312–1327

1325

Fig. 15. The first IMF of the signal shown in Fig. 14 obtained by OEMD.

Fig. 16. The envelope spectrum of the IMF shown in Fig. 15.

Fig. 17. The envelope spectrum obtained by the conventional Hilbert envelope analysis method with band filtering (the pass band is [1000, 2000] Hz).

lines can be found at the characteristic frequency fr (15 Hz) and its double frequency. Therefore, for the conventional envelope analysis method, the selection of central frequency and bandwidth of band-pass filter

brings great subjectivity that may influence the accuracy of diagnosis result, whereas OEMD method can overcome the limitation which occurs in the conventional envelope analysis method.

Fig. 18. The envelope spectrum obtained by the conventional Hilbert envelope analysis method with band filtering (the pass band is [2000, 3000] Hz).

1326

Y. Qin et al. / Measurement 44 (2011) 1312–1327

According to the process of ISS (IBS), we can write mj 6 mk1. Hence, we have

5. Conclusions OEMD is a new method for decomposing the signal into IMFs adaptively, which can be applied to multicomponent signal decomposition. The proposed method has strict orthogonality and completeness and solves some major disadvantages of EMD. The basic principle of this method is to obtain each IMF and the residue with self-adaptive filtering. Since the orthogonal basis function is not compactly supported and the filtered results should preserve ideal band-pass characteristics in the frequency domain for OEMD, analytical band-pass filtering is presented. Then its perfect filtering performance in the frequency domain gets proved, whereas the spectrum aliasing may bring error to the filtered result. Based on analytical band-pass filtering, this paper proposes two fast algorithms for implementing OEMD, i.e. IMF sequential searching (ISS) and IMF binary searching (IBS), and the two algorithms are compared. Taking advantage of these algorithms, not only the IMFs but also their Hilbert spectra can be obtained fast, which is very beneficial for the later calculation of Hilbert spectrum. Applications of this proposed method to synthetic signals and mechanical vibration signals show its effectiveness for various typical signals. Compared with EMD, OEMD well solves some problems in EMD, however, it still has bad performance when the spectra of some components lap over in a wide frequency range like other methods. Moreover, if the frequencies of two signals are too close to each other, the compositive signal can be regarded as an IMF, therefore OEMD cannot correctly separate the two components from such signal. In despite of these limitations, it can be seen from the above application examples that OEMD has a high performance for some typical mechanical vibration signals.

N1 X

C j ðmÞC k ðmÞ ¼ 0

m¼0

Using Parseval theorem, we obtain N1 X

cj ðnÞck ðnÞ ¼

n¼0

N1 X

C j ðmÞC k ðmÞ ¼ 0

ðA:1Þ

m¼0

Eq. (A.1) shows that two arbitrary IMFs obtained by ISS (IBS) are orthogonal. Then let us examine the orthogonality between the residue r(t) and any IMF cj(t). Assuming L IMFs are obtained after performing ISS (IBS) on x(n). From (24), the discrete spectrum of r(n) is given by

8 > < XðmÞ 0 6 m < mL RðmÞ ¼ 0 mL 6 m 6 N  mL > : XðmÞ N  mL < m < N Since mj P mL, we get N1 X

RðmÞC j ðmÞ ¼ 0

m¼0

Similarly, by means of Parseval theorem, we can write N1 X

rðnÞcj ðnÞ ¼

n¼0

N1 X

RðmÞC j ðmÞ ¼ 0

ðA:2Þ

m¼0

It follows from Eqs. (A.1) and (A.2) that both ISS and IBS are orthogonal decomposition methods. Moreover, from the decomposition rule of ISS (IBS) and the spectra of IMFs and residue, we can obtain

XðmÞ ¼

L X

C j ðmÞ þ RðmÞ

j¼1

Acknowledgment This work was supported by the Fundamental Research Funds for the Central Universities (No. CDJZR10 28 00 07).

Equivalently,

xðnÞ ¼

L X

cj ðnÞ þ rðnÞ

ðA:3Þ

j¼1

Appendix A Let x(n) (n = 0, 1, . . . , N  1) be a discrete signal. cj(n) and ck(n), (j < k) are two arbitrary IMFS obtained by performing ISS (IBS) on the signal. From Eq. (23), the discrete spectra of cj(n) and ck(n) can be given by

C j ðmÞ ¼

8 0 0 6 m < mj > > > > > > < XðmÞ mj 6 m < mj1

0 mj1 6 m 6 N  mj1 > > > XðmÞ N  mj1 < m 6 N  mj > > > : 0 N  mj < m < N

8 0 > > > > > > < XðmÞ C k ðmÞ ¼ 0 > > > > XðmÞ > > : 0

0 6 m < mk mk 6 m < mk1 mk1 6 m 6 N  mk1 N  mk1 < m 6 N  mk N  mk < m < N

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] N.E. Huang, M.L. Wu, S.R. Long, A confidence limit for the empirical mode decomposition and Hilbert spectral analysis, Proceedings of the Royal Society of London Series A 459 (2003) 2317–2345. [3] 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. [4] M.B. Velasco, B. Weng, K.E. Barner, ECG signal denoising and baseline wander correction based on the empirical mode decomposition, Computers in Biology and Medicine 38 (1) (2008) 1–13. [5] J.W. Essex, A.P. Wiley, Application of the Hilbert–Huang transform to the analysis of molecular dynamics simulations, Journal of Physical Chemistry 107 (2003) 4869–4876. [6] J. Chen, Y.L. Xu, R.C. Zhang, Modal parameter identification of Tsing Ma suspension bridge under Typhoon Victor: EMD-HT method, Journal of Wind Engineering and Industrial Aerodynamics 92 (10) (2004) 805–827.

Y. Qin et al. / Measurement 44 (2011) 1312–1327 [7] D. Ghosh, A.R. Chowdhury, P. Saha, On the various kinds of synchronization in delayed Duffing–Van der Pol system, Communications in Nonlinear Science and Numerical Simulation 13 (4) (2008) 790–803. [8] T.Y. Li, Y. Zhao, N. Li, G. Fen, H.H. Gao, A new method for power quality detection based on HHT, Proceedings of the CSEE 25 (17) (2005) 52–56. [9] Z.K. Peng, P.W. Tse, F.L. Chu, A comparison study of improved Hilbert–Huang transform and wavelet transform: application to fault diagnosis for rolling bearing, Mechanical System and Signal Processing 19 (5) (2005) 974–988. [10] S.R. Qin, Y.M. Zhong, A new envelope algorithm of Hilbert–Huang transform, Mechanical System and Signal Processing 20 (8) (2006) 1941–1952. [11] J.S. Cheng, D.J. Yu, Y. Yang, Application of support vector regression machines to the processing of end effects of Hilbert–Huang transform, Mechanical Systems and Signal Processing 21 (2007) 1197–1211.

1327

[12] H.L. Li, L.H. Yang, D.R. Huang, The study of the intermittency test filtering character of Hilbert–Huang transform, Mathematics and Computers in Simulation 70 (1) (2005) 22–32. [13] Y. Qin, S.R. Qin, Y.F. Mao, Research on iterated Hilbert transform and its application in mechanical fault diagnosis, Mechanical Systems and Signal Processing 22 (8) (2008) 1967–1980. [14] P. Flandrin, G. Rilling, P. Goncalves, Empirical mode decomposition as a filter bank, IEEE Signal Processing Letters 11 (2) (2004) 112– 114. [15] Q.S. Cheng, L.W. Wu, The empirical mode frequency decomposition of time series analysis, Mathematics in Practice and Theory 36 (5) (2006) 151–153. [16] S.L. Hahn, Hilbert Transforms in Signal Processing, Artech House, Boston, 1996. [17] G. Rilling, P. Flandrin, On empirical mode decomposition and its algorithms, IEEE-EURASIP Workshop Nonlinear Signal Image Process, Grado, Italy, 2003.