Self-Adaptive Fault Feature Extraction of Rolling Bearings Based on Enhancing Mode Characteristic of Complete Ensemble Empirical Mode Decomposition with Adaptive Noise

Fang Ma 1,†, Liwei Zhan 1,*,† , Chengwei Li 2, Zhenghui Li 1 and Tingjian Wang 3 1 Aero Engine Corporation of China Harbin bearing Co., Ltd., Harbin 150500, China; mafang521184@163.com (F.M.); zhenghuili90@163.com (Z.L.) 2 School of Electrical Engineering and Automation, Harbin Institute of Technology, Harbin 150001, China; ChengweiLI@hit.edu.cn 3 College of Mechanical Engineering, Tianjin University of Technology and Education, Tianjin 300222, China; wangtingjian1231@163.com * Correspondence: zhanliwei333@163.com † The author Liwei Zhan and Fang Ma are equal contribution for this paper.


Introduction
Rolling bearings are a key part of rotary machines and are directly related to machine health status.Thus, it is very essential to introduce an effective method to accurately detect potential abnormalities in bearing.Failure to take such measures risks the occurrence of unexpected and potentially disastrous events [1][2][3].
Currently, vibration analysis technology plays important effect on fault diagnosis domain of rolling bearing.However, traditional signal analysis method (such as Fast Fourier Transform (FFT)) is not suitable for processing vibration signals because of its non-stationary and non-linear characteristics [4].Wavelet transform is another typical application in the fault diagnosis domain [5,6].However, it is difficult to select a wavelet base function.Afterwards, empirical mode decomposition (EMD) was introduced by Huang to process the signal which has complex non-linear and non-stationary characteristic [7].This method can automatically decompose a complicated signal into sets of intrinsic mode function (IMFs).Each IMF is a mono-component and represents a characteristic of modulated amplitude signal.However, EMD presents drawbacks in terms of end effects and mode mixing.For phase-amplitude coupling (PAC) assessment [8], the drawbacks can lead to in-nacuracy in frequency resolution; for signals with shape edges and spikes, it can lead to potential spurious amplitude-amplitude coupling (AAC) [9].Thus, when it is used to analyze fault signals, it is also limited in terms of extracting characteristic information for rolling bearings.To overcome these drawbacks, Deering et al. [10] and Yeh et al. [11] used a masking signal to improve the decomposed performance of EMD.Reference [11] indicates that the essential relation was illustrated quantifiably between low-and high-frequency bands using the masking phase and amplitude coupling method.However, the authors did not adequately check the choosing procedure of a proper masking signal.Afterwards, Wu et al. [12] and Yeh et al. [13] proposed the ensemble EMD (EEMD) and Complementary EEMD (CEEMD) methods, respectively.However, EEMD and CEEMD bring new issues of spurious modes, because of interactions of added noise and signal in signal decomposition.Then, Torres [14] introduced a complete EEMD with an adaptive noise (CEEMDAN) method to alleviate the drawbacks of EEMD and CEEMD.However, when CEEMDAN is used to analyze the vibration signal, two issues were observed: (1) how to obtain mode information of the original signal with physical meanings; (2) how to identify the IIM of the vibration signal.
(1) On the issue of decomposing an intrinsic mode function with physical meaning When the vibration signal is decomposed by CEEMDAN, it is necessary to determine two critical parameters: the first is the amplitude of added white noise, and the other is ensemble trails.Lei et al. [15,16] and Zhou et al. [17] believe that no specific expression is used to guide the choice of the two parameters.Chang et al. [18], Zhang et al. [19] and Zvokelj et al. [20] use the signal-noise ratio (SNR) to choose a corresponding amplitude.However, it is assumed that the signal can be analyzed.For vibration signals, Zvokelj et al. [21] still apply experimental methods to judge the amplitude; this may result in added noise not being fully eliminated, and in the decomposed IMFs containing the interference component because of interactions between the added noise and the signal.Finally, although CEEMDAN can avoid mode mixing, it also produces over-decomposition.When IMFs are used to extract fault features, they contain interference frequency components that influence the extraction of the actual bearing fault.To avoid the drawbacks associated with the use of spurious mode, Marcelo et al. [22] proposed a method to extract the k-th IMF by adding an EMD decomposed mode.Although this method may reduce residual noise and obtain more modes with physical meanings, it still has interference in the spurious mode.Yang et al. [23] introduced ensemble local mean decomposition (ELMD) to address the drawbacks of mode mixing in local mean decomposition (LMD).However, the ELMD still has the problem of residue noise and spurious mode.Afterwards, Wang et al. [24] developed a complete ensemble local mean decomposition with adaptive noise (CELMDAN) to avoid the drawbacks in ELMD.However, the CELMDAN still suffers from the same problems as CEEMDAN.Thus, the LMD improved versions (ELMD and CELMDAN) are not also suitable for decomposing vibration signals.Gao et al. [25] proposed combining the methods of adjacent IMFs to obtain a minimum number IMFs with physical meanings; however, this method only relies on IMFs obtaining combined IMF, and the results may introduce distorted information, thereby impeding accurate fault feature extraction.
(2) Issue of identifying IIM among IMFs Generally, the main fault information of bearing rolling concentrates on one CIMF.A proper method must be developped to find this mode.Traditionally, Kurtosis [26] is sensitive to vibration Symmetry 2019, 11, 513 3 of 20 signals.However, when a fault become severe, it is not consistent with the development of a fault.The correlation coefficient [27,28] and Kullback-Leibler (K-L) divergence [29][30][31] do not work well for identifying faults.Recently, sample entropy (SE) [32,33], approximate entropy (AE) [34,35] and fuzzy entropy (FE) [36,37] were introduced into the fault diagnosis domain.However, at the boundary, SE is changed suddenly and not continuously.AE has the drawback of being heavily reliant upon the data length.FE originates in membership functions, and the function is difficult to determine.
Generally, there were two fault diagnosis methods: one is the discriminant method based on the model, and the other is based on the feature frequency.The former combines artificial intelligence and statistical approaches, such as artificial neural networks (ANN), principal component analysis (PCA), fuzzy logic, Kalman filters and so on, to set up a discriminant model to extract fault [38,39].For example, Nguyen et al. [40] proposed percentile measurements and a gamma process model to estimate the remaining useful life (RUL) in batch manufacturing processes.Samir et al. [41] proposed a dynamic model to generate fault indicators to identify a cluster for each operation state.This method estimates RUL without prior knowledge.For the second method, Wang et al. [24] first employed ELMD and CELMDAN to decompose the vibration signal and to obtain mode components of the vibration signal, respectively.Then, the kurtosis maximum principle was used to identify IIM and for analyses of the envelope spectrum for IIM.The method can extract fault features, but there are some interference spectra.Li et al. [42] introduced a method based on frequency band entropy to solve the selected problem of IIM.However, the final results showed that the different identified method leads to big differences in the extracted fault frequency.
To overcome the above problem, firstly, this paper proposes a method automatically combining IMFs into a minimum number of IMFs with physical meanings to enhance CEEMDAN decomposition.Here, this method produces new mode function (CIMFs) by combining adjacent IMFs; then, the FFT is done for each CIMF to enhance the difference features of the frequency component in IMFs; all information can then be reflected by PDF estimation.The difference of CIMFs is qualified by MHD [43] to obtain new IMFs with a minimum number and physical meanings.Finally, the comprehensive evaluated index based on Kurtosis and de-trended fluctuation analysis (DFA) [44][45][46] was proposed for the extraction of vibration IIMs.The experiment results show that the proposed method outperforms the other methods for rolling bearing fault diagnoses.
The remainder of this paper is organized as follows: Section 2 introduces the related works to illustrate the principles of existing algorithms; Section 3 describes proposed works (enhancing decomposed algorithm of CEEMDAN and identifying IIM algorithm); Section 4 evaluates the performance of the proposed method for fault diagnosis; Section 5 is the conclusion.

CEEMDAN Algorithm
It is well known that EMD [7] has the drawback of mode mixing and end effect.EEMD [12] was introduced to solve this problem.However, it introduces residual noise into the IMFs and produces a spurious mode.To alleviate these phenomena, CEEMDAN [14] was introduced to decompose non-liner and non-stationary signals.Currently, this method is also widely applied in the fault diagnosis domain to obtain feature information on vibration signals [47,48].The following is the basic principle of CEEMDAN: (1) Apply EMD to decompose signal (x(t) + w 0 ε i (t)) to extract the first IMF; Here, ε(t) is white noise with unit variance, and w 0 is the corresponding amplitude.(2) Obtain the differential signal by Equation (2); (3) Extract the first mode by decompose r 1 (t) + w 1 E 1 ε i (t) and appoint the second IMF using Equation (3); (4) Decompose the k-th residue and extract the first mode, and the (k + 1)-th mode can be obtained using the following equation (k = 2, . . ., K): Here, E k (•) indicates the symbol function of the k-th IMF.(5) Obtain last mode by Equation ( 5) when the residue has fewer than (or equal to) two extrema by repeating step, Signal x(t) is decomposed using equation,

De-Trended Fluctuation Algorithm (DFA)
Recently, DFA [44] has been known as a sales tool, and has been used to analyze non-stationary signals, especially for the volatility of time series.The DFA principle is as follows: (1) Divide the following series into integrated time: Here, x indicates the mean of the series x(i).
(2) Slice y(k) into n length sub-section; (3) Apply least square method fitting to obtain the local trend y n (k); (4) Extract y n (k) from the y(k) to obtain fluctuation function F(n); (5) Obtain different F(n) by different length segments; (6) Calculate the slope α (fractal scaling index) between log(F(n)) and log(k); the bigger slope α of signal, the smoother time series (as shown in Equation ( 9)): Symmetry 2019, 11, 513 5 of 20

Modified Hausdorff Distance (MHD)
Hausdorff distance (HD) [49] was first used in image processing.HD computes the similarity between two objects [50,51].The quantity principle is based on object shape.The following expression is the distance of two point (a and b): Then, the distance between a point a and a set of points B = b 1 , . . ., b N b is defined as: The Equation ( 11) is defined as direct distance between two series A = {a 1 , . . . ,a N a } and The modified directed distance (MHD) can been obtained using the following equation: The Equation ( 13) is the undirected distance: Thus, the modified Hausdorff distance (MHD) is the MHD = f(d 2 (A, B), d 2 (B, A)).When the signal is influenced by noise, MHD can obtain the real distance.

Identifying of IMF with Minimum Number and Physical Meaning
In this section, a method is proposed to obtain the minimum number of IMFs and physical meanings based on frequency characteristics of IMFs.Firstly, it combines the adjacent IMFs into new modes (sign as CIMFs).Then, a FFT analysis is performed for each CIMF (obtained waveform is sign as FC).In this way, the FFT can organize the IMF containing the same or similar frequency bands into a class.In addition, because a PDF can reflect a state of knowledge of the desired signal rather than only a frequency based on Bayesian theory, a PDF contains complete information about the FC, and is expected to capture the feature information of FCs (sign as PFC).Finally, the similarity is quantized by the MHD of the adjacent PFC to adaptively obtain the final IMFs (FIMFs).
Here, the following simulated signal is introduced to illustrate the identified process.
In Equation ( 15), the amplitude A 1 and A 2 are 120 and 50, and frequency f 1 and f 2 are 10 Hz and 75 Hz, respectively.n(t) obey distribution N 0, 11 2 .The sample frequency is 1 kHz, and the data length is 20148 (as shown Figure 1).Theoretically, the four IMFs should be obtained through the decomposition of signal ( ) x t (the first is the noise component, the second is the IMF with a frequency component of 75 Hz, the third is that of 10 Hz, and last one is residue).However, it was found that twelve IMFs were obtained by CEEMDAN decomposition (as shown Figure 2).Here, the amplitude of white noise is 0.2, and the ensemble size is 100; this means that there are some spurious modes because of the interactions between the signal and the added noise.It is difficult to identify real frequency information of the original signal after it had been subjected to the interference of spurious mode.To solve this problem, the proposed method is introduced to automatically distinguish the frequency components of the simulated signal.Firstly, the adjacent IMFs are combined by Equation ( 16) to obtain the CIMFs (as shown Figure 3).Theoretically, the four IMFs should be obtained through the decomposition of signal x(t) (the first is the noise component, the second is the IMF with a frequency component of 75 Hz, the third is that of 10 Hz, and last one is residue).However, it was found that twelve IMFs were obtained by CEEMDAN decomposition (as shown Figure 2).Here, the amplitude of white noise is 0.2, and the ensemble size is 100; this means that there are some spurious modes because of the interactions between the signal and the added noise.It is difficult to identify real frequency information of the original signal after it had been subjected to the interference of spurious mode.To solve this problem, the proposed method is introduced to automatically distinguish the frequency components of the simulated signal.Firstly, the adjacent IMFs are combined by Equation ( 16) to obtain the CIMFs (as shown Figure 3).Theoretically, the four IMFs should be obtained through the decomposition of signal ( ) x t (the first is the noise component, the second is the IMF with a frequency component of 75 Hz, the third is that of 10 Hz, and last one is residue).However, it was found that twelve IMFs were obtained by CEEMDAN decomposition (as shown Figure 2).Here, the amplitude of white noise is 0.2, and the ensemble size is 100; this means that there are some spurious modes because of the interactions between the signal and the added noise.It is difficult to identify real frequency information of the original signal after it had been subjected to the interference of spurious mode.To solve this problem, the proposed method is introduced to automatically distinguish the frequency components of the simulated signal.Firstly, the adjacent IMFs are combined by Equation ( 16) to obtain the CIMFs (as shown Figure 3).Next, FFT is used to extract the frequency component of the CIMF.Here, the waveform, which is obtained by FFT of CIMF, is referred to as FC.To distinguish the feature difference among FCs, the main frequency band (the dotted box in Figure 4), which contains from 0 Hz to 100 Hz in all FCs, is enlarged (except the frequency band from 100 Hz to 500 Hz for FC1 and FC2). Figure 5 illustrates the frequency spectrum and amplitude information of the original and enlarged FC.It was found that the frequency spectrum of FC1 is same as that of FC2, and that the maximum amplitude of FC1 is less than that of FC2.Among FC3, FC4 and FC5, the frequency spectra were similar.However, the maximum amplitude of FC3 is less than that of FC4 and FC5; the residue FCs (FC6, FC7, FC8, FC9, FC10 and FC11) are identical; this means that the IMF1 and IMF2 can been divided into a class, as can IMF6, IMF 7, IMF 8, IMF 9, IMF 10 and IMF 11.However, it is impossible to determine whether IMF3, IMF4 and IMF5 may also organized into a class.Next, the PDF and MHD were used to enhance their feature difference and quantify classification results.Next, FFT is used to extract the frequency component of the CIMF.Here, the waveform, which is obtained by FFT of CIMF, is referred to as FC.To distinguish the feature difference among FCs, the main frequency band (the dotted box in Figure 4), which contains from 0 Hz to 100 Hz in all FCs, is enlarged (except the frequency band from 100 Hz to 500 Hz for FC1 and FC2). Figure 5 illustrates the frequency spectrum and amplitude information of the original and enlarged FC.It was found that the frequency spectrum of FC1 is same as that of FC2, and that the maximum amplitude of FC1 is less than that of FC2.Among FC3, FC4 and FC5, the frequency spectra were similar.However, the maximum amplitude of FC3 is less than that of FC4 and FC5; the residue FCs (FC6, FC7, FC8, FC9, FC10 and FC11) are identical; this means that the IMF1 and IMF2 can been divided into a class, as can IMF6, IMF 7, IMF 8, IMF 9, IMF 10 and IMF 11.However, it is impossible to determine whether IMF3, IMF4 and IMF5 may also organized into a class.Next, the PDF and MHD were used to enhance their feature difference and quantify classification results.Here, the PDF is used to estimate each FC (the obtained waveform is designated as PFC) to enhance discernment of the difference features among FCs (especially for FC3, FC4 and FC5) and subsequent quantitative ability.Figure 5 shows the estimated PDF and local zoom of the interest components in each PFC.It was found that the geometry feature of each PFC was more obvious than that of the corresponding FC.Next, FFT is used to extract the frequency component of the CIMF.Here, the waveform, which is obtained by FFT of CIMF, is referred to as FC.To distinguish the feature difference among FCs, the main frequency band (the dotted box in Figure 4), which contains from 0 Hz to 100 Hz in all FCs, is enlarged (except the frequency band from 100 Hz to 500 Hz for FC1 and FC2). Figure 5 illustrates the frequency spectrum and amplitude information of the original and enlarged FC.It was found that the frequency spectrum of FC1 is same as that of FC2, and that the maximum amplitude of FC1 is less than that of FC2.Among FC3, FC4 and FC5, the frequency spectra were similar.However, the maximum amplitude of FC3 is less than that of FC4 and FC5; the residue FCs (FC6, FC7, FC8, FC9, FC10 and FC11) are identical; this means that the IMF1 and IMF2 can been divided into a class, as can IMF6, IMF 7, IMF 8, IMF 9, IMF 10 and IMF 11.However, it is impossible to determine whether IMF3, IMF4 and IMF5 may also organized into a class.Next, the PDF and MHD were used to enhance their feature difference and quantify classification results.Here, the PDF is used to estimate each FC (the obtained waveform is designated as PFC) to enhance discernment of the difference features among FCs (especially for FC3, FC4 and FC5) and subsequent quantitative ability.Figure 5 shows the estimated PDF and local zoom of the interest components in each PFC.It was found that the geometry feature of each PFC was more obvious than that of the corresponding FC.Finally, modified Hausdorff distance (MHD) is used to calculate the geometry distance between adjacent PFCs (as shown Figure 6).Here, the PDF is used to estimate each FC (the obtained waveform is designated as PFC) to enhance discernment of the difference features among FCs (especially for FC3, FC4 and FC5) and subsequent quantitative ability.Figure 5 shows the estimated PDF and local zoom of the interest components in each PFC.It was found that the geometry feature of each PFC was more obvious than that of the corresponding FC.
Finally, modified Hausdorff distance (MHD) is used to calculate the geometry distance between adjacent PFCs (as shown Figure 6).Finally, modified Hausdorff distance (MHD) is used to calculate the geometry distance between adjacent PFCs (as shown Figure 6).The second and fifth indices of MHD are the local maxima; this means that CIMF2 and CIMF3, CIMF5 and CIMF6 are not similar.CIMF2 combines IMF1 and IMF2, and CIMF3 combines IMF1, IMF2 and IMF3; this illustrates that IMF3 is not similar to IMF1 and IMF2.Likewise, IMFs after IMF6 are not similar to IMF3, IMF4 and IMF5; this means that the IMF3, IMF4 and IMF5 can been organized into a class.Figure 7 shows a relational graph of CIMFs ('Y' means that there is a similarity among adjacent CIMFs, while 'N' means non-similarity).Thus, the CIMFs are obtained using following equation:   The second and fifth indices of MHD are the local maxima; this means that CIMF 2 and CIMF 3 , CIMF 5 and CIMF 6 are not similar.CIMF 2 combines IMF 1 and IMF 2 , and CIMF 3 combines IMF 1 , IMF 2 and IMF 3 ; this illustrates that IMF 3 is not similar to IMF 1 and IMF 2 .Likewise, IMFs after IMF 6 are not similar to IMF 3 , IMF 4 and IMF 5 ; this means that the IMF 3 , IMF 4 and IMF 5 can been organized into a class.Figure 7 shows a relational graph of CIMFs ('Y' means that there is a similarity among adjacent CIMFs, while 'N' means non-similarity).Thus, the CIMFs are obtained using following equation: Finally, modified Hausdorff distance (MHD) is used to calculate the geometry distance between adjacent PFCs (as shown Figure 6).The second and fifth indices of MHD are the local maxima; this means that CIMF2 and CIMF3, CIMF5 and CIMF6 are not similar.CIMF2 combines IMF1 and IMF2, and CIMF3 combines IMF1, IMF2 and IMF3; this illustrates that IMF3 is not similar to IMF1 and IMF2.Likewise, IMFs after IMF6 are not similar to IMF3, IMF4 and IMF5; this means that the IMF3, IMF4 and IMF5 can been organized into a class.Figure 7 shows a relational graph of CIMFs ('Y' means that there is a similarity among adjacent CIMFs, while 'N' means non-similarity).Thus, the CIMFs are obtained using following equation:    residue, respectively; this proves that the proposed method can be used to distinguish the frequency components of a complicated signal and avoid the interference of spurious mode.

Selection of Intrinsic Information Mode
Generally, only one FCMF contains the main intrinsic information of bearing faults in FCIMFs.Thus, it was necessary to develop a proper method to identify the FCMF.Traditionally, Kurtosis [26] is often used to indicate the features of an impulse signal and to extract early faults in rolling bearings.However, Kurtosis does not indicate the fault level when it becomes more serious.The expression of Kurtosis is as follows: ( ) Here, μ indicates the mean, σ is the standard deviation and ( ) Recently, DFA [44] has been used to reflect fluctuation of time series, especially for non-stationary signals such as vibration signals; the more seriousness the bearing fault, the larger fluctuation of the vibration signal, and the larger the DFA (see Section 2.2 for more information).This means that DFA can been used to evaluate the degree of seriousness of the bearing fault.Thus, this paper proposes a comprehensive evaluation index, i S , which combines Kurtosis and DFA, to select IIM in FIMFs. Here,

Selection of Intrinsic Information Mode
Generally, only one FCMF contains the main intrinsic information of bearing faults in FCIMFs.Thus, it was necessary to develop a proper method to identify the FCMF.Traditionally, Kurtosis [26] is often used to indicate the features of an impulse signal and to extract early faults in rolling bearings.However, Kurtosis does not indicate the fault level when it becomes more serious.The expression of Kurtosis is as follows: Here, µ indicates the mean, σ is the standard deviation and E(•) is expectation operator.Recently, DFA [44] has been used to reflect fluctuation of time series, especially for non-stationary signals such as vibration signals; the more seriousness the bearing fault, the larger fluctuation of the vibration signal, and the larger the DFA (see Section 2.2 for more information).This means that DFA can been used to evaluate the degree of seriousness of the bearing fault.Thus, this paper proposes a comprehensive evaluation index, S i , which combines Kurtosis and DFA, to select IIM in FIMFs.
Here, ku i and α i denote the i-th values of Kurtosis and DFA, respectively, N F indicates the number of FIMFs.
The proposed method is presented as follows, and a flow chart is shown in Figure 9:

Diagnose Inner Raceway Fault of Rolling Bearing
In this section, the proposed method is applied to identify the feature frequencies of fault signal,

Diagnose Inner Raceway Fault of Rolling Bearing
In this section, the proposed method is applied to identify the feature frequencies of fault signal, which draws upon the website of the Case Western Reserve University [52], to illustrate its performance.The experiment rig consisted of a 2 HP motor, a torque transducer, a dynamometer and control electronics (as shown Figure 10).Table 1 presents the parameters of the rolling bearing (bearings type: 6205-2RS SKF deep groove ball).

Diagnose Inner Raceway Fault of Rolling Bearing
In this section, the proposed method is applied to identify the feature frequencies of fault signal, which draws upon the website of the Case Western Reserve University [52], to illustrate its performance.The experiment rig consisted of a 2 HP motor, a torque transducer, a dynamometer and control electronics (as shown Figure 10).Table 1 presents the parameters of the rolling bearing (bearings type: 6205-2RS SKF deep groove ball).The following equation is the theoretical feature frequency of inner raceway fault ( in f ).
Here, r f is the rotation frequency.
The test bearing presented a single point fault with fault diameters of 0.018 mm by electrodischarge machining.The sample rate was 12 kHz.The experiment data was selected from the drive end and motor load 0 (HP), and the motor speed was 1797 RPM.The data length was 2048 (as shown in Figure 11).The following equation is the theoretical feature frequency of inner raceway fault ( f in ).
Here, f r is the rotation frequency.
The test bearing presented a single point fault with fault diameters of 0.018 mm by electro-discharge machining.The sample rate was 12 kHz.The experiment data was selected from the drive end and motor load 0 (HP), and the motor speed was 1797 RPM.The data length was 2048 (as shown in Figure 11).Here, CEEMDAN is used to decompose the vibration signal (ensemble size is 100, and the amplitude of the added noise is 0.2).The eleven IMFs and one residue were obtained (as shown Figure 12).It was found that it does not reflect the real component of original signal for each IMF completely, and that it contains some spurious information (such as IMF10).To extract the IIM of the vibration signal, the CIMFs were first obtained using Equation ( 16), and then FFT and the corresponding local zoom (main information) were conducted for each CIMF (as shown in Figure 13).It was found by comparing FC1 with FC2 that there were additional frequency bands.Similarly, the frequency component in FC6 is different from that of FC7 (the detail difference is shown Figure 14).That in FC8 is also different from that in FC9; this means that the feature difference is easy to distinguished among IMFs using the proposed method.Here, CEEMDAN is used to decompose the vibration signal (ensemble size is 100, and the amplitude of the added noise is 0.2).The eleven IMFs and one residue were obtained (as shown Figure 12).It was found that it does not reflect the real component of original signal for each IMF completely, and that it contains some spurious information (such as IMF 10 ).To extract the IIM of the vibration signal, the CIMFs were first obtained using Equation ( 16), and then FFT and the corresponding local zoom (main information) were conducted for each CIMF (as shown in Figure 13).
It was found by comparing FC1 with FC2 that there were additional frequency bands.Similarly, the frequency component in FC6 is different from that of FC7 (the detail difference is shown Figure 14).That in FC8 is also different from that in FC9; this means that the feature difference is easy to distinguished among IMFs using the proposed method.
Here, CEEMDAN is used to decompose the vibration signal (ensemble size is 100, and the amplitude of the added noise is 0.2).The eleven IMFs and one residue were obtained (as shown Figure 12).It was found that it does not reflect the real component of original signal for each IMF completely, and that it contains some spurious information (such as IMF10).To extract the IIM of the vibration signal, the CIMFs were first obtained using Equation ( 16), and then FFT and the corresponding local zoom (main information) were conducted for each CIMF (as shown in Figure 13).It was found by comparing FC1 with FC2 that there were additional frequency bands.Similarly, the frequency component in FC6 is different from that of FC7 (the detail difference is shown Figure 14).That in FC8 is also different from that in FC9; this means that the feature difference is easy to distinguished among IMFs using the proposed method.Here, CEEMDAN is used to decompose the vibration signal (ensemble size is 100, and the amplitude of the added noise is 0.2).The eleven IMFs and one residue were obtained (as shown Figure 12).It was found that it does not reflect the real component of original signal for each IMF completely, and that it contains some spurious information (such as IMF10).To extract the IIM of the vibration signal, the CIMFs were first obtained using Equation ( 16), and then FFT and the corresponding local zoom (main information) were conducted for each CIMF (as shown in Figure 13).It was found by comparing FC1 with FC2 that there were additional frequency bands.Similarly, the frequency component in FC6 is different from that of FC7 (the detail difference is shown Figure 14).That in FC8 is also different from that in FC9; this means that the feature difference is easy to distinguished among IMFs using the proposed method.The MHD is used to calculate the similar of adjacent PCs.If they are similar, the MHD is relatively small and vice versa.Figure 15 is the MHD of adjacent FCs.It is found that the theoretical cut-off point in Figure 13 is originally the third and sixth point.However, the MHD does not reflect the local maximum; it illustrates that when the MHD calculates the similarity of adjacent FCs directly, the difference of frequency features is not identified.The MHD is used to calculate the similar of adjacent PCs.If they are similar, the MHD is relatively small and vice versa.Figure 15 is the MHD of adjacent FCs.It is found that the theoretical cut-off point in Figure 13 is originally the third and sixth point.However, the MHD does not reflect the local maximum; it illustrates that when the MHD calculates the similarity of adjacent FCs directly, the difference of frequency features is not identified.
The MHD is used to calculate the similar of adjacent PCs.If they are similar, the MHD is relatively small and vice versa.Figure 15 is the MHD of adjacent FCs.It is found that the theoretical cut-off point in Figure 13 is originally the third and sixth point.However, the MHD does not reflect the local maximum; it illustrates that when the MHD calculates the similarity of adjacent FCs directly, the difference of frequency features is not identified.To identify the cut-off point of FCs, the PDF of each FC is estimated by the normal Kernel density function.Meanwhile, the MHD is also used to calculate similar of adjacent PFCs (as shown in Figure 16).It was found that there are four local maxima points.The Figure 17 is the relational graph among the CIMFs; this illustrates that CIMF1 is not similar to CIMF2, CIMF3 is not similar to CIMF4, CIMF6 is not similar to CIMF7, CIMF8 is not similar with CIMF9, and the residues are similar.To identify the cut-off point of FCs, the PDF of each FC is estimated by the normal Kernel density function.Meanwhile, the MHD is also used to calculate similar of adjacent PFCs (as shown in Figure 16).
The MHD is used to calculate the similar of adjacent PCs.If they are similar, the MHD is relatively small and vice versa.Figure 15 is the MHD of adjacent FCs.It is found that the theoretical cut-off point in Figure 13 is originally the third and sixth point.However, the MHD does not reflect the local maximum; it illustrates that when the MHD calculates the similarity of adjacent FCs directly, the difference of frequency features is not identified.To identify the cut-off point of FCs, the PDF of each FC is estimated by the normal Kernel density function.Meanwhile, the MHD is also used to calculate similar of adjacent PFCs (as shown in Figure 16).It was found that there are four local maxima points.The Figure 17 is the relational graph among the CIMFs; this illustrates that CIMF1 is not similar to CIMF2, CIMF3 is not similar to CIMF4, CIMF6 is not similar to CIMF7, CIMF8 is not similar with CIMF9, and the residues are similar.It was found that there are four local maxima points.The Figure 17 is the relational graph among the CIMFs; this illustrates that CIMF 1 is not similar to CIMF 2 , CIMF 3 is not similar to CIMF 4 , CIMF 6 is not similar to CIMF 7 , CIMF 8 is not similar with CIMF 9 , and the residues are similar.
The MHD is used to calculate the similar of adjacent PCs.If they are similar, the MHD is relatively small and vice versa.Figure 15 is the MHD of adjacent FCs.It is found that the theoretical cut-off point in Figure 13 is originally the third and sixth point.However, the MHD does not reflect the local maximum; it illustrates that when the MHD calculates the similarity of adjacent FCs directly, the difference of frequency features is not identified.To identify the cut-off point of FCs, the PDF of each FC is estimated by the normal Kernel density function.Meanwhile, the MHD is also used to calculate similar of adjacent PFCs (as shown in Figure 16).It was found that there are four local maxima points.The Figure 17 is the relational graph among the CIMFs; this illustrates that CIMF1 is not similar to CIMF2, CIMF3 is not similar to CIMF4, CIMF6 is not similar to CIMF7, CIMF8 is not similar with CIMF9, and the residues are similar.Thus, five FIMFs are obtained using the proposed method (as shown in Equation ( 21)). Figure 18 presents the FIMFs and residue.
Next, the method in Section 3.2 was used to identify IIM (scales range of DFA is from 10 to 100 with step 5). Figure 19 is the comprehensive evaluated index.It was found that the biggest value is the S 2 of the FIMF 2 ; this illustrates that FIMF 2 should be used to extract fault information from the vibration signal.Meanwhile, to compare performance of proposed method, the original CEEMDAN, combined with comprehensive evaluated index, is also used to extract fault information for rolling bearings (as shown in Figure 20).It was found that the second point was the maximum point; this means that the IMF 2 contains the main fault features of vibration signal.Figure 21 shows the identified IMF2 and FIMF2.
Thus, five FIMFs are obtained using the proposed method (as shown in Equation ( 21)). Figure 18 presents the FIMFs and residue.Next, the method in Section 3.2 was used to identify IIM (scales range of DFA is from 10 to 100 with step 5). Figure 19 is the comprehensive evaluated index.It was found that the biggest value is the 2 S of the FIMF2; this illustrates that FIMF2 should be used to extract fault information from the vibration signal.Meanwhile, to compare performance of proposed method, the original CEEMDAN, combined with comprehensive evaluated index, is also used to extract fault information for rolling bearings (as shown in Figure 20).It was found that the second point was the maximum point; this means that the IMF2 contains the main fault features of vibration signal.Figure 21 shows the identified IMF2 and FIMF2.Next, the method in Section 3.2 was used to identify IIM (scales range of DFA is from 10 to 100 with step 5). Figure 19 is the comprehensive evaluated index.It was found that the biggest value is the 2 S of the FIMF2; this illustrates that FIMF2 should be used to extract fault information from the vibration signal.Meanwhile, to compare performance of proposed method, the original CEEMDAN, combined with comprehensive evaluated index, is also used to extract fault information for rolling bearings (as shown in Figure 20).It was found that the second point was the maximum point; this means that the IMF2 contains the main fault features of vibration signal.Figure 21 shows the identified IMF2 and FIMF2.Next, the method in Section 3.2 was used to identify IIM (scales range of DFA is from 10 to 100 with step 5). Figure 19 is the comprehensive evaluated index.It was found that the biggest value is the 2 S of the FIMF2; this illustrates that FIMF2 should be used to extract fault information from the vibration signal.Meanwhile, to compare performance of proposed method, the original CEEMDAN, combined with comprehensive evaluated index, is also used to extract fault information for rolling bearings (as shown in Figure 20).It was found that the second point was the maximum point; this means that the IMF2 contains the main fault features of vibration signal.Figure 21 shows the identified IMF2 and FIMF2.Finally, the envelope spectrum was used to identify the frequency components of IMF2 and FIMF2 (as shown in Figure 22a,b).It was found that 2 r f , i f , 2 i f , 3 i f are all identified for IMF2 and FIMF2.The envelope spectrum of FIMF2 does not contain any interference frequency.However, that of IMF2 abundantly contains the interference frequency; this means that the proposed method can effectively identify fault information of rolling bearings.Meanwhile, the original signal was also used for the envelope spectrum (as shown in Figure 23).It was found that there is a greater interference component compared to IMF2.In addition, the computational load of the proposed method is compared with CEEMDAN; they are 20.24 s and 19.26 s, respectively (the amplitude of added white noise and ensemble trails are 0.1 and 100, respectively).There was only a difference of 0.98 s (the used laptop has Core i5-3230M CPU, 2.60 Ghz, RAM is 6 GB).However, the decomposed effectiveness for Finally, the envelope spectrum was used to identify the frequency components of IMF 2 and FIMF 2 (as shown in Figure 22a,b).It was found that 2 f r , f i ,2 f i ,3 f i are all identified for IMF 2 and FIMF 2 .The envelope spectrum of FIMF 2 does not contain any interference frequency.However, that of IMF 2 abundantly contains the interference frequency; this means that the proposed method can effectively identify fault information of rolling bearings.Meanwhile, the original signal was also used for the envelope spectrum (as shown in Figure 23).It was found that there is a greater interference component compared to IMF2.In addition, the computational load of the proposed method is compared with CEEMDAN; they are 20.24 s and 19.26 s, respectively (the amplitude of added white noise and ensemble trails are 0.1 and 100, respectively).There was only a difference of 0.98 s (the used laptop has Core i5-3230M CPU, 2.60 Ghz, RAM is 6 GB).However, the decomposed effectiveness for proposed method outperformed that of CEEMDAN.Thus, the proposed method is more suitable for fault diagnoses of rolling bearings.
Finally, the envelope spectrum was used to identify the frequency components of IMF2 and FIMF2 (as shown in Figure 22a,b).It was found that 2 r f , i f , 2 i f , 3 i f are all identified for IMF2 and FIMF2.The envelope spectrum of FIMF2 does not contain any interference frequency.However, that of IMF2 abundantly contains the interference frequency; this means that the proposed method can effectively identify fault information of rolling bearings.Meanwhile, the original signal was also used for the envelope spectrum (as shown in Figure 23).It was found that there is a greater interference component compared to IMF2.In addition, the computational load of the proposed method is compared with CEEMDAN; they are 20.24 s and 19.26 s, respectively (the amplitude of added white noise and ensemble trails are 0.1 and 100, respectively).There was only a difference of 0.98 s (the used laptop has Core i5-3230M CPU, 2.60 Ghz, RAM is 6 GB).However, the decomposed effectiveness for proposed method outperformed that of CEEMDAN.Thus, the proposed method is more suitable for fault diagnoses of rolling bearings.

Diagnose Outer Raceway Fault of Rolling Bearing
In this experiment, the following setup (as shown in Figure 24) was used to extract the fault features of rolling bearings.It consisted of four rolling bearings (symbols 1-4), an electric spindle FIMF2 (as shown in Figure 22a,b).It was found that 2 r f , i f , 2 i f , 3 i f are all identified for IMF2 and FIMF2.The envelope spectrum of FIMF2 does not contain any interference frequency.However, that of IMF2 abundantly contains the interference frequency; this means that the proposed method can effectively identify fault information of rolling bearings.Meanwhile, the original signal was also used for the envelope spectrum (as shown in Figure 23).It was found that there is a greater interference component compared to IMF2.In addition, the computational load of the proposed method is compared with CEEMDAN; they are 20.24 s and 19.26 s, respectively (the amplitude of added white noise and ensemble trails are 0.1 and 100, respectively).There was only a difference of 0.98 s (the used laptop has Core i5-3230M CPU, 2.60 Ghz, RAM is 6 GB).However, the decomposed effectiveness for proposed method outperformed that of CEEMDAN.Thus, the proposed method is more suitable for fault diagnoses of rolling bearings.

Diagnose Outer Raceway Fault of Rolling Bearing
In this experiment, the following setup (as shown in Figure 24) was used to extract the fault features of rolling bearings.It consisted of four rolling bearings (symbols 1-4), an electric spindle

Diagnose Outer Raceway Fault of Rolling Bearing
In this experiment, the following setup (as shown in Figure 24) was used to extract the fault features of rolling bearings.It consisted of four rolling bearings (symbols 1-4), an electric spindle (symbol 5), four temperature sensors (for monitoring the temperature of the bearing out ring; symbol 6), two non-contact temperature sensors (for monitoring the temperature of the inner raceway; symbol 7), and a load cell (symbol 8).
The test bearing presented a single point fault which was caused by electro-discharge machining (as shown in Figure 25).Table 2 shows the bearing parameters.The test bearing presented a single point fault which was caused by electro-discharge machining (as shown in Figure 25).Table 2 shows the bearing parameters.Equation ( 22) is the theoretical feature frequency of the outer raceway fault ( o f ): The sample rate was 4 kHz.The selected data length was 2000 (as shown in Figure 26).Equation ( 22) is the theoretical feature frequency of the outer raceway fault ( f o ): The sample rate was 4 kHz.The selected data length was 2000 (as shown in Figure 26).The test bearing presented a single point fault which was caused by electro-discharge machining (as shown in Figure 25).Table 2 shows the bearing parameters.Equation ( 22) is the theoretical feature frequency of the outer raceway fault ( o f ): The sample rate was 4 kHz.The selected data length was 2000 (as shown in Figure 26).Firstly, the proposed method was used to obtain FIMFs.The two critical parameters in CEEMDAN were 0.2 (the amplitude of added white noise) and 100 (ensemble size).FIMFs are shown in Figure 27.Then, the comprehensive evaluated index was used to identify the IIM (scales range of DFA was from 10 to 100, with step 5 in this paper).Meanwhile, the AE, SE and FE were also used to identify IIM to compare the identified performance with the proposed method.Figure 28 is the identified result.It was found that the second FIMF is IIM for AE, FE and the proposed method.However, the difference in values between the first and second FIMFs for the proposed method is more obvious than in the other methods (the corresponding ratios are 1.01, 1.02 and 1.21 for the three kinds of method); this proves that the proposed method has a strong ability to identify faults, relative to the AE and FE methods.The envelope spectrum was conducted for the second FIMF to extract fault information.Meanwhile, to compare the extracted ability of the rolling bearing, the first FIMF identified by the SE method, VMD [53] and CELMDAN combining with kurtosis [24] were employed to extract fault information.Figure 29a,b shows the kurtosis for the VMD and CELMDAN methods.It was found that the largest value was the index of the forth IMF and the second PF. Figure 30 is the corresponding envelope spectrum.It shows that the feature frequency (rotational frequency ( f r = 82 Hz) and the outer raceway fault frequency ( f o = 344 Hz)) were more obvious than those obtained using the other methods.Here, the theoretical f r and f o are 83.33Hz and 344.76 Hz, as obtained by Equation ( 22) (the rotational speed is 9000R/min), respectively.The corresponding amplitude is bigger than the others methods.However, for the AE, VMD and CELMDAN methods, it was very difficult to extract rotational speeds and fault frequencies; this is because the extra information interferes with the normal vibration signal, and the rotational frequency and outer raceway fault frequency were all submerged.Thus, it also indicates that the proposed method is more effective at extracting faults in the outer raceways of rolling bearings.corresponding envelope spectrum.It shows that the feature frequency (rotational frequency ( =82Hz r f ) and the outer raceway fault frequency ( o =344Hz f )) were more obvious than those obtained using the other methods.Here, the theoretical r f and o f are 83.33Hz and 344.76 Hz, as obtained by Equation ( 22) (the rotational speed is 9000R/min), respectively.The corresponding amplitude is bigger than the others methods.However, for the AE, VMD and CELMDAN methods, it was very difficult to extract rotational speeds and fault frequencies; this is because the extra information interferes with the normal vibration signal, and the rotational frequency and outer raceway fault frequency were all submerged.Thus, it also indicates that the proposed method is more effective at extracting faults in the outer raceways of rolling bearings.)) were more obvious than those obtained using the other methods.Here, the theoretical r f and o f are 83.33Hz and 344.76 Hz, as obtained by Equation ( 22) (the rotational speed is 9000R/min), respectively.The corresponding amplitude is bigger than the others methods.However, for the AE, VMD and CELMDAN methods, it was very difficult to extract rotational speeds and fault frequencies; this is because the extra information interferes with the normal vibration signal, and the rotational frequency and outer raceway fault frequency were all submerged.Thus, it also indicates that the proposed method is more effective at extracting faults in the outer raceways of rolling bearings.

Conclusion
This paper develops a novel method to extract the feature frequencies of rolling bearing.Figure 31 is a flowchart of diagnosis process.The following conclusions may be drawn: (1) Our research introduces a method which combines adjacent IMFs into a novel mode function (CIMF) to enhance difference characteristics among IMFs (as shown in step 2).( 2) It proposes a method of combing FFT, PDF and MHD to obtain final intrinsic information mode (FIMF) of minimum number and with physical meanings (as shown in step 3).(3) It introduces a comprehensive evaluate index (DFA and Kurtosis) to identify intrinsic information mode (IIM) to extract the feature frequencies of rolling bearings (as shown in step 4).
The experiment results indicate that this method can minimize interference of spurious mode in original CEEMDAN decomposition, and that it can accurately extract information.Furthermore, the method is also compared with some typical methods.The results indicate that the proposed method has remarkable advantages for the extraction of fault features of rolling bearings.

Conclusions
This paper develops a novel method to extract the feature frequencies of rolling bearing.Figure 31 is a flowchart of diagnosis process.The following conclusions may be drawn: The experiment results indicate that this method can minimize interference of spurious mode in original CEEMDAN decomposition, and that it can accurately extract information.Furthermore, the method is also compared with some typical methods.The results indicate that the proposed method has remarkable advantages for the extraction of fault features of rolling bearings.

Figure 8
Figure 8 is the final combined FIMF and corresponding Fast Fourier transform (FFT).It was found that the identified frequency component of FIMF2 and FIMF3 are almost in accordance with theoretic frequencies in Equation (15) (75 Hz and 10 Hz).The FIMF1 and R are the noisy IMF and

Figure 8
Figure 8 is the final combined FIMF and corresponding Fast Fourier transform (FFT).It was found that the identified frequency component of FIMF2 and FIMF3 are almost in accordance with theoretic frequencies in Equation (15) (75 Hz and 10 Hz).The FIMF1 and R are the noisy IMF and

Figure 8
Figure8is the final combined FIMF and corresponding Fast Fourier transform (FFT).It was found that the identified frequency component of FIMF 2 and FIMF 3 are almost in accordance with theoretic frequencies in Equation (15) (75 Hz and 10 Hz).The FIMF 1 and R are the noisy IMF and residue, respectively; this proves that the proposed method can be used to distinguish the frequency components of a complicated signal and avoid the interference of spurious mode.

Figure 9 .
Figure 9. Flow chart of proposed method.

Figure 9 .
Figure 9. Flow chart of proposed method.

Figure 14 .
Figure 14.Detail difference of frequency component between FC6 and FC7.

Figure 23 .
Figure 23.Envelope spectrum of original signal.

Figure 23 .
Figure 23.Envelope spectrum of original signal.

Figure 23 .
Figure 23.Envelope spectrum of original signal.

Symmetry 2019 ,
11, x FOR PEER REVIEW 15 of 21 (symbol 5), four temperature sensors (for monitoring the temperature of the bearing out ring; symbol 6), two non-contact temperature sensors (for monitoring the temperature of the inner raceway; symbol 7), and a load cell (symbol 8).

Figure 25 .
Figure 25.Bearings with outer raceway fault and corresponding housing.

Figure 26 .
Figure 26.Fault signal of outer raceway fault of rolling bearing.

Figure 25 .
Figure 25.Bearings with outer raceway fault and corresponding housing.

Figure 25 .
Figure 25.Bearings with outer raceway fault and corresponding housing.

Figure 26 .
Figure 26.Fault signal of outer raceway fault of rolling bearing.Figure 26.Fault signal of outer raceway fault of rolling bearing.

Figure 26 .
Figure 26.Fault signal of outer raceway fault of rolling bearing.Figure 26.Fault signal of outer raceway fault of rolling bearing.

Figure 28 .
Figure 28.Comparison of identified IIM with AE, SE, FE and the proposed method.

Figure 28 .
Figure 28.Comparison of identified IIM with AE, SE, FE and the proposed method.

Figure 28 . 21 Figure 29 .
Figure 28.Comparison of identified IIM with AE, SE, FE and the proposed method.Symmetry 2019, 11, x FOR PEER REVIEW 17 of 21

Figure 30 .
Figure 30.Comparison of feature frequencies for different methods.

Figure 29 .
Figure 29.(a) Identifying of IIM with kurtosis method for IMF; (b) Identifying of IIM with kurtosis method for PF.

Figure 29 .
Figure 29.(a) Identifying of IIM with kurtosis method for IMF; (b) Identifying of IIM with kurtosis method for PF.

Figure 30 .
Figure 30.Comparison of feature frequencies for different methods.

Figure 30 .
Figure 30.Comparison of feature frequencies for different methods.

Table 2 .
Bearing parameters in this experiment.

Table 2 .
Bearing parameters in this experiment.

Table 2 .
Bearing parameters in this experiment.

Table 2 .
Bearing parameters in this experiment.