Fault diagnosis of rolling element bearing using ACYCBD based cross-correlation spectrum

Rolling element bearings are crucial components in all kinds of rotating machinery. Its fault detection is of great importance, as it ensure the performance of the whole machine. Periodic transient impulses caused by bearing defects are usually submerged in strong background noise which poses a challenge for effective fault feature extraction. To detect bearing faults reliably, a new fault feature extraction method is presented. First, the adaptive maximum second-order cyclostationary blind deconvolution (ACYCBD) is utilized to recover bearing fault related impulses, while the optimal filter length is chosen based on the harmonic significance index (HSI) which quantifies the diagnostic information contained in a deconvoluted signal. Second, cross-correlation is calculated between the teager energy operator (TEO) and the envelope of the deconvoluted signal to further eliminate the irrelevant noise. Finally, fast fourier transform (FFT) is employed to acquire the cross-correlation spectrum and the fault features can be extracted successfully. The performance of the proposed method is verified on both simulation signals and experimental signals acquired from a test rig. The superior abilities of noise reduction and fault detection are shown clearly when compared with some traditional method.


Introduction
Rolling element bearings (REBs) are widely used in rotating machinery and play an important role in modern industries. The working condition of the whole machine is usually closely related to the performance of the REB. Hence, detecting bearing defects at an early stage is of great importance and has drawn considerable attention. Since the transient impulses caused by fault occurance are weak which often masked by background noise and interference components [1], along with the unknown transmission path, the extraction of bearing fault features could be quite difficult.
Vibration based signal processing method is most commonly applied to bearing fault diagnosis as the measured vibration signals contain abundant information concerning bearing healthy conditions. In the past few decades, many effective approaches were proposed, such as the spectral analysis [2], the wavelet transform [3], the spectral kurtosis [4], the singular value decomposition (SVD) [5][6], the empirical mode decomposition (EMD) [7], the variation mode decomposition (VMD) [8], the morphological filter [9][10] and so on. Although the effectiveness of the above methods had been verified, some problems still exist. For example, both proper center frequency and filtering bandwidth are needed for the envelope analysis. The SVD is an efficient denoising method which may still be useless if the original signal has low signal-to-noise ratio (SNR). The EMD is usually used for noise reduction due to its excellent signal decomposition ability, however, the performance of EMD suffers from end effects and mode mixing greatly. The VMD is a newly proposed non-recursively adaptive time-frequency analysis method that can decompose a complex multi-components signal into several intrinsic mode components [11], which has been introduced into bearing fault detection [12][13]. However, the results obtained by VMD may be limited by the selection of mode number K and balancing parameter a, hence, many attentions have been drawn to solve the problem [14][15]. The morphology filter is a nonlinear analysis method which is often utilized to eliminate the noise and enhance the impulse components, but proper selection of morphology operators and structure elements are essential to achieve ideal filtered results.
In real application, the vibration signals are often measured far from the defective points, the complicated transmission path may devote a great impact to the fault detection. As a result, the blind deconvolution (BD) method has been widely used which can recover the impulsive components from the noisy observations. As a typical BD technique, the minimum entropy deconvolution (MED) was originally pioneered by Wiggins [16] to process the seismic signal, which aims to enhance the impulses in raw vibration signals based on the maximization of kurtosis. The MED has been widely used for bearing fault diagnosis and is usually combined with other signal processing techniques, for example, the spectral kurtosis [17], the auto-regressive model [18] and the morphological filter [19]. However, since the noise with single impulse components may also have high kurtosis values, which may influence the filtered results and the failure related impacts can be lost. Moreover, the selection of filter coefficients is also of great importance, which can be solved by optimization algorithms, such as the shuffled frog leaping algorithm (SFLA) [20] and the particle swarm optimization (PSO) [21].
McDonald et al proposed the maximum correlated kurtosis deconvolution (MCKD) to enhance the impulse components in the vibration signals related to a certain fault period [22]. However, the MCKD entails rigorous requires for input parameters, only when all the parameters are set properly, can the superiority of MCKD be highlighted. To overcome this limitation, the cuckoo search algorithm (CSA) that aims to achieve the maximum feature energy ratio is employed to select the optimal filter length and number of shift, which was proved to be effective compared with some traditional methods [23]. In addition, a non-iterative deconvolution approach called multipoint optimal minimum entropy deconvolution adjusted (MOMEDA) was then proposed [24], which shows some advantages over MED and MCKD.
Recently, a BD technique based on the cyclostationary maximization was proposed by Buzzoni et al [25], which is presented for both single-input-single-output (SISO) system and single-input-multi-output (SIMO) system. The BD method based on the maximization of the second-order indicators of cyclostationary (ICS2), called CYCBD, is suitable for bearing fault detection in both time and angle domain due to characteristic of bearing fault signals. The CYCBD was proved to provide excellent diagnostic performance when compared with other BD methods in early fault detection. Although the CYCBD seems to be robustness with different filter length and defect size for the values of the indicator of cyclostationary (ICS), an arbitrary selection of filter length may still be improper. Therefore, to extract the weak fault features from the vibration signals masked by strong noise, a new method named ACYCBD-based cross-correlation spectrum is proposed. The main contribution of this paper is the use of harmonics significance index (HSI) to choose the optimal filter length for the CYCBD and the propose of cross-correlation spectrum to further enhance fault features of the deconvoluted signals. The proposed method uses the HSI to determine the optimal filter length for the CYCBD. With detected optimal filter length, the more impulsive signal can be obtained. Then, the cross-correlation analysis is employed to further reduce the irrelevant noise, which is calculated between the teager energy operator (TEO) and the envelope of the filtered signal. Finally, fast Fourier transform (FFT) is applied to the result of cross-correlation analysis to acquire the cross-correlation spectrum, and the bearing fault features can be extracted successfully.
The remainder of this paper is organized as follows. Section 2 briefly reviews the theoretical background of relevant techniques employs in this study. Simulation verification is conducted in Section 3. Section 4 shows the performance of the proposed method using signals acquired from a bearing test rig.
Conclusion are presented in Section 5.

Basic theory of CYCBD
CYCBD is able to identify impulsive signals exhibiting cyclostationary behavior due to localized defects. In fact, the BD method aims to recover the source signals s0 from a noisy observed signal x.
where L and N represents the length of the discrete signal s and h.
The ICS2 can be defined as follows [25].  (4) where RXWX and RXX denote the weighted correlation matrix and the correlation matrix, respectively and the weighting matrix W reads: By solving the following generalized eigenvalue problem, the maximum ICS2 can be reached [25]: where the maximum eigenvalue λ corresponds to the maximum ICS2. As the initialization of weighting matrix W need to be carried out by guess, hence, an iterative algorithm is required to acquire the maximum ICS2, and the detailed steps are shown as follows: step 1: an auto-regressive model is adopted to initialize the whitening filter h, and the coefficients of the filter are obtained. step 2: the weighting matrix W is estimated using observed signal X and filter h； step 3: the maximum λ and its corresponding h are obtained by solving Eq. (9); step 4: revert back to step 2 using the newly estimated h until convergence.
The blind deconvolution method above is suitable for SISO system. Different from MED and MCKD, the CYCBD can be easily extended to SIMO system, which can enhance ability of the bearing fault detection [25]. Assuming that there are Q number of observed signals xq, Rqq represents the auto-correlation matrix of xq and q qR represents the cross-correlation matrix between xq and q x ( Thus, RXX can be changed as [25]: where Rqq and q qR are placed at the diagonal and off-diagonal position of the matrix RXX, respectively. RXWX can be defined similarly. The deconvolution problems can still be solved by Eq. (9) which returns a filter h consists of filter coefficients related to qth response. The maximum eigenvalue λ still represents maximum criterion.
The qth contribution to the source s can be calculated as: where Xq is composed as X in Eq. (3) which is related to xq, finally, the estimated source signal s can be calculated as:

Adaptive filter length selection method based on harmonic significant index
As the analyzed result of CYCBD for a vibration signal with bearing defect is influenced by the filter length and the longer filter length may lead to more time during convergence, the proper determine of filter length is the key point in the BD process. Just like the traditional method MCKD and MED, the filter length can selected by experience or adaptively chosen within a certain range based on a predefined criterion. In this paper, the second method is chosen and the remaining problem is to determine an appropriate index to select the optimal filter length.
According to Ref. [25], the increment of ICS2 seems to be consistent with the increase of the filter length, so choosing the optimal filter length based on the maximum ICS2 may result in a larger filter length, which may be improper as large filter length can lead to more calculation time.
Harmonic product spectrum (HPS) was originally developed for speech and voice analysis [26][27], which is able to detect the maximum coincidence of harmonics in the frequency domain and can be expressed as: where F(ω) represents the magnitude at frequency ω in the envelope spectrum and K is highest order of harmonics to be considered. However, the original HPS still entails some limitations when using as a criterion to quantify certain frequency components in the vibration signal [28]. To deal with the problem, some modifications are made on the original HPS and an improved harmonic product spectrum (IHPS) is proposed that can be expressed as [28]: (14) where N(ω) is the level of background noise around the frequency ω, which can be calculated using a moving average filter (which is the mean value of ten points' magnitude around the frequency ω in this paper, and 5 points are selected on both sides of ω) and the magnitude of frequency ω should be removed in advanced to achieve an accurate result. P(ω)=F(ω)/N(ω) is the magnitude ratio which can reflect the significant level of a certain frequency component. Thus, H(ω) can represent the significant level of the ω-spaced harmonics above the local background noise, while for bearing fault detection, ω would be set as defect characteristic frequency and its harmonics, and the value of H(ω), which can also be called harmonic significant index, will increase if the bearing fault features are prominent in the envelope spectrum. Hence, HSI is chosen as the criterion to determine the optimal filter length for CYCBD and the highest harmonic order K is chosen as 5 in this paper.

Cross-correlation spectrum
As the estimated source signal s cannot achieve the ideal result that still contains some noise components, which may affect the effectiveness of fault feature extraction. Based on the result shown in Ref. [29], the cross-correlation may amplify the frequency components contained in both signals while the uncorrelated components are attenuated. As a result, the cross-correlation spectrum may provide more prominent harmonics for desired frequency.
Since teager energy operator (TEO) is conveniently calculated which is able to enhance the transient impact components, so as to extract the early weak fault characteristics of the rolling bearings [25]. As a result, in this paper, the cross-correlation is calculated between TEO and envelope of the deconvoluted signal, due to both two signals contain the desired frequency components caused by bearing defects, the cross-correlation may amplify the components related to bearing fault characteristic frequency and attenuate unwanted noise and uncorrelated peaks.
The envelope of filtered signal s is acquired by Hilbert transform: where As s is a discrete-time signal, the TEO can be expressed as: Analogously, the zero-mean TEO can be expressed as: Subsequently, the normalized cross-correlation between these two kinds of signal can be calculated as:  (19) where ϕ and R denotes the conjugate complex number of ϕ and R, respectively.

Summary for the proposed method
To availably diagnose bearing faults from vibration signal, a new fault feature extraction method is proposed with the combination of the ACYCBD and the cross-correlation analysis. Fig.1 depicts the flowchart of the proposed method and the detailed process of the proposed method is described as follows.  Step 1: The vibration signals of REBs with localized defects are acquired by multi-sensors.
Step 2: Calculate the theoretical fault characteristic frequency f using the geometrical parameters of the test bearings and the cyclic frequency is set based on the value of f. Since large filter length would result in more calculation time, especially for a long signal. So, the searching range of filter length is chosen between 40 and 200 in this paper.
Step 3: Normalize the original signals, then some filtered results are acquired using CYCBD with different filter length and HSI values of all filtered signals are calculated. The optimal filter length is determined adaptively according to the maximum HSI value.
Step 4: The original vibration signal is processed by CYCBD with optimal filter length so as to acquire the final filtered result.
Step 5: Cross-correlation is conducted between the TEO and the envelope of the optimal filtered result.
Step 6: Cross-correlation signal is subject to fast fourier transform (FTT) to obtain the cross-correlation spectrum, and the defective features of a REB are acquired.

Simulation verification
A simulation example for bearing with inner-race defect is built in this part to verify the effectiveness of the proposed method for the extraction of fault features. When a localized defect is existed on the rolling element bearing, its vibration signals are presented as periodic transient impulses. However, these impulses are usually masked by interference components, such as the harmonic components and white noise. Hence, the model of simulation signal can be expressed as follows [30]:  (20) where Ai is is adopted to simulate the amplitude modulation which cycle is 1/fr, and fr=42Hz; fn=3200 and r=0.05 denote the natural frequency of the system and the damping coefficient, respectively; The repetition period of the periodic impulses is T=1/185s which means the fault characteristic frequency is f=185Hz; ti represents uncertainties on the arrival time (jitters) of the ith impact (where E{ti}=0 and E{titj}=δijσt, E{·} represents the expected value, δij represents the Kronecker symbol and σt represents the standard deviation of the sequence ti, in this simulation, σt is set to be 0.04T). The harmonic B(t) is employed to simulate the presence of interferences in the signal while f1 and f2 are 50Hz and 100Hz, respectively. The sampling frequency in simulation is 32768Hz and the signal is simulated for 1 second. zero-mean random noise n1(t) with standard deviation of 1 is simulated by MATLAB based on the function randn (1, n) in advance, while n(t)=3*n1(t) is used in the simulation. spectrum. It can be seen from Fig. 2(a) that periodic impulses are submerged by strong background noise which can hardly be observed, while the envelope spectrum shown in Fig. 2(b) is also greatly influenced by the noise and only the fault characteristic frequency is able to be detected. Therefore, the proposed method is applied to the simulation signal. The HSI values of the filtered signals over the selected range of filter length are calculated and presented in Fig. 3(a). The maximum HSI value can be achieved at the position L=175 in the HSI curve. The filtered signal using CYCBD with the optimal filter length is shown in Fig. 3(b), the cross-correlation spectrum described in section 2.3 is displayed in Fig. 3(c) and the fault characteristic frequency along with four of its higher harmonics are extracted for the fault judgement on the inner-race. Moreover, the interference of background noise is almost eliminated.
For comparison, the envelope spectrum of the optimal filtered signal is shown in Fig, 3 To further demonstrate the validity of the proposed method, first, two other blind deconvolution approaches named MED and MCKD are used. In this paper, for both the two comparison methods, the filter length is selected equal to the optimal length of the CYCBD, and for the MCKD, the period of deconvolution is determined according to the fault characteristic frequency and sample frequency, the number of shift is chosen to be M=3. Analyzed results using this two method are presented in Fig. 4 and Fig.   5. Fig. 4(a) and Fig. 5(a) shows the filtered signal using the MED and the MCKD respectively, where no obvious impact components can be recognized, the corresponding envelope spectra are exhibited in Fig. 4(b) and Fig. 5(b). In the envelope spectrum, the fault characteristic frequency and its harmonics are unable to be observed due to the influence of heavy noise. Moreover, inappropriate selection of the parameters in the two methods may also account for the failure of fault feature extraction. Based on the analysis above, it can conclude that the proposed method is effective in fault feature extraction from heavy noise and the cross-correlation analysis is very important for noise reduction, meanwhile, proper selection of filter length and other parameters for the denconvolution approaches is also of great essential. Then the classic Fast Kurtogram proposed in Ref. [31] is applied to the simulation signal and results are displayed in Fig. 6. As can be seen from Fig. 6(a), the central frequency and bandwidth of the optimal analysis frequency band can be determined as 3072Hz and 682Hz (red dotted line) whose square envelope spectrum is exhibited in Fig. 6(b) where only the fault characteristic frequency is able to be observed, thus, the Fast Kurtogram cannot achieve the ideal results when analyzing the simulation signal, which shows the superiority of proposed method for the extraction of impulsive components from the vibration signal with heavy noise.

Application to experimental signal
The experimental signals acquired from a test rig is employed to further demonstrate the validity of the proposed method. The test rig shown in Fig. 7 consists several main components which are bearing support structure, main shaft, experimental bearing, lubricating oil system, servo-driven motor, radial loading device, axial loading device and control system.

Fig. 7. The test rig
As sensors may be impossible to be mounted near the defective bearing in practical applications, the fault diagnosis using the signals acquired at a long distance may be more meaningful. Hence, two measuring points (shown in Fig .7) are selected away from the experimental bearing to collect the vibration signals.
These signals may contain various noise and interference components which may make it more difficult for the fault diagnosis of REB.
The type of test bearing is 6010, the outer and inner diameter of the bearing are 80mm and 50mm respectively, thus the pitch diameter is D=65mm; the diameter and number of rolling elements are d=9mm and Z=13; the contact angle is α=0°. The experimental signals are collected by the B&K vibration test system, the models of the acceleration sensor and the signal conditioning and acquisition module are B&K4354B-004 and B&K3053-B-120, respectively. During the experiment, the rotating speed is set as fr=3000r/min. The sampling frequency of the data acquisition device is 32768Hz. The fault characteristic frequencies for different localized defects can be calculated by the following equations: where fi and fo denote the ball pass frequency of the inner race (BPFI) and the ball pass frequency of the outer race (BPFO). Hence, the theoretical BPFI and BPFO are 370Hz and 280Hz, respectively.

Bearing inner-race defect identification
Since the the fault detection ability of CYCBD can be improved using vibration signals collected by multi-sensors [25], in the experimental verification, vibration signals collected at point 1 and 2 are used together for the estimation of source signals. Fig. 8(a) shows the time domain waveform of signal with inner-race defect measured at point 1, the repetitive transient impulses are not obvious due to the existence of noise. The envelope spectrum in Fig .8(b) is also masked by noise, where only two harmonics related to BPFI are able to be observed. Fig. 8 the time domain waveform of signal measured at point 2 and its envelope spectrum is shown in Fig. 8(d), where the harmonics are even harder to identify. Therefore, the interference is heavier for signals collect at point 2, and the envelope analysis isn't enough if applied to the experimental signals directly. Then, the proposed method is applied to the experimental signals. As mentioned in Fig. 9(a), the optimal filter length should chosen to be 177 and the filtered signal is displayed in Fig. 9(b) where the presence of noise is reduced obviously. The cross-correlation spectrum achieved by ACYCBD is provided in Fig. 9(c). As can be seen from Fig. 9(c), the fault characteristic frequency and four of its harmonics are prominent, while the interference of noise is reduced significantly.
Moreover, the side-bands modulated by rotating frequency around the BPFI and its harmonics are able to be observed. All the features above clearly indicate the existence of a defect on the inner race.
Similar to simulation analysis, the envelope spectrum of the optimal deconvoluted signal is presented in Fig. 9(d) which also shows the fault features clearly. However, noise and some other irrelevant frequency components are more obvious than those in the cross-correlation spectrum. (d) the envelope spectrum of (c).
What's more, filter length corresponds to lower HSI value (N=114) is employed for analysis and results are shown in Fig. 9(e) and (f). Fig. 9(e) shows the filtered signal and the noise is obviously reduced, but only the fault characteristic frequency and two of its harmonics are successfully extracted in the envelope spectrum, which means a proper selection of filter length is of great significance.  inner-race fault can be found in the envelope spectrum presented in Fig. 11(b), but some are not prominent.
When MCKD is used to analysis the measured signal at point 2, the filter signal and its corresponding envelope spectrum are presented in Fig. 11(c) and (d), only four harmonics related to inner-race fault are able to be detected and the noise components are more obvious than those in Fig, 11(b). Comparing Fig.   10(b), Fig. 10(d), Fig. 11(b) and Fig. 11(d) to Fig. 9(c), it reveals that the proposed technique can largely extract the fault related components from the raw vibration signal and eliminate the influence of noise to enable a more accurate detection of bearing fault, which demonstrate the effectiveness of proposed method.
What's more, the Fast Kurtogram [31] is applied to the measured signal, Fig. 12

Bearing outer-race defect identification
Analogously, bearing vibration signals with outer-race defect are also collected. Then, the proposed method is introduced to the experimental signals and results are provided in Fig. 14.
Based on the HSI curve presented in Fig, 14(a), the maximum HSI value can be achieved at the filter length 109 and the deconvoluted signal is displayed in Fig. 14(b) where the noise is removed significantly. Fig. 14(c) exhibits the cross-correlation spectrum of the filtered signal where the fault characteristic frequency and its higher harmonics are very prominent, which indicates the existence of a defect on the outer race. Fig. 14(d) provides the envelope spectrum of the optimal filtered signal, and the fault features are also extracted successfully. This may result from that the outer-race fault detection is easier than inner-race fault detection in the same situation. However, the presence of noise can still be observed in Fig. 14(d), which cannot achieve the optimal result. What's more, another filter length N=110 with lower HSI value is employed and the results are displayed in Fig. 14(e) and (f). The 3rd to 5th order harmonics of the fault characteristic frequency aren't very prominent in the cross-correlation spectrum, which reveals that a improper selection of filter length is conducted. (d) the envelope spectrum of (c).
For further comparison, the MED and MCKD is also applied to the measured signals. Fig. 15(a) displays the deconvoluted result of signal measured at point 1 using the MED, and its envelope spectrum is provided in Fig. 15(b), which extracts the fault characteristic frequency and its harmonics, the noise has only a little influence. Fig. 15(c) provides the filtered results of signal measured at point 2 using the MED, and the corresponding envelope spectrum is presented in Fig. 15(d) where the influence of noise is more obvious and the features are not as obvious as those in Fig. 15(b). Fig. 16 displays the analyzed results of the measured signals with outer-race fault using the MCKD. When signal at point 1 is used for analysis, the envelope spectrum presented in Fig, 16(b) shows prominent fault characteristic frequency components but the 5th order harmonic are unable to be observed. Fig. 16(c) and (d) display the analyzed results of signal at point 2. As signal at point 2 contains more noise than signal at point 1, in the envelope spectrum, only 4 spectral peaks related to outer race fault are obvious, the 5th order harmonic is not prominent while the 6th order and 7th order harmonics are unable to be detected.
The comparison above indicates that the MED and MCKD may result in improper diagnosis when a bearing defect signal is contaminated by noise, and a randomly selection of parameters can also lead to the decrease of the ability for fault feature extraction, which further demonstrates the validity of proposed method for bearing fault detection from vibration signals with strong background noise.
Moreover, the Fast Kurtogram [31] is also employed. Fig. 17(a) displays the Kurtogram of signal measured at point 1, which indicates an optimal band with central frequency and bandwidth of 13653Hz and fault characteristic frequency and three of its higher harmonics are able to be observed, the 5-order harmonic is not prominent while the 3-order and 4-order harmonics are unable to be detected. Fig. 17(c) displays the Kurtogram of signal measured at point 2, the optimal band is determined with central frequency and bandwidth of 15680Hz and 128Hz. The square envelope spectrum displayed in Fig .17(d) surely reveals the wrongly selection of optimal frequency band. The comparison above indicates that the Fast Kurtogram may result in improper diagnosis when a bearing defect signal is contaminated by noise, which also demonstrates the validity of proposed method for bearing fault detection from vibration signals with strong background noise.
(a) (b) (c) (d) Fig. 16. Analyzed results of experimental signals with outer-race defect using the MCKD; (a) the filtered result using signal measured at point 1; (b) the envelope spectrum of (a); (c) the filtered result using signal measured at point 2; (d) the envelope spectrum of (c).

conclusions
This paper presents a new fault diagnosis method for REB based on ACYCBD and cross-correlation spectrum. In this approach, the ACYCBD is employed to recover impulsive patterns related to bearing defects, and the optimal filter length is determined adaptively using the maximum value of HSI as a criterion.
The cross-correlation analysis is then employed to further attenuate the uncorrelated background noise, and bearing fault features can be extracted by the cross-correlation spectrum. The effectiveness of the proposed technique is verified using both simulation and experimental signals. Results show that the proposed method can largely remove the presence of noise and provide prominent bearing fault characteristic frequency components for the fault detection even when the signals are masked by strong background noise, which also shows advantages over some other traditional fault detection methods. However, the vibration signals employed for verification is limited in this paper, hence, the proposed method should be applied to analyzing more vibration signals with bearing defects in future work, especially the signals acquired from real machine.