A Comparison of Signal Analysis Techniques for the Diagnostics of IMS Rolling Element Bearing Dataset

1. Introduction

Rolling element bearings (REBs) play a fundamental role within industrial rotating machineries to support moving parts, transmit motion and power and reduce friction among them. As such, they experience constant stress and are seriously prone to suffer from wear which ultimately results in faults. If neglected, a damaged bearing is likely to worsen and affect adjacent components due to physical contact, affecting the whole system functioning and turning into a catastrophic failure, imposing an immediate stop of production or service with deleterious economical consequence.

It’s therefore necessary to prevent such a harmful situation, identifying in advance characteristic signs of its incoming occurrence, without waiting for defects to become visible, i.e. already at an advanced stage, as done in the past with the run-to-break approach, or relying on periodic replacements carried out blindly, as in the time-based maintenance.

Thanks to the use of more powerful sensors and technologies, in recent years predictive maintenance has established itself as an effective discipline to real-time trace the healthy state of mechanical parts and find out the best economical time to execute maintenance intervention. Symptoms of different nature can be taken into consideration, from thermal and pressure fields, noise, IR emissions to current consumption and lubricant-based analysis, but the most common and effective approach is that based on the study of bearing vibrational signals.

Vibration-based condition monitoring is a dynamic research field that is performed with tools susceptible to constant progress; nowadays, besides the mathematical signal-based techniques, new data-driven approaches emerge as promising methods. These latters are particularly appreciated for fault prognosis, which is the third and last phase involved in the process, aiming to predict the future evolution of the vibration signal beyond the present state, and therefore the remaining useful time of the mechanical component; the two phases that precede prognosis are fault detection, which consists in defining the component’s level of healthiness, namely the potential presence of faults, followed by fault diagnosis, which allows the specific sub-component affected by faults and thus their particular nature to be properly identified.

The main disadvantage regarding Machine Learning techniques is that they require huge amount of data to be trained on, are not always generalizable as they would like to be and may have onerous computational demands.

We refer to dedicated books [1,2] and articles [3,4] for an exhaustive description of the principal vibration-based techniques used for diagnostic purposes, including most of those applied in this work.

The goal of this paper is to compare different signal-based analysis methods applying them in the study of the first IMS bearing dataset, which represents a proper paradigmatic test bench to establish their potential. Unlike the other two datasets, involving each the presence of a single fault and so often regarded as an ideal case study to assess the feasibility of new techniques (see for example [5,6,7]), dataset 1 is well-known in the Scientific Community as a demanding challenge to analyze due to the occurrence of three different faults, two of them on the same bearing. However its analysis is not blind because the Center itself provided researchers and industrialists with the diagnostic report of the faults detected at the end of the experiment and since its publication the dataset has been intensively studied.

In [8], all the three datasets are successfully analyzed using the Short Time Fourier Transform (STFT), the Squared Envelope Spectrum (SES) based on the Kurtogram and the Cyclic Spectral Coherence (CSC), but, especially for bearing 4, faults signatures remained barely detectable. [9] draws a comparison between the performance of wavelet decomposition-based denoising and Morlet wavelet filter-based denoising methods applied on dataset 1, attesting the latters as more suitable and reliable to detect a weak signature of mechanical impulse-like defect signals. The study carried out in [10] proposes a hybrid scheme which combines EWT- (Empirical Wavelet Transform) threshold filtering and the Teager energy spectrum method to approach dataset 1, showing better modal aliasing suppression effect and more explicit frequency division segments than the traditional EMD- (Empirical Mode Decomposition) threshold filtering. In [11], Authors test the potential of Hjorth’s parameters as diagnostic tools. In [12] a data-driven method based on Gath-Geva clustering operated on energy spectrum and RMS value is used to identify the degraded state of IMS bearings, proving good classification performance resulting in fault diagnostic accuracy, while in [13] a LAD (Logical Analysis of Data )-based decision model, an innovative supervised learning data mining methodology, is designed and tested to detect faults on the very same bearings. [14] deals with the phase of fault prognosis presenting a Bayesian approach for the RUL (Remaining Useful Life) prediction of rolling bearings: a modified Paris crack growth model is developed and RUL distribution is estimated applying the Metropolis-Hasting algorithm to the Markov chain Monte Carlo method. In [15], the so-called feature engineering, aiming at Prognostics and Health Management (PHM), is conducted applying three metrics (Correlation, Monotonicity and Robustness) to assess the suitability of statistical parameters as health indexes, but only considering the third bearing of the test rig.

The state-of-the-art on the analysis of IMS dataset 1 just summed up could be further extended, but it’s a fact that, probably due to the complexity of the case, these signals have been so far much more frequently studied with data-driven methods (primarily Neural Networks, as in [16,17]) or hybrid approaches ([18]) than with signal-based techniques.

The main contributions the present paper aims to bring are:

• a practical comparison about the efficacy of various signal-based condition monitoring techniques, from the most elementary and general-purpose statistical parameters to complex cyclostationary methods, touching all the three phases involved in predictive maintenance, included signal pre-treatment strategies;
• a detailed study of IMS bearing dataset 1 which echoes the results of prolific research works on the theme, but proves to be innovative for the ample amount of techniques applied successfully, highlighting fault signatures with remarkable clarity, without resorting to data-driven techniques;
• a concrete proof of the validity of two of the most fruitful tools used for the study, which are Cepstrum Pre-Whitening (CPW) and the cyclostationary technique of the Improved Envelope Spectrum (IES), never applied before for the analysis of IMS dataset 1, to the best of the Authors’ knowledge.

The paper is structured as follows: Section 2 contains the description of the experimental test rig from which the dataset is derived and a brief introduction of the techniques used for its analysis, Section 3 gives the most relevant results coming from their application, Section 4 is a recap on the level of efficacy proved by each method and on the major achievements they’ve been able to offer, finally in Section 5 an overall conclusion is drawn.

2. IMS bearing dataset and signal analysis techniques

2.1. Description of IMS dataset 1

As already mentioned, the signals which will be analyzed come from an open access dataset released in 2014 by the NASA Center for Intelligent Maintenance Systems (IMS) [19]; the focus is on the results of the first of three run-to-failure experimental campaigns carried out at the University of Cincinnati, Ohio, USA, with support from Rexnord Corp., Milwaukee, Wisconsin.

The test rig included four rolling bearings installed on a common shaft, set in rotation at a constant speed of 2000 RPM by an AC motor coupled to the shaft via rub belts (Figure 1). The four double row bearings were of type Rexnord ZA-2115, force lubricated by a circulation system that regulated the flow and the temperature and charged with a 26690 N (6000 lbs) constant radial load applied by a spring mechanism. Eight PCB 353B33High Sensitivity Quartz ICP accelerometers were installed on the bearing housing (two sensors for each bearing, mounted orthogonally along the $x$ and $y$ axes), the sampling frequency was set to 20480 Hz. The dataset consists of individual files that are 1-second vibration signal snapshots recorded at specific intervals, every 5 (for the first 54 acquisitions) or 10 minutes, but sometimes subjected to a series of interruptions which makes the time history not continuous (Figure 2).

The experiment was stopped conventionally when the accumulation of debris on a magnetic plug exceeded a certain level, indicating the possibility of an impending failure. The resulting endurance duration was equal to 49 680 minutes (i.e. 34 days and 12 hours), exceeding bearing designed life time, which was more than 100 million revolutions.

The most relevant experimental set-up characteristics and findings of IMS dataset 1 are reported in Table 1. It should be stressed that this dataset is particularly valuable because the bearing degradation was left evolve naturally and was not artificially induced, as it may often happen to speed up the experimental test.

Vibrational data collection was conducted with the National Instruments LabVIEW program thanks to a NI DAQ Card-6062E data acquisition and the analysis that will be presented in this paper was developed with the software platform MATLAB 2021b by The MathWorks [20].

Table 2. Rolling bearing mechanical characteristics and working conditions (from [19]).

Model	Rexnord ZA-2115
Pitch diameter D	71.5 mm (2.815 inch)
Rolling element diameter d	8.4 mm (0.331 inch)
Number of rolling elements per row n	16
Load contact angle ϕ	15.17°
Static load Q	26690 N (6000 lbs.)
Shaft angular velocity ω	2000 rpm

Table 2. Rolling bearing mechanical characteristics and working conditions (from [19]).

Model	Rexnord ZA-2115
Pitch diameter D	71.5 mm (2.815 inch)
Rolling element diameter d	8.4 mm (0.331 inch)
Number of rolling elements per row n	16
Load contact angle ϕ	15.17°
Static load Q	26690 N (6000 lbs.)
Shaft angular velocity ω	2000 rpm

Considering the rolling bearing mechanical characteristics and working conditions of Table 2, it is possible to compute the value of the characteristic frequencies of the faults they’re typically subjected to during their operating life, i.e., the fault signatures, by using the following formulas:

$$\mathrm{B}\mathrm{P}\mathrm{F}\mathrm{O}=\frac{n{f}_r}{2}\left(1-\frac{d}{D}\mathrm{c}\mathrm{o}\mathrm{s}\varphi \right)$$

(1)

$$\mathrm{B}\mathrm{P}\mathrm{F}\mathrm{I}=\frac{n{f}_r}{2}\left(1+\frac{d}{D}\mathrm{c}\mathrm{o}\mathrm{s}\varphi \right)$$

(2)

$$\mathrm{F}\mathrm{T}\mathrm{F}=\frac{f_r}{2}\left(1-\frac{d}{D}\mathrm{c}\mathrm{o}\mathrm{s}\varphi \right)$$

(3)

$$\mathrm{B}\mathrm{S}\mathrm{F}=\frac{f_{r}D}{2d}\left[1-{\left(\frac{d}{D}\mathrm{c}\mathrm{o}\mathrm{s}\varphi \right)}^2\right]$$

(4)

where:

• $\mathrm{B}\mathrm{P}\mathrm{F}\mathrm{O}\left(BallPassFrequency,OuterRace\right)$ is the frequency representing the fault signature on the outer race,
• $\mathrm{B}\mathrm{P}\mathrm{F}\mathrm{I}\left(BallPassFrequency,InnerRace\right)$ is the frequency representing the fault signature on the inner race,
• $\mathrm{F}\mathrm{T}\mathrm{F}\left(FundamentalTrainFrequency\right)$ is the frequency representing the fault signature on the bearing cage,
• $\mathrm{B}\mathrm{S}\mathrm{F}\left(BallSpinFrequency\right)$ is the frequency representing the fault signature on the rolling element,

see the rolling element bearing of Figure 3.

The nominal rotation frequency of the transmission shaft $f_r$ can be directly obtained from the working conditions:

$$f_r=\frac{\omega }{2\pi }=\frac{2000\mathrm{r}\mathrm{p}\mathrm{m}}{2\pi ·60\mathrm{s}/\mathrm{m}\mathrm{i}\mathrm{n}}=33.3\mathrm{H}\mathrm{z}$$

(5)

where the corresponding characteristic frequencies of the faults are reported in Table 3.

2.2. Signal analysis techniques

In this Section, a list and a brief mathematical introduction of the diagnostic techniques used to analyze vibration signals are provided.

Techniques are divided into the three main phases of the condition monitoring process, i.e., fault detection, fault diagnosis and fault prognosis.

Moreover, within the fault diagnosis techniques, signal denoising techniques, filters and cyclostationary techniques are considered.

Fault detection techniques:

• Mean;
• Variance;
• Standard Deviation;
• Peak;
• Root Mean Squared value (RMS);
• Skewness;
• Crest Factor;
• Clearance Factor;
• Shape Factor;
• Impulse Factor;
• Peak-to-Peak;
• Hjorth’s Parameters;
• Detectivity.

Fault diagnosis techniques:

• Short-Time Fourier Transform (STFT);
• Power Spectral Density (PSD);
• Squared Envelope Spectrum (SES);
• Autoregressive Linear Prediction (ALP);
• Time-Synchronous Averaging (TSA);
• Kurtosis (time domain);
• Kurtogram.

Signal denoising techniques:

• Daubechies’ Wavelets;
• Cepstrum Pre-Whitening (CPW);

Filters:

• FIR;
• Adaptive;

Cyclostationary techniques:

• Spectral Correlation (SC);
• Cyclic Spectral Coherence (CSC);
• Improved Envelope Spectrum (IES);
• Wigner-Ville Distribution (WVD).

Fault prognosis techniques:

• Correlation;
• Monotonicity;
• Robustness.

The use of numerous statistical parameters is due to their limited analytical and computational complexity and their general and intuitive meaning for detection purposes, although they just play the role of high-level indicators and need to be followed by more in-depth, specific and elaborate techniques. Table 4 shows the mathematical definition of the most common statistical parameters.

Hjorth’s parameters and Detectivity are statistical time-domain quantities, but their application in the context of signal diagnosis is still unusual. The formers (Activity, Mobility, Complexity) were introduced by Hjorth and Elema-Schönander in 1970 and are commonly used in the analysis of electroencephalography (EEG) signals for feature extraction [21] and in the tactile signal analysis in robotic area [22]. They’re related to the variance of the signal and of its subsequent derivatives, and can be computed with minimum effort and resources, suggesting their use in real-time application for condition monitoring of ball bearings or as a feature array for machine learning techniques.

“Activity” of a vibration signal ($\mathrm{A}\mathrm{c}\mathrm{t}\left(x\right)$) is defined as the variance of its amplitude, therefore it has the same dimension of the signal (the one of the accelerations in the current case) and it’s directly related to its power.

“Mobility” ($\mathrm{M}\mathrm{o}\mathrm{b}\left(x\right)$) is related to the variance of the vibration signal velocity, so it’s associated to the jerk (first time derivative of the signal) and its dimension is expressed as a ratio per time unit.

“Complexity” ($\mathrm{C}\mathrm{o}\mathrm{m}\left(x\right)$) is a dimensionless parameter, connected to the variance of the acceleration of the signal, hence to the snap (second time derivative of the signal and four-time derivative of space).

Moreover, these three parameters can be also derived in the frequency domain, starting from the spectral moments $m_n$ of the time signal $x\left(t\right)$:

$$\mathrm{A}\mathrm{c}\mathrm{t}\left(x\right)={\sigma }^2\left(x\right)=\underset{T⟶\infty }{\lim }\frac{1}{T}{\int }_0^T{x}^2\left(t\right)dt=m_0={\int }_{-\infty }^{+\infty }S\left(\omega \right)d\omega {\left[\mathrm{m}/{\mathrm{s}}^2\right]}^2$$

(6)

$$\mathrm{M}\mathrm{o}\mathrm{b}\left(x\right)=\sqrt{\frac{\mathrm{A}\mathrm{c}\mathrm{t}\left(\dot{x}\right)}{\mathrm{A}\mathrm{c}\mathrm{t}\left(x\right)}}=\sqrt{\frac{m_2}{m_0}}\left[{\mathrm{t}}^{-1}\right]$$

(7)

$$\mathrm{C}\mathrm{o}\mathrm{m}\left(x\right)=\sqrt{\frac{\mathrm{M}\mathrm{o}\mathrm{b}\left(\dot{x}\right)}{\mathrm{M}\mathrm{o}\mathrm{b}\left(x\right)}}=\sqrt{\frac{m_4/m_2}{m_2/m_0}}\left[\right]$$

(8)

where $S\left(\omega \right)$ is the power spectral function, defined as the product of the Fourier transform $F\left(\omega \right)$ and its complex conjugate function $F^{\mathrm{*}}\left(\omega \right)$:

$$S\left(\omega \right)=F\left(\omega \right)·F^{\mathrm{*}}\left(\omega \right)={\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}x\left(t\right)·x^{\mathrm{*}}\left(t\right)dt$$

(9)

and the generic spectral moment of order $n$ can be computed in the form:

$$m_n=F\left(\omega \right)·F^{\mathrm{*}}\left(\omega \right)={\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}{\omega }^n·S\left(\omega \right)d\omega$$

(10)

As for their interpretation in the field of predictive maintenance, Activity and Complexity are proportional to the amplitude of the fault frequency components, while Mobility is expected to reduce with fault progression.

In order to sum up the informative contribution of Hjorth’s parameters into a single value, another quantity, called Detectivity ($\mathrm{D}\mathrm{e}\mathrm{t}\mathrm{e}\mathrm{c}\mathrm{t}$), was introduced in [23] and defined as:

$$\mathrm{D}\mathrm{e}\mathrm{t}\mathrm{e}\mathrm{c}\mathrm{t}={\mathrm{A}\mathrm{c}\mathrm{t}}_{\left(dB\right)}-{\mathrm{M}\mathrm{o}\mathrm{b}}_{\left(dB\right)}+{\mathrm{C}\mathrm{o}\mathrm{m}}_{\left(dB\right)}$$

(11)

where:

$${\mathrm{A}\mathrm{c}\mathrm{t}}_{\left(dB\right)}=10·{\log }_{10}\left(\frac{\mathrm{A}\mathrm{c}\mathrm{t}}{{\mathrm{A}\mathrm{c}\mathrm{t}}_{\mathrm{r}\mathrm{e}\mathrm{f}}}\right)$$

(12)

$${\mathrm{M}\mathrm{o}\mathrm{b}}_{\left(dB\right)}=10·{\log }_{10}\left(\frac{\mathrm{M}\mathrm{o}\mathrm{b}}{{\mathrm{M}\mathrm{o}\mathrm{b}}_{\mathrm{r}\mathrm{e}\mathrm{f}}}\right)$$

(13)

$${\mathrm{C}\mathrm{o}\mathrm{m}}_{\left(dB\right)}=10·{\log }_{10}\left(\frac{\mathrm{C}\mathrm{o}\mathrm{m}}{{\mathrm{C}\mathrm{o}\mathrm{m}}_{\mathrm{r}\mathrm{e}\mathrm{f}}}\right)$$

(14)

The three terms in the right member of (11) are expressed in Decibel, after being normalized computing the ratio between their value on the signal to be diagnosed and a reference value coming from their application to the data collected on the faultless component. Mobility behaves as an attenuator due to its inverse trend with respect to the evolution of the fault, while Activity and Complexity behave as amplifiers.

Due to the particular nature of the vibration signals produced by damaged rolling bearings, which are classified as quasi-pseudo-ciclostationary of the second order, it’s convenient to use dedicated techniques, specifically formulated for their treatment. As previously anticipated, four powerful and complex cyclostationary techniques have been employed for the analysis of the dataset: i) Spectral Correlation $S_{2x}\left(\alpha ,f\right)$, defined in the spectral frequency ($f$) – cyclic frequency ($\alpha$) domain as the bidimensional Fourier transform of the Instantaneous Autocorrelation Function$R_{2x}\left(t,\tau \right)$; ii) Cyclic Spectral Coherence ${\mathrm{C}\mathrm{S}\mathrm{C}\mathrm{o}\mathrm{h}}_{2x}\left(\alpha ,f\right)$, which is a normalized version of the previous technique; iii) Improved Envelope Spectrum ${\mathrm{I}\mathrm{E}\mathrm{S}}_{\mathrm{C}\mathrm{S}\mathrm{C}\mathrm{o}\mathrm{h}}\left(\alpha \right)$, obtained as the integral of ${\mathrm{C}\mathrm{S}\mathrm{C}\mathrm{o}\mathrm{h}}_{2x}$ computed on a limited spectral frequency band $\left[F_1,F_2\right]$, providing information similar to that given by the traditional envelope analysis, but from a different and unconventional perspective, and iv) Wigner-Ville Distribution ${\mathrm{W}\mathrm{V}\mathrm{D}}_{2x}(t,f,\varphi )$, belonging to the class of the time-frequency domain techniques [24]:

$$S_{2x}\left(\alpha ,f\right)={FT}_{t\to \alpha }{FT}_{\tau \to f}\left\{R_{2x}\left(t,\tau \right)\right\}$$

(15)

$${\mathrm{C}\mathrm{S}\mathrm{C}\mathrm{o}\mathrm{h}}_{2x}\left(\alpha ,f\right)=\frac{S_{2x}\left(\alpha ,f\right)}{\sqrt{S_{2x}\left(0,f\right)S_{2x}\left(0,f-\alpha \right)}}$$

(16)

$${\mathrm{I}\mathrm{E}\mathrm{S}}_{\mathrm{C}\mathrm{S}\mathrm{C}\mathrm{o}\mathrm{h}}\left(\alpha \right)=\frac{1}{F_2-F_1}{\int }_{F_1}^{F_2}\left|{\mathrm{C}\mathrm{S}\mathrm{C}\mathrm{o}\mathrm{h}}_{2\mathbf{x}}\left(\alpha ,f\right)\right|df$$

(17)

$${\mathrm{W}\mathrm{V}\mathrm{D}}_{2x}\left(t,f,\varphi \right)={FT}_{\tau \to f}\left\{{Rs}_{2x}\left(t,\tau \right)\right\}$$

(18)

where:

$${Rs}_{2x}\left(t,\tau \right)={\int }_{-\infty }^{+\infty }x(u+\tau /2)x^{*}(u-\tau /2)\varphi (\left(t-u,\tau \right)du$$

(19)

As for fault prognosis, the suitability of the statistical parameters (excluded Hjorth’s ones) as metrics exploitable to predict signal future trend, namely the so called “health indicators” (or “health indexes” – HI), was assessed with the computation of the prognostic values of Correlation ($\mathrm{C}\mathrm{o}\mathrm{r}\mathrm{r}$), Monotonicity ($\mathrm{M}\mathrm{o}\mathrm{n}$) and Robustness ($\mathrm{R}\mathrm{o}\mathrm{b}$), defined as:

$$\mathrm{C}\mathrm{o}\mathrm{r}\mathrm{r}\left(T,X\right)=\frac{\mathrm{C}\mathrm{o}\mathrm{v}\left(T,X\right)}{{\sigma }_T{\sigma }_X}$$

(20)

$$\mathrm{M}\mathrm{o}\mathrm{n}\left(T,X\right)=\frac{\mathrm{C}\mathrm{o}\mathrm{v}\left(r_T,r_X\right)}{{\sigma }_{r_T}{\sigma }_{r_X}}$$

(21)

$$\mathrm{R}\mathrm{o}\mathrm{b}\left(T,X\right)=\frac{1}{k}\sum _{k}exp\left(-\left|\frac{X\left(k\right)-\tilde{X}\left(k\right)}{X\left(k\right)}\right|\right)$$

(22)

in which $X$ is the evaluated quantity (as a statistical parameter), $T$ is the time domain, $\mathrm{C}\mathrm{o}\mathrm{v}$ is the covariance, $r$ represents the rank of a vector, composed by the positions of its elements according to a given sorting, $\stackrel{~}{}$ is the symbol of the smooth value, obtained with the use of a moving-average filter applied on all the $k-1$ elements before the $k$-th one.

Specifically, Correlation (aka Pearson’s Correlation) measures the linearity between $X$ and $T$ and is expressed by a number included in the interval [-1,1] (where 1 states perfect linearity), Monotonicity (aka Spearman Correlation) is defined on the same range and gives information about the continuity of increasing or decreasing of the quantity over time, and Robustness is related to the amount of noise which affects the signal, thereby its value is desirably as low as possible not to compromise the prognostic ability of the parameter.

3. Results

In this Sections, the diagnostic results derived from the application of the aforementioned techniques to the considered dataset will be presented. To avoid redundancy, only the more meaningful elaborated signal between those coming from the couple of orthogonal sensors mounted on each bearing will be considered.

3.1. Fault detection

The computation of statistical parameters is of paramount importance to reveal the presence of an incoming fault, accomplishing the first phase of condition monitoring.

If parameter trend is examined along the time duration of the whole registration (made continuous with the exclusion of the interruptions of acquisition occurred in real time), it can be observed that Mean, Variance and Standard Deviation of bearings 1 and 2 (Figure 4(a) and Figure 5(a)) remain almost constant, with a slight increase just at the very end of the experiment, probably due to the propagation of vibrations produced by other damaged bearings and detected by sensors 1 and 3 because of the physical continuity of the testing, in which all the components are installed on a single shaft.

RMS value, Skewness and Kurtosis (Figure 4(b-d) and Figure 5(b-d)), albeit oscillating, don’t exhibit a global monotonicity too, except for in the last measurements for the reason just mentioned.

On the contrary, the presence of faults on the bearings 3 and 4 can be clearly deduced from the RMS value (Figure 6(b) and Figure 7(b)) and even more from the Kurtosis (Figure 6(d) and Figure 7(d)), which starting from a value around 3 typical of a Gaussian distribution, raises over 70 from the 1831^st second (31^st day of test) for the third bearing and from the 1486^th second (28^th day of test) for the fourth one. Consistent information is provided by other parameters such as Clearance Factor (Figure 10(a) and Figure 11(a)) and Impulse Factor (Figure 10(b) and Figure 11(b)), while Mean, Variance and Standard Deviation (Figure 6(a) and Figure 7(a)) remain mostly stable and Skewness (Figure 6(c) and Figure 7(c)) doesn’t show a clear reduction, as expected in case of fault.

Figure 6. First group of statistical parameters applied to the signal from sensor 5 on bearing 3: (a) Mean (red), Variance (green), Standard Deviation (blue); (b) RMS value; (c) Skewness; (d) Kurtosis.

Figure 6. First group of statistical parameters applied to the signal from sensor 5 on bearing 3: (a) Mean (red), Variance (green), Standard Deviation (blue); (b) RMS value; (c) Skewness; (d) Kurtosis.

Figure 7. First group of statistical parameters applied to the signal from sensor 7 on bearing 4: (a) Mean (red), Variance (green), Standard Deviation (blue); (b) RMS value; (c) Skewness; (d) Kurtosis.

Figure 7. First group of statistical parameters applied to the signal from sensor 7 on bearing 4: (a) Mean (red), Variance (green), Standard Deviation (blue); (b) RMS value; (c) Skewness; (d) Kurtosis.

Figure 8. Second group of statistical parameters applied to the signal from sensor 1 on bearing 1: (a) Peak (red), Crest Factor (blue), Clearance Factor (green); (b) Shape Factor (red), Impulse Factor (blue), Peak-to-Peak (green).

Figure 8. Second group of statistical parameters applied to the signal from sensor 1 on bearing 1: (a) Peak (red), Crest Factor (blue), Clearance Factor (green); (b) Shape Factor (red), Impulse Factor (blue), Peak-to-Peak (green).

Figure 9. Second group of statistical parameters applied to the signal from sensor 3 on bearing 2: (a) Peak (red), Crest Factor (blue), Clearance Factor (green); (b) Shape Factor (red), Impulse Factor (blue), Peak-to-Peak (green).

Figure 9. Second group of statistical parameters applied to the signal from sensor 3 on bearing 2: (a) Peak (red), Crest Factor (blue), Clearance Factor (green); (b) Shape Factor (red), Impulse Factor (blue), Peak-to-Peak (green).

Figure 10. Second group of statistical parameters applied to the signal from sensor 5 on bearing 3: (a) Peak (red), Crest Factor (blue), Clearance Factor (green); (b) Shape Factor (red), Impulse Factor (blue), Peak-to-Peak (green).

Figure 10. Second group of statistical parameters applied to the signal from sensor 5 on bearing 3: (a) Peak (red), Crest Factor (blue), Clearance Factor (green); (b) Shape Factor (red), Impulse Factor (blue), Peak-to-Peak (green).

Figure 11. Second group of statistical parameters applied to the signal from sensor 7 on bearing 4: (a) Peak (red), Crest Factor (blue), Clearance Factor (green); (b) Shape Factor (red), Impulse Factor (blue), Peak-to-Peak (green).

Figure 11. Second group of statistical parameters applied to the signal from sensor 7 on bearing 4: (a) Peak (red), Crest Factor (blue), Clearance Factor (green); (b) Shape Factor (red), Impulse Factor (blue), Peak-to-Peak (green).

Figure 12, Figure 13, Figure 14 and Figure 15 display the trend of Hjorth’s parameters (omitting Activity, which is equivalent to Variance) and Detectivity; the latter has been computed using the first 5 days files (first 51 seconds of acquisition) to obtain the reference values involved in its formulation (11).

Their behaviour agrees with that of the previous parameters: Mobility (Figure 12(a) and Figure 13(a)) and Complexity (Figure 12(b) and Figure 13(b)) of the bearings 1 and 2 have opposite trends than those which would be associated to a fault, while for bearings 3 and 4 Mobility slightly decreases (Figure 14(a) and Figure 15(a)) and Complexity increases (Figure 14(b) and Figure 15(b)).

However, Detectivity is the parameter that offers the clearest signature of fault with a significant increase of its value (Figure 14(c) and Figure 15(c)), proving its efficacy in the field of predictive maintenance as a combination (and an informative summa) of Hjorth’s parameters. The same temporal occurrence of faults previously deduced via statistical parameters is also observed considering Detectivity, with a quite late fault on bearing 3 and an earlier one on bearing 4.

3.2. Fault diagnosis

Beside pointing out the presence of faults, the phase of fault detection is essential to locate temporally their occurrence and then to identify the limited portion of signal to analyze in the second step of the study, fault diagnosis, which is based on the application of more complex and with higher computational effort techniques to determine the exact nature of each fault, i.e., the particular bearing sub-component affected by it.

Since the statistical parameters previously adopted for the fault detection revealed the presence of faults only on the bearings 3 and 4, and not on the bearings 1 and 2, then these last two bearings are not further considered, and the fault diagnosis is carried out only on the bearings 3 and 4.

The signal from the last 49 seconds of acquisition, belonging to the last part of the 33^rd and all the 34^th day of experiment, will be considered for bearing 3 as well as the signal from the 26^th day until the end of the campaign for bearing 4, in the light of the earlier occurrence of faults on it, which tend to gradually reduce their impact on the signal and their impulsive nature due to progressive smoothing of the damaged surfaces.

3.2.1. Short-Time Fourier Transform (STFT)

Before filtering the signal, it’s advisable to study it in its raw version to prevent the loss of important information on account of a premature and inaccurate pre-treatment.

The first applied technique is the Short-Time Fourier Transform, defined on the time-frequency domain, assuming a Hamming window with a length equal to 1024 and without overlap, as suggested in [8].

It’s possible to detect a general amplitude increase of the spectral components of bearing 3 starting from the 12^th second (of the last 49 considered), even more evident beyond the 45^th (Figure 16), with the appearance of the fourth and third harmonics of the BPFI (297 Hz), respectively, still completely invisible in the first days of screening (Figure 17); in the chromatic graph related to bearing 4 signal (Figure 18), a clear display of the fourth harmonic of the BPFO is detectable (although considering 238 Hz instead of 236 Hz, probably because of slipping among parts) as well as the presence of the second and sixth harmonics of the BSF (278 Hz). It’s worth noting that traces of the same harmonic series are visible even in the signal acquired by the third accelerometer on bearing 2 (Figure 19), especially those related to the fault occurring on bearing 3 which is contiguous to it on the shaft.

3.2.2. Power Spectral Density (PSD)

Starting again from the raw signal, the Power Spectral Density technique has been applied, getting information about the frequency distribution of the signal power; its formulation implies the computation of the signal Fourier transform, the corresponding square (power), the division of the result by the sampling frequency $f_s$ and the number of sampled values $N$ (density), a chain of operations which may be equivalently replaced by the evaluation of the periodogram:

$$PSD\left(f\right)=\frac{1}{f_s·N}{\left|\sum _{n=0}^{N-1}{x}_n{e}^{-j2\pi f∆tn}\right|}^{2}with-\frac{1}{2}∆t\le f\le \frac{1}{2}∆t,∆t=\frac{1}{f_s}$$

(23)

In the graphs, shown in linear scale for purposes of clarity, some harmonics of the BPFI for bearing 3 (Figure 20, with harmonics highlighted in magenta) and of the BPFO for bearing 4 (Figure 22, with harmonics highlighted in red as multiples of 238 Hz instead of 236 Hz, as stated before) are already visible, but their detection is strongly affected by the emergence of other frequency lines related to the test rig mechanical structure and hence still present at the beginning of the experiment (Figure 21); specifically, they are located at 160 Hz and their harmonic frequencies (320 Hz and 480 Hz) and at 505 Hz and its striking double (1011 Hz). There are also peaks at very high frequencies, around 4000-5000 Hz and 9000 Hz, probably related to the system natural resonances. It’s curious to observe that the first harmonic of the BPFO is more clearly visible in the spectrum of bearing 3 than in that of bearing 4 due to the cross-talk among contiguous sensors.

3.2.3. Squared Envelope Spectrum (SES)

The Squared Envelope Spectrum is considered a benchmark among the diagnostic techniques and it’s attained through signal demodulation, firstly filtering the signal (in the case in exam, with a low-pass filter with cut-off frequency equal to $f_s$/16 = 1280 Hz), then getting the analytic signal with the Hilbert transform and finally computing the envelope signal as the absolute value of the analytic one and applying the Fourier transform to return to the frequency domain.

It’s known from literature that squaring the signal at the very beginning of this process (i.e., correlating the signal with itself) is beneficial because it emphasizes the modulation and the clearness of information detectable in the final spectrum (while a squaring in the frequency domain would harmfully blur it). This operation is conveniently introduced for the considered dataset which is characterized by inner race and rolling element faults, typically modulated by shaft rotation frequency $f_r$ and cage frequency FTF, respectively (Figure 23); as a matter of fact, in Figure 24 and Figure 25 the BPFI and BSF harmonics are accompanied by the characteristic modulation sidebands.

Specifically, in Figure 24, regarding bearing 3, the even harmonics of the shaft rotation frequency appear distinctly, while the odd ones emerge from noise only when near to the BPFI harmonics, as their sidebands.

Truthfully, the BPFI harmonics are not present directly, but indirectly, through the excitation exerted on those of shaft rotation, namely the spectrum shows frequency lines at 301 Hz and 602 Hz, multiples of $f_r$ = 33.3 Hz, rather than at 297 Hz and 594 Hz, fundamental and second harmonic frequencies of the BPFI; in spite of the closeness of these values, their difference is not negligible or imputable to a mere approximation.

In Figure 25 about bearing 4, at least the first three harmonics of the BPFO can be detected, not subjected to modulation effects, and those of the BSF (centered on the multiple frequencies of 282 Hz instead of the nominal 278 Hz) with sidebands distant from the fault frequencies by the cage frequency, with the left band higher than the right one, unlike the theoretical case (Figure 23).

3.2.4. Time-Synchronous Average (TSA)

While fault frequency components have been already recognized thanks to the techniques applied so far, in order to accomplish a comprehensive and detailed study it’s necessary to consider pre-processing techniques with the aim of denoising and filtering the raw signal to remove the presence of noise and masking components unrelated to faults, such as the harmonic series of the frequencies at 160 Hz and 505 Hz, making fault signatures even more easy to detect.

The following techniques will be applied re-starting each time from the original signal and will be followed by the use of the SES to make a consistent comparison among their effect.

The Time-Synchronous Average is a common time-domain diagnostic tool used as a pre-treatment technique to remove the deterministic components from the signal, bringing to light in the residual part those connected to bearing faults, which are cyclostationary; for the signal under test, the average has been made on time intervals which are synchronous to the shaft rotation for bearing 3 and to the cage rotation for bearing 4. The resulting graphs (Figure 26 and Figure 27) are similar to those already obtained (Figure 24 and Figure 25, respectively), with no added information, making the technique not useful.

3.2.5. Autoregressive Linear Prediction (ALP)

According to the common practice, the Autoregressive Linear Prediction technique is used as a pre-whitening action and it implies at first to get the weighting coefficients of the fourth order autoregressive model (with a form similar to that of a IIR filter), applying Yule-Walker equations and the Levinson-Durban recursive algorithm, then to compute the prediction and finally to derive the residual signal as the difference between the initial and the predicted ones.

In the present case, ALP and SES techniques remove part of the background noise, but the improvement on the information content remains poor (Figure 28 and Figure 29).

3.2.6. Daubechies’ Wavelets

A more advanced technique for denoising purposes is that based on the wavelets, a family of functions with a fixed shape, but which can be shifted and dilated in the time starting from the definition of a root function, i.e., the mother wavelet. Specifically, Daubechies’ orthogonal wavelets of the fifth level have been chosen, determining the extent and threshold values of the blocks with James-Stein method.

The signal processed via Wavelets and SES techniques appears denoised and cleaned up from masking components (Figure 30 and Figure 31); however, the fault spectral lines, albeit visible, are not as emerging as desirable.

3.2.7. Kurtogram and Filters

With the intention of selecting the most impulsive frequency bands, Kurtogram has been used to identify the centre frequency and bandwidth on which to configure the intervention of a filter (Figure 32).

A 10^th order band-pass FIR filter has been chosen, which is able to slightly arise the signal Kurtosis (Figure 33 and Figure 35); nevertheless, after the computation of SES technique (Figure 34 and Figure 36), it appears clear that the filter brought out a frequency band in which the fault components were not predominant, probably because the impulsiveness detected by the Kurtogram was related to phenomena different from the impacts caused by faults on the surfaces or ascribable to the signal transmission path.

The use of an adaptive filter applied in the very same way is not beneficial too.

3.2.8. Cepstrum Pre-Whitening (CPW)

Unlike the mostly unsuccessful trials just reported, the use of the Cepstrum Pre-Whitening proves to be a powerful resource for the application in exam; it’s a relatively new technique which gained appreciation both in literature [25] and in industrial practice [26], but it’s not included into the most traditional and ever invoked set of diagnostics techniques yet. It can be employed in two ways: the first one exploits real cepstrum to perform the “editing”, removing through liftering the peaks of the periodic components which correspond to impulses at particular rahmonics; the second method, far more radical, sets the whole real cepstrum to zero except for the amplitude at the null quefrency, removing both resonances and discrete frequencies, and then resumes to time-domain, after recombination with the unwrapped phase of the original signal, neglected in the real cepstrum. This latter approach is actually easily implemented by dividing the signal Fourier transform by its absolute value and after inverse transforming it back to the time-domain:

$$x_{CPW}\left(t\right)=F^{-1}\left\{\frac{F\left(x\left(t\right)\right)}{\left|F\left(x\right(t)\right|}\right\}$$

(24)

The conceptual and operative simplicity of this technique and its limited computational request are as appreciable as the results it provides. The only drawback is represented by the reduction of various orders of magnitude of all the amplitudes, even those that are useful; however this aspect does not affect its efficacy since the desired components clearly stand out from the others and from noise, which are both significantly mitigated.

Envelope spectra are indeed dominated by the frequency lines of the fault harmonic series: as for bearing 3, Figure 37 and Figure 38 confirm the previously stated conclusion that the inner race fault harmonics are not visible directly, but indirectly through the excitation of those of the shaft rotation frequency which are near to them and of its even ones in general thanks to modulation effects; in the spectrum of bearing 4 (Figure 40), the cage frequency harmonics are undoubtedly present as well as those of the BPFO until the 11^th one and those of the BSF to a lower extent (until the 6^th one, barely visible), accompanied by the modulation sidebands at a distance equal to the FTF; it’s important to note that the strong presence of the harmonics of the shaft and cage rotation frequencies should not be considered a masking factor to be removed, on the contrary it’s just related to faults and ought to be seen as their indirect sign, not present at the beginning of the experiment (Figure 39).

CPW has therefore been capable of deleting all the components unrelated to faults, included the harmonic series based on the fundamental frequencies at 160 Hz and 505 Hz and the peaks at very high frequency coming from structural resonances, leaving the three fault signatures thoroughly visible and unmistakably detectable.

3.2.9. Cyclostationary Analysis

Given that the analysis made so far has allowed to perform the second phase of condition monitoring, namely the fault diagnosis, without taking into consideration the particular nature of the signal in question, for the sake of completeness it’s worth trying to apply some techniques specifically formulated for the study of second order cyclostationary signals, which are the ones actually produced by the damaged bearings. The purpose is twofold: i) validating the previous conclusions by taking a substantially different point of view, ii) testing the efficacy of these tools, which are potentially more powerful but also more complex to understand, interpret and handle.

It’s well known that there is a point of connection between stationary and cyclostationary techniques owed to the mathematical identity between the integral of the Spectral Correlation function made along all the spectral frequencies and the Squared Envelope Spectrum: it seems then reasonable to compute the Improved Envelope Spectrum (IES) in order to verify that it offers information similar to that got with the SES and to understand the differences between them, if any.

Three-dimensional chromatic maps displaying Spectral Correlation and Cyclic Spectral Coherence are not immediate to figure out due to the abundance of information, but seen from above, i.e. projected on the frequency plane, they show a preferential arrangement of the peaks at cyclic frequencies $\alpha$ which correspond to those of faults and their modulating ones (even harmonics of the shaft rotation frequency and odd harmonics near to those of the BPFI for bearing 3 in Figure 41 and Figure 42 and harmonics of the BPFO, FTF and BSF for bearing 4 in Figure 43 and Figure 44) and at spectral frequencies $f$ which are always exact multiples of 160 Hz or 80 Hz (320 Hz, 3120 Hz, 3760 Hz, 3680 Hz, 3840 Hz, 5280 Hz); these last frequencies were already observed as related to system periodicities, although it’s not possible to state surely their physical origin with the information provided on the dataset.

Coming back to the dependence from a single variable, the IES (Figure 45 and Figure 46), computed as the integral of the Cyclic Spectral Coherence on the spectral frequency band [0 10000 Hz] (nearly more like an EES than an IES), puts the stress on the same spectral frequencies which represent fault signatures, even if with a bit less emphasis with respect to Figure 37, Figure 38 and Figure 40: this makes a lot of sense because the IES was applied on the raw signal while, in the previous paragraph, the signal was pre-processed via the radical technique of the CPW.

Applying the previous cyclostationary techniques to excerpts of the original signal all along the experiment time duration (omitting the resulting plots for brevity), the fault history can be traced more precisely than what done before: the inner race fault on bearing 3 occurs almost at the end of the test (33^rd-34^th day), while the outer race fault on bearing 4 can be detected since day 28 and that on the rolling element is still present from around the 26^th day and dominates on the other, but it’s already less visible in Figure 46, probably due to the progressive smoothing of the fault pit caused by friction.

3.2.10. Wigner-Ville Distribution (WVD)

Strictly speaking, WVD has to be ranked among the time-frequency domain techniques, even if it’s characterized by a mathematical formulation which is rather similar to that of the cyclostationary techniques and it’s suitable for the analysis of the same type of signals.

Looking at the respective plots, the presence of cyclic phenomena related to fault development, totally absent in the initial phase, can be detected, albeit more vaguely than before. Figure 47 shows a series of peaks at the spectral frequency of 5280 Hz (still an exact multiple of 160 Hz), spaced from each other by a period (cycle) equal to 30 ms, i.e. the inverse of the shaft rotation frequency f_r. A marked impulsiveness features Figure 48(a) at the spectral frequency of 4000 Hz (25*160 Hz), and its zoom (Figure 48(b)) is dominated by a train of alternated chromatic bands at the spectral frequencies of 240 Hz (1.5*160 Hz) and 480 Hz (3*160 Hz), spaced apart by a time interval equal to the inverse of the BPFO and its double, respectively (about 4 ms and 2 ms).

3.3. Fault prognosis

After fault detection and diagnosis, the third phase required to accomplish the whole condition monitoring of a signal is the fault prognosis, which implies to make a prediction about the signal future trend, beyond the end of real data, with a certain level of probabilistic confidence, with the goal of estimating the component’s Remaining Useful Life (RUL). This operation represents the key point to identify the optimal time to replace a mechanical part, in the trade-off between the cost minimization due to the intervention on the system and the cost of the new component, fully exploiting the service life of the one to remove, on one hand, and the major cost to bear in case of a catastrophic failure occurs on the other, leading to a harmful stop of production, which is absolutely to be prevented.

Prognosis issues are often carried out making use of advanced data-driven techniques belonging to the domain of AI and specifically of Machine Learning, which are programmatically beyond the scope and the intention of this paper. This section is meant for providing with a first glance at that field, identifying the most recommended statistical parameters for prognostic considerations on the IMS dataset.

As stated before, bearings’ end of life was conventionally determined through a lubricant-based analysis; however it’s still feasible to determine how reliable information extrapolated with statistical operators may be about the behavior of the damaged components if brought to the extreme of degradation.

For this reason, the previously mentioned parameters of Correlation, Monotonicity and Robustness, typical indicators used in health management, have been computed on the statistical tools, similarly to what was done in [15], at first on every second of acquisition, then averaging the results on the entire test duration and finally getting a single quantity for each statistical parameter and each bearing as the average of the three indexes, in the same way adopted for Hjorth’s parameters and Detectivity.

Histograms depicted in Figure 49 clearly suggest that Mean and Skewness are not proper instruments for prognostic operations, as they were not so useful for fault detection too, while RMS value, Variance, Standard Deviation and Peak-to-Peak prove to be the most reliable quantities to be used as a solid base for future estimations. It’s interesting to notice that, as may be expected, the parameters evaluated on the signals coming from the two damaged bearings perform better as prognostic tools than when applied on the healthy components.

4. Discussion

From the analysis conducted, the following considerations can be drawn about the efficacy of the involved diagnostic techniques for the study of IMS dataset 1.

• The first statistical approach allowed to execute the fault detection phase, bringing to light the occurrence of different faults on bearings 3 and 4 and suggesting the time excerpts of signal to focus on in the following analysis: in particular, RMS value, Kurtosis and Detectivity proved to be the most informative quantities.
• Diagnostic techniques such as STFT, PSD and SES, even individually, gave a fundamental contribution to identify the exact faults’ nature and their temporal location throughout the experiment, also taking advantage from comparison with faultless signals, at the beginning of the test or coming from healthy bearings; it was then possible to confirm the diagnostic report diffused by the IMS Center itself which attests the presence of:

  -

  a rather late inner race fault on bearing 3, occurring since day 32, which does not show directly, but through the excitation of the spectral harmonics of the shaft rotation frequency; the latter has also modulation effects on the BPFI (even though the interpretation of this fault as located at multiple frequencies of 301 Hz instead of the nominal 297 Hz is still potentially acceptable),

  -

  an outer race fault on bearing 4 detectable for certain from day 28 due to the emergence of various harmonics of the BPFO, without modulating phenomena,

  -

  a rolling element fault on bearing 4, detectable from day 26 with the appearance of the harmonics of the BSF, modulated by those of the cage rotation (FTF), raised themselves by the fault.

• The CPW technique turned out to be a powerful and efficient tool to whiten signal, removing all the masking components and noise, and highlighting fault signatures.
• TSA, ALP and Daubechies’ Wavelets didn’t prove to be sufficiently useful as signal pre-treatment techniques.
• FIR and adaptive filters based on Kurtogram were not able to select the frequency band dominated by fault signs.
• The second order cyclostationary techniques, directly applied on the raw signal, enabled to further confirm the previous conclusions, providing detailed and articulate plots from a different point of view, deeply related to the true nature of the signal under test, qualifying themselves as useful means to use alongside the more traditional ones and justifying their higher complexity.
• The brief insight made into the prognosis domain, to quantify the suitability of statistical parameters to make predictions on the signal future trends, represents an interesting starting point which may be worthy to be further developed.

5. Conclusions

The successful application of traditional as well as more advanced and non-conventional vibration signal-based condition monitoring techniques to analyze the first of the three datasets made available by NASA IMS Center can be considered as a positive example of the efficacy of these mathematical tools to accomplish predictive maintenance of rotating machineries. Besides the possibility to deepen the knowledge of this dataset, the study allowed to make a comparison among techniques’ performance, showing their diagnostic potential, which we believe may be strengthen by the combined use of data-driven techniques, but not completely replaced by them.