A hybrid method combining Teager Kaiser energy operator, empirical mode decomposition and minimum entropy deconvolution for monitoring gears damages

– Amplitude demodulation is an essential tool for diagnosing gears faults. The quality of the demodulation determines the eﬃciency of the spectrum analysis in detecting the defect. A signal analysis technique combining minimum entropy deconvolution (MED), empirical mode decomposition (EMD) and Teager Kaiser energy operator (TKEO) is presented. The proposed method consists in enhancing the vibratory signal by using MED, decomposing this signal in intrinsic mode functions (IMF) and selects only the IMF which presents the highest correlation coeﬃcient with the original signal. After that, TKEO is used to track the modulation energy. The spectrum is applied to the instantaneous amplitude. The simulation and experimental results show that an envelope spectrum analysis based on MED-EMD and TKEO provides a consistent signal analysis tool.


Introduction
Gearboxes play an important role in industrial applications, and unexpected failures often result in significant economic losses. Numerous papers considering gear condition monitoring through vibration measurements were published over the years. Good reviews of most techniques used for diagnosis applied to gears are presented in [1,2]. The most important components in gear vibration spectra are the tooth-meshing frequency and its harmonics, together with sidebands due to modulation phenomena. The increment in the number and amplitude of such sidebands may indicate a fault condition [3]. Demodulation amplitude using high-frequency resonance technique is one of most techniques used to diagnosis gears. On other hand, Teager energy operator (TEO) [4] originally proposed for nonlinear signal processing provides an effective approach to extract the instantaneous amplitudes and frequencies of a modulated signal by energy operator demodulation approach and the demodulation effect is good compared with Hilbert demodulation method; meanwhile the computation time greatly decreases [5]. TEO was first used for the frequency and amplitude demodulation of the speech signals [4][5][6][7] and was adopted a Corresponding author: marc.thomas@etsmtl.ca for vibration signals of rotating machines [7][8][9]. Recently, TKEO was used to diagnosis gears defects; Li [10][11][12] used TKEO mixed with EEMD to diagnosis a gear fault and a pre-process of the signal was used to perform the TKEO. Soltani et al. [13] proposed a combined method based on wavelet and TKEO to improve the diagnosis of gear faults. Zhong et al. [14] decomposed gear faults of a wind turbine gearbox on a multiple component using EMD; after that, TKEO was applied to the selected IMF. Zhipeng et al. [15] proposed a method for fault diagnosis for wind turbine planetary gearboxes via demodulation analysis based on ensemble empirical mode decomposition (EEMD) and energy separation. The TKEO incorporates both amplitude and frequency demodulation and improves the fault signal to interference ratio and thus, could enhance detection efficiency. The energy separation algorithm is completely data driven without the need to construct any basis functions. It is therefore adaptive to the local structure of a signal. It can track the energy and identify the instantaneous frequency and instantaneous amplitude of mono-component signals.
However, it is well-known that the vibration signal from gears is a multi-component signal. For this reason, before computing TKEO, we propose to decompose the signal in multi-components by using the empirical Article published by EDP Sciences mode decomposition method (EMD), since it can adaptively decompose the signal into intrinsic mode functions (IMF) [16,17]. The TKEO can be computed for the IMF which is the most correlated with the original signal. We thus propose to use the coefficient of Pearson's correlation as an indicator to select the most effective IMF. So, the selected IMF is the one which presents the highest correlation coefficient with the original signal. However, in the early stages of gear failure, the signal amplitude is somewhat weak and often hidden by large noises and other signals from other sources. The selected IMF may thus be submerged by noise and cannot effectively discriminate the real defect signal from noise. Furthermore the TKEO is highly sensitive to the noise level. A preprocessing method is thus needed for signal decomposition by using the EMD method. We thus propose to use the Minimum Entropy Deconvolution (MED) which has been proven to be efficient in detecting incipient faults buried in large noise, and allows for an optimal filter to extract impulsive signals [18,19]. We propose to combine all these three methods for the diagnosis of gears faults. The simulation and experimental results show that this hybrid method can effectively diagnose the faults of the bearing faults. The paper is organized as follows: Section 2 presents a theoretical background of the three methods to be used. Section 3 gives the different steps of the proposed method. Sections 4 and 5 give the results from simulation and experimental data of gear faults. Finally, the conclusions are provided in Section 6.

Empirical mode decomposition (EMD)
The EMD method decomposes the time signal into a finite set of oscillatory functions called intrinsic mode functions (IMF). An IMF is a function that meets the following conditions: the number of extrema and the number of zero crossings must either equal or differ at most by one; the value of the moving average envelope defined by local maxima and the envelope defined by local minima is zero.
The decomposition method in IMF may be summarized as follows [17]: If r i+1 still has least 2 extrema then go to step (2) with i = i + 1, else, the decomposition process is finished and r i+1 is the residue of the signal.

Ensemble empirical mode decomposition (EEMD)
EMD is a well-known and widely used method. However, it suffers from the major drawback of mode mixing. To overcome the mode mixing separation problem, a new noise assisted data analysis method is proposed by Wu and Huang [20]. EEMD defines the true IMF components as the mean of an ensemble of trials, each consisting of the signal plus a white noise of finite amplitude. The proposed EEMD is defined as follows: (1) add a white noise series to the targeted data x(t); (2) decompose the data x(t) with added white noise into IMFs; (3) repeat step 1 and step 2 again and again, but with different white noise series each time; (4) obtain the (ensemble) means of corresponding IMFs of the decompositions as the final result (Eq. (1)) The critical parameters that directly affect the performance of the EEMD method are: the amplitude of the added white noise and the ensemble number needed. Lei et al. [16,[21][22][23] and Zhou et al. [24] convinced that there were no specific equations reported in the literature to guide the choice of EEMD parameters, especially the noise amplitude. In most of cases, they tried different noise amplitudes and made their selection from them. Wu and Huang [20] gave a relationship among the ensemble number N E , the amplitude of the added white noise a and the standard deviation of error e by using the equation: ln (e) + a 2 ln (N E ) = 0. They suggested that the amplitude of the added white noise is approximately 0.2 of a standard deviation of the original signal and the value of ensemble is a few hundred.

Teager Kaiser energy operator (TKEO)
The Teager-Kaiser energy operator (TKEO) is a nonlinear operator considered as a high-resolution energy estimator quantifying a product of both frequency and amplitude. TKEO can track the modulation energy and identify the instantaneous amplitude and frequency. This operator was developed by Teager and Kaiser [4,5]. TKEO, applied on a continuous signal x (t), is defined as: The discrete function equivalent to the energy operator is given by: where x(n) is a discrete-time signal.

610-page 2
This operator can extract the amplitude and frequency modulations from the signal. The TKEO has a great time resolution since the operator only needs three samples into the signal to be computed. This excellent time resolution provides the ability to capture the energy fluctuations. Furthermore, this operator is very easy to implement efficiently. The Energy Separation Algorithm (ESA) developed by Maragos et al. [6,7] uses the TKEO to separate the signal x(n) into its amplitude envelope |a(t)| and the mono-component AM-FM demodulation f (t): The ESA is a very simple demodulating technique for AM-FM demodulation. It presents a good time resolution but is very sensitive to noise. For this reason, Minimum Entropy deconvolution (MED) is proposed to extract the fault impulses while minimizing the noise before computing EMD and the TKEO.

Minimum entropy deconvolution (MED)
MED was originally proposed for applications on seismic recordings by Wiggins since 1978 [25]. The MED is a deconvolution aimed for extracting the fault impulses while minimizing the noise [26]. Higher entropy corresponds to a tendency to become random and the aim of MED is to enhance the structured information into the signal by searching for an optimum set of filter coefficients that recover the output signal with the maximum value of Kurtosis. The Kurtosis is an indicator that reflects the "peakiness" of a signal, and therefore the property of impulses.
The detailed steps of this implementation can be recalled from Sawalhi et al. [18].

Proposed method
The choice of the IMFs to be analyzed is usually realized by visual or experience criteria by the user. However the process is not automatic in this way and an interaction with the user is required. We propose to use the coefficient of Pearson's correlation as an indicator to select the most efficient IMF. Consequently the selected IMF is the IMF which presents the highest correlation coefficient with the original signal. The procedure of the proposed method based on MED-EMD and TKEO is given as follows: 1. enhance the peakiness of the signal by using the MED method; 2. decompose the time domain the signal into IMFs using EMD; 3. select the IMF which has the highest correlation coefficient; 4. compute the instantaneous amplitude using TKEO; 5. analyze the spectrum of the envelope computed.
In order to demonstrate the effectiveness of the proposed method for fault detection of ball bearings, we proposed to compare between three methods.
the first method consists in applying the TKEO to the original signal; the second method consists in decomposing the signal using EMD into IMFs selecting the IMF with higher correlation coefficient and computing the TKEO (see Fig. 1); the third method consists in decomposing the signal using EEMD into IMFs selecting the IMF with higher correlation coefficient and computing the TKEO (see Fig. 1).

Simulation
As illustrated in Figure 2, the vibration signal x(t) of a gearbox can be described as a convolution between the impulse response function (IRF) of the transmission path h and the combined effect of an anomaly caused by a localized gear tooth fault (fault impulses) y(t), the deterministic signals d(t) inherent in operating gears and the noise. The model of vibration signal from a gear can be written as [19]: where " * " represents convolution operator. In general, vibration signals generated from a healthy gearbox are dominated by gear meshing vibration accompanied by some modulation effects caused by geometric, assembly errors, speed and load fluctuations of the gears. The modulation effects include both amplitude and phase modulations. In a complete revolution of the gear, the gearbox vibration signal can be described by the following equation [3,27]: where m (0, 1, . . . , M) is the meshing harmonic number; A m is the amplitude at mth harmonic frequency f m . f m = m · N f s , where N is the gear tooth number, f s is the gear shaft rotation frequency. The functions a m (t) and b m (t) are the amplitude and phase modulation functions, respectively, and β m the initial phase.
In the case of a localised fault, an impact will be produced by the fault. This impact excites the resonance of the gear which may be represented as a series of damped oscillations, superposed on system vibration [3,13,19,[28][29][30][31]. The response is influenced by the  transmission path (h) between the gears meshing area. For simplicity, the effects of transfer function on gear meshing vibration are assumed to be relatively insignificant. The impact produced by the gear can be represented as: where ζ reflects the damping characteristic of the system; w 0 is the excited resonance frequency; A i is the amplitude of the ith impulse and T s = 1/f s is the time period associated with the gear rotational speed.
So the vibration of the faulty gear signal can be modeled as: We simulate vibrations of a single-stage healthy and defective gearbox. The parameters of the simulation are given in Table 1.
The pinion and output gear have 20 and 21 teeth, respectively. The instantaneous input shaft speed is f s = 15 Hz. Three harmonics of the meshing frequency are included into the signal. A normally distributed random signal with 0 mean and standard deviation of 0.5 is added into the simulated signal. Figure 3 shows the amplitude spectra of simulated vibration signal for (A) a normal gear and (B) a defective gear, respectively. The major difference between the two spectra in Figure 3 is in the region around 4000 Hz where a damped resonance is assumed excited.    The Blue Color (full line) represents the healthy gear and the red color (dotted line) represents the effect of the defect (with repetitive impacts). The figure shows clearly the gear mesh and its harmonics. In the very low frequency range, we note the presence of the frequency rotation and its harmonic (15,30 Hz). It must be noticed that except the gear variation in amplitude, there is no significant variation in frequency between both cases. This may be caused by the fact that the TKEO becomes effective when the signal is mono-component. So, it must be concluded that it may be difficult to achieve a diagnosis by applying TKEO directly to the signal.

Method 2: Applying EMD and TKEO
The simulated defective gear was processed using the EMD technique. We have obtained 18 IMF. Figure 5 shows the IMFs of the decomposition and Table 2 lists the Pearson's correlation coefficients between the signal and their corresponding IMFs.
It may be noticed that most of the noise imbedded into the raw signal resides in the high frequency range. Figure 5 shows that the repetitive impulses (at 15 Hz) reside in IMF 4 . Its spectrum is presented in Figure 6B. It can be noticed that it exhibits the modulations and its harmonics. This is due to the fact that the EMD decomposes the signal from high to low frequencies and the impulses due to the defect are located around the damped resonance 4000 Hz.
However, the aim of the method is only to use the IMF which presents the highest coefficient value with the raw signal in order to automatically select it. According to Table 2, the IMF 9 (which represents the low frequency components) is, in this case, the component with the highest correlation value. The spectrum of TKEO is presented in Figure 6A. We note the presence of the fundamental gear mesh (300 Hz) and its harmonic, but the information related to the defect is not revealed. We must conclude that the method 2 is not enough efficient compared to the first one.

Method 3: Applying EEMD and TKEO
Applying TKEO after EEMD was developed by Zhipeng et al. [15]. Method 3 applies this method. The EEMD algorithm is applied with an ensemble of, N e = 100, and white noise amplitude of 0.2 times standard deviation of the data. Figure 7 displays the ensemble empirical mode decomposition in ten IMFs. The second IMF represents the frequency components excited by the defect. Table 3 lists the Pearson's correlation coefficients between the signal and their corresponding IMFs.
According to Table 3, the C 4 is the component which presents the highest correlation value. The spectrum of TKEO after applying EEMD is presented in Figure 8A. We note the presence of the fundamental gear mesh (300 Hz) and its harmonics. At low frequencies, we also detect the modulation related to the defect (15 Hz). We must conclude that the method 3 is more efficient compared to the second method (Fig. 6A).
The repetitive impulses reside in C 2 . Its spectrum is presented in Figure 8B. It can be noticed that it exhibits the modulations and its harmonics until 800 Hz which is better than the result observed in Figure 6B. But in this case the choice of the IMF is not selected according to the Pearson's correlation coefficients. The selection is visual. We must conclude that the method 3 is not enough efficient compared to the method 2.

Method 4: Applying MED-EMD and TKEO
Before computing EMD, the MED is used. The MED depends on the number of iterations to perform and the choice of the length of the filter (N ). Figure 9 shows the evolution of the kurtosis obtained after applying MED to the signal x(t). The result obtained from MED filter output is denotedŷ(t). Figure 9 shows the convergence characteristics of the Kurtosis ofŷ(t) over the number of iterations performed. The convergence properties of the MED filters with different filter lengths (16, 32, 64, 128, 256. 512 and 1024 samples) are compared. For the filters (16 and 32), the convergence is achieved within 7 iterative steps, while, for the other filter lengths, the convergence is achieved within 12 iterative steps regardless of  the lengths of the MED filter. The length of the MED filter has an important influence on the value of the kurtosis. We can see that when the length of the filter increases, the kurtosis values increase also. We may also notice that there is no much difference on the kurtosis values when the length of the filter is greater than 128. Consequently, we may consider that N = 128 is a good choice. However, the comparison of the time signals of the MED outputsŷ(t), performed with a filter length equal to 128 in Figure 10, shows that there is no considerable improvement from having more than 5 iterations. As cautioned by Endo et al. [19], a greater number of iterations of the MED filter can erroneously enhance small impulses, unrelated to gear faults, mainly in case of real vibration signal.
When the EMD is applied to the results obtained by MEDŷ(t), the estimated correlation is improved. Table 4 lists the Pearson's correlation coefficients between the signal and their corresponding IMFs.
The first three IMFs exhibit higher correlation coefficients than the others, and thus Figure 11A only shows the first three IMFs. Most of the noise is distributed into IMF 2 and IMF 3 , and only the periodic impulses due to the impact reside into IMF 1 . As indicated in Table 4, IMF 1 has a stronger correlation coefficient (0.77) than all the other signal components and contains the main components in the filtered signal. Hence, IMF 1 was selected as the final resultant signal recovered from the raw signal. The spectrum of TKEO for IMF 1 (obtained by MED-EMD) is presented in Figure 11B. We clearly detect   the modulation due the impact and its harmonics until 800 Hz. Consequently, the hybrid method gives improved results than the three others presented previously.

Experimental study
The recordings of vibration signal were carried out at CETIM 1 , on a gear system with a train of gearing with a ratio of 20/21 functioning continuously until its destruction. Table 5 gives the details of the gear test rig 1 Centre desÉtudes Techniques des Industries Mécaniques de Senlis, France. parameters. The test duration was 13 days with a daily mechanical appraisal; measurements were collected every 24 h except at the first day. Table 6 gives a description of  These signals have already been used on several occasions to demonstrate diagnostic procedures [28,31]. The EMD (method 2) and the MED-EMD (method 4) algorithms were compared. The length of the filter used is 1024 and the number of iteration is 5. To simplify the procedure, the IMF 1 was automatically selected. Table 7 lists the coefficient correlation between IMF 1 (obtained by applying EMD) and original signal, and between the IMF 1 (obtained by applying MED-EMD) and the filtered signal. Figure 12 shows the original signal (Figs. 12A and 12B) and IMF 1 (Figs. 12C and 12D) obtained after applying MED-EMD, for the day 2 and 3, respectively. A decrease of noise may be noticed and the periodic shocks appear clearly. The period of the shocks is equal to the shaft rotation period. The same observations are noted for the other days (see Fig. 13).
The spectrum of TKEO for the selected IMF (obtained by MED-EMD) is presented in Figures 14 and 15. Figure 14 shows all spectrums from the days 2 to 11 and Figure 15 shows all spectrums from the days 11 to 13. For the first 10 days, we clearly detect the modulation due to the impact and its harmonics. This shows that the TKEO can effectively extract the AM/FM modulations.
The change in the spectrums is marked from day 5 and day 6. We note an increase of the amplitude of the peak related to the frequency rotation of the shaft compared to the days 2, 3 and 4. According to Table 6, the chipping has only visually been observed after the day 6. According to the spectrum of TKEO of the selected IMF, the beginning of the chipping was in fact initiated at the day 5. From day 7 to 11, we note a little decrease of the amplitude, which stays still higher than the first 3 days. But we observe an increase of the number of the harmonics. The number of the harmonics can be used as an indicator to characterise the damage. The last two days (12 and 13) are noticeable by the growth of amplitude of the peak of the frequency rotation of the shaft and its harmonics. At this stage, the gear is completely damaged. The Kurtosis values for the experimental signal calculated from days 2 to 13 (day1: no signal was taken) are shown in Figure 16. The Kurtosis computed for the raw signal (blue color) indicates a change after the day 11. It can be seen that the Kurtosis increases greatly only after the day 11. But, when the signal is pre-processed by MED and decomposed by EMD, the Kurtosis values computed for the selected IMF (red color) start increasing after day 5. Consequently, the Kurtosis analysis of IMF could be considered as an advantageous indicator for early detection and characterization of faults.

Conclusion
In this paper, a new approach for fault diagnosis of gears faults is presented. The algorithm involves a simple transformation and a spectral analysis step. Four methods are compared. The selected method is hybrid and consists in first, filtered the signal by using MED. After that, the obtained signal is decomposed into several IMFs. The study of the correlation coefficient between each IMF and the original signal allows for automatically selecting the better IMF. In our study, the first component IMF 1 , which represents the highest frequency components, was considered to accomplish the diagnosis by using the TKEO metric in order to compute the instantaneous amplitude. We illustrated the principle and the effectiveness of the proposed method by analyzing a simulated signal and experimental signals of a gearbox. The analysis results show that the proposed method (method 4) is potentially useful to detect and locate gear faults.