An algorithm for the simulation of faulted bearings in non-stationary conditions

In the field of condition monitoring the availability of a real test-bench is not so common. Furthermore, the early validation of a new diagnostic technique on a proper simulated signal is crucial and a fundamental step in order to provide a feedback to the researcher and to increase the chances of getting a positive result in the real case. In this context, the aim of this paper is to detail a step-by-step analytical model of faulted bearing that the reader could freely and immediately use to simulate different faults and different operating conditions. The vision of the project is a set of tools accepted by the community of researchers on condition monitoring, for the preliminary validation of new diagnostics techniques. The tool proposed in this paper is focused on ball bearing, and it is based on the well-known model published by Antoni in 2007. The features available are the following: selection of the location of the fault, stage of the fault, cyclostationarity of the signal, random contributions, deterministic contributions, effects of resonances in the machine and working conditions (stationary and non-stationary). The script is provided for the open-source Octave environment. The output signal is finally analysed to prove the expected features.

Abstract In the field of condition monitoring the availability of a real test-bench is not so common. Furthermore, the early validation of a new diagnostic technique on a proper simulated signal is crucial and a fundamental step in order to provide a feedback to the researcher and to increase the chances of getting a positive result in the real case. In this context, the aim of this paper is to detail a step-by-step analytical model of faulted bearing that the reader could freely and immediately use to simulate different faults and different operating conditions. The vision of the project is a set of tools accepted by the community of researchers on condition monitoring, for the preliminary validation of new diagnostics techniques. The tool proposed in this paper is focused on ball bearing, and it is based on the well-known model published by Antoni   Rolling bearings, together with gears, are one of the most studied components. They are common components in mechanical design and they allow the relative motion between two or more elements of the machine. Unfortunately, the continuous movement between the parts of the bearing leads to wear phenomena and subsequent failure. The degradation of the bearing conditions can be revealed and monitored analysing the vibration signal produced by the contact among the bearing elements. There are other types of techniques to determine the state of health of the bearings, such as monitoring the temperature or analysing the chemical content of the lubricant; however, the vibration analysis is, de facto, the main technique used in condition monitoring, despite the ease the noise and disturbances may enter into the measurement. So far, thousand of algorithms have been published in the literature trying to reject disturbances and to obtain a clear and telltale signal to assess the health status of the bearing [1]. All these publications usually provide results on both simulated signals and real measurements, more rarely on only one of those. It is a matter of fact that the availability of a real test-bench is not so common, and this is proven by the number of scientific papers validated on few on-line available data centers (e.g. the Case Western University) providing real measurement data. On the contrary, simulated signals are always available, since they are created on the same software for scientific computing used in the post processing. The main advantage of a simulated signal is to avoid the complexity of a real environment, focusing only on the main contributions the developer decided to include. The main drawback is that a too simple model may be too far from reality, making the proposed algorithm not useful. The foundation of a faulted bearing simulation signal is the model proposed by McFadden and Smith [2][3][4]. The bearing is modelled as an epicyclic gear, where the inner ring is the sun gear, rolling elements are the planet gears, the outer ring is the annular gear and the cage is the planet carrier. This simple but powerful model allows the computation of characteristics fault frequencies which are the fingerprints of a damage on the bearing [5]. Moreover, the model takes into account also the modulation effects due cyclic passage of the rolling elements on the load zone. Su and Lin [6] studied the models under variable load due to shaft and roller errors. The ''gearbox'' model for the bearings has a main limitation: the contact among the bearing components is supposed to be a pure rolling contact, while some slippery effect is always present due to the presence of the cage. Ho and Randall [7] proposed to model the bearing fault vibrations as a series of impulse responses of a single-degree-of-freedom system, where the timing between the impulses has a random component simulating the slippery effect. The next fundamental contribution to the modelling of bearings came from the works of Antoni and Randall [8][9][10]. Starting from the work of Gardner [11], Antoni and Randll proposed to model the vibration signal from a ball bearing as a cyclostationary signal, i.e. a random process with a periodic autocorrelation function. Cyclostationarity better describes the effect of slippery and has paved the way for later development.
Most recent developments regard the modelling of the vibration signal in non-stationary conditions [12, 13], i.e. taking into account the speed or load variations in the working conditions of the machine. Unfortunately, as the proposed models have become more detailed, the implementation of the algorithms has become more complex. If the model of McFadden could be easily taught in an introductory course at an engineering school, concepts like cyclostationarity and nonstationary conditions are hardly present in advanced courses at engineering faculties. As a consequence, it could be a gap between the theoretical description of a vibration signal and the algorithm implemented to generate that vibration signal on a computer. A wrong implementation leads to wrong simulated signals used to test diagnostics procedures. In this scenario, the aim of this paper is to provide a detailed step-by-step algorithm for the simulation of the vibration signal provided by a faulted ball bearing. The script is developed in Octave environment, an open source high-level interpreted language, primarily intended for numerical computations and quite similar to Matlab.
The base of this model is the one proposed by Antoni [14] with some improvements. In particular, the model of incipient faults at constant speed has been extended to variable speed applications. In the distributed fault model, the mathematical formulation is completely original and developed by the authors of this paper. Details on the characteristics that the model takes into account will be explained in the next sections. The final goal is to start a discussion with the readers to define a bearing model that can be used as a benchmark, recognized by the scientific community. The paper is structured as follows: Sect. 2 covers the theoretical background of the bearing model and the numerical implementation. Section 3 focuses a numerical example, showing the output results of the proposed algorithm. In Sect. 4 experimental end simulated data are compared to validate the signal model. After the conclusions in Sect. 5, ''Appendix'' lists the script of the algorithm as Octave code.

Theoretical background
At first glance, the vibration signal model of a localized fault in a rolling element bearing could be considered as the repetition of impact forces when a defect in one bearing surface strikes a mating surface, which may excite resonances in the bearing and in the machine. The repetition frequency of these impacts uniquely depends on the defect location, being the defect on the inner race, outer race, or in one of the rolling elements. Even if several resonances can be present in the actual response, for simplicity, it will be assumed in the remaining discussion that only one resonance occurs.
The vibration signal of a localized fault in a rolling element bearing can be reasonably modelled as [9,14]: where h(t) is the impulse response to a single impact as measured by the sensor; q(t) takes into account the periodic modulation due to the load distribution, possible bearing unbalance or misalignment, as well as the periodic changes in the impulse response due to the movement of the faults towards and backwards with respect to the sensor; T is the inter-arrival time between two consecutive impacts; s i accounts for the uncertainties on the inter-arrival time (jitters) of the ith impact due to the necessary random slip of the rolling elements; n(t) gathers the background noise.
Since this work focus the attention on the numerical implementation of Eq. (1), instead of taking into account uncorrelated (white) jitters s i [14], an uncorrelated (white) inter-arrival time difference s iþ1 À s i is used [9]: where r s is the standard deviation and d ij is the Kronecker's symbol. Even if Eq. (1) embodies a well defined harmonic structure, the presence of very slight random fluctuations of the inter-arrival time of consecutive impulses causes the rapidly turns of the signal into a random one. Therefore, weak harmonic components can be located in the lower-frequency range, and a dominating random cyclostationary component can be located in the higher-frequency range (pseudocyclostationary). A detailed theoretical explanation of the frequency content of Eq. (1) can be found in [9,14].
When a localized fault propagates on the surface where it was initiated, a larger area of the bearing becomes involved in the genesis of the vibration signature. In this scenario, no sharp impulses are generated, but the fault signature becomes purely cyclostationary (as opposed to pseudo-cyclostationary) [8,15]. This pure cyclostationary content is the result of a randomly distributed phase, caused by the different position on the rough surface of the rolling elements for every revolution. However, strong periodic components are generated at the shaft periodicity, when the fault only extends over a limited sector of the race. Moreover, if the bearing is highly loaded, a periodic component can be initiated by the bearing stiffness variation due to the changing numbers and positions of the rolling elements in the load zone. The distributed fault vibration signature may be written [8]: where p(t) accounts for the periodic component such as shaft and stiffness variation periodicities and B(t) for the purely cyclostationary content with EfBðtÞg ¼ 0.

Numerical implementation
This work focuses the attention on the numerical implementation of the vibration signal models of Eqs. (1) and (3). In particular, these models are extended to cover generic speed profile of the bearing shaft. In order to include a speed variation, the vibration signal is firstly defined in the angle domain and then transformed back to the time domain according with the chosen speed profile. Let hðtÞ be the rotation angle of a bearing moving race (inner and/or outer). Without loss of generality, in the following the bearing outer race is considered fixed whilst the inner race is rotating. A generic speed profile in the angle domain can be constructed as: where f c is the carrier component of the rotation frequency, f d is the frequency deviation and f m is the frequency modulation. The main terms (f c , f d and f m ) of Eq. (4) can or cannot be angle dependent. Figure 1 depicts an example of Eq. (4) for a case of sinusoidally speed varying profile. Without loss of generality, hereafter it is assumed that at time t ¼ 0 the defect is located at the position h ¼ 0 and it is in contact with a rolling element.
Concerning localized fault in ball bearing, the angle between two consecutive impulses can be easily obtained from the ''gearbox'' model of the rolling element bearing (see Table 1 for the usual bearing fault frequencies), for a inner-race fault: Equation (5) can be used to obtain the angular position of a series of equispaced impulses, i.e. a purely deterministic signal. As stated before, in order to take into account the necessary random slip of the rolling elements a random contribution must be added to Eq. (5). The angle between two consecutive impulses is strictly positive, and so the gamma law is the best candidate; nonetheless when the variance is low with respect to the mean value, the gamma distribution is well approximated by a normal distribution with the Fig. 1 Example of sinusoidally speed varying profile same mean and variance. In this work, the random contribution is taken into account by generating normally distributed random numbers with mean Dh imp and variance r 2 Dh . As the speed profile is defined in terms of rotation angle h [Eq. (4)], the inter-arrival time among the impulses can be obtained by the generated random numbers as: where DT i is the ith inter-arrival time, Dh i is the ith angle between two consecutive impulses randomly generated with mean Dh imp and variance r 2 Dh and f r ðhÞ is the angular dependent rotation frequency.
The results of Eq. (6) are the inter-arrival times of each impulse with the speed profile defined in Eq. (4). These times define the beginning of each impulse response hðt À iT À s i Þ in the time signal itself; such a signal can be obtained in a Matlab/Octave environment as follows: 1. Generate a L point vector filled with zeros, corresponding at times t ¼ l=f s , where f s is the sample frequency in Hz and l is a index ranging from 0 to L À 1, 2. Place 1 at index values obtained by dividing each inter-arrival time DT i by the chosen sample frequency f s , 3. Weight the so generated vector with the weighting function q(iT), 4. Filter the weighted vector with the FFT-based method of overlap-add by choosing as filter coefficients the impulse response function of a SDOF system in terms of acceleration.
Several methods can be found in the literature in order to obtain the impulse response of a SDOF system [7,16]; they deal with the implementation of such a response in the frequency domain and then transform it back in the time domain via the Inverse Fourier Transform. However, this procedure involves the generation of a low pass filter as well as a phase correction [7]. In this work the authors decided to generate the response of the SDOF system to unit impulse in the time domain as: where F is the amplitude of the force exciting the SDOF system, m the system mass, f the damping coefficient, x n the natural frequency in [rad/s] and . The response in terms of acceleration can be simply obtained by a double derivative with respect to time. In this scenario, the numerical derivative does not add high frequency noise inside the signal, because no noise is present in the generated x SDOF ðtÞ.
From the procedure heretofore described, the key point is to find inside the time signal, the index l corresponding to the beginning of the impulse. De facto, l must be an integer number. However by dividing DT i by the selected sample frequency f s , a rational number is usually obtained. Instead of using an interpolation procedure on the time signal itself, the authors decide to rounding the rational numbers to the nearest integers. With this operation, an error is introduced that depends on the selected sample frequency f s (the greater f s , the lower is the error), which affects both mean and variance of the theoretical DT i . Let DT i the value of DT i obtained via the rounding procedure, the error term is: with mean and variance: Finally, the last term of Eq. (1) which deals with the noise component, can be added to the signal by generating randomly distributed number with a given power. The power of the noise can be set with a desired signal-to-noise ratio (SNR), which is a measure that compares the level of a desired signal to the level of background noise. The SNR is defined as: where P signal is the power of the signal without noise and P noise is the noise power. Figure 2 depicts the schema of the proposed procedure. Moreover, in ''Appendix'' an Octave function called bearingSignalModelLocal has been inserted in order to easily implement equation (1).
The same procedure can be efficiently extended for the case of distributed faults in rolling element bearing. This vibration signal model is a mixture of two terms, one deterministic and one purely cyclostationary. Once the speed profile has been defined with respect to rotation angle h, the deterministic part can be described in the angular domain as: where q rot and q stiff are two positive numbers which weigh the amplitude of the deterministic components, whilst s stiff is a geometrical bearing parameter which can be obtained by the ''gearbox'' bearing model as: De facto, in rolling element bearings, the frequency of the stiffness variation is equal to the frequency of an outer-race fault. As done before, by the knowledge of the speed profile, the angle signals can be transformed by a simple interpolation in the time domain. The purely cyclostationary component (B(t)) is a random modulated noise, where the modulation frequency is the fault frequency. Once the speed profile is selected, the modulating function for an inner-race fault (see Table 1 for other types of fault), can be expressed in the angle domain as: where q Fault is a number governing the amplitude of the modulating function and s Fault is a geometrical parameter which can be obtained by the ''gearbox'' model of the rolling element bearing as: Equation (13) can be transformed in the time domain by a simple interpolation and the purely cyclostationary component can be obtained by modulating normally distributed random number with the time domain q Fault ðtÞ. Finally, stationary noise can be added to the signal with a given SNR via the use of Eq. (10). Figure 3 depicts the schema of the proposed procedure. Moreover, in ''Appendix'' an Octave function called bearingSignalModelDist has been inserted in order to easily implement equation (3).  Table 2 shows the vibration signal model parameters.

Numerical example
As stated beforehand, a speed profile has to be generated. The selected speed profile used hereafter in the numerical examples is depicted in Fig. 4, and it deals with a constant rotation frequency of 10Hz modulated at 1Hz with an amplitude of 0.8Hz (see Table 2). From now on both localized and distributed faults in the inner-race of a rolling element bearing are taken into account (see ''Appendix'' for Octave scripts). The mean and variance of the random contribution related to the rolling element slips are set in the angle domain as Dh imp and 0:04Dh imp respectively, that lead to 7:8981E À 3 (the inverse of the fault frequency) and 3:0224E À 07 due to the selected speed profile. As stated in the previous section, the impulse locations in the time domain signal are approximated by a neighbour interpolation that introduces an error term in both the selected mean and variance, which is related to the sample frequency of the time signal itself. In particular, the final mean and variance are 7:8980E À 3 and 3:0245E À 07 showing that the error is negligible in the generation of the vibration signal for the usual sample frequencies. Figure 5 depicts the simulated time signal in case of inner-race localized fault following the data of Table 2. At a first glance, the signal seems strictly deterministic, showing a series of impulse responses (Fig. 5a, b). However, the random slips of the rolling elements turn the signal to strictly random. This effect can be easily seen from the PSD signal. Figure 5c plots the PSD computed with the Welch's method, by using an Hanning window with a 75% of overlap. It is clearly visible from the PSD signal that the noise signal is purely random in nature, in particular the harmonic series related to the repetition of the impulses are strongly masked by the background noise and die quickly.
De facto, Antoni and Randall [9,14] proved that the decay of the harmonic structure strictly depends on the selected variance due to the low-pass filter nature of Eq. 1. In order to highlight the fault frequency cyclostationary analysis has to be carried out. The main signal processing technique in the cyclostationary field is Spectral Correlation Density function (SCD), which depicts the cyclostationary content with respect to the frequency content of the signal. This technique has to be used in case of constant speed, however when the speed is changing a cyclo-nonstationary signal is generated. D'Elia et al. [13,15] were the firsts to explore the order-frequency approach extending the SCD to speed varying signals. Abboud et al. [17] proposed a more rigorous approach to the analysis of cyclo-non-stationary signals. Figure 7a depicts the Order-Frequency Spectral Correlation function (OFSC) for the synthesized signal in case of localized fault. It is possible to see how the order related to the inner-race fault (see Table 2) is highlighted around a frequency region of 6kHz, which is the resonance frequency excited by the bearing impulses. Moreover, the OFSC also highlights the amplitude modulation due to the periodic variation of the load distribution. Figure 6a, b depict the time signal for a inner-race distributed fault with and without noise addiction. It is possible to see how the signal seems strictly random. In particular, even without noise the deterministic component related to the stiffness variation as well as shaft rotation are hidden. Figure 6c highlights the PSD of such a signal, where the random contribution is clearly visible in the medium/high frequency range, whilst the deterministic components are depicted in the low frequency region. Moreover, due to the speed variation, modulation around the bearing stiffness variation frequency can be easily detected. As done before, in order to highlight the fault frequency the OFSC function is evaluated on the simulated signal. Figure 7b plots the result of these operations. The

Experimental validation
In this section, the proposed algorithm is validated on experimental data of faulted bearing. The data are provided by Prof. Gareth Forbes at Curtain University, by Creative Commons Attribution 4.0 International License, through the Data-acoustics.com Database [18]. The provided Matlab files contain radial vibration measurements on the bearing housing of the SpectraQuest Machinery Fault Simulator test rig. The set of measurements contain two files: a known inner and outer race bearing fault, respectively. The validation of the algorithm focuses on the outer race bearing fault case. The measured parameters are: • Radial bearing housing acceleration (m/s 2 ) • Tacho-once per revolution pulse (V) Bearing dimensions and setup characteristics are listed in Table 3. Figure 8 shows the raw data loaded from file and the corresponding spectrum.
First, the raw data is analyzed to characterize the frequency content. Since the bearing data is a cyclostationary signal [14], the cyclic modulation spectrum of the raw data is shown in Fig. 9.
The spectrum is characterized by three main components at Ball Pass Frequency of Outer ring (BPFO) and harmonics (a-axis). Moreover, there is a relevant component at rotational frequency of the shaft (29 Hz) and successive harmonics. Probably there is an imbalance on the shaft, although not reported on the test description. The signal has a resonance band around 2800 Hz (f-axis) and a secondary one around 10,400 Hz. It must be noted that the fault components  10 10 are present at the first resonance only, while the imbalance of the shaft is present on both resonances. The simulation of the faulted bearing will focus on the outer race fault components only, not covering the imbalance effects. From the tacho signal the instantaneous rotational speed in angle is computed (Fig. 10), verifying that the test was done at constant speed with a small fluctuation of the rotational speed.
The speed profile is given as input to the signal model, providing also the bearing information of Table 3 and resonance frequency highlighted in the cyclic modulation spectrum (Fig. 9). In addition to the parameters listed in Table 3, the used model variables are collected in Table 4. It is worth noting that the carrier component of the shaft speed (f c ), the frequency deviation (f d ) and the modulation frequency (f m ) are not necessary, since the instantaneous rotational frequency is directly computed from the tacho signal. The output data are compared with experimental raw data in Fig. 11. The time domain comparison highlights that the signal periodicity is captured by the simulated signal, albeit differences in terms of signal amplitude occurs. However, it has to be underlined that the primary goal of a signal model is to correctly represent the frequency content of the experimental signal, less its amplitude. Finally, Fig. 12 shows the cyclic modulation spectrum of the faulted bearing simulated signal. The characteristic fault frequency and its harmonics are evident, like in the experimental cyclic modulation spectrum in Fig. 9. The spectrum components at the rotational frequency and harmonics are not present since the model focuses on the fault   Table 4.