Shock and Vibration
Volume 2015 (2015), Article ID 425989, 9 pages
http://dx.doi.org/10.1155/2015/425989

Research Article

Feature Extraction of Faulty Rolling Element Bearing under Variable Rotational Speed and Gear Interferences Conditions

Dezun Zhao, Jianyong Li, and Weidong Cheng

School of Mechanical, Electronic and Control Engineering, Beijng Jiaotong University, Beijing 100044, China

Received 6 May 2015; Revised 15 July 2015; Accepted 26 July 2015

Academic Editor: Wahyu Caesarendra

Copyright © 2015 Dezun Zhao et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

In the field of rolling element bearing fault diagnosis, variable rotational speed and gear noise are main obstacles. Even though some effective algorithms have been proposed to solve the problems, their process is complicated and they may not work well without auxiliary equipment. So we proposed a method of faulty bearing feature extraction based on Instantaneous Dominant Meshing Multiply (IDMM) and Empirical Mode Decomposition (EMD). The new method mainly consists of three parts. Firstly, IDMM is extracted from time-frequency representation of original signal by peak searching algorithm, which can be used to substitute the bearing rotational frequency. Secondly, resampled signal is obtained by an IDMM-based resampling algorithm; then it is decomposed into a number of Intrinsic Mode Functions (IMFs) based on the EMD algorithm. Calculate kurtosis values of IMFs and an appropriate IMF with biggest kurtosis value is selected. Thirdly, the selected IMF is analyzed with envelope demodulation method which can describe the fault type of bearing. The effectiveness of the proposed method has been demonstrated by both simulated and experimental mixed signals which contain bearing and gear vibration signal.

1. Introduction

There are two main obstacles for rolling bearing diagnosis; they are variable rotational speed and gear vibration noise [1]. On one hand, time-varying speed will directly influence frequency and amplitude of rolling element bearing vibration signal. In terms of frequency, the fault impulse will not recur at a constant interval; hence demodulated resonance technique and any other methods based on the premise of constant speed are no longer applicable. Furthermore, the amplitude of the bearing vibration signal changes with the rotational speed in time domain. On the other hand, the amplitude of gear meshing noise is higher than bearing in rotating machinery. Even though the gear does not have fault, impulses of faulty bearing will be contaminated by gear noise, and the extraction of resonance frequency excited by bearing fault also will be influenced.

Many efforts have focused on the study of time-varying rotational speed and gear vibration noise. Order tracking [2, 3] is one of the most effective methods to deal with variable speed signal. The essence of the order tracking is resampling the original vibration signal at a constant angle increment which converts the nonstationary signal (in time domain) into the stationary one (in angular domain) to remove the effects of speed fluctuation. Order tracking algorithm can be divided into three types: hardware-based order tracking, computed order tracking [2–4], and order tracking based on instantaneous frequency estimation [5, 6]. Compared with hardware-based order tracking and computed order tracking, the order tracking method based on instantaneous frequency estimation can overcome the limits of auxiliary equipment and costs. Hence it has become a new focus for many researchers, and its principle is extracting rotational frequency from original vibration signal then angular domain resampling. However, because of gear noise interference, it is difficult to directly extract the bearing rotational frequency from the mixed vibration signal. Traditional methods to eliminate the interference of the gear noise mainly include Time Synchronous Average (TSA) [7], Linear Prediction [8], Adaptive Noise Cancellation (ANC) [9], Self-Adaptive Noise Cancellation (SANC) [10], Discrete/Random Separation (DRS) [11], and Edited Cepstrum [12]. These algorithms require rotational speed to keep stable except the ANC algorithm. ANC algorithm requires auxiliary equipment to obtain reference signals of gear, and the installation position of the auxiliary equipment is very important.

For faulty rolling element bearing feature extraction under time-varying rotational speed and gear noise, Borghesani et al. [13] proposed a method. This method consists of three steps: removing the effects of speed fluctuation based on the order tracking algorithm, eliminating gear noise by traditional denoising method in the gear angular domain, returning the denoised signal to time domain based on the reverse order tracking and angular resampling by bearing rotational frequency, and, at last, getting envelope order spectrum to determine bearing fault type. On the other hand, Wang [14] presented an autoregressive model-based method to diagnose faulty bearing under time-varying speed and gear noise conditions. The methods proposed by Borghesani and Wang effectively solve the problem of bearing fault diagnosis under time-varying rotational speed and gear noise, but these methods need speed information of bearing and gear which is obtained by auxiliary equipment. Hence it may not work in some conditions that have restrictions of installation space and costs. Afterwards, Wang [15, 16] proposed a new concept, named Instantaneous Dominant Meshing Multiply (IDMM), which can be directly extracted from mixed signal. Based on IDMM, a function is constructed as the reference signal for the ANC algorithm which is used to remove gear noise. Signal resampling based on instantaneous fault characteristic frequency (IFCF) which is extracted from time-frequency representation of remaining bearing signal component, then, obtains fault characteristic order (FCO) spectrum to determine bearing fault type. Although the effectiveness of the improved ANC algorithm has been validated, the process is very complex, because we need to try many times to obtain optimal parameters of the ANC algorithm, and the process must rely on FCO spectrum. Even though we have optimal parameters of the ANC algorithm, the multiples of the gear meshing frequency except for the IDMM cannot be removed. On the other hand, IFCF has errors which will influence bearing fault diagnosis.

In view of above analysis, a new method based on IDMM and EMD is proposed. The proposed method does not need auxiliary equipment to measure speed information, avoids complex parameters selecting, and simplifies the process of rolling element bearing fault diagnosis. It has three main steps: (a) extracting IDMM trend from time-frequency representation of mixed signal and signal resampling via the IDMM-based resampling algorithm, (b) obtaining IMFs by the EMD algorithm and selecting suitable IMF by kurtosis values, and (c) obtaining envelope order spectrum to diagnose rolling element bearing.

2. Theoretical Background

2.1. IDMM Trend

In the previous studies, gear noise has been generally classified as narrowband periodic signal. Due to this characteristic, we can easily find gear meshing frequency and its multiples. The amplitude of the dominant meshing multiple has fixed proportion with gear rotational frequency, and it can be represented as

where is a number of the collection and () is the meshing multiple number; the value of is different in different working conditions; in the most general case, it is 2 or 3; is the gear meshing frequency; is the number of gear teeth; is the rotational speed of gear whose unit is r/min.

If gear rotational frequency changes with time, Instantaneous Dominant Meshing Multiple will change in proportion to it. The IDMM is defined as an instantaneous meshing frequency harmonic that has the highest amplitude among all harmonics of the instantaneous meshing frequency [15], and the IDMM trend is a collection of IDMM. IDMM can be easily extracted from time-frequency representation of the mixed signal, because it has the highest amplitude among all harmonics of the meshing frequency and the rolling element bearing signal component has relatively weaker one. Short-Time Fourier Transform (STFT) is used to compute time-frequency representation of mixed signal in this paper, and the IDMM is calculated as follows:

where is the time; is the frequency; is a Gaussian window function; represents instantaneous frequency spectrum; is the sequence number of instantaneous frequency spectrum whose value is from 1 to ; the number of instantaneous frequency spectrum is determined by the length of Gaussian window function; arg max function is used to find frequency domain abscissa of when has the largest value.

2.2. EMD Algorithm

EMD is a self-adaptive analysis method proposed by Huang et al. [17] in 1998. It is based on the local characteristic time scale of signal and decomposes the signal into a number of narrowband components, called IMFs, that satisfies two conditions:

(a)In the whole data set, the number of extrema and the number of zero-crossings must either be equal or differ at most by one.(b)At any point, the mean value of the envelope defined by the local maxima and the envelope defined by the local minima is zero.

The process of EMD is shown as follows:

(1)Obtain the local maxima and minima of the signal . Interpolate, using natural cubic splines, along the points of local maxima and minima in order to form an upper and a lower envelope.(2)The mean of the two envelopes is designated as . , the difference between and , is calculated as

 If satisfies all the requirements of IMF, it is the first IMF of .(3)If not, is treated as the original signal and processed by previous steps. The mean of the two envelopes is designated by . Let . If is not an IMF, we can repeat the cycle until is an IMF; that is,

 Then we define ; it is the first IMF of the original signal.(4)Separate from as follows:

The residue, , is treated as the original signal and subjected to the same cyclic process as described above. Repeat cycle times and get IMFs as follows:

Stop cyclic process until becomes a monotonic function from which no more IMF can be extracted. Via EMD decomposition, original signal can be represented as

where is a residue, which represents the signal’s central tendency.

2.3. IDMM-Based Resampling Algorithm

We obtain resampled signal by the IDMM-based resampling algorithm. Traditional resampling algorithm is firstly to fit speed curve, secondly to calculate function relationship between the angle of reference shaft and the time , and thirdly to obtain resampled signal via interpolation algorithm. However, calculated burden of the traditional algorithm is too heavy, and its efficiency is quite low, as it needs to solve a number of quadratic equations. We use the IDMM-based resampling algorithm to resample the signal. The essence of the algorithm is to divide signal into a number of parts, larger speed parts with larger sampling rate and lower speed parts with lower one, and the ratio of sampling rate of two different time segments is equal to the ratio of rotational frequency. The process of IDMM-based resampling algorithm is shown as follows:

(1)Divide mixed signal into ( is 200 in this paper) parts, named , and every part of signal has a rotational frequency: .(2)Define the initial sample rate, namely,

 where is sampling rate of the original signal. The initial IDMM is calculated as

(3)Based on the ratios of rotational frequency of , reset sampling rate of is determined as

(4)Obtain resampled signal based on the computed sampling rates, . Firstly, determine the new sampling point based on new sampling rate, ; secondly, based on original signal , determine resampled signal by polynomial interpolation; lastly, the s are arranged by the value of , and the resampled signal is shown as .

3. Feature Extraction of Faulty Bearing Based on IDMM and EMD

Instantaneous Dominant Meshing Multiply (IDMM) is easily extracted from time-frequency representation of the mixed signal by peak searching algorithm. In addition, IDMM trend is proportional to bearing rotational frequency. Hence we proposed resampling the signal based on IDMM trend. EMD is a self-adaptive algorithm, which decomposes signal by local characteristic time scale of the signal. Furthermore, IMFs are distributed from high to low in turn according to their frequency. Bearing signal is high frequency resonance signal, and gear meshing signal generally has lower frequency. Hence we use the EMD algorithm to remove interference of gear noise. In the proposed method, the IMF with the biggest Kurtosis value is selected for further analysis, which means that it includes much more fault information than other IMFs. The Kurtosis value of a random variable is defined as

where represents mathematical expectation of a random variable and and are standard deviation and mean value of , respectively.

The steps of the proposed method can be given as follows. A detailed flowchart is presented in Figure 1.

(1)Extract IDMM trend from time-frequency representation (TFR) of mixed signal by peak searching algorithm. Then obtain resampled signal by the IDMM-based resampling algorithm.(2)Decompose signal by EMD algorithm and obtain IMFs. Compute Kurtosis value of IMFs based on (11) and select an appropriate IMF.(3)Obtain envelope order spectrum of the selected IMF to determine bearing fault type.

Figure 1: Flowchart of rolling element bearing faulty feature extraction based on IDMM and EMD.

4. Simulation and Experimental Validation

4.1. Simulation Study

In this subsection, a synthetic speed-up signal, which includes faulty rolling element bearing, gear, and the white Gaussian noise, is constructed to demonstrate the effectiveness of the proposed method.

The model representing the faulty rolling element bearing signal under speed-up, , is given by

where is the amplitude of the th fault impulse, is the structural damping frequency, is the fault exiting resonance frequency, is a unit step function, and is the occurrence time of the th impulse and it can be calculated as

where , = 0.01~0.02, is the errors caused by the random slippage of the rolling elements; is the rotational frequency function of rolling element bearing and ; and is the number of fault impacts per shaft revolution.

The model representing the gear signal under speed-up can be written as

where () is the meshing multiple number and is the amplitude of th harmonic frequency. is the tooth number. is integral function of the gear rotational frequency, . is the ratio of gear rotational frequency to bearing; namely, . The amplitudes of gear meshing frequency and amplitudes of its 2nd and 3rd multiples are set as 5 and 0.6, which means that the 2nd multiple is the dominant meshing multiple. The other parameters in the model are listed in Table 1.

Table 1: Parameters of mixed bearing and gear simulation model under time-varying rotational speed.

Parameters
Value

Length of the signal
8 s

Length of
5 Hz–25 Hz

Length of
6 Hz–30 Hz

Sampling rate
20000 Hz

Fault characteristic coefficient
3

Bearing resonance frequency
5000 Hz

Number of teeth L
18

Rotating frequency ratio
1.2

The simulation mixed signal can be represented as

Four seconds of the simulation signal is selected for analysis, and the selected simulation signal is shown in Figure 2(a). The amplitude of the mixed signal goes up in proportion to the rotational frequency and the fault impulses are not obvious. Time-frequency representation of the mixed signal is shown in Figure 2(b). From Figure 2(b), the resonance frequency can be found at around 5000 Hz. The three linearly rising lines in the lower part of Figure 2(b) are the three multiples of the gear meshing frequency, and 2nd meshing frequency is quiet obvious.

(a)

(b)

Figure 2: (a) Mixed simulation signal in time domain. (b) Time-frequency representation of the mixed signal.

IDMM trend is extracted using peak searching algorithm from time-frequency representation of the mixed simulation signal, and it is shown in blue solid line in Figure 3, together with the preset gear meshing frequency trend (red solid line) and its 2nd multiple (red dash line). The extracted IDMM trend and the 2nd multiple of the preset gear meshing frequency are almost identical. Hence we can conclude that IDMM trend can be easily extracted. Because IDMM trend has invariable proportion (the proportion coefficient is 43.2, which is calculated by (1) and parameters of the simulation model) to bearing rotational frequency, it can be used to substitute bearing rotational frequency.

Figure 3: Extracted IDMM trend and preset gear meshing frequency.

We obtain the resampled signal by the IDMM-based resampling algorithm. The angular resampled signal is shown in Figure 4(a). Figure 4(b) is the time-frequency representation of the angular resampled signal. The EMD algorithm is used to process the resampled signal, and the decomposed first six IMFs are shown in Figure 5. We respectively calculate Kurtosis values of IMFs and Kurtosis values of IMFs are listed in Table 2. IMF1 is selected for further analysis as the biggest kurtosis value.

Table 2: Kurtosis values of IMFs.

IMF
1
2
3
4
5
6

Kurtosis value
16.91
1.84
6.37
5.18
1.76
1.45

(a)

(b)

Figure 4: (a) Angular resampled signal. (b) Time-frequency representation of angular resampled signal.

Figure 5: Decomposition result of the resampled signal.

Figure 6 is time-frequency representation of the IMF1.Compared to Figure 4(b), the results revealed that the three multiples of the gear meshing frequency in the mixed signal have been removed by EMD algorithm. Apply envelope demodulation method to process the IMF1, and the corresponding envelope order spectrum is shown in Figure 7(a). In this picture, fault characteristic order (FCO) and its 2nd and 3rd harmonics are easily identified whose values are 0.0696, 0.139, and 0.2086, respectively. Bearing Rotational Order (BRO) and its 2nd harmonic are also easily identified whose values are 0.0232 and 0.0462. The ratio of the FCO and the BRO is 3, which is equal to the preset fault characteristic coefficient. Hence the fault type of rolling element bearing can be identified.

Figure 6: Time-frequency representation of IMF1.

(a)

(b)

Figure 7: (a) Envelope order spectrum of IMF1. (b) Envelope order spectrum of resampled signal.

For the purpose of comparison, the resampled signal is directly analyzed with envelope demodulation method, and the corresponding envelope order spectrum is shown in Figure 7(b), from which the gear interference peaks are clearly revealed. Figures 7(a) and 7(b) indicated that the proposed method can effectively remove interference of gear meshing frequency. Furthermore, the amplitude of each harmonic is obviously increased. The changed ratio of the FCO amplitude between the resampled signal and the denoising signal is 113%; other amplitudes of harmonics (see column 2 and 3 of Table 3) and changed ratios (see column 4 of Table 3) are given in Table 3, respectively.

Table 3: Comparison of the envelope order spectrum.

Harmonic
order
Amplitudes of the harmonics
Changed ratio

Resampled signal
Denoising signal

BRO
0.035
0.036
2.9%

2BRO
0.019
0.027
42.1%

FCO
0.116
0.247
113%

2FCO
0.119
0.251
111%

3FCO
0.100
0.186
86%

4FCO
0.0787
0.125
58.8%

4.2. Experimental Analyses

The effectiveness of the proposed method is further verified by the mixed signal, which is measured on the vibration test rig. The layout of the test rig is illustrated in Figure 8. The gearbox is driven by an electrical motor, and the speed of the motor is controlled by a speed controller. A tachometer is fixed on the output shaft to measure speed information. The sampling rate is 24,000 Hz. The signal is collected by an acceleration sensor which is mounted on the top surface of the faulty bearing, YE6231 acquisition card, and its corresponding software. Gears do not have fault, and their reduction ratio is 3/11. A simulated outer race faulty bearing, 6206, is assembled in the forepart of the input shaft. The parameters of faulty rolling element bearing and fault characteristic coefficient of outer race () are shown in Table 4.

Table 4: Parameters of bearing in experiment.

Bearing type
Number of balls
Ball diameter d/mm
Pitch diameter D/mm
Contact angle

6206
9
9.52
51
0
3.66

Figure 8: Layout of the test rig.

The measured signal in time domain is shown in Figure 9(a), and the length of the signal is 1 second. Figure 9(b) is time-frequency representation of the measured signal. From the picture we can easily find the IDMM trend. The IDMM trend is extracted from Figure 9(b) using the peak searching algorithm. Figure 10 displays the IDMM trend (blue solid line) and the measured gear meshing frequency (red solid line) and its 3rd multiple trend (red dash line) for comparison. The measured gear meshing frequency is calculated by the measured speed information and (1). As shown in Figure 10, the IDMM trend and the 3rd multiple of gear meshing frequency almost overlap with each other. The proportional relationship between the IDMM trend and the gear meshing frequency is invariable, which is calculated by the parameters of the gearbox. Therefore, the IDMM trend can be reliably used for signal resampling.

(a)

(b)

Figure 9: (a) Measured signal in time domain. (b) Time-frequency representation of the measured signal.

Figure 10: Extracted IDMM trend and measured gear meshing frequency.

Resampled signal is obtained by an IDMM-based resampling algorithm, and the EMD algorithm is used to decompose the resampled signal. The IMFs are shown in Figure 11, respectively. Kurtosis values of IMFs are 25.34, 1.98, 4.18, 4.6, 3.31, and 1.12, respectively. Because IMF1 has the biggest Kurtosis value, we select IMF1 for analysis. The envelope order spectrum of IMF1 is shown in Figure 12(a). As shown in Figure 12(a), the FCO and its harmonics are obvious; we can also easily find the BRO and its 2nd harmonic. Values of the FCO and the BRO are 0.0786 and 0.0213, respectively. The calculated fault characteristic coefficient is 3.62, and the actual value is 3.66; hence we can determine that the fault is on the outer race of the rolling element bearing.

Figure 11: Decomposition result of the resampled signal.

(a)

(b)

Figure 12: (a) Envelope order spectrum of IMF1. (b) Envelope order spectrum of resampled signal.

We directly analyze resampled signal by envelope demodulation method, and the envelope order spectrum is shown in Figure 12(b). There are obvious interference peaks, and the BRO cannot be found in the picture. Compared to Figure 12(b), the BRO and its 2nd harmonic are easily recognized in the envelope order spectrum of IMF1. In addition, the FCO and the BRO have higher amplitude. For easy comparison, amplitude of each harmonic and the changed ratios of harmonics between the resampled signal and the denoising signal are listed in Table 5.

Table 5: Comparison of the envelope order spectrum.

Harmonic
order
Amplitudes of the harmonics
Changed ratio

Resampled signal
Denoising signal

BRO

0.017

2BRO

0.034

FCO
0.065
0.092
41.5%

2FCO
0.054
0.072
33.3%

3FCO
0.045
0.064
42.2%

4FCO
0.035
0.049
40%

5. Conclusions

A novel method of rolling element bearing fault feature extraction based on IDMM and EMD under variable rotational speed and gear interferences conditions is proposed. The major innovative features of the proposed method are reflected by the following:

(a)IDMM trend is used to substitute bearing rotational frequency, and the proposed method does not need auxiliary equipment to measure speed information.(b)The EMD algorithm is applied to remove gear noise, and the results are satisfying.(c)The proposed method avoids selecting complex parameters and simplifies the process of rolling element bearing diagnosis under time-varying rotational speed and gear interferences conditions.

Simulation and experiment analysis have demonstrated that the proposed method is effective in rolling element bearing fault diagnosis under time-varying rotational speed and gear noise. In addition, compared with the envelope order spectrum of the resampled signal, the BRO and the FCO in the envelope order spectrum obtained by the proposed method have higher amplitude.

It should be pointed out that the proposed method is only suitable for some mechanism whose bearing rotational frequency is proportional to gear rotational frequency, such as gearbox.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgment

This study is supported by the National Natural Science Foundation of China under Grant no. 51275030. The support is greatly appreciated.

References

T.-Y. Wang, J.-Y. Li, and W.-D. Cheng, “Fault diagnosis of rolling bearing under a variable rotational speed and gear vibration noise based on revised ANC algorithm and FCO spectrum,” Journal of Vibration and Shock, vol. 33, no. 18, pp. 7–13, 2014.

K. R. Fyfe and E. D. S. Munck, “Analysis of computed order tracking,” Mechanical Systems and Signal Processing, vol. 11, no. 2, pp. 187–202, 1997.

W. Cheng, R. X. Gao, J. Wang, T. Wang, W. Wen, and J. Li, “Envelope deformation in computed order tracking and error in order analysis,” Mechanical Systems and Signal Processing, vol. 48, no. 1-2, pp. 92–102, 2014.

P. N. Saavedra and C. G. Rodriguez, “Accurate assessment of computed order tracking,” Shock and Vibration, vol. 13, no. 1, pp. 13–32, 2006.

J. Urbanek, T. Barszcz, and J. Antoni, “A two-step procedure for estimation of instantaneous rotational speed with large fluctuations,” Mechanical Systems and Signal Processing, vol. 38, no. 1, pp. 96–102, 2013.

G. Yu, Q. Shuren, T. Baoping, and et al, “Order tracking of rotating machinery based on instantaneous frequency estimation,” Chinese Journal of Mechanical Engineering, vol. 39, no. 3, pp. 32–36, 2003.

P. D. McFadden and M. M. Toozhy, “Application of synchronous averaging to vibration monitoring of rolling element bearings,” Mechanical Systems and Signal Processing, vol. 14, no. 6, pp. 891–906, 2000.

R. B. Randall and J. Antoni, “Rolling element bearing diagnostics—a tutorial,” Mechanical Systems and Signal Processing, vol. 25, no. 2, pp. 485–520, 2011.

Y. Wang R, Bearing fault detection and oil debris monitoring by adaptive noise cancellation [M.S. thesis], University of Ottawa, Ottawa, Canada, 2008.

D. Ho and R. B. Randall, “Effects of time delay, order of FIR filter and convergence factor on self-adaptive noise cancellation,” in Proceedings of the 5th International Congress on Sound and Vibration, Adelaide, Australia, 1997.

P. Borghesani, P. Pennacchi, R. B. Randall, and R. Ricci, “Order tracking for discrete-random separation in variable speed conditions,” Mechanical Systems and Signal Processing, vol. 30, pp. 1–22, 2012.

R. B. Randall and N. Sawalhi, “A new method for separating discrete components from a signal,” The Sound and Vibration, vol. 45, no. 5, pp. 6–9, 2011.

P. Borghesani, R. Ricci, S. Chatterton, and P. Pennacchi, “A new procedure for using envelope analysis for rolling element bearing diagnostics in variable operating conditions,” Mechanical Systems and Signal Processing, vol. 38, no. 1, pp. 23–35, 2013.

W. Wang, “Autoregressive model-based diagnostics for gears and bearings,” Insight: Non-Destructive Testing and Condition Monitoring, vol. 50, no. 8, pp. 414–418, 2008.

T. Wang, M. Liang, J. Li, W. Cheng, and C. Li, “Bearing fault diagnosis under unknown variable speed via gear noise cancellation and rotational order sideband identification,” Mechanical System and Signal Processing, vol. 62-63, pp. 30–53, 2015.

T. Wang, M. Liang, J. Li, and W. Cheng, “Rolling element bearing fault diagnosis via fault characteristic order (FCO) analysis,” Mechanical Systems and Signal Processing, vol. 45, no. 1, pp. 139–153, 2014.

N. E. Huang, Z. Shen, S. R. Long 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 Series A, vol. 454, no. 1971, pp. 903–995, 1998.