<figure><img src="image_1.png" /></figure>

HAL
open science

# A Novel Feature Extraction Method for the Condition Monitoring of Bearings

Abdenour Soualhi, Bilal El Yousfi, Hubert Razik, Tianzhen Wang

► To cite this version:

Abdenour Soualhi, Bilal El Yousfi, Hubert Razik, Tianzhen Wang. A Novel Feature Extraction Method for the Condition Monitoring of Bearings. Energies, 2021, 14 (8), pp.2322. 10.3390/en14082322. [hal-03682964]

HAL Id: hal-03682964

https://hal.science/hal-03682964v1

Submitted on 11 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Article

# A Novel Feature Extraction Method for the Condition Monitoring of Bearings

Abdenour Soualhi <fnref n="1" />*, Bilal El Yousfi ¹ID, Hubert Razik <fnref n="2" />,³ID and Tianzhen Wang ³ID

**Abstract:** This paper presents an innovative approach to the extraction of an indicator for the monitoring of bearing degradation. This approach is based on the principles of the empirical mode decomposition (EMD) and the Hilbert transform (HT). The proposed approach extracts the temporal components of oscillating vibration signals called intrinsic mode functions (IMFs). These components are classified locally from the highest frequencies to the lowest frequencies. By selecting the appropriate components, it is possible to construct a bank of self-adaptive and automatic filters. Combined with the HT, the EMD allows an estimate of the instantaneous frequency of each IMF. A health indicator called the Hilbert marginal spectrum density is then extracted in order to detect and diagnose the degradation of bearings. This approach was validated on two test benches with variable speeds and loads. The obtained results demonstrated the effectiveness of this approach for the monitoring of ball and roller bearings.

**Keywords:** signal processing; empirical mode decomposition; Hilbert transform; bearing fault; vibration analysis; feature extraction

<figure><img src="image_3.png" /></figure>

check for
updates

**Citation:** Soualhi, A.; El Yousfi, B.; Razik, H.; Wang, T. A Novel Feature Extraction Method for the Condition Monitoring of Bearings. *Energies* **2021**, *14*, 2322. https://doi.org/10.3390/en14082322

Academic Editor: Maciej Sułowicz

Received: 20 March 2021

Accepted: 15 April 2021

Published: 20 April 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_4.png" /></figure>

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

## 1. Introduction to Bearings

Bearings are mechanical components used to provide a mobile connection between two rotating elements. This system increases the rotation speed while reducing the friction and the transmission of forces by replacing sliding with rotational motion. Bearings are therefore found where rotation is necessary (cars, electric motors, tools, etc.).

Studies conducted on induction motors have classified the most common fault mechanisms according to bearing faults, stator faults, rotor faults, and other faults (see Figure 1). These studies were carried out by the engine reliability group of the IEEE-Industry Applications Society (IAS), which studied 1141 motors, and the Electrical Energy Research Institute (EERI), which studied 6312 motors rated above 148 KW.

<figure><img src="image_5.png" /><figcaption>**Figure 1.** Distribution of faults in induction motors [1,2].</figcaption></figure>

Additionally, as shown in Figure 2, studies dealing with faults located in power transmissions based on gearboxes have demonstrated that one of the most susceptible

<footnote n="1">¹ Laboratory LASPI EA-3059, University of Jean Monnet, 42100 Saint Etienne, France; bilal.elyousfi@univ-st-etienne.fr</footnote>

<footnote n="2">² Laboratory Ampère UMR 5005, University of Lyon, 69007 Lyon, France; hubert.razik@univ-lyon1.fr</footnote>

<footnote n="3">³ Logistics Engineering College, Shanghai Maritime University, Shanghai 201306, China; tzwang@shmtu.edu.cn</footnote>

<footnote>* Correspondence: abdenour.soualhi@univ-st-etienne.fr</footnote>
components to fail is a bearing [3]. Beyond manufacturing or assembly faults, the causes of damage are numerous and give rise to anomalies with several degrees of severity. It is therefore logical, from an industrial and scientific point of view, to focus diagnostic efforts mainly on bearings.

<figure><img src="image_6.png" /><figcaption>**Figure 2.** Distribution of faults in gearboxes.</figcaption></figure>

There exists in the literature a wide variety of signal processing tools used in bearing diagnostics. These methods can be classified into three categories based on one, two, and three dimensional indicators, respectively.

- One-dimensional indicators

These methods are based on statistical indicators that quantify the information contained in the signals. Monitoring the evolution of these indicators provides more or less precise and complete information about the health of bearings. These indicators are mainly based on the calculation of the rms value, the crest factor, and kurtosis. The rms value was one of the first indicators used in the industry. Excessive variation in the rms level usually means a change in health status and therefore possible failure. One of the major drawbacks of using the rms is that it can give a late alarm, especially in the case of bearing faults where the variation in the signal due to the appearance of the fault is masked by other components of high amplitude (see Figure 21). The crest factor enables fault detection by measuring the ratio between the maximum value of the signal modulus (crest value) and its rms value. Other indicators have been developed based on the crest factor, such as the K-factor, by multiplying the peak value by the rms value or the peak-to-peak value, measuring the difference between the amplitudes of the upper and lower peaks. Kurtosis is derived from the statistical moment of order 4 and is defined as the ratio between the mean value of the signal raised to the power of 4 and the square of its variance. It is a statistical quantity used to analyze the flattening of a distribution and therefore to observe the shape of the signal. For normal bearing operation, kurtosis is close to 3 and increases dramatically as soon as pulses occur due to the occurrence of a fault.

- Two-dimensional indicators

Two-dimensional methods are based on the signal transformation in the frequency domain or other auxiliary domain. We can cite the Fourier transform, cepstrum, and spectral kurtosis [4,5]. The Fourier transform analyzes the vibration signals of a bearing into a multitude of elementary components of "amplitude-frequency" characterizing the frequencies of a bearing to be identified and therefore its fault frequencies or amplitudes. Cepstral analysis is defined as the inverse Fourier transform of the decimal logarithm of the Fourier transform. Like the Fourier transform, cepstral analysis makes it possible to highlight the periodicities of a time signal. However, the proper use of this tool always requires combining it with other techniques such as the Fourier transform. This is in particular due to its high sensitivity to background noise. Spectral kurtosis (SK) has been proposed by [6]. This tool makes it possible to identify the best demodulation band from the signal and without the need for a reference signal. In a complementary study, authors in [7] discussed the application of SK to the vibration maintenance of rotating machines, providing a robust
means of detecting incipient faults and designing optimal filters to improve the mechanical
signature of these faults.

- Three-dimensional indicators

Examples of three-dimensional methods are wavelet analysis, spectrograms [8], scalo-
grams, Wigner-Ville distribution [9], and kurtograms [10], which aim to display a signal in
time-frequency domains. These methods generally suffer from a need for strong expertise
on the part of the user, both for their application and for their interpretation. The use of
these tools in the diagnosis of rotating machines is hampered in practice by their high
computational cost, especially in real-time applications.

Authors in [11] listed the main challenges reported by past researchers of bearing
condition monitoring, such as the signal background noise, the interference of the bearing
signal with other components of the studied system (gears, shafts, etc.), the speed and load
fluctuations, etc. To overcome these issues, several approaches were proposed by later
works, including wavelet transform-based techniques [12], the Hilbert-Huang transform
technique (HHT) [13], empirical mode decomposition (EMD), and ensemble empirical
mode decomposition (EEMD) [14,15]. Additional signal analysis techniques have also been
proposed for the condition monitoring of bearings, such as order tracking methods [16],
cyclostationary based methods [17], and morphological signal processing techniques [18].

As we can see, the EMD is well-known in the literature [19–25]. The novelty of the
proposed paper concerns the improvement made to the EMD and its use in the filtration of
the vibration signal in order to extract a new health indicator called the Hilbert marginal.
Indeed, the contributions of the paper can be listed as follows:

- Resolution of the problem of the edge effect. The classic method for the determination of envelopes is the use of a cubic spline function. In this method, monotony between successive points on the same line is not respected, and the tangential criterion is also not respected. If the approximate curve passes through an extrema (max or min), this curve does not change its meaning (tangential criterion), but continues in its direction up to a precise level and then changes meaning. The problem here is that the first and last samples of the signal are not minima and maxima. So, there will be two intervals left with only one point to calculate the interpolation function, which is not possible. To solve this problem, we proposed extending the extrema sets by adding data at the start and end of the given interval.

- Need for a stop criterion before a perfect IMF. The objective of the sifting process is to extract an IMF from the residue left by the previous iteration of the EMD algorithm. We will see in Section 2 that the definition of an IMF will reside in two points: the difference from at most one between the number of zeros and the number of extrema, and the null local average. The first condition stems from the second, so it is necessary to focus on it to find an adequate stopping criterion. Indeed, if the sifting is processed too much, low frequencies will be captured in the first IMFs, while if the sifting is not processed enough, the high frequencies will spread in all IMFs. A stopping criterion has been proposed in this paper to solve the problem.

- The creation of an adaptive tool for the decomposition and automatic filtration in vibration signals. Indeed, one of the limits of the wavelet approach is the need to predefine the basic functions to decompose the signal. The use of the EMD as a sub-band, local, and self-adaptive decomposition method for the analysis of non-stationary signals therefore seems to be the best approach for the adaptive denoising of the vibration signals. Indeed, the EMD is fully driven by the data and, unlike the Fourier or wavelet transform, this decomposition is not based on any family of functions (mother wavelet) defined as a priori. The main advantage of our approach using the EMD is that the basic functions are derived from the signal itself. Therefore, the analysis is adaptive, unlike traditional methods where the basic functions are fixed. Thus, the first step in our approach is to decompose the vibration signals into a finished sum of components called intrinsic function modes (IMFs). The construction
of IMFs requires a method for defining the envelopes to the edges of the signal that must limit the error at the edges.

- New health indicators based on the calculation of signal energy from the instantaneous frequency and amplitude.

After a presentation of the motivations that led to the choice of the EMD for the signal processing of vibration signals, we will explain its principle, stop criterion, and estimation, and compare it to spectral and multi-resolution analysis approaches to highlight its advantages and limitations. Finally, we will show how the EMD is associated with the Hilbert transform (HT) to detect the envelope of vibration signals for the determination of their instantaneous frequencies and will apply it for the denoising of vibration signals in order to extract our health indicator.

## 2. Feature Extraction for the Health Monitoring of Bearings

Studies conducted on the health monitoring of bearings indicate that it is possible to extract health indicators (features) by using two approaches: an experience-based approach given by the manufacturer and a data-driven approach [26–33]. There is a third one called the model-based approach. This latter approach is not adapted to the case of bearings. Another objective of this paper is to demonstrate that it is more interesting to monitor bearings by analyzing sensor data rather than using statistical data (experience-based approach). The proposed approach allows extracting a robust indicator able to monitor the degradation of bearings. To do this, our study is carried out on ball and roller bearings. In the proposed experiments, vibration signals will be analyzed because any mechanical component in operation generates vibrations. The degradation of a bearing generates a modification in the distribution of the vibration energy, most often leading to an increase in the level of the vibrations. By observing the evolution of this level, it is therefore possible to obtain very useful information about the state of the bearings.

The vibration analysis is based on the idea that the structures of the bearing, excited by dynamic forces, emit vibrations where the frequency is equal to those of the forces that caused them. The overall measurement taken at a point is the sum of the vibration responses of the structure compared to the various excitatory forces. We can therefore, thanks to sensors, record the vibrations transmitted by the structure of a bearing and, thanks to their analysis, identify the origin of the forces to which it is subjected. In addition, if we have the vibration “signature” of the bearing when it is considered in a good state, through comparison we can observe the evolution of its condition or detect the appearance of new dynamic forces resulting from a degradation in progress.

In general, signal processing techniques, such as spectrograms or the Wigner-Ville distribution (WVD), are based on the Fourier transform (FT), which provides a global analysis rather than a local analysis of signals. Therefore, these methods inherit certain limitations from the FT. For example, in the case of the WVD, the presence of interference terms interferes with the extraction of information in the representations made and the interpretation of the signal studied. The FT or the wavelets analysis require the use of a basic function such as the mother wavelet: hence, the advantage of choosing a decomposition adapted to the signal, not involving external functions, and allowing us to obtain a frequency representation of the signal studied. Based on these constraints, the empirical mode decomposition (EMD was proposed. The EMD adaptively breaks down a signal as a sum of oscillating components. Combined with the HT, the EMD allows an estimate of the instantaneous frequency of each component (IMF). So, knowing these modes, it is possible to know the frequency content of the signal.

The proposed approach is performed in two steps. In the first step, the EMD decomposes the signal into intrinsic mode functions (IMFs) representing the mean trend of the signal. These IMFs allow the elimination of noises in the signal. The second step applies the Hilbert transform to the IMFs to obtain the instantaneous frequencies and instantaneous amplitudes of the analyzed signal. As noted previously, the EMD decomposes a signal into a set of IMFs representing oscillating modes. Usually, the oscillations with the shortest
period (high frequency) are identified and decomposed in the first IMF. The oscillations with long periods (low frequencies) are then identified and decomposed in the following IMFs. The EMD is based on the assumption that any signal is the sum of the different IMFs. For each intrinsic mode function (IMF), the same number of extrema and zero crossings will be assigned. There is only one extremum between two successive zero crossings. Each IMF must be independent from the others. In this way, it is possible to eliminate the noise in the signal by removing the IMF generating noises. Each IMF must satisfy the following constraints:

- The number of extrema and the number of zero crossings must be greater than or equal to 1;

- The average value of the envelope defined by the local maxima and the envelope defined by the local minima must be close to zero.

The EMD decomposes a signal denoted $\delta(t)$ into four steps:

1. Step 1: The local extrema of $\delta(t)$ are identified and connected by a line to form an upper envelope (green line in Figure 3).

2. Step 2: The local minima of $\delta(t)$ are identified and connected by a line to form a lower envelope (blue line in Figure 3).

3. Step 3: The average value of the upper and lower envelope is calculated and denoted $m_1(t)$ (purple line in Figure 3). The subtraction of $m_1(t)$ from $\delta(t)$ allows creating a component denoted $h_1(t)$ as shown in Figure 4. The expression of $h_1(t)$ is given as follows:

$$ \delta(t) - m_1(t) = h_1(t) \qquad (1) $$

4. If $h_1(t)$ satisfies the constraints of the IMF, then $h_1(t)$ can be considered as the first IMF$_1$ of $\delta(t)$. If not, steps 1 to 3 will be applied to $h_1(t)$ as shown in Equation (2):

$$ h_1(t) - m_{11}(t) = h_{11}(t) \qquad (2) $$

where $m_{11}(t)$ is the average value of the upper and lower envelope of $h_1(t)$. This procedure is called the sifting process and is repeated k-times on $h_{1k}(t)$ until the average line between the upper and lower envelopes is close to zero. This procedure is represented by the following equation:

$$ h_{1(k-1)}(t) - m_{1k}(t) = h_{1k}(t) \qquad (3) $$

<figure><img src="image_7.png" /><figcaption>Figure 3. Sifting process: original signal $\delta(t)$ with the upper (green line) and lower envelope (blue line) and middle line $m_1(t)$ (purple line).</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 4. First component of the signal δ(t).</figcaption></figure>

The IMF$_{1}$, denoted $c_{1}(t)$ with $c_{1}(t) = h_{1k}(t)$, represents the high frequency component of $\delta(t)$. $c_{1}(t)$ is then extracted from $\delta(t)$ in order to repeat the sifting process in a new component denoted $r_{1}(t)$ (see Equation (4)):

$$ \left\{ \begin{array}{l} r_{1}(t) = \delta(t) - c_{1}(t) \\ r_{n}(t) = r_{n-1}(t) - c_{n}(t) \end{array} \right. \qquad (4) $$

A stop criterion must be defined to ensure that the signal obtained verifies the properties of an IMF, while limiting the number of iterations. The proposed criterion is based on the calculation of the relative variation of the signal between two consecutive iterations. It is given as follows:

$$ CR = \sum_{t=0}^{T} \left[ \frac{\left| h_{1(k-1)}(t) - h_{1k}(t) \right|^2}{h_{1(k-1)}^2(t)} \right] \qquad (5) $$

The upper and lower envelopes are obtained by an interpolation function. The interpolation function works on the minima (or maxima) values of the signal and on the corresponding abscissas. The problem here is that the first and last samples of the signal are not minima and maxima. So, there will be two intervals left with only one point to calculate the interpolation function, which is not possible. One solution would be to keep the computed interpolation polynomial over the adjacent interval. But as Figure 5 shows, the envelopes calculated with this solution (broken lines) are not satisfactory. To solve this problem, we added data at the start and end of the given interval, using the method proposed by Boustane and Quellec [34].

<figure><img src="image_9.png" /><figcaption>Figure 5. Interpolation of the extrema using a cubic spline.</figcaption></figure>
Let $\delta(t)$ be a time series of length T whose maximum and minimum envelopes must be determined and $\delta_{max} = \{\delta_{max2}, \dots, \delta_{max T}\}$ and $\delta_{min} = \{\delta_{min2}, \dots, \delta_{min T}\}$ is the set of maxima and minima.

- We will fix the first point by applying the following rule:

♦ If: $\delta_{min2} \le \delta_1 \le \delta_{max2}$

• $\delta_{max1} = \delta_{max2}$

• $\delta_{min1} = \delta_{min2}$

♦ If not

• If $\delta$ increasing: $\delta_{min1} = \delta_1$ and $\delta_{max1} = \delta_{max2}$

• If $\delta$ decreasing: $\delta_{max1} = \delta_1$ and $\delta_{min1} = \delta_{min2}$

- For the last point, we will apply:

♦ If: $\delta_{min T} \le \delta_T \le \delta_{max T}$

• $\delta_{max T} = \delta_{max T-1}$

• $\delta_{min T} = \delta_{min T-1}$

♦ If not

• If $\delta$ increasing $\delta_{max T} = \delta_T$ and $\delta_{min T} = \delta_{min T-1}$

• If $\delta$ decreasing $\delta_{min T} = \delta_T$ and $\delta_{max T} = \delta_{max T-1}$

Figure 5 shows the interpolation of the extrema using a cubic spline. This figure shows that the method of Boustane and Quellec makes it possible to minimize the errors on the bounds of the interpolation interval (the envelopes in solid lines).

The implementation of the EMD method improved by the stop criterion and the avoidance of the edge effect is given by the following Algorithm 1:

<table><caption>**Algorithm 1.** Presentation of the different stages of the improved EMD</caption><thead><tr><th></th><th></th></tr></thead><tbody><tr><td>1.</td><td>Step 1: Define a threshold μ: i &lt; 1 (i-th IMF) (μ = 0.25 in our case).</td></tr><tr><td>2.</td><td>Step 2: r<sub>i-1</sub>(t) &lt; δ(t) (residue).</td></tr><tr><td>3.</td><td>Step 3: Extract the i-th IMF.<br/>a. h<sub>i,k-1</sub>(t) &lt; r<sub>i-1</sub>(t); k &lt; 1 (k, iteration of the sifting process)<br/>b. Extraction of the maxima and minima points of h<sub>i,k-1</sub>(t)<br/>c. Calculation of the upper envelope E<sup>up</sup><sub>i,k-1</sub>(t) and lower envelope E<sup>low</sup><sub>i,k-1</sub>(t) by interpolation (Boustane and Quellec method)<br/>d. Calculation of the average envelope: m<sub>i,k-1</sub>(t) = <sup>E<sup>up</sup><sub>i,k-1</sub>(t) + E<sup>low</sup><sub>i,k-1</sub>(t)</sup>&frasl;<sub>2</sub><br/>e. Update: h<sub>i,k</sub>(t) &lt; (h<sub>i,k-1</sub>(t)-m<sub>i,k-1</sub>(t)), k &lt; k + 1<br/>f. Criterion stop calculation: CR(k) = &sum;<sub>t=0</sub><sup>T</sup> [<sup>|h<sub>i(k-1)</sub>(t)-h<sub>ik</sub>(t)|<sup>2</sup></sup>&frasl;<sub>h<sub>i(k-1)</sub><sup>2</sup>(t)</sub>] (T: length of the signal)<br/>g. Repeat steps (b) to (f) while CR(k) &ge; μ and affect IMF<sub>i</sub>(t) &lt; h<sub>i,k</sub>(t) (i-th IMF).</td></tr><tr><td>4.</td><td>Step 4: Update the residue r<sub>i</sub>(t) = r<sub>i-1</sub>(t)-IMFi(t) until the number of extrema in r<sub>i</sub>(t) is less than 2. Each extracted IMF is denoted c<sub>i</sub>(t) &lt; IMFi(t).</td></tr><tr><td>5.</td><td>Step 5: Repeat step 3 with i &lt; i + 1, until the number of extrema in r<sub>i</sub>(t) is less than 2.</td></tr></tbody></table>

The sifting process is repeated on $r_1(t)$ n-times to obtain successive components with decreasing frequency bands. This procedure allows obtaining (n-1)-IMFs from the original signal $\delta(t)$. When $r_n(t)$ becomes a monotonic function, the decomposition process can be stopped. The original signal $\delta(t)$ is composed by the (n-1)-IMFs in addition to the last component $r_n(t)$ designated as a residue, as shown in the following equation:

$$ \delta(t) = \sum_{i=1}^{n-1} c_i(t) + r_n(t) \qquad (6) $$
The IMFs $c_1(t)$, $c_2(t)$, ..., $c_{n-1}(t)$ include different frequency bands ordered from top to bottom. The frequency band contained in each component is different from another component. These frequency bands change according to the variation of the signal $\delta(t)$.

After the decomposition of $\delta(t)$ into IMFs, a Hilbert transform is applied to these latter in order to extract the instantaneous amplitude $a_i(t)$ and frequency $f_i(t)$. With the instantaneous amplitude $a_i(t)$ and frequency $f_i(t)$, it is possible to express the signal $\delta(t)$ by:

$$ \delta(t) = PR \sum_{i=1}^{n} a_i(t) \exp \left[ j2\pi \int_{0}^{T} f_i(t) dt \right] \quad (7) $$

where $PR$ is real part and $T$ is the length of the signal $\delta(t)$. This approach is analogous to the Fourier transform; it allows amplitude and frequency modulation simultaneously.

The instantaneous amplitude $a_i(t)$ is obtained by the following equation:

$$ a_i(t) = \sqrt{(c_i(t))^2 + (H[c_i(t)])^2} \quad 1 \le i \le n-1 \quad (8) $$

With $H[c_i(t)]$: the Hilbert transform of $c_i(t)$ expressed as follows:

$$ H[c_i(t)] = \frac{1}{\pi} \int_{-\infty}^{+\infty} \frac{c_i(\tau)}{t-\tau} d\tau = c_i(t) * \frac{1}{\pi t} \quad (9) $$

where the symbol "*" denotes convolution operation.

We can define the instantaneous frequency $f_i(t)$ as follows:

$$ f_i(t) = \frac{1}{2\pi} \frac{d\theta_i(t)}{dt} \quad (10) $$

$\theta_i(t)$ is the instantaneous phase defined by the following equation:

$$ \theta_i(t) = \tan^{-1} \left( \frac{H[c_i(t)]}{c_i(t)} \right) \quad (11) $$

For each IMF, the Hilbert spectrum density is defined as the square of the instantaneous amplitude:

$$ H(f_i(t)) = (a_i(t))^2 \quad (12) $$

$H(f_i(t))$ gives a local value of the energy in the time-frequency representation, and it is possible to extract the joint probability density $p(f_i(t), a_i(t))$ of the frequency and the amplitude from the set of the original series $a_i$ and $f_i$ for the IMFs set $i = 1, \dots, n-1$. This allows us to estimate the Hilbert marginal spectrum density:

$$ h_i(f) = \int_{0}^{\infty} H(f_i(t)) dt \quad (13) $$

The Hilbert marginal spectrum density is an energy density for a frequency $f$, which can be compared to the energy spectral density estimated by the Fourier transform. The objective of the proposed approach is to carry out the decomposition at the local level, in a time-frequency space.

The implementation of the feature extraction approach is given by the following Algorithm 2:
<caption>**Algorithm 2.** Presentation of the health indicator algorithm</caption>

1. Apply the EMD to the vibration signal $\delta(t)$

2. Obtain the IMF $C_i(t)$ $1 \le i \le n-1$.

3. Apply the HT to each $C_i(t)$ $1 \le i \le n-1$.

4. Obtain instantaneous amplitudes of each $C_i(t)$ $1 \le i \le n-1$.

5. Obtain instantaneous frequencies of each $C_i(t)$ $1 \le i \le n-1$.

6. Calculate the Hilbert marginal spectrum density for each $C_i(t)$ $1 \le i \le n-1$.

7. Select the Hilbert marginal spectrum corresponding to the frequency bands of bearings failures.

The flowchart given in Figure 6 summarizes the steps of an automatic extraction of the Hilbert spectrum density related to the three characteristic frequencies $f_{or}$, $f_{ir}$, and $f_b$ corresponding, respectively, to the frequencies of the outer race, inner race, and balls:

<figure><img src="image_10.png" /><figcaption>**Figure 6.** Health indicators extraction process.</figcaption></figure>

**3. Comparison of the Proposed Approach with Wavelet Analysis**

The EMD decomposition of a signal $\delta(t)$ shows that the components are classified locally from the highest frequencies to the lowest frequencies. The first modes have the fastest oscillations and the last modes have the slowest oscillations. The EMD can be considered as a self-adapting and automatic filter. The EMD can therefore be used as a signal denoising tool. It is therefore natural to compare it with the wavelet analysis. The wavelet decomposition is an iterative process that decomposes the signal $\delta(t)$ using a high-pass filter and a low-pass filter in quadrature. The high-pass filter is achieved by the convolution of the signal and the chosen wavelet. With a mean value equal to zero, the chosen wavelet captures the rapid variations of the signal $h(t)$. On the other hand, the second signal $g(t)$, coming from the low-pass filter, is downsampled, and the process is applied recursively, as shown in Figure 7. From a theoretical point of view, this downsampling consists of dividing by 2 the frequency of the wavelet and therefore capturing lower frequencies. In this sense, the EMD and discrete wavelet transform algorithms both analyze the signal starting with the extraction of details and ending with the low-frequency components. However, and unlike wavelets, the EMD is a non-parametric approach that does not require a projection basis. This is a non-linear temporal decomposition not involving an external function (mother wavelet) with a predetermined basis, while the results with the wavelet transform depend on the choice of the wavelet. Moreover, the EMD is a local, iterative, adaptive approach that describes the signal as local oscillatory components going from the high to low frequencies while the wavelet method is a parametric approach that describes the signal in two categories, the details and approximations, by ordering them from lower to higher frequencies.

One of the advantages of the EMD method is the simplicity and the rapidity of denoising a signal by reconstituting the original signal as demonstrated by formula (6). This principle is noted reciprocity. This reciprocal is a transformation often analogous to the decomposition as it is the case in the reconstruction of the wavelet where the components must be treated by the mirror filters in the reverse order of the decomposition, as illustrated in Figure 8. On the other hand, the EMD carries out extractions of signals by successive
subtraction of IMF, which makes the reconstruction fast and visual because it consists of a simple summation. In this sense, the EMD has an asymmetric character atypical of signal decompositions.

<figure><img src="image_11.png" /><figcaption>**Figure 7.** Discrete wavelet decomposition and associated filter bank analogy with EMD.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 8.** Decomposition and then reconstruction of the signal using wavelets: the two operations are symmetrical.</figcaption></figure>

Another advantage of the EMD is the stop criterion and its adaptability. Indeed, in the wavelet transform, a sub-sampling is carried out. The algorithm terminates due to finite precision of the data, but the decomposition is theoretically infinite. This stop criterion has practical interest since one can know in advance the number of components that it will be necessary to calculate. In contrast, the EMD is based on the variations of the signal for extraction of the IMFs. The number of IMFs cannot be predicted in advance. This is both a drawback and strength of the algorithm because one might think that there is a better correspondence between an IMF and the origins of its oscillations. As noted previously, the EMD therefore adapts better to the input data, and this without the choice of parameters for the decomposition (type of wavelet, number of sub-samples to be performed, etc.).

**4. Experimental Results**

First, tests were conducted on NSK 6804RS ball bearings to show the efficiency of the feature extraction approach and to show also the limits of the statistical method for monitoring ball bearings (Bbs). The characteristics of the studied Bbs are given in Table 1. Tests on Bbs corresponded to accelerated aging tests obtained by the application of a radial overload of 4000 Newtons. The used aging test bench was called "PRONOSTIA" [35] (see Figure 9). The chronogram of the aging tests is represented in Figure 8. In these tests, Bbs were subjected to constant stresses in speed and load. As shown in Figure 10, the chronogram is composed of a time axis showing the periodicities of vibration signal records (see arrows). The acquisition steps were measurement steps intended to record and verify the state of Bb degradation. Vibration signals were recorded every 10 s with a sampling frequency of 25.6 kHz.
<table><caption>**Table 1.** NSK 6804RS ball bearing characteristics.</caption><tr><td>Ball diameter (mm)</td><td>3.5</td></tr><tr><td>Number of rollers</td><td>13</td></tr><tr><td>Diameter of outer race (mm)</td><td>29.1</td></tr><tr><td>Diameter of inner race (mm)</td><td>22.1</td></tr><tr><td>Bearing mean diameter (mm)</td><td>25.6</td></tr></table>

<figure><img src="image_13.png" /><figcaption>**Figure 9.** Aging test bench "PRONOSTIA".</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 10.** Aging with constant load and speed.</figcaption></figure>

Probability calculations allowed the estimation of bearing life for a given application, but these estimations only took into account the phenomenon of fatigue (a natural phenomenon of aging in a bearing). The nominal life corresponded to the number of revolutions that at least 90% of bearings can perform under the same conditions before the appearance of the first signs of fatigue. It is denoted $L_{10}$:

$$L_{10} = \left(\frac{C}{P}\right)^{\beta} \qquad (14)$$

with:

* P: equivalent radial load applied on the bearing in Newtons (N).

* C: dynamic load supported by a bearing to achieve a nominal life of 1 million revolutions, in N.

* $\beta$ is a coefficient where: $\begin{cases} \beta = 3 \\ \beta = 10/3 \end{cases}$ for ball bearings for roller bearings.

If the speed of the Bb is constant, we can express the bearing life in hours by:

$$L_{10h} = \frac{10^6}{60n} \left(\frac{C}{P}\right)^{\beta} \qquad (15)$$

with:
* *n*: rotation speed in rpm.

Formula (15) gives a life estimation of the NSK 6804RS ball bearing for the degradation conditions (CD₁ and CD₂) described as follows:

- CD₁: rotation speed of 1800 rpm and a radial load of 4000 N → L10 h = 9 h 15 min

- CD₂: rotation speed of 1650 rpm and a radial load of 4200 N → L10 h = 8 h 43 min

Two ball bearings were tested for each degradation condition (CD). For CD₁, the bearings were respectively numbered Bb11 and Bb12. The first number corresponded to the CD, and the second number corresponded to the number of the tested bearing. For CD₂, bearings were respectively numbered Bb21 and Bb22. As shown in Figure 11, the Bbs (11 and 12), in CD₁, were degraded after 7 h 25 min and 2 h 25min, respectively.

<figure><img src="image_15.png" /><figcaption>Figure 11. Evolution of the degradation of Bbs according to the operating conditions (CDs).</figcaption></figure>

These results indicated a mismatch between the expected life given by the manufacturer (see formula (15)) and those observed for the tested bearings 11 and 12. These results also revealed a large gap between bearings 11 and 12, of 5 h 50 min. Figures 12 and 13 show the vibration signals of bearings 11 and 12, respectively. These signals corresponded to defective bearings when the acceleration exceeded 20 m/s².

<figure><img src="image_16.png" /><figcaption>Figure 12. Vibration signal of the degraded Bb11.</figcaption></figure>

We can notice from these two figures that the acceleration peaks exceeded 20 m/s² for bearing 11 and 13 m/s² for bearing 12. If we consider bearing 11 as a defective bearing, this corresponded to an rms (Root Mean Square) equal to 6.7 m/s² (see Figure 11). Th rms of bearing 12 was equal to 2.1 m/s². Even if an rms of 2.1 m/s² was a high value compared to a healthy bearing, which equals 0.45 m/s², does it mean that bearing 12 was defective. In addition, there was no visible periodicity between the acceleration peaks giving an indication of its fault frequency. The empirical mode decomposition (see Figure A1 in
Appendix A) followed by a Hilbert transform on the vibration signal of the Bb11 makes it possible to obtain the Hilbert spectrum of each IMF (see Figure 14).

<figure><img src="image_17.png" /><figcaption>Figure 13. Vibration signal of the degraded Bb12.</figcaption></figure>

Our approach was based on filtering the vibration signals around the characteristic frequencies of bearing faults. Formulas for calculating the characteristic fault frequencies of bearings are given as follows:

$$f_{or} = \frac{n_b}{2} \cdot f_r \cdot \left[ 1 - \frac{BD}{PD} \cdot \cos \psi \right] \quad (16)$$

$$f_{ir} = \frac{n_b}{2} \cdot f_r \cdot \left[ 1 + \frac{BD}{PD} \cdot \cos \psi \right] \quad (17)$$

$$f_b = \frac{PD}{BD} \cdot f_r \cdot \left[ 1 - \frac{BD^2}{PD^2} \cdot \cos^2 \psi \right] \quad (18)$$

with:

- $f_{ir}$: characteristic frequency of the inner ring;

- $f_{or}$: characteristic frequency of the outer ring;

- $f_b$: characteristic frequency of the ball;

- $f_r$: frequency of rotation of the rotor;

- $\Psi$: contact angle;

- $n_b$: number of balls;

- $BD$: ball diameter;

- $PD$: pitch diameter of the bearing.

For a rotation speed of 1800 rpm ($f_r = 30$ Hz), and for the geometry of the NSK 6804RS ball bearing, we obtain:

- $f_{or} = 168$ Hz;

- $f_{ir} = 221$ Hz;

- $f_b = 215$ Hz.

The Hilbert spectra of IMFs 4–8 allow us to identify the frequency corresponding to the rotation speed (30 Hz) and the frequency corresponding to the inner race fault ($f_{ir}$). However, IMFs 1–3 do not show significant anomalies (the amplitude of the defect is embedded in the noise). These spectral anomalies may be due to the following:

- Transmission belt in poor condition: The poor condition of a "V" belt (variation in width, deformation) creates variations in tension able to induce vibrations of frequency equal to the rotation of the belt. If the pulleys are not properly aligned, there will be an axial component inherent in this vibration.

- Chipping of bearings: Chipping of a bearing railway causes shocks and bearing resonances that are easily identified with a shock-wave measuring device. In spectral
analysis, this phenomenon appears at high frequencies with a spectral density that increases when the bearing deteriorates. If the bearing damage only included a single point, a peak could appear at a frequency linked to the speed of the bearing and to the dimensions of the bearing (diameter of inner and outer race, number of balls).

Friction: The friction of surfaces with asperities generates vibrations generally located in the high-frequency range. The state of the surface and the nature of the materials in contact have an influence on the intensity and frequency of the vibrations thus created. These phenomena are often sporadic and therefore difficult to analyze and monitor.

Given that the vibration signal $\delta(t)$ of Bb11 is the sum of the IMFs, the filtration of the high frequencies is obtained by summing the IMFs, which allows the low frequencies to appear. We thus obtain Figure 15. This figure shows a filtered vibration signal. We can observe a periodicity of 0.033 s corresponding to the rotation speed of 1800 rpm (30 Hz). The rms of this signal is equal to 3 m/s². We can also observe from Figures 16–18 that the spectrum of the 6th-IMF indicates strong amplitudes around the characteristic frequencies of the defects of Bb NSK 6804RS.

<figure><img src="image_18.png" /><figcaption>Figure 14. Hilbert spectrum of IMFs in the vibration signal of Bb11.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>Figure 15. Vibration signal of defective Bb11 after filtering.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Figure 16. Hilbert spectrum of IMF6 in the vibration signal of Bb11.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Figure 17. Hilbert spectrum of IMF6 in the vibration signal of Bb12.</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>**Figure 18.** Hilbert spectrum of IMF6 in the vibration signal of Bb21 considered as defective.</figcaption></figure>

By repeating the same steps on the vibration signal of bearing 12, we obtain Figure 19:

<figure><img src="image_23.png" /><figcaption>**Figure 19.** Vibration signal of defective Bb12 after filtering.</figcaption></figure>

In the Bb12, the Hilbert spectrum of IMF 6 was equal to 4.3 m/s² around 222 Hz and therefore around the characteristic frequencies of the inner race (see Figure 17). The same point may be made for the Bb11, where the amplitudes in the IMF6 spectrum were equal to 9.9 m/s² (see Figure 16). For Bb 21, we can see from Figure 16 that the maximum amplitudes in the IMF6 spectrum were around 204 Hz. Using the Formulas (16)–(18) for a rotation speed of 1650 rpm ($f_r$ = 27.5 Hz) and the geometry of the NSK 6804RS ball bearing, we obtain: $f_r$ = 154 Hz, $f_{ir}$ = 203 Hz, $f_b$ = 197 Hz. It is clear from Figure 18 that the increase of the amplitude around 204 Hz was due to an inner race and ball defect.

Each Bb operated under very specific constraints. We can see from Figures 16–18 that Hilbert spectrum amplitudes located around the fault frequencies of Bbs 11, 12, and 21 were clearly distinct and quantifiable in terms of severity. We can conclude that this method allows detecting, locating, and quantifying the severity of degradation for a given speed and load. In fact, if we take the case of Bb21, we can see from Figures 20 and 21 that the degradation of Bb12 was reflected by the Hilbert marginal spectrum and the rms. In Figure 21, the classical rms of the vibration signal of Bb12 is represented. We chose three measurement points. The first point corresponded to the beginning of the test. This measurement corresponded to a Bb in good health. The third measurement was taken at the end of the test and represented a faulty Bb. The second measurement point was taken when the first signs of degradation were detected. By comparing the two figures, it can be said that it is easier to follow the evolution of the degradation of a bearing by filtering and extracting its Hilbert marginal spectral density in order to detect the first signs of degradation. The first signs of degradation were detected at $t = 122$ min for our
indicator while for the rms, the detection was validated only after 140 min into the test. It was explained in [36] that classical spectral analysis cannot detect rolling defects due to the presence of background noise and the randomness of the vibrations. By applying the proposed approach, it is possible to follow the degradation of the Bb at the beginning, as shown in Figure 22. Figure 20 shows the Hilbert spectrum of the vibration signal recorded after 122 min into the test (point 2 in Figure 20). With this figure, we can see that the beginning of the degradation was due to an anomaly in the inner race of the bearing.

<figure><img src="image_24.png" /><figcaption>**Figure 20.** Evolution of the degradation of Bb21 shown by the Hilbert Marginal spectrum density.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>**Figure 21.** Evolution of the degradation of the Bb21 shown by the rms.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>**Figure 22.** Hilbert spectrum of IMF6 in the vibration signal of Bb21 after 122 min of test.</figcaption></figure>

In order to verify the validity of this method for other types of bearings, we decided to use it to detect the appearance of defects in roller bearings. A database corresponding to the accelerated aging of four roller bearings was used. These bearings were placed on a test bench designed by Rexnord Corp. The test bench was composed of four ZA-2115 roller
bearings placed on a shaft, as shown in Figure 23. The characteristics of these bearings are detailed in Table 2.

<figure><img src="image_27.png" /><figcaption>Figure 23. Roller bearing aging bench.</figcaption></figure>

<table><caption>**Table 2.** Roller bearing characteristics.</caption><tr><td>Roller diameter (mm)</td><td>8.4</td></tr><tr><td>Number of rollers</td><td>16</td></tr><tr><td>Pitch diameter (mm)</td><td>71.62</td></tr><tr><td>Contact angle (°)</td><td>15.17</td></tr></table>

The rotation speed was fixed at 2000 rpm and a radial load of 26,689 N was placed on the shaft and bearings by a spring mechanism. Vibration signals were recorded every 10 min with a sampling frequency of 20 kHz. The vibration signal of bearing number 3 recorded on the last day of testing, was used to assess the performance of the proposed approach. This signal corresponded to an inner race fault. The vibration signal of bearing 3 is shown in Figure 24:

<figure><img src="image_28.png" /><figcaption>Figure 24. Vibration signal of an inner race defect in bearing 3.</figcaption></figure>

From this signal, the EMD and the HT were performed in order to obtain the marginal spectral density of each IMF. Fifteen IMFs we obtained. Knowing the rotational speed of the bearings and their characteristics, we were able to deduce the frequency of faults. They are given as follows:
- Outer ring fault: 236.4 Hz;

- Inner ring fault: 296.8 Hz;

- Roller fault: 280.8 Hz.

The Hilbert spectrum amplitudes of IMFs 1 to 6 indicated that the low-frequency spectrum showed no significant anomalies. On the other hand, IMF 7 indicated high peaks (0.74 m/s²) around 296 Hz and therefore around the characteristic frequencies of the inner ring (see Figures 25 and 26).

<figure><img src="image_29.png" /><figcaption>**Figure 25.** Spectrum of the IMF 7 for bearing 3 at $t = 0$ (red line) and $t = 35$ days (blue line).</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>**Figure 26.** Defect in the inner ring of bearing 3.</figcaption></figure>

**5. Conclusions**

This paper presented an approach for the condition monitoring of bearings based on the EMD and the HT. As presented previously, the EMD decomposed the studied signal into IMFs. These components were classified locally from the highest frequencies to the lowest frequencies. This reflected a behavior in a bank of self-adaptive and automatic filters of the EMD and its interest in signal denoising. However, some errors related to the principle of EMD affect the quality of signal decomposition. Indeed, EMD decomposition has certain limitations linked to its basic principle, which makes use of the sifting process, interpolation, to determine the upper and lower envelopes, or to the nature of the studied signal, its sampling. The most important errors are due to the boundary effects and accumulations of errors introduced by the sifting algorithm and interpolation (Gibbs phenomenon). A solution was proposed to solve the problem of the boundary effects. The obtained results demonstrated the effectiveness of the proposed approach for the detection and diagnosis of bearing defects (ball or roller) for different speeds and loads. Implementing the EMD requires making choices, especially with regard to the criteria for stopping the sifting process, which is deeply linked to the free interpretation of the zero local mean value in the definition of an IMF. Another negative point is the fact that the EMD is very sensitive to boundary errors. The final drawback of the presented EMD was its lack of analytical expression. Because the EMD is an algorithm, there is no underlying mathematical foundation for this decomposition, such as Fourier or wavelets transforms.
## 6. Perspectives

In addition to detecting the presence of a defect and locating it, it would also be interesting to identify the different degradation stages before they appear. We therefore recommend retesting the NSK 6804RS ball bearings under less severe operating conditions than those presented previously. Based on formula (15), the life of the NSK 6804RS ball bearing for the new degradation conditions (CD₁, CD₂ and CD₃) is given as follows:

- CD₁: rotation speed of 1800 rpm and a radial load of 2000 N → L10 h = 74 h 04 min

- CD₂: rotation speed of 1650 rpm and a radial load of 2500 N → L10 h = 41 h 22 min

- CD₃: rotation speed of 1500 rpm and a radial load of 3000 N → L10 h = 26 h 20 min

In addition to the steady state, two other types of aging tests need to be considered:

* Non-stationary tests (NS1): ball bearings are subjected to constant stresses in speed and variable stresses under load;

* Non-stationary tests (NS2): ball bearings are subjected to constant stresses under load and variable stresses in speed.

### Tests in non-stationary state (NS1):

The proposed timing diagram used for aging in the non-stationary state (NS1) is presented in Figure 27. The Bbs will be subjected to a constant speed of 2000 rpm and a radial load varying between 2.5 and 3.5 kN. When a Bb is subjected to a variable load and a constant speed, an equivalent load can be used to determine the life of the Bb. It is given by the following formula:

$$C_{eq} = (t_1 \cdot C_1^{\alpha} + t_2 \cdot C_2^{\alpha} + \dots + t_z \cdot C_z^{\alpha})^{\alpha^{-1}} \quad (19)$$

<figure><img src="image_31.png" /><figcaption>Figure 27. Chronogram used for monitoring the aging of Bbs in non-stationary test (NS1).</figcaption></figure>

- $C_{eq}$: Equivalent load applied to the bearing during the test;

1. $t_i$: Duration of load use with: $\sum_{i=1}^{z} t_i = 1$;

2. $C_i$: Applied load in the Bb for a period $t_i$;

3. $\alpha$: Constant equal to 3 for a ball bearing and 10/3 for a roller bearing.

If a load of 3.5 kN is applied to the Bb for 40% of the test duration and then 2.5 kN in the remaining duration (60%), then the equivalent load will be 3 kN. The lifetime of the Bb will be equal to 19 h 45 min. In the end, we would apply a radial load of 3.5 kN for 7 h 54 min, then a radial load of 2.5 kN for 11 h 51 min.

### Tests in non-stationary state (NS2):

The chronogram used for aging in the non-stationary test (NS2) is presented in Figure 28. The Bbs will be subjected to a constant radial load of 3 kN and a speed which varies between 1650 and 2000 rpm. When a Bb is subjected to a variable speed and a constant load, an equivalent speed can be used to determine the life of the Bb. It is given by the following formula:

$$V_{eq} = t_1 \cdot V_1 + t_2 \cdot V_2 + \dots + t_z \cdot V_z \quad (20)$$
with:

- $V_{eq}$: equivalent speed applied to the Bb during the test;

- $V_t$: rotation speed applied to the Bb for duration $t_i$.

<figure><img src="image_32.png" /><figcaption>**Figure 28.** The chronogram used for monitoring the aging of Bbs in non-stationary test (NS2).</figcaption></figure>

If a rotation speed of 1650 rpm is applied to the Bb for 40% of the test duration and then 2000 rpm for the remaining duration (60%), then the equivalent speed will be equal to 1860 rpm. The lifetime of the Bb will be equal to 21 h 23 min. In the end, we would apply a speed of 1650 rpm for 8 h 33 min, then a speed of 2000 rpm for 12 h 49 min.

**Author Contributions:** All authors contributed to the research in this paper. A.S. conceived the methodology and model of the research system of this article. B.E.Y. performed experimentation and analyzed the data. H.R. finished the original draft of the article. T.W. edited and reviewed the article. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

## Appendix A

<figure><img src="image_33.png" /><figcaption>**Figure A1.** *Cont.*</figcaption></figure>
<figure><img src="image_34.png" /><figcaption>Figure A1. IMFs of the vibration signal of Bb11.</figcaption></figure>

References

1. Albrecht, P.F.; Appiarius, J.C.; Sharma, K.D. Assessment of the reliability of motors in utility applications-Updated. *IEEE Trans. Energy Convers.* **1986**, *1*, 39-46. [CrossRef]

2. IEEE Committee Report. Report of large motor reliability survey of industrial and commercial installation, Part I and Part II. *IEEE Trans. Ind. Appl.* **1985**, *21*, 853-872.

3. Stewart, R. The Specification and Development of a Standard for Gearbox Monitoring. In *Proceedings of the 2nd International Conference on Vibration in Rotating Machines, Cambridge, UK, 1-4 September 1980; Mechanical Engineering Publications: London, UK, 1980; pp. 353-358*.

4. Randall, R.; Antoni, J.; Chobsaard, S. The Relationship between spectral correlation and envelope analysis in the diagnostics of bearing faults and other cyclostationary machine signals. *Mech. Syst. Signal Process.* **2001**, *15*, 945-962. [CrossRef]

5. Choi, Y.-C.; Kim, Y.-H. Fault detection in a ball bearing system using minimum variance cepstrum. *Meas. Sci. Technol.* **2007**, *18*, 1433-1440. [CrossRef]

6. Antoni, J. The spectral kurtosis: A useful tool for characterising non-stationary signals. *Mech. Syst. Signal Process.* **2006**, *20*, 282-307. [CrossRef]

7. Antoni, J.; Randall, R. The spectral kurtosis: Application to the vibratory surveillance and diagnostics of rotating machines. *Mech. Syst. Signal Process.* **2006**, *20*, 308-331. [CrossRef]

8. Wang, Y.; He, Z.; Zi, Y. A demodulation method based on improved local mean decomposition and its application in rub-impact fault diagnosis. *Meas. Sci. Technol.* **2009**, *20*, 025704. [CrossRef]

9. Meng, Q.; Qu, L. Rotating machinery fault diagnosis using Wigner distribution. *Mech. Syst. Signal Process.* **1991**, *5*, 155-166. [CrossRef]

10. Antoni, J. Fast computation of the kurtogram for the detection of transient faults. *Mech. Syst. Signal Process.* **2007**, *21*, 108-124. [CrossRef]

11. Rai, A.; Upadhyay, S. A review on signal processing techniques utilized in the fault diagnosis of rolling element bearings. *Tribol. Int.* **2016**, *96*, 289-306. [CrossRef]

12. Chen, B.; Shen, B.; Chen, F.; Tian, H.; Xiao, W.; Zhang, F.; Zhao, C. Fault diagnosis method based on integration of RSSD and wavelet transform to rolling bearing. *Measurement* **2019**, *131*, 400-411. [CrossRef]

13. Yu, X.; Ding, E.; Chen, C.; Liu, X.; Li, L. A Novel Characteristic Frequency Bands Extraction Method for Automatic Bearing Fault Diagnosis Based on Hilbert Huang Transform. *Sensors* **2015**, *15*, 27869-27893. [CrossRef] [PubMed]

14. Guo, T.; Deng, Z. An improved EMD method based on the multi-objective optimization and its application to fault feature extraction of rolling bearing. *Appl. Acoust.* **2017**, *127*, 46-62. [CrossRef]

15. Lee, D.-H.; Ahn, J.-H.; Koh, B.-H. Fault Detection of Bearing Systems through EEMD and Optimization Algorithm. *Sensors* **2017**, *17*, 2477. [CrossRef] [PubMed]
16. Yang, T.; Guo, Y.; Wu, X.; Na, J.; Fung, R.-F. Fault feature extraction based on combination of envelope order tracking and cICA for rolling element bearings. *Mech. Syst. Signal Process.* **2018**, *113*, 131–144. [CrossRef]

17. Abboud, D.; Antoni, J.; Eltabach, M.; Sieg-Zieba, S. Angle\time cyclostationarity for the analysis of rolling element bearing vibrations. *Measurement* **2015**, *75*, 29–39. [CrossRef]

18. Yan, X.; Liu, Y.; Jia, M. A Feature Selection Framework-Based Multiscale Morphological Analysis Algorithm for Fault Diagnosis of Rolling Element Bearing. *IEEE Access* **2019**, *7*, 123436–123452. [CrossRef]

19. Jiang, F.; Zhu, Z.; Li, W.; Ren, Y.; Zhou, G.; Chang, Y. A Fusion Feature Extraction Method Using EEMD and Correlation Coefficient Analysis for Bearing Fault Diagnosis. *Appl. Sci.* **2018**, *8*, 1621. [CrossRef]

20. Yang, Z.; Zhong, J. A Hybrid EEMD-Based SampEn and SVD for Acoustic Signal Processing and Fault Diagnosis. *Entropy* **2016**, *18*, 112. [CrossRef]

21. Liang, J.; Zhong, J.-H.; Yang, Z.-X. Correlated EEMD and Effective Feature Extraction for both Periodic and Irregular Faults Diagnosis in Rotating Machinery. *Energies* **2017**, *10*, 1652. [CrossRef]

22. Yi, C.; Lin, J.; Zhang, W.; Ding, J. Faults Diagnostics of Railway Axle Bearings Based on IMF's Confidence Index Algorithm for Ensemble EMD. *Sensors* **2015**, *15*, 10991–11011. [CrossRef]

23. Zhao, H.; Sun, M.; Deng, W.; Yang, X. A New Feature Extraction Method Based on EEMD and Multi-Scale Fuzzy Entropy for Motor Bearing. *Entropy* **2016**, *19*, 14. [CrossRef]

24. Li, R.; Ran, C.; Zhang, B.; Han, L.; Feng, S. Rolling Bearings Fault Diagnosis Based on Improved Complete Ensemble Empirical Mode Decomposition with Adaptive Noise, Nonlinear Entropy, and Ensemble SVM. *Appl. Sci.* **2020**, *10*, 5542. [CrossRef]

25. Zhan, L.; Ma, F.; Zhang, J.; Li, C.; Li, Z.; Wang, T. Fault Feature Extraction and Diagnosis of Rolling Bearings Based on Enhanced Complementary Empirical Mode Decomposition with Adaptive Noise and Statistical Time-Domain Features. *Sensors* **2019**, *19*, 4047. [CrossRef]

26. Khlaief, A.; Nguyen, K.; Medjaher, K.; Picot, A.; Maussion, P.; Tobon, D.; Chauchat, B.; Cheron, R. Feature Engineering for Ball Bearing Combined-Fault Detection and Diagnostic. In Proceedings of the 2019 IEEE 12th International Symposium on Diagnostics for Electrical Machines, Power Electronics and Drives (SDEMPED), Toulouse, France, 27–30 August 2019; pp. 384–390.

27. He, L.; Hao, L.; Qiao, W. Remote Monitoring and Diagnostics of Pitch Bearing Defects in a MW-Scale Wind Turbine Using Pitch Symmetrical-component Analysis. In Proceedings of the 2019 IEEE Energy Conversion Congress and Exposition (ECCE), Baltimore, MD, USA, 29 September–3 October 2019; pp. 1–6.

28. Xiaozhi, L.; Ganggang, S.; Yinghua, Y. Fault Diagnosis of Rolling Bearing Based on Wavelet Packet Transform and GA-Elman Neural Network. In Proceedings of the 2019 Chinese Control and Decision Conference (CCDC), Nanchang, China, 3–5 June 2019; pp. 462–466.

29. Wang, H.; Xu, J.; Yan, R.; Gao, R.X. A New Intelligent Bearing Fault Diagnosis Method Using SDP Representation and SE-CNN. *IEEE Trans. Instrum. Meas.* **2020**, *69*, 2377–2389. [CrossRef]

30. Tian, Y.; Liu, X. A deep adaptive learning method for rolling bearing fault diagnosis using immunity. *Tsinghua Sci. Technol.* **2019**, *24*, 750–762. [CrossRef]

31. Qiao, M.; Yan, S.; Tang, X.; Xu, C. Deep Convolutional and LSTM Recurrent Neural Networks for Rolling Bearing Fault Diagnosis under Strong Noises and Variable Loads. *IEEE Access* **2020**, *8*, 66257–66269. [CrossRef]

32. Liu, Z.; Tang, X.; Wang, X.; Mugica, J.E.; Zhang, L. Wind Turbine Blade Bearing Fault Diagnosis Under Fluctuating Speed Operations via Bayesian Augmented Lagrangian Analysis. *IEEE Trans. Ind. Inform.* **2020**, *1*. [CrossRef]

33. Mao, W.; Liu, Y.; Ding, L.; Li, Y. Imbalanced Fault Diagnosis of Rolling Bearing Based on Generative Adversarial Network: A Comparative Study. *IEEE Access* **2019**, *7*, 9515–9530. [CrossRef]

34. Boustane, T.; Quellec, G. *Implantation de la Méthode EMD en C avec Interface Matlab*; Institut Supérieur d'Informatique de Modélisation et de leurs Applications: Aubiere, France, 2004.

35. Medjaher, K.; Zerhouni, N.; Baklouti, J. Data-driven prognostics based on health indicator construction: Application to PRONOSTIA's data. In Proceedings of the 2013 European Control Conference (ECC), Zurich, Switzerland, 17–19 July 2013; pp. 1451–1456.

36. Antoni, J. Cyclic spectral analysis of rolling-element bearing signals: Facts and fictions. *J. Sound Vib.* **2007**, *304*, 497–529. [CrossRef]