Fault Diagnosis of Rolling Element Bearings with a Spectrum Searching Method

Rolling element bearing faults in rotating systems are observed as impulses in the vibration signals, which are usually buried in noises. In order to effectively detect the fault of bearings, a novel spectrum searching method is proposed. The structural information of spectrum (SIOS) on a predefined basis is constructed through a searching algorithm, such that the harmonics of impulses generated by faults can be clearly identified and analyzed. Local peaks of the spectrum are located on a certain bin of the basis, and then the SIOS can interpret the spectrum via the number and energy of harmonics related to frequency bins of the basis. Finally bearings can be diagnosed based on the SIOS by identifying its dominant components. Mathematical formulation is developed to guarantee the correct construction of the SISO through searching. The effectiveness of the proposed method is verified with a simulation signal and a benchmark study of bearings.


Introduction
Rotating machinery is widely used in many industrial fields. Fault diagnosis in rotating machinery is of great importance for system maintenance and process automation. In practice, faulty bearings contribute to most of the failures in rotating machinery [1,2]. It is reported that about 40-90% of failures are related to rolling element bearing damage in large and small machines [3]. In particular, inner-race or outer-race flaws dominate most bearing failures. Periodic sharp impulses characterize these faults, and the characteristic frequencies can be computed theoretically. However, those impulses are of low power and are usually buried in noise [4,5]. Moreover, these signals are also usually modulated by some high-frequency harmonic components, which results in a series of harmonics of characteristic frequencies [6]. The above characteristics have presented great challenges in performing fault diagnosis.
In order to realize bearing fault diagnosis, different approaches have been proposed in the time domain, frequency domain, and time-frequency domain over the past few years. In time-domain approaches the time-domain waveform is analyzed directly to extract statistical indices such as the rootmean-square amplitude, skewness, and kurtosis [7]. Once a monitoring feature triggers its predefined damage threshold, the bearing can be considered damaged. However, it is difficult to determine the threshold value of the damage, in particular in different machines. Frequency-domain approaches are usually employed to find the fault's characteristic frequencies via frequency analysis, such as the Fourier spectrum, cepstrum analysis, and the envelope spectrum [8][9][10]. This approach is characterized by its simplicity and intuitive nature for locating the components corresponding to shaft frequency in Measurement Science and Technology Fault diagnosis of rolling element bearings with a spectrum searching method Wei Li 1,2 , Mingquan Qiu 1,2 , Zhencai Zhu 1,2 , Fan Jiang 1,2 and Gongbo Zhou 1,2 the spectrum. Time-frequency-domain approaches including wavelet analysis, the fast Fourier transform (FFT), Wigner-Ville distribution, and Hilbert-Huang transform, etc, which investigate waveform signals in both the time and frequency domain, and can provide more information about the fault signature [11][12][13][14]. However, these time-frequency methods are commonly more complicated than time-or frequency-domain methods in real applications.
In the process of mechanical fault diagnosis, demodulation or filtering techniques are often employed to extract bearing fault features [15,16]. Spectral kurtosis is developed to identify the characteristic frequencies of the bearing, where a filter is designed to obtain the signal with the maximum kurtosis in the spectrum, and then the envelope analysis is usually applied to show the characteristic frequencies [11,17]. The wavelet techniques are widely used to decompose the vibration signals in order to find the most useful filter for fault diagnosis [18][19][20]. Conventional band-pass filters are also applied, whose parameters are optimized through genetic algorithms or adaptive algorithms [21,22]. The above-mentioned methods are relatively complex due to the complicated computations involved. Indeed, spectral kurtosis based methods are to find the resonant frequency band of the vibration signals which contains a train of high-energy harmonics of characteristic frequency, and then this resonant frequency band can be transformed into a low-frequency band through envelope analysis. It is suggested that the characteristic frequency could be identified by finding harmonics with high energies in the spectrum. It was also pointed out that the impulses of bearing faults could be detected when a series of harmonics of the characteristic frequency are identified in the spectrum [6].
This paper proposes a simple method to detect bearings faults by searching the harmonics in the spectrum generated with a FFT algorithm. As discussed above, the impulses generated by bearing faults are usually modulated, and harmonics of the bearing's characteristic frequency exist in the spectrum of vibration signals. Hence we simply search the local peaks in the spectrum of vibration signals. Then the local peaks related to harmonics of a certain frequency are projected onto a predefined frequency grid, such that the so-called structural information of the spectrum (SIOS) is constructed. The SIOS includes two defined indices, which provide the information about the number and the power of harmonics of certain frequency. The dominant and significant components in the SIOS simply correspond to the characteristic frequency of the bearings and therefore the bearing faults can be diagnosed.
The rest of the paper is organized as follows. Section 2 introduces the spectrum searching method and the SIOS. In section 3, this method is applied to diagnose a bearing with a simulated signal, and then this is compared with a benchmark study. The discussion and concluding remarks are provided in sections 4 and 5, respectively.

Spectrum searching method
As illustrated in [6], the faults of rolling element bearings generate impulses and excite frequency resonances of the whole structure between the bearing and the transducer. The low harmonics of the bearing characteristic frequencies are usually masked by other vibration components. The harmonics can be found more easily in a higher frequency range, but higher harmonics may smear over one another. In the case of heavy noise the harmonic series can usually not be directly recognized in the spectrum. The harmonics do exist in the spectrum, but it is not possible to determine the bearing's characteristic frequencies by measuring the spacing of the harmonic series. Hence we propose a spectrum searching method to identify the bearing's characteristic frequencies. The main steps are as follows: • find the local peaks with locally larger amplitudes in the spectrum by searching over the whole frequency range; • construct the SIOS on a predefined frequency grid by projecting components of the spectrum with local peaks onto components of the frequency grid; and • identify the dominant or significant components on the frequency grid according to the SIOS.
In this section we will describe the first two steps. For convenience in the rest of the paper, P(k) is defined as the singlesided power amplitude of the kth frequency component in the spectrum; F(k) is defined as the frequency in hertz corresponding to the kth frequency component; and I(k) is defined to indicate whether the kth frequency component has a local peak or not. The resolution of the spectrum is denoted as ∆ s . The sampling rate of the vibration signal is denoted as F s .

Find the peaks in the spectrum
The spectrum consists of the power amplitude of each frequency component. In order to construct the SIOS, we define the local peaks of the spectrum as follows.
Definition of local peaks. Given three frequency components of the spectrum, i.e. F(k − 1), F(k), and F(k + 1), P(k) is called a local peak if P(k) > P(k − 1) and P(k) > P(k + 1). (1) By searching the spectrum, all frequency components satisfying inequality (1) can be found. Figure 1 gives an example of local peaks. Such a definition may lead to a large number of local peaks. Obviously we are not interested in the frequency components with small amplitudes which may be related to noise. The harmonics of bearing characteristic frequencies usually have relatively larger amplitudes. Hence a threshold is proposed to suppress the influence of noise to some extent as follows: where l and δ are nonnegative constants. If P(k) satisfies inequality (1) and then I(k) = 1, and otherwise I(k) = 0. The threshold in equation (2) varies in terms of k. The first part of the threshold is the moving average of the power amplitudes, and the second part, i.e. δ, is used to control the total number of identified local peaks.

Constructing the SIOS
Assume the bearing characteristic frequencies are within the range from F l to F h in hertz (F l and F h are both frequency comp onents of the spectrum), and the frequency grid is defined as where ∆ G is a selectable positive constant. Therefore the number of frequency components of G is (F h − F l )/∆ G + 1. We use G(i) to represent the ith frequency component of G. ∆ G is the interval between two adjacent components of the frequency grid, and it is suggested to be the resolution of the spectrum or its demultiplier, and we denote it as θ is a positive integer and θ 1.
If a local peak is found on the kth frequency component of the spectrum, then it is projected onto the ith component of G if is an integer with i = 1, 2, · · · , (F h − F l )/∆ G + 1. In other word, we try to project all frequency components of the spectrum onto a frequency grid according to equation (5). Figure 2 illustrates the relation between the spectrum and the frequency grid.
Based on the frequency grid, we define two indexes to represent the SIOS. The first is the number of local peaks projected onto the ith component of G, i.e. N(i), i = 1, 2, · · · , (F h − F l )/∆ G + 1, and where G(i) and F(k) satisfy equation (5).
If F(k) is projected onto a component of G according to equation (5), then the harmonics of F(k) will also be projected onto the same component of G.
The second index is the total power of local peaks pro- This index is used to distinguish the useful signal from the noise. If G(i) corresponds to a characteristic frequency of bearing, E(i) should be relatively large.
In the case that the characteristic frequency of the bearing and its harmonics are multipliers of the ith frequency component of G, all of them will be projected onto the ith component of G. Then the corresponding N(i) and E(i) could be a dominant component in the SIOS. Nevertheless ∆ G and ∆ s cannot be infinitely small. Hence the characteristic frequency of the bearing and its harmonics may not be the frequency components of the spectrum, and they may also not be projected onto the component of G according to equation (5). For instance, a given frequency, i.e. f m , within the range of G described by The jth harmonic of f m denoted by β j f m , may also not be projected onto G, where β j = j. In order to overcome this problem, we reformulate equation (5) as where σ(i, j) > 0 and • is the flooring operator. When inequality (7) is satisfied, we say the jth harmonic of f m is projected onto the ith component of G. Now the key question is how to select σ(i, j). Without loss of generality, f m can be assumed as We would like to find σ(i, j) in inequality (7), such that f m and its harmonics are all projected onto the ith component of G when inequalities (1) and (7) are satisfied.
Considering the limitation of interval of G, we have Assume β j is smaller than α G(i) (assumption 1), then we have according to inequality (8). Clearly f m and its harmonics can be projected onto the same component on G(i) when inequality (7) and equation (9) are satisfied.
Next, we will shown that assumption 1 made on β j and α G(i) can generally be satisfied. According to inequality (6), α G(i) becomes larger when θ is larger, as Since F l f m F h , β j in equation (9) is constrained by We could take a sufficient large sampling length of the vibration signal to obtain a fine spectrum resolution, such that Then based on inequality (10) and inequality (12), we have and consequently β j θ/α G(i) < 1. This means that assumption 1 can be fulfilled if inequality (12) is satisfied. With a defined F l , inequality (12) can always be satisfied with a sufficient large sampling length of vibration signals. Hence assumption 1 can generally be satisfied.
For instance, we are interested in the characteristic frequency in the range from 100 Hz to 200 Hz, and therefore the frequency grid is defined by F l = 100 Hz and F h = 200 Hz. The sampling rate is 12 kHz. According to in inequality (12), assumption 1 can be satisfied when We take 2 15 as the sampling length, and then ∆ s = 0.3662 which satisfies inequality (12). According to inequalities (10) and (11), we have and therefore max i,j (β j /α G(i) ) < 1, which means that assumption 1 is satisfied. The SIOS, i.e. N and E on G, gives the information about the harmonics on the frequency components of G. Index N(i) represents the number of harmonics of G(i) found in the spectrum and index E(i) represents the total power of the harmonics of G(i). In fact the spectrum is interpreted in terms of G(i), N(i), and E(i).
The frequency component in SIOS is treated as a dominant component, when it is significant in N and E at the same time. If E(i) is relatively large and N(i) is small, then G(i) is very likely to be a discrete component of the spectrum without harmonics. If N(i) is relatively large and E(i) is small, then G(i) is likely to correspond to noise.
A flowchart for constructing the SIOS is given in figure 3. Moreover, let us suppose that the sampling frequency of a discrete time series x(i), i = 1, 2, . . . , n, satisfies the limiting condition in inequality (12), then the pseudo-code of the SIOS algorithm for x(i) is given in table 1. Once N and E are derived, they can be applied to perform bearing fault diagnosis. Next, a simulation and two experimental studies are used to describe how to realize bearing diagnosis using the SIOS.

Simulation analysis
In order to demonstrate the proposed method, we use the similar simulated bearing fault signals according to the vibration model in [23]: +n(t).
(13) As described in equation (13), x bear (t) is used to simulate the impulsive signals excited by a fault bearing, where T a is the time interval between two adjacent impulses, and τ i is a random variable to simulate the bearing slip. x ext (t) represents random impulses caused by external knocks on the bearing seat. x har (t) is applied to simulate the harmonic interferences coming from shaft imbalance, electrical components, etc, and n(t) is the white Gaussian noise to simulate the background disturbance.
Considered one resonant frequency, the bearing fault signal is simulated similarly according to that presented in [19] as follows: where β is equal to 900, f m is the fault characteristic frequency ( f m = 74 Hz), F s is the sampling frequency (F s = 12 000 Hz), τ r is a uniformly distributed random number which is used to simulate the randomness caused by the slip, and f is the resonant frequency ( f = 3900 Hz). k − r × F s /f m − τ r 0 is used to ensure the causality of the exponential function. It is known that slipping of the bearing could cause smearing of harmonics. In this case the power amplitudes of the harmonics will be distributed to adjacent components. Such effects will certainly influence the diagnosis results, and therefore we consider τ r ∈ [−0.1, 0.1]. Furthermore, the interferences from shaft imbalance and electrical component are considered in this simulation analysis, and x har (t) is expressed as: where f r represents the rotating frequency of the shaft ( f r = 15 Hz), f e is the power supply frequency ( f e = 50 Hz), and the phases of these two interferences are set as zero, i.e. β 1 = 0, β 2 = 0.
The simulated bearing fault signal, the harmonics of the shaft and electrical component, the random impulse signal, the noise with −10 dB, and the mixed signal are shown in figure 4, where only samples of 0.3 s are displayed. The averaged spectrum of the mixed signal is given in figure 5.
The SIOS is given in figures 6 and 7, where l = 6000, δ = 0.0002, and the frequency grid is selected as [10Hz, 170Hz) with ∆ G = 0.1∆ s . The power amplitudes of the harmonics are significantly reduced due to smearing, as shown in figure 5. However, the proposed method is still effective as shown in figures 6 and 7, where f m is clearly identified. The main reason is that the SIOS can provide abundant information about harmonics even with small amplitudes through the proposed algorithm.
In addition, five kinds of frequency grids are employed to demonstrate the computational efficiency of the proposed method, where experiments were performed on an Intel Core i3-4130 desktop computer with a 3.40 GHz CPU and 4 GB RAM running Matlab R2014a. As shown in table 2, the second column is the average time (50 trials) spent in the SIOS, excluding the time spent in the FFT algorithm. It can be seen that the time consumed increases with the range of the frequency grid. Fortunately, our approach can be of high efficiency when the range of the predefined frequency grid is not large (e.g. the time consumed is 0.03577 s when the frequency grid is [10 Hz, 170 Hz)). However, the frequency grid can be selected in a narrow range because the characteristic frequencies of the bearing can be estimated in advance. Therefore, the proposed approach seems to be a promising tool for fault diagnosis in bearings from the point of effectiveness and efficiency. Next, two kinds of bearing data are applied to further validate its performance. Determine the required sampling length according to inequality (12) and perform FFT.
Select δ and find the peaks according to Eq. (2) and inequality (3). (See section 4.1 for the selection of ).
Project the components of the spectrum onto the frequency grid according to inequality (7) where is determined via Eq. (9).
Compute N(i) and E(i) of the SIOS. added to the shaft and bearing. Four Rexnord ZA-2115 double row bearings were installed on one shaft. An accelerometer was mounted on the housing of each bearing. Vibration data of the bearings were collected every 10 min. The data sampling rate is 20 kHz. For more detailed information about this experiment, please refer to [24], and the data can be downloaded from Prognostics Center Excellence (PCoE) through the prognostic data repository contributed by Intelligent Maintenance System (IMS), University of Cincinnati.
The data set collected from 12 February 2004 10:32:39 to 19 February 2004 06:22:39 is used for further analysis. At the end of the test-to-failure experiment, outer-race failure occurred in bearing 1. Here we take record 510, which is from 15 February 2004 23:22:39, i.e. the very early stage of the fault. Figures 8 and 9 show the waveform and the spectrum of record 510.
The frequency grid is selected as [200 Hz, 300 Hz), and ∆ G = 0.1∆ s . Figures 10 and 11 depict the SIOS of record 510 with the outer-race fault. Although 246 Hz is dominant in figure 11, the number of peaks found around 246 Hz is quite small as shown in figure 10. Hence 246 Hz is not treated as the dominant component on the frequency grid. The situation is similar for 272. 5 Hz, which has a smaller power than 230.4 Hz. Only 230.4 Hz is significant in both figures, and therefore it is identified as the dominant component on the frequency grid. It clearly corresponds to the outer-race fault [25].
As demonstrated by the results, the outer-race fault is detected after around 3.54 running days. We also apply a spectral kurtosis based method to find the characteristic frequency, and the envelope spectrum of the filtered signal is shown in figure 12, in which no evident characteristic frequency can be observed. The proposed method also detects the fault earlier than the method in [25], where the fault is detected after around 3.8 running days.

Case 2: a rolling element bearing benchmark.
The vibration data from the Case Western Reserve University (CRWU) Bearing Data Center [26] are analyzed and compared with the published benchmark results in [27]. The test stand consists of a 2 hp motor (left), a torque transducer/ encoder (center), a dynamometer (right), and control electronics. The test bearings support the motor shaft. With the help of electrostatic discharge machining, inner-race and outerrace faults of different sizes are made. The vibration data are collected using accelerometers attached to the housing with magn etic bases. In this study, the driver end (DE) data, with the sampling frequency 12 000 Hz, are analyzed. The characteristic frequencies of bearings, i.e. the ball-pass frequency of the outer-race (BPFO), the ball-pass frequency of the innerrace (BPFI), the fundamental train frequency (FTF), and the ball spin frequency (BSF), are shown in table 3. And we use f r to denote the rotating frequency.
The frequency grid is selected as [100 Hz, 180 Hz). ∆ G = 0.1∆ s , and l = 10 000 are chosen for all records. δ is set as 0.0002 for most records, where δ is set as 0.002 for records 3005-3008.
We firstly take records 105, 130, and 118 as examples to demonstrate the results, which demonstrate a inner-race fault, outer-race fault, and ball fault, respectively. It is interesting that, all significant components in figures 13 and 15 seem to be harmonics of 0.2 × f r . Similar observations have been found in [27], where some records with a ball fault were analyzed with the envelope spectrum. As stated in [27], it is quite likely that the amount of mean slip in the bearing has adjusted itself to lock onto an exact subharmonic of a dominant frequency such as shaft speed. In  = 1, 2, . . . , L G ; Search local peaks: k+l i=k−l P(i) + δ; 6 if P(k) > P(k − 1) and P(k) > P(k + 1) and P(k) > J th then 7 P(k) is a local peak: I(k) ← 1; 8 end 9 end 10 F I ← (the frequency component of I whose element equals 1); Construct the SIOS: this case BPFI seems to lock onto 5.4 × f r and BPFO seems to lock onto 3.6 × f r , which are both multiples of 0.2 × f r . This is also confirmed by most of the records with an innerrace fault or outer-race fault, whose harmonics of 0.2 × f r are significant on the SIOS and the component corresponding to 5.4 × f r or 3.6 × f r is the dominant one in N(i) and E(i) of the SIOS. Figures 17 and 18 depict the SIOS of record 118 with the ball fault. An interesting feature is that the ball fault records often present evidence of outer-and inner-race faults, which was also discussed in [27]. Most of the marked components in figure 17 are harmonics of 0.2 × f r , among which 107.68 Hz is related to 3. If one of the above rules is fully satisfied, the diagnosis result is marked by Y; and if one of the above rules is partly satisfied, the diagnosis result is marked by P. If no rules are satisfied, the diagnosis result is marked by N. The diagnosis results are compared with the benchmark study in [27], where envelope analysis, spectral kurtosis, and cepstrum techniques were applied. As listed in table B2 given by [27], records 198, 199, and 200 with an outer-race fault cannot be diagnosed, or can only be partly diagnosed, through the benchmark method. We take records 198 and 200 to illustrate the effectiveness of the proposed method, as shown in figures [19][20][21][22]. It is clear that BPFO is the dominant   In addition, most records with a ball fault were marked by N1 in [27], because it was difficult to explain why some BSF components in the envelope spectrum were much stronger than the others. In fact the stronger BSF components in the envelop spectrum often coincide with the BPFI or BPFO comp onents. In other words, once the BSF component conforms with a BPFI or BPFO component, it could be a strong one. That is also the reason why we observe significant components corresponding to BPFI and BPFO in the SIOS of records with ball faults (see figure 17). In this sense the diagnosis results of the records with ball faults can also be marked by Y or P according to the rules defined in [27]. We will not compare those records due to the ambiguous results of the benchmark. Tables 4-6 give the full diagnosis results with the proposed method and a comparison with the benchmark. Records 121 and 188 with a ball fault cannot be fully diagnosed by the proposed method, where the BFP can not be fully identified and 2× BSF is not dominant in the SIOS. With the proposed method satisfactory results also cannot be obtained for records 3001-3004, and the reason is similar to that demonstrated in [27]. (2) According to equation (2), the number of local peaks is determined by the selection of l and δ. Since the average power in different frequency bands may vary, we use l in equation (2) to build a frequency-dependent baseline in order to suppress frequency components with low power amplitudes. The selection of l is not strict, and it can be simply chosen such that 2l + 1 components represent any desired frequency bandwidth, e.g. 50 Hz, 100 Hz. Then the number of local peaks       will be controlled through the selection of δ. If δ is too large, then only a few local peaks can be found and therefore some harmonics may be lost. If δ is too small, the number of local peaks could be large and the searching effort could be very high. Since the harmonics of bearing characteristic frequencies usually have relatively larger amplitudes, in practice we suggest selecting δ such that 0.5% ∼ 3% of amplitudes are treated as local peaks.

On the selection of l and δ in equation
In the benchmark study l is set as 10 000, such that the bandwidth of 114 Hz is used to compute the moving average. And δ is set as 0.0002 for most records except records 3005-3008. Since the spectra of different loads and different faults (with different sizes) are quite different, 0.5% ∼ 2.1% of amplitudes are treated as local peaks in the SIOS of different spectra with the same l and δ. From this point of view,      Y -, -, P the use of the same parameters for those different spectra have already demonstrated that l and δ are not sensitive to the construction of the SIOS if they are within a reasonable range.

On the determination of the frequency grid
If the characteristic frequencies of a given type of bearing are roughly estimated according to the geometrical parameters, the frequency grid can be selected according to the range of estimated characteristic frequencies. If the characteristic frequencies are completely unknown, the range of the frequency grid could be the same as that of the spectrum. In this way the harmonics of all components of the original spectrum are considered in the SISO. Indeed there is no any restriction on the selection of G. The purpose of selecting F l and F h in equation (4) is to reduce the computation effort of the search.

On the influence of slip
Slip always exists in the running process of bearings, and the higher harmonics smear over one another with even a small amount of slip. However the harmonics of the characteristic frequency would present some peaks over the whole spectrum. Due to the slip and noise, some peaks of the spectrum are not related to the fault; and some harmonics of the characteristic frequency may not present peaks.
With the proposed method the influence of the slip and noise is reduced, because local peaks of the whole spectrum are taken into consideration and the sample length of signals is sufficiently large. On one hand the peaks caused by noises can be distinguished by considering N(i) and E(i) together. On the other hand once part of the harmonics of the characteristic frequency present peaks, N(i) and E(i) could give information about the fault. The simulation and experimental results all demonstrate that the proposed method can handle the slip and noise to some extend.

Advantages and limitations
The proposed method is based on a simple searching algorithm, and it is effective in finding the harmonics of the frequency range of interest. Even the harmonics with small amplitudes could be found and projected onto the frequency grid. Although noises and other random impulses may introduce some local peaks unrelated to the faults, the significant components in the SIOS can still be clearly recognized and related to the characteristic frequencies of the bearings.
The proposed method is robust against noise and random impulses. Similar methods, e.g. cepstrum, usually could not find enough periodic components of the signal in the case of heavy noise and other interference. As illustrated in the benchmark study of bearings, the proposed method could provide more information about the harmonics of impulses than the envelope analysis and better diagnosis results can be achieved.
Limitations of the proposed method are as follows: (a) the SIOS can only provide qualitative information, where N(i) and E(i) are not the true values due to the finite spectrum resolution and noise; (b) the computational effort could be high when the range of the frequency grid or the number of local peaks is large; and (c) the method cannot work if the fault has no signatures made of peaks.

Conclusion
A simple and effective method for detecting bearing faults has been proposed based on a searching algorithm. The SIOS of the vibration signals is constructed such that the information of the train of harmonics is clearly represented with N(i) and E(i) on a defined frequency grid. The dominant or significant components of the SIOS correspond to the characteristic frequencies of the bearing. Based on the SIOS the faults bearings are successfully diagnosed. Its effectiveness is verified with simulated and experimental bearing signals as well as benchmark results.