Fault Feature Extraction and Diagnosis of Rolling Bearings Based on Enhanced Complementary Empirical Mode Decomposition with Adaptive Noise and Statistical Time-Domain Features

In this paper, a novel method is proposed to enhance the accuracy of fault diagnosis for rolling bearings. First, an enhanced complementary empirical mode decomposition with adaptive noise (ECEEMDAN) method is proposed by determining two critical parameters, namely the amplitude of added white noise (AAWN) and the ensemble trails (ET). By introducing the concept of decomposition level, the optimal AAWN can be determined by judging the mutation of mutual information (MI) between adjacent intrinsic mode functions (IMFs). Furthermore, the ET is fixed at two to reduce the computational cost. This method can avoid disturbance of the spurious mode in the signal decomposition and increase computational speed. Enhanced CEEMDAN demonstrates a more significant improvement than that of the traditional CEEMDAN. Vibration signals can be decomposed into a set of IMFs using enhanced CEEMDAN. Some IMFs, which are named intrinsic information modes (IIMs), effectively reflect the vibration characteristic. The evaluated comprehensive factor (CF), which combines the shape, crest and impulse factors, as well as the kurtosis, skewness, and latitude factor, is developed to identify the IIM. CF can retain the advantage of a single factor and make up corresponding drawbacks. Experiment results, especially for the extraction of bearing fault under variable speed, illustrate the superiority of the proposed method for the fault diagnosis of rolling bearings over other methods.


Introduction
The rolling bearing is an important part of rotating machinery. Fault diagnosis is significant to ensure normal machinery operation [1][2][3]. A diagnosis method which concentrates on the time domain, frequency domain, and time-frequency domain has been proposed to minimize the interference and extract the fault feature due to interferences from heavy background noise and other unsteady operation states [4,5]. However, the time and frequency domain methods are unsuitable in fault diagnosis in the case of non-stationary and non-linear vibration signals [6]. Time-frequency domain methods, such as wavelet transform [7] and Wigner-Ville distribution [8], have been used to extract the fault feature of rolling bearings. However, these methods cannot obtain the ideal time-frequency resolution subject to 1. Generate the noisy signal x i , x i = x + ω 0 ξ i (i = 1, · · · , N), where ξ i is white noise with unit variance, ω 0 is coefficient of added white noise, and x is the discrete signal.

2.
Extract the first IMF (c i 1 ) in each noisy signal x i decomposition by EMD.

3.
Obtain the first IMF (c 1 ) by taking the average of each c i 1 .

5.
Apply EMD to decompose r 1 + w 1 E 1 ε i and extract the first IMF to obtain the second decomposed IMF (c 2 ). w 1 is the coefficient of added white noise for this stage, and operator E m (·) is the m-th IMF by EMD decomposition. 6.
Repeat the above steps until the residual R contains less than two extrema.
The discrete signal x can be written as: where K is the number of IMFs.

Mutual Information
Information theory, which originally deals with communication problems, was developed by Shannon [29]. In this study, the key measure was entropy. Afterward, information theory was also applied to feature extraction [30][31][32].
Let x i be a discrete series with n data length. p(x i ) is the probability mass function for discrete series x i . The entropy H(X) is defined as: The joint entropy can be expressed as: where p x i , y j is the joint probability distribution of discrete series x i and y j . The conditional entropy H(X|Y ) for discrete series x i is defined as: m j=1 p x i , y j log p x i y j (8) where p x i y j is the posterior probabilities of X given Y, which can be rewritten as: Thus, the MI [33] is introduced to quantify the shared information by discrete series x j and y j : Moreover, the entropy and MI can be related (Equation (11)): I(X; Y) = H(X) − H(X|Y ) (11) By combining Equations (9) and (11), the MI is re-expressed as:

ECEEMDAN
Using the traditional CEEMDAN method, a discrete signal x can be decomposed into the following modes c i (i = 1, . . . , N): Equation (13) shows that the decomposition effect is related to two parameters, namely AAWN and ET. Traditionally, the AAWN is 0.2-0.5 times the strand deviation of the signal, and the ET is fixed at several hundred ensemble trails [17]. However, the specified guidelines for the selection of the two critical parameters are lacking. To solve this problem, first, ET is fixed at two ensemble trails. If the ET Sensors 2019, 19, 4047 5 of 24 is fixed at one trail, then the process of IMF decomposition using CEEMDAN is the same as the EMD method, preventing the elimination of mode mixing. If the ET is bigger than the two trails, then the computational speed will decrease. Hence, the optimal ET is two trails. Then, the different AAWN is defined using the following expression: where σ is the strand deviation of signal and L i is the decomposition level. Thus, Equation (13) can be rewritten as: If the decomposition level L i is determined, then the two critical parameters are also identified. MI is introduced in this study to determine the optimal decomposition level (L opt ), and (L 1 , . . . , L N ) is the defined decomposition level. A set of IMFs (c ij ( j = 1, . . . , N i )) can be obtained for each decomposition level L i using CEEMDAN. Then, the MI between adjacent IMFs (c ij and c i( j+1) ) is calculated. The MI for each L i can be expressed as: 11 MI 12 · · · MI 1(N 1 −2) MI 21 MI 22 · · · MI 2(N 2 −2) . . .
If no interference of spurious modes is observed, then the MI of corresponding IMF is almost the same under each decomposed level (Equation (17)): . . .
where M is the order of the decomposition level. Once the spurious modes are available in IMFs, the relationship between adjacent IMFs will be broken, indicating that the unsuitability of Equation (17). Hence, the optimal L opt can be determined by judging the change in MI. The Figure 1 presents the flowchart of ECEEMDAN algorithm.

Evaluation Methodology of Decomposing Effect
A signal x can be decomposed into a set of IMFs by ECEEMDAN. Ideally, each IMF has an independent frequency component, and any two IMFs are mutually orthogonal [34]. This means that the scalar product ( should be zero. However, the orthogonality, which means that the product is not zero, will be broken due to the interaction between the added white noise and the signal. Thus, this paper introduces the concept of an orthogonality index (OI) to evaluate the decomposing effect.
(3) Equation (19) is reorganized and quantified. The OI is defined as: A small OI results in an improved decomposition effect, and vice versa. Moreover, the decomposition number and computational cost are considered to be the evaluation methodology. The two indexes can reflect the interference situation of spurious mode and computational speed.

Evaluation Methodology of Decomposing Effect
A signal x can be decomposed into a set of IMFs by ECEEMDAN. Ideally, each IMF has an independent frequency component, and any two IMFs are mutually orthogonal [34]. This means that the scalar product ( c i · c j (i j)) should be zero. However, the orthogonality, which means that the product c i · c j (i j) is not zero, will be broken due to the interaction between the added white noise and the signal. Thus, this paper introduces the concept of an orthogonality index (OI) to evaluate the decomposing effect.
(3) Equation (19) is reorganized and quantified. The OI is defined as: A small OI results in an improved decomposition effect, and vice versa. Moreover, the decomposition number and computational cost are considered to be the evaluation methodology. The two indexes can reflect the interference situation of spurious mode and computational speed.

Simulation Signal Analysis
Some typical signals, which include pure (noise-free) signals, noisy signals, and simulated vibration signals, were decomposed to verify the effectiveness of ECEEMDAN. Meanwhile, the

Free-Noise Signal
In this section, ECEEMDAN is used to decompose the simulation signal, which contains three different frequency components: where the amplitudes a 1 , a 2 , and a 3 were 4, 2.5, and 5, respectively, and the frequencies f 1 , f 2 , and f 3 were 20 Hz, 10 Hz, and 3 Hz, respectively. The sample frequency ( f s ) was 1 kHz. The data length was 1024. Figure 2 presents the simulation signal.

Simulation Signal Analysis
Some typical signals, which include pure (noise-free) signals, noisy signals, and simulated vibration signals, were decomposed to verify the effectiveness of ECEEMDAN. Meanwhile, the decomposition effect of traditional methods (EEMD, CEEMD and CEEMDAN) was compared with the ECEEMDAN method.

Free-Noise Signal
In this section, ECEEMDAN is used to decompose the simulation signal, which contains three different frequency components: 1 2 3 x t x t x t x t = + + (22) where the amplitudes 1 a , 2 a , and 3 a were 4, 2.5, and 5, respectively, and the frequencies 1 f , 2 f , and 3 f were 20 Hz, 10 Hz, and 3 Hz, respectively. The sample frequency ( s f ) was 1 kHz. The data length was 1024. Figure 2 presents the simulation signal. First, ECEEMDAN was used to decompose the simulation signal (SSC). Figure 3 shows the MI of adjacent IMFs under each decomposition level. The decomposition level in this study was set from 10 −12 to 10 −1 , with a step of 10 1 . The results showed that the MI and decomposition number were the same when the decomposed levels were between 10 −12 and 10 −6 . However, when decomposition levels were higher than 10 −6 , MI mutation occurred (shown as the red ellipse in Figure 3). This finding indicates that a low decomposition level results in an improved decomposition effect. The decomposed level (L = 10 −6 ) was considered as the optimal level in this study. First, ECEEMDAN was used to decompose the simulation signal (SSC). Figure 3 shows the MI of adjacent IMFs under each decomposition level. The decomposition level in this study was set from 10 −12 to 10 −1 , with a step of 10 1 . The results showed that the MI and decomposition number were the same when the decomposed levels were between 10 −12 and 10 −6 . However, when decomposition levels were higher than 10 −6 , MI mutation occurred (shown as the red ellipse in Figure 3). This finding indicates that a low decomposition level results in an improved decomposition effect. The decomposed level (L = 10 −6 ) was considered as the optimal level in this study.  Meanwhile, improved EMD versions (EEMD, CEEMD and CEEMDAN) were also used to decompose the SSC (here, the ET and AAWN were fixed at 200 and 0.2, respectively). Figure 4 shows the comparison of ECEEMDAN with the EMD improved versions. The results show three IMFs and one residual in the ECEEMDAN method, whereas the number of IMFs using the EMD improved versions are all more than four.   Meanwhile, improved EMD versions (EEMD, CEEMD and CEEMDAN) were also used to decompose the SSC (here, the ET and AAWN were fixed at 200 and 0.2, respectively). Figure 4 shows the comparison of ECEEMDAN with the EMD improved versions. The results show three IMFs and one residual in the ECEEMDAN method, whereas the number of IMFs using the EMD improved versions are all more than four. Meanwhile, improved EMD versions (EEMD, CEEMD and CEEMDAN) were also used to decompose the SSC (here, the ET and AAWN were fixed at 200 and 0.2, respectively). Figure 4 shows the comparison of ECEEMDAN with the EMD improved versions. The results show three IMFs and one residual in the ECEEMDAN method, whereas the number of IMFs using the EMD improved versions are all more than four. and and x t x t ,x t x t ) in the original signal was also compared with the decomposed IMFs (IMF1 and IMF2, IMF2 and IMF3) using the ECEEMDAN method ( Figure 5), and the relative errors were 13% and 2%, respectively. These results illustrate that the IMFs, including the information component, were also almost the same as that of the original signal. Moreover, the spurious modes exist in IMFs using improved EMD versions. Meanwhile, the MI of the adjacent sub-signal (x 1 (t) and x 2 (t), x 2 (t) and x 3 (t)) in the original signal was also compared with the decomposed IMFs (IMF1 and IMF2, IMF2 and IMF3) using the ECEEMDAN method ( Figure 5), and the relative errors were 13% and 2%, respectively. These results illustrate that the IMFs, including the information component, were also almost the same as that of the original signal. Moreover, the spurious modes exist in IMFs using improved EMD versions. Moreover, the computational cost and OI were compared when using the improved EMD versions ( Table 1). The results show that the computational cost and OI using the ECEEMDAN method are smaller than that of the improved EMD versions. Hence, the decomposition effect for the proposed method is superior than that of the traditional methods (EEMD, CEEMD, and CEEMDAN).  Moreover, the computational cost and OI were compared when using the improved EMD versions ( Table 1). The results show that the computational cost and OI using the ECEEMDAN method are smaller than that of the improved EMD versions. Hence, the decomposition effect for the proposed method is superior than that of the traditional methods (EEMD, CEEMD, and CEEMDAN).

Noisy Signal
In this section, the noisy signals, which were 'blocks', 'bumps', 'heavy sine waves' and 'quadchirps', were decomposed using ECEEMDAN to prove the performance of the decomposition effect under the effect of noise interference. Here, the data length was 2048. Figure 6 plots the noisy signals.

Noisy Signal
In this section, the noisy signals, which were 'blocks', 'bumps', 'heavy sine waves' and 'quadchirps', were decomposed using ECEEMDAN to prove the performance of the decomposition effect under the effect of noise interference. Here, the data length was 2048. Figure 6 plots the noisy signals.
The Gaussian white noise, which had an input signal noise ratio (SNR) from −5 to 5 dB with step of 1 dB, was added into the four noisy signals. Moreover, the improved EMD versions were used to decompose the noisy signal. Evaluation methodologies, which include the computational cost, OI and decomposition number, were employed to analyze the decomposition effect under each SNR. The compared results are shown in Figures 7-10. The results show that the computational cost and OI using ECEEMDAN is evidently smaller than that of the other methods, and the decomposing number is almost smaller than that of the improved EMD versions. These results indicate that the proposed method has a strong decomposition capability when the signal is subject to noise interference.  The Gaussian white noise, which had an input signal noise ratio (SNR) from −5 to 5 dB with step of 1 dB, was added into the four noisy signals. Moreover, the improved EMD versions were used to decompose the noisy signal. Evaluation methodologies, which include the computational cost, OI and decomposition number, were employed to analyze the decomposition effect under each SNR. The compared results are shown in Figures 7-10. The results show that the computational cost and OI using ECEEMDAN is evidently smaller than that of the other methods, and the decomposing number is almost smaller than that of the improved EMD versions. These results indicate that the proposed method has a strong decomposition capability when the signal is subject to noise interference.

Comparison of computional cost
Input signal noise ratio (dB) Input signal noise ratio (dB) Input signal noise ratio (dB)

Comparison of computional cost
Input signal noise ratio (dB) Input signal noise ratio (dB) Input signal noise ratio (dB) Figure 9. Comparison of the decomposition effect for the 'heavy sine wave' signal.

Simulated Vibration Signal
In this section, ECEEMDAN is used to decompose the following simulated vibration signal (SVS):

Comparison of computional cost
Input signal noise ratio (dB) Input signal noise ratio (dB) Input signal noise ratio (dB) Figure 10. Comparison of the decomposition effect for the 'quadchirp' signal.

Simulated Vibration Signal
In this section, ECEEMDAN is used to decompose the following simulated vibration signal (SVS): where G refers to the amplitude of the vibration signal, i refers to the number of oscillation periods of vibration signal, T refers to the period of vibration oscillation and reciprocal of characteristic frequency of bearing fault ( f c ), and n 1 (t) refers to the noise component. u(·) is the operator of impulse response: where d refers to the attenuation constant and f 0 refers to the resonant frequency.
The following parameters were used in the simulated signal: G at 1.4 m/s 2 , d at 620, f 0 at 1 kHz, f c at 117 Hz, and a sample frequency ( f s ) of 20 kHz. Here, the data length was 1197. White Gaussian noise, with SNR of 6 dB, was added to the simulated signal. Figure 11 presents the noise-free signal, the noise, and the noisy signal.

Simulated Vibration Signal
In this section, ECEEMDAN is used to decompose the following simulated vibration signal (SVS): where d refers to the attenuation constant and 0 f refers to the resonant frequency.
The following parameters were used in the simulated signal: G at 1. White Gaussian noise, with SNR of 6 dB, was added to the simulated signal. Figure 11 presents the noise-free signal, the noise, and the noisy signal.  Input signal noise ratio (dB) Figure 11. Simulated vibration signal.
The decomposition level was specified from 10 −12 to 10 −1 , with a step of 10 1 . Figure 12 plots the MI of adjacent IMFs under each decomposition level, and Table 2 shows the decomposition number under each decomposition level. When the decomposition level was higher than 10 −7 , the MI demonstrated a marked change (shown as a red ellipse in Figure 12). The decomposition number changes from nine to ten, illustrating that the optimal decomposing level is smaller than 10 −7 . The decomposition level was specified from 10 −12 to 10 −1 , with a step of 10 1 . Figure 12 plots the MI of adjacent IMFs under each decomposition level, and Table 2 shows the decomposition number under each decomposition level. When the decomposition level was higher than 10 −7 , the MI demonstrated a marked change (shown as a red ellipse in Figure 12). The decomposition number changes from nine to ten, illustrating that the optimal decomposing level is smaller than 10 −7 . Decomposition Level (L) 10 −12 10 −11 10 −10 10 −9 10 −8 10 −7 10 −6 10 −5 10 −4 10 −3 10 −2 10 −1 IMF number 9 9 9 9 9 9 10 11 12 13 12 12 Moreover, the improved EMD versions were employed to decompose the SVS (ET and AAWN were 200 and 0.2, respectively). Figure 13 shows the decomposition result. The results indicate that the IMF numbers obtained using EEMD, CEEMD, and CEEMDAN were 10, 10, and 11, respectively, whereas that obtained using ECEEMDAN was only nine. These results illustrate that ECEEMDAN can restrain the spurious mode.  Decomposition Level (L) 10 −12 10 −11 10 −10 10 −9 10 −8 10 −7 10 −6 10 −5 10 −4 10 −3 10 −2 10 −1 IMF number 9 9 9 9 9 9 10 11 12 13 12 12 Moreover, the improved EMD versions were employed to decompose the SVS (ET and AAWN were 200 and 0.2, respectively). Figure 13 shows the decomposition result. The results indicate that the IMF numbers obtained using EEMD, CEEMD, and CEEMDAN were 10, 10, and 11, respectively, whereas that obtained using ECEEMDAN was only nine. These results illustrate that ECEEMDAN can restrain the spurious mode. Moreover, the improved EMD versions were employed to decompose the SVS (ET and AAWN were 200 and 0.2, respectively). Figure 13 shows the decomposition result. The results indicate that the IMF numbers obtained using EEMD, CEEMD, and CEEMDAN were 10, 10, and 11, respectively, whereas that obtained using ECEEMDAN was only nine. These results illustrate that ECEEMDAN can restrain the spurious mode. Then, the OI and computational cost were utilized to compare the decomposition effect. Table 3 shows the compared results, which indicate that the computational cost and OI were smaller than MI of adjacent IMF Then, the OI and computational cost were utilized to compare the decomposition effect. Table 3 shows the compared results, which indicate that the computational cost and OI were smaller than those of the three others decomposition methods. This finding proves that the proposed method is more suitable for vibration signal decomposition than the traditional methods.

Identification of IIM
Although ECEEMDAN can decompose a signal into a set of IMFs with physical meaning, each IMF does not completely and effectively reflect vibrational characteristics (intrinsic information). One of the IMFs mainly concentrates on the intrinsic information of the vibration signal. Traditionally, some statistical methods, such as the shape factor, crest factor, and impulse factor, as well as the kurtosis, skewness, and latitude factor, are used to identify IIMs. Table 4 shows the expressions of the statistical methods. However, these methods have different sensibilities and stabilities for various fault patterns. The latitude factor, impulse factor, and kurtosis are sensitive for initial issuance faults.
Sensitivity significantly decreases a when fault is developed to a certain degree. The shape factor method is sensitive to fault features with low frequency. Kurtosis is suitable for diagnosing fault components from medium to high frequency. These findings indicate that a single method does not effectively identify IIMs when the fault pattern development is from slight to severe. Thus, this paper combines these statistical methods into the comprehensive factor (CF) to avoid the drawbacks of a single statistical method. Assuming that ECEEMDAN decomposes signals into a set of IMFs (IMF1, IMF2, . . . , IMFN). The normalized shape factor, crest factor, and impulse factor, as well as the kurtosis, skewness, and atitude factor can represented as Sh i , Cr i , Im i , Ku i , Sk i and La i , respectively, for the i-th IMF. The CF can be expressed as: Note: x and x σ refer to the mean and standard deviation, respectively, of time series x(i).
Here, the following selected criterion is used: where k is the index of the corresponding IIM. This is the adaptive threshold identified method. Table 5 shows the framework of the identified IIM.
La i ;

End for
Find maximum of CF(CF max ); Find index corresponding CF max : k = arg max(CF); Output: the intrinsic information mode (IMF k ). In this study, the decomposed IMFs by ECEEMDAN in Figure 13 are used to select the IIM to prove the effectiveness of the CF. Figure 14 presents the CF and the single statistical method. The results show that the CF of IMF3 (CF 3 ) is the largest. The envelope spectrum of IIM (IMF3) was determined to extract the feature frequency. Moreover, the original signal was analyzed to compare the identified effectiveness of feature frequency. Figure 15a,b shows the identified result of the feature frequency. Figure 15c,d comprehensively shows the identified result. The fault frequency ( f c ) and second harmonic frequency (2 f c ) are completely identified using the IIM and the original signal. However, reflecting the third harmonic frequency (3 f c ) is difficult for the original signal, thus, 3 f c is reflected in the identified IIM. The peak frequencies were 334.2Hz and 367.6Hz for the original signal and 350.9Hz for IIM (here, 3 f c is consistent with the real signal, as shown in Figure 16 and Equation (23)). In addition, Figure 14 shows that the shape factor slightly changes and is unsuitable for IIM selection. The maximum value for the skewness factor is reflected in IMF4. The corresponding envelope spectrum for IMF4 is plotted in Figure 17. The result shows that 3 f c is difficult to identify.
Find maximum of CF ( CF max ); Find index corresponding CF max : Output: the intrinsic information mode ( IMF k ).
In this study, the decomposed IMFs by ECEEMDAN in Figure 13 are used to select the IIM to prove the effectiveness of the CF. Figure 14 presents the CF and the single statistical method. The results show that the CF of IMF3 (CF3) is the largest. The envelope spectrum of IIM (IMF3) was determined to extract the feature frequency. Moreover, the original signal was analyzed to compare the identified effectiveness of feature frequency. Figure 15a Figure 16 and Equation (23)). In addition, Figure 14 shows that the shape factor slightly changes and is unsuitable for IIM selection. The maximum value for the skewness factor is reflected in IMF4. The corresponding envelope spectrum for IMF4 is plotted in Figure 17. The result shows that 3 c f is difficult to identify.  . Figure 16. Octave frequency of the real signal.     The different input SNR and fault frequency (as shown in the input parameters in Table 6) were added to Equation (23) to further evaluate the discriminant capability of IIM. First, ECEEMDAN was used to decompose the SVS of corresponding cases, such as cases 1, 2, 3, 4, and 5 in Table 5. Then, the CF was used to identify the IIM for each case. Finally, the envelope spectrum of IIM was determined to identify the feature frequency. Meanwhile, the original signal is also used in the envelope spectrum to extract feature frequency and compare the identified effects. The identified result is shown in     The different input SNR and fault frequency (as shown in the input parameters in Table 6) were added to Equation (23) to further evaluate the discriminant capability of IIM. First, ECEEMDAN was used to decompose the SVS of corresponding cases, such as cases 1, 2, 3, 4, and 5 in Table 5. Then, the CF was used to identify the IIM for each case. Finally, the envelope spectrum of IIM was determined to identify the feature frequency. Meanwhile, the original signal is also used in the envelope spectrum to extract feature frequency and compare the identified effects. The identified result is shown in  The different input SNR and fault frequency (as shown in the input parameters in Table 6) were added to Equation (23) to further evaluate the discriminant capability of IIM. First, ECEEMDAN was used to decompose the SVS of corresponding cases, such as cases 1, 2, 3, 4, and 5 in Table 5. Then, the CF was used to identify the IIM for each case. Finally, the envelope spectrum of IIM was determined to identify the feature frequency. Meanwhile, the original signal is also used in the envelope spectrum to extract feature frequency and compare the identified effects. The identified result is shown in Table 7. The 'yes' means that the feature frequency can been correctly identified, whereas 'no' indicates the opposite. The fault frequency ( f c ) and corresponding second harmonic frequency (2 f c ) are both identified by IIM and the original signal. However, the third harmonic frequency (3 f c ) can only be identified by the IIM. Thus, these results prove that the proposed CF method is suitable for the effective selection of IIM.  Thus, this paper proposes a novel identification method, namely ECEEMDAN with the addition of CF (ECEEMDAN-CF) for the diagnosis of fault features concerning rolling bearings: (1) Obtain the IMF without interference of a spurious mode via ECEEMDAN.
(2) Find the maximum index of CF to obtain the IIM.   Thus, this paper proposes a novel identification method, namely ECEEMDAN with the addition of CF (ECEEMDAN-CF) for the diagnosis of fault features concerning rolling bearings: (1) Obtain the IMF without interference of a spurious mode via ECEEMDAN.
(2) Find the maximum index of CF to obtain the IIM.

Performance Analysis
In this section, the proposed method (ECEEMDAN-CF) is used to extract fault features from rolling bearings. This method is comprised of two parts: (i) The data obtained from the Case Western Reserve University (CWRU) [35] and the (ii) real measurements from the experimental rig.

Experiment Data from the Case Western Reserve University
The test rig was comprised of a 2 HP motor, a torque transducer, a dynamometer, and control electronics ( Figure 19). The bearing type was 6205-2RS SKF, and the parameters are shown in Table  7.

Performance Analysis
In this section, the proposed method (ECEEMDAN-CF) is used to extract fault features from rolling bearings. This method is comprised of two parts: (i) The data obtained from the Case Western Reserve University (CWRU) [35] and the (ii) real measurements from the experimental rig.

Experiment Data from the Case Western Reserve University
The test rig was comprised of a 2 HP motor, a torque transducer, a dynamometer, and control electronics ( Figure 19). The bearing type was 6205-2RS SKF, and the parameters are shown in Table 7.

Case 1 Bearing with an Inner Raceway Fault
The following expression is the theoretical characteristic frequency of the inner raceway fault ( in f ): First, ECEEMDAN is employed to decompose the fault signal. The decomposed level (L) is set from 10 −14 to 10 0 , with a step of 10 1 , and the ET is fixed to two. The MI of the adjacent IMFs corresponding to each decomposition level are shown in Figure 21. The MI suddenly changed (shown as red ellipse in Figure 21) when the decomposed level was higher than 10 −6 , which is considered to be the optimal decomposition level. Figure 22 presents the decomposition number of the IMFs for each decomposition level.

Case 1 Bearing with an Inner Raceway Fault
The following expression is the theoretical characteristic frequency of the inner raceway fault ( f in ): First, ECEEMDAN is employed to decompose the fault signal. The decomposed level (L) is set from 10 −14 to 10 0 , with a step of 10 1 , and the ET is fixed to two. The MI of the adjacent IMFs corresponding to each decomposition level are shown in Figure 21. The MI suddenly changed (shown as red ellipse in Figure 21) when the decomposed level was higher than 10 −6 , which is considered to be the optimal decomposition level. Figure 22 presents the decomposition number of the IMFs for each decomposition level.      Figure 23 shows the CF of each IMF. Here, the first CF has the maximum value. This illustrates that the first IMF is the IIM. Meanwhile, the traditional CEEMDAN method was used to decompose the inner raceway fault signal, and the first IMF was selected to identify the envelope spectrum ( Figure 24). The results show that the proposed method can identify the inner fault frequency and its corresponding harmonics (from second to seventh harmonics). The fault frequency is consistent with the theoretical frequency (Equation (26)). However, identifying the corresponding harmonics using CEEMDAN is difficult, because interference components disturb the extraction of fault information. Moreover, the computational cost is compared for the two methods. The computational costs of ECCEMDAN and CEEMDAN are 0.48 s and 9.16 s, respectively ( Figure 25). Thus, the proposed method outperforms the traditional method in terms of fault frequency extraction for rolling bearings.  Figure 23 shows the CF of each IMF. Here, the first CF has the maximum value. This illustrates that the first IMF is the IIM. Meanwhile, the traditional CEEMDAN method was used to decompose the inner raceway fault signal, and the first IMF was selected to identify the envelope spectrum ( Figure 24). The results show that the proposed method can identify the inner fault frequency and its corresponding harmonics (from second to seventh harmonics). The fault frequency is consistent with the theoretical frequency (Equation (26)). However, identifying the corresponding harmonics using CEEMDAN is difficult, because interference components disturb the extraction of fault information. Moreover, the computational cost is compared for the two methods. The computational costs of ECCEMDAN and CEEMDAN are 0.48 s and 9.16 s, respectively ( Figure 25). Thus, the proposed method outperforms the traditional method in terms of fault frequency extraction for rolling bearings.

Case 2 Bearing with an Outer Raceway Fault
In this section, ECEEMDAN is used to decompose the fault signal. A total of 14 IMFs were obtained under the condition of the optimal decomposition level (10 −6 ). Then, the CF was employed to identify the IIM (Figure 26). The results show that CF1 is the maximum, meaning that the first IMF contains the main information of the vibration signal. Meanwhile, complete ensemble local mean decomposition with adaptive noise (CELMDAN) -kurtosis and variational mode decomposition (VMD) -kurtosis [36] are also utilized to extract the IIM. Figure 27 shows the identified result of IIM for the two methods. Moreover, minimum entropy deconvolution (MED) [37] combined with autoregressive filter (AR) (MED-AR) was used to identify the fault frequency. The first and the forth IMF were the identified IIMs for CELMDAN kurtosis and VMD kurtosis, respectively.

Case 2 Bearing with an Outer Raceway Fault
In this section, ECEEMDAN is used to decompose the fault signal. A total of 14 IMFs were obtained under the condition of the optimal decomposition level (10 −6 ). Then, the CF was employed to identify the IIM (Figure 26). The results show that CF 1 is the maximum, meaning that the first IMF contains the main information of the vibration signal. Meanwhile, complete ensemble local mean decomposition with adaptive noise (CELMDAN) -kurtosis and variational mode decomposition (VMD)-kurtosis [36] are also utilized to extract the IIM. Figure 27 shows the identified result of IIM for the two methods. Moreover, minimum entropy deconvolution (MED) [37] combined with autoregressive filter (AR) (MED-AR) was used to identify the fault frequency. The first and the forth IMF were the identified IIMs for CELMDAN kurtosis and VMD kurtosis, respectively.    Figure 28d).
Whereas, auto-regression-minimum entropy deconvolution (AR-MED), VMD kurtosis and CELMDAN kurtosis could only identify the second harmonics. From the results, more obvious disturbance components existed with the VMD kurtosis method than that with the proposed method (shown as red circle in Figure 28b). The amplitude energy of the fault characteristic for the AR-MED method was smaller than that of the proposed method (Figure 28a). This means that the proposed method had a better performance for the identification of fault features than the other methods.    Figure 28d).
Whereas, auto-regression-minimum entropy deconvolution (AR-MED), VMD kurtosis and CELMDAN kurtosis could only identify the second harmonics. From the results, more obvious disturbance components existed with the VMD kurtosis method than that with the proposed method (shown as red circle in Figure 28b). The amplitude energy of the fault characteristic for the AR-MED method was smaller than that of the proposed method (Figure 28a). This means that the proposed method had a better performance for the identification of fault features than the other methods.  Figure 28 shows the identified result of the fault frequency. It was found that the proposed method could identify the tenth fault harmonic (10 f out , shown as the red point in Figure 28d). Whereas, auto-regression-minimum entropy deconvolution (AR-MED), VMD kurtosis and CELMDAN kurtosis could only identify the second harmonics. From the results, more obvious disturbance components existed with the VMD kurtosis method than that with the proposed method (shown as red circle in Figure 28b). The amplitude energy of the fault characteristic for the AR-MED method was smaller than that of the proposed method (Figure 28a). This means that the proposed method had a better performance for the identification of fault features than the other methods. Figure 28 shows the identified result of the fault frequency. It was found that the proposed method could identify the tenth fault harmonic (10 out f , shown as the red point in Figure 28d).
Whereas, auto-regression-minimum entropy deconvolution (AR-MED), VMD kurtosis and CELMDAN kurtosis could only identify the second harmonics. From the results, more obvious disturbance components existed with the VMD kurtosis method than that with the proposed method (shown as red circle in Figure 28b). The amplitude energy of the fault characteristic for the AR-MED method was smaller than that of the proposed method (Figure 28a). This means that the proposed method had a better performance for the identification of fault features than the other methods.

Data from Real Measurement
In this section, the fault signal is acquired on the following test rigs ( Figure 29). The defective bearing with an inner raceway fault and a healthy bearing were installed on the experiment rig. The temperatures of the inner and outer raceways were monitored by two wireless temperature sensors and two thermocouple sensors, respectively. Accelerometer sensors were installed on bearing house in the vertical direction.

Data from Real Measurement
In this section, the fault signal is acquired on the following test rigs ( Figure 29). The defective bearing with an inner raceway fault and a healthy bearing were installed on the experiment rig. The temperatures of the inner and outer raceways were monitored by two wireless temperature sensors and two thermocouple sensors, respectively. Accelerometer sensors were installed on bearing house in the vertical direction. Below, Table 8 presents the bearing parameters. The sample frequency ( s f ) was 4 kHz. The rotating speeds of motor were set to 5000, 7500 and 9000 RPM. Figure 30 shows the acquired vibration signal, temperature signal, and rotational speed. At increased rotation speeds, the temperatures of the inner and outer raceways for fault bearing also increased, and the temperature of the outer raceway was higher than that of inner raceway for faulty and healthy bearings. Below, Table 8 presents the bearing parameters. The sample frequency ( f s ) was 4 kHz. The rotating speeds of motor were set to 5000, 7500 and 9000 RPM. Figure 30 shows the acquired vibration signal, temperature signal, and rotational speed. At increased rotation speeds, the temperatures of the inner and outer raceways for fault bearing also increased, and the temperature of the outer raceway was higher than that of inner raceway for faulty and healthy bearings.
Combining Table 8 and Equation (27), the theoretical feature frequency for faulty bearings under each test condition was 488.5 Hz at 5000 rpm, 732.7 Hz at 7500 rpm and 879.32 Hz at 9000 rpm.
The proposed method was used to extract fault features under different rotation speeds. Meanwhile, the CELMDAN kurtosis, final intrinsic mode functions-de-trended fluctuation analysis (FIMF-DFA) [38], VMD kurtosis and AR-MED methods were used to extract the fault feature. The data length was 2048 for each test condition. Figures 31-33 show the identified results at three different speeds. This finding indicates that the CELMDAN kurtosis and VMD kurtosis methods did not completely identify fault features for the three different test conditions. For the rotation speeds of 5000 and 9000 RPM, evident interference components ( f int ) are observed around the fault frequency. When the rotation speed was 7500RPM, the amplitude of the fault frequency was considerably small, such that it was difficult to extract the fault frequency using the AR-MED and FIMF-DFA methods. The identified fault frequencies at 5000, 7500, and 9000 RPM, using the proposed method, were 488 Hz, 716.4 Hz and 878.4 Hz, respectively, with relative errors of 0.1%, 2.22% and 0.1%, respectively. Evidently, the fault frequency ( f in ) can be identified by the proposed method, proving the suitability of the proposed method for fault component extraction under variable speed conditions. Below, Table 8 presents the bearing parameters. The sample frequency ( s f ) was 4 kHz. The rotating speeds of motor were set to 5000, 7500 and 9000 RPM. Figure 30 shows the acquired vibration signal, temperature signal, and rotational speed. At increased rotation speeds, the temperatures of the inner and outer raceways for fault bearing also increased, and the temperature of the outer raceway was higher than that of inner raceway for faulty and healthy bearings. Combining Table 8 and Equation (27), the theoretical feature frequency for faulty bearings under each test condition was 488.5 Hz at 5000 rpm, 732.7 Hz at 7500 rpm and 879.32 Hz at 9000 rpm.
The proposed method was used to extract fault features under different rotation speeds. Meanwhile, the CELMDAN kurtosis, final intrinsic mode functions-de-trended fluctuation analysis (FIMF-DFA) [38], VMD kurtosis and AR-MED methods were used to extract the fault feature. The data length was 2048 for each test condition. Figures 31, 32 and 33 show the identified results at three different speeds. This finding indicates that the CELMDAN kurtosis and VMD kurtosis methods did not completely identify fault features for the three different test conditions. For the rotation speeds of 5000 and 9000 RPM, evident interference components ( int f ) are observed around the fault frequency. When the rotation speed was 7500RPM, the amplitude of the fault frequency was considerably small, such that it was difficult to extract the fault frequency using the AR-MED and FIMF-DFA methods. The identified fault frequencies at 5000, 7500, and 9000 RPM, using the proposed method, were 488 Hz, 716.4 Hz and 878.4 Hz, respectively, with relative errors of 0.1%, 2.22% and     Meanwhile, the feature frequency at different rotation speeds was also compared ( Figure 34). The results show that a high rotation speed results in the large energy amplitude for the fault frequency.      Meanwhile, the feature frequency at different rotation speeds was also compared ( Figure 34). The results show that a high rotation speed results in the large energy amplitude for the fault frequency.

Conclusion
This paper proposes a novel method for fault feature identification concerning rolling bearings. The following conclusions are drawn: (1) The ECEEMDAN method, which sets ET as two, was introduced and used to determine the AAWN by the MI of adjacent IMFs under different decomposition levels. This method can improve the computational speed and avoid the interference of the spurious mode under interactions of added white noise and signal. ECEEMDAN also enhances the decomposition effectiveness of the original CEEMDAN method.
(2) The CF, which is established through combination with statistical methods, is proposed to identify the IIM of vibration signals. This method can compensate for the drawbacks of single statistical methods, such as insensitivity and unsteadiness, improving the identified effectiveness of

Conclusions
This paper proposes a novel method for fault feature identification concerning rolling bearings. The following conclusions are drawn: (1) The ECEEMDAN method, which sets ET as two, was introduced and used to determine the AAWN by the MI of adjacent IMFs under different decomposition levels. This method can improve the computational speed and avoid the interference of the spurious mode under interactions of added white noise and signal. ECEEMDAN also enhances the decomposition effectiveness of the original CEEMDAN method. (2) The CF, which is established through combination with statistical methods, is proposed to identify the IIM of vibration signals. This method can compensate for the drawbacks of single statistical methods, such as insensitivity and unsteadiness, improving the identified effectiveness of the IIM.
In comparison with typical methods, the results indicate that the proposed method outperforms typical methods for the fault feature extraction of rolling bearings.
Notably, when the proposed method is used in various applications, the decomposition level should be investigated in future works.