A Fast Signal Estimation Method Based on Probability Density Functions for Fault Feature Extraction of Rolling Bearings

Abstract

Fault diagnosis of rolling bearings is essential to ensure the efficient and safe operation of mechanical equipment. The extraction of fault features of the repetitive transient component from noisy vibration signals is key to bearing fault diagnosis. However, the bearing fault-induced transients are often submerged by strong background noise and interference. To effectively detect such fault-related transient components, this paper proposes a probability- and statistics-based method. The maximum-a-posteriori (MAP) estimator combined with probability density functions (pdfs) of the repetitive transient component, which is modeled by a mixture of two Laplace pdfs and noise, were used to derive the fast estimation model of the transient component. Subsequently, the LapGauss pdf was adopted to model the noisy coefficients. The parameters of the model derived could then be estimated quickly using the iterative expectation–maximization (EM) algorithm. The main contributions of the proposed statistic-based method are that: (1) transients and their wavelet coefficients are modeled as mixed Laplace pdfs; (2) LapGauss pdf is used to model noisy signals and their wavelet coefficients, facilitating the computation of the proposed method; and (3) computational complexity changes linearly with the size of the dataset and thus contributing to the fast estimation, indicated by analysis of the computational performance of the proposed method. The simulation and experimental vibration signals of faulty bearings were applied to test the effectiveness of the proposed method for fast fault feature extraction. Comparisons of computational complexity between the proposed method and other transient extraction methods were also conducted, showing that the computational complexity of the proposed method is proportional to the size of the dataset, leading to a high computational efficiency.

1. Introduction

Rolling bearings are key parts of much of the modern mechanical equipment which plays an essential role in social production and living, such as engines, generators, machine tools, and so on. As rolling bearings usually operate under harsh working conditions and need to work long hours, they are prone to failures [1]. So, it is crucial to diagnose faults of the rolling bearings as soon as possible to avoid safety accidents and maintain good working efficiency.

Vibration signals generated by faulty bearings directly reflect the dynamic characteristics of faulty bearings [2]. Therefore, it is effective and reasonable to extract the features of the repetitive transient components by analyzing the vibration signals of faulty bearings. However, the vibration signals measured are always drowned out by strong noise, which makes the features of the repetitive transient component difficult to extract [3]. Therefore, methods to effectively extract the features of the repetitive transient component from the noisy vibration signals are key to bearing fault diagnosis.

To deal with this issue, various effective methods, such as time–frequency analysis [4], spectral kurtosis [5], empirical mode decomposition [6], wavelet analysis [7], stochastic resonance [8], deep learning [9], sparse representation [10] and bispectrum analysis [11] have been proposed over the past decade. The time-frequency analysis method with high time-frequency aggregation has achieved good results in the extraction of frequency conversion curve and the extraction of fault feature frequency [12]. The empirical mode decomposition method has good time-frequency resolution and adaptability and is able to separate the noise and effective component by reasonably selecting the intrinsic mode functions without giving a given transform base [13]. The spectral kurtosis filtering method, based on the non-Gaussian property exhibited by the fault component, designs the filter to select the resonant frequency band of the fault component according to the principle of maximizing the statistic such as kurtosis, and then determines the fault feature frequency according to envelope spectrum analysis [14]. The deep learning method can adaptively extract features based on the model constructed, which provides a new idea for improving the more accurate and effective fault diagnosis of rolling bearings [15]. Sparse representation has been successfully applied to mechanical fault diagnosis and achieved significant results. As an optimization method, the sparse representation method needs to optimize an objective function, which consists of the fidelity term, penalty term and regularization parameter [16]. In the field of mechanical fault diagnosis, research on sparse representation has focused on the construction of penalty functions [17]. At the same time, in the sparse representation method an over-complete dictionary according to the characteristics of the components needs to be constructed to be estimated in advance [18]. Modulation signal bispectrum (MSB) is a representative of bispectrum analysis-based methods. The MSB method has inherent noise suppression capabilities and can help enhance the modulation effects of a bearing fault and achieve the goal of frequency band optimization for fault detection [19]. The aforementioned methods have achieved good results in fault feature extraction. However, most methods described above cannot satisfy the requirement of timeliness very well. From the point of view of the timeliness of fault feature extraction, a method to quickly estimate the transient components of noisy bearing fault signals and extract fault features is needed.

Many methods based on probability and statistics have been proposed and successfully applied to many fields, including fault diagnosis and detection [20]. These methods can mainly be divided into three categories: maximum-likelihood (ML) [21], Bayesian [22], and maximum-a-posteriori (MAP) approaches [23]. The methods based on ML are common and use a simple approach, which needs to maximize a log-likelihood function. However, these approaches might not work for denoising problems [24]. The Bayesian-oriented methods make the optimal estimation based on the known probability distribution and observed data, which is often used to solve the optimization and approximation problems of nonlinear signal models or processing systems [25]. MAP approaches use an assumed prior distribution of the related signal to regularize the estimation process, and the basis of MAP is Bayesian formula. These methods usually have deduced probabilistic models and the pdfs of the related signals play an important role in these methods. Furthermore, most of these methods need to iteratively estimate the parameters of the models deduced. But, generally speaking, these methods have high computational efficiency.

In this paper, we propose an effective wavelet-based method with low computational complexity to fast extract the fault feature components from noisy bearing fault signals. This method only needs one prior, the pdfs of related signals. The analysis of simulation and experimental noisy fault vibration signals prove that our proposed method can efficiently estimate the repetitive transient component and extract fault feature frequency. The proposed method uses the MAP estimator and derives the general model of estimating the transient components based on the pdfs of the transient components of faulty bearings and noise, which means that our method has a scientific explanation of probability theory. From the point of view of timeliness, the proposed method in this paper only requires a small amount of iterative parameter estimation of the model derived using the expectation-maximization (EM) algorithm. At the same time, the computational cost of the discrete wavelet transform is also very low. Therefore, compared with some other popular fault feature extraction methods, the computational complexity of our proposed method is much lower and the running time also proves the high computation efficiency. Our method can not only get the result at a fast speed, but also ensure the accuracy of fault feature extraction.

This paper is organized as follows. In Section 2, the derivation of the model to estimate the transient components and the algorithm to estimate parameters are presented in detail. In Section 3, the proposed method is used to analyze a simulated signal. In Section 4, three experimental bearing fault signals are analyzed using the proposed method respectively. Lastly, Section 5 presents the conclusions of the study.

2. Theoretical Description

2.1. Model Derivation

In the time domain, a bearing fault signal corrupted by white Gaussian noise can be represented as

$$g=x+n_t,$$

(1)

where $g\in R^N$ denotes the observed noisy bearing fault signal, $x\in R^N$ denotes the noise-free transient components and $n_t\in R^N$ denotes noise in time domain. It is necessary to estimate the transient components $x$ as accurately as possible based on the noisy signal $g$ and their characteristics. In order to make the most of the statistical distribution characteristics of related signals, a linear transformation is needed. Wavelet transform is an efficient linear transformation and it transforms a signal in time domain into time-frequency domain, which can better observe the local characteristics of the signal and simultaneously observe the time and frequency information of the signal [26]. At the same time, the computational complexity of wavelet transform is also low, which can ensure the fast extraction of fault features. When the wavelet transform with orthogonality is performed, the model can be represented as

$$y=w+n,$$

(2)

where $y$ denotes the wavelet coefficients of the observed noisy bearing fault signal $g$, $w$ denotes the wavelet coefficients of the noise-free transient components $x$ and $n$ denotes noise.

Through the analysis of bearing fault mechanism and statistical analysis of simulation signals, we found that the pdfs of the noise-free transient components and its wavelet coefficients are mixture of two Laplace pdfs [27]:

$$p\left(w\right)=k_1\mathrm{lap}\left({\sigma }_1\right)+k_2\mathrm{lap}\left({\sigma }_2\right),$$

(3)

where

$$\mathrm{lap}(\sigma )=\frac{1}{\sqrt{2}\sigma }{e}^{-\frac{\sqrt{2}}{\sigma }|w|},$$

(4)

and $k_1+k_2=1$. And noise follows Gaussian distribution and is independent of transient components [28].

In order to minimize the noise contained in the estimated components $w$, the MAP estimator is adopted to estimate $w$ according to noisy wavelet coefficients $y$ [23]. The MAP estimator for $w$ in Equation (2) is defined as

$$\widehat{w}\left(y\right)=\arg \underset{w}{\max }p(w|y).$$

(5)

Based on the Bayes’ theorem, the maximum a posteriori estimator can be expressed as

$$\widehat{w}\left(y\right)=\arg \underset{w}{\max }\frac{p(y|w)p\left(w\right)}{p\left(y\right)}.$$

(6)

As the pdf $p(y)$ on the denominator is irrelevant with $w$, the MAP estimator of $w$ is transformed into the following form

$$\widehat{w}\left(y\right)=\arg \underset{w}{\max }p(y|w)\cdot p\left(w\right).$$

(7)

The Equation (7) is equivalent to solving the following formula

$$\widehat{w}\left(y\right)=\arg \underset{w}{\max }\ln \left[p(y|w)\cdot p\left(w\right)\right].$$

(8)

Consequently, the MAP estimate of $w$ can be obtained by solving the following formula

$$\frac{\partial }{\partial w}\ln p(y|w)+\frac{\partial }{\partial w}\ln p\left(w\right)=0.$$

(9)

According to the form of the above formula, it is necessary to know the conditional pdf $p(y|w)$. As for the model in Equation (2), $p(w|y)$ denotes the pdf of $n$ when the parameter $w$ is given

$$p\left(y|w\right)=\frac{1}{\sqrt{2\pi }{\sigma }_n}{e}^{-\frac{{\left(y-w\right)}^2}{2{\sigma }_n^2}},$$

(10)

where ${\sigma }_n^2$ represents the variance of the white zero-mean Gaussian noise. The standard variance of the independent noise ${\sigma }_n$ also needs to be estimated. The parameter can be estimated through a frequently used method when implementing denoising in wavelet domain

$${\sigma }_n=\frac{\mathrm{MAD}}{0.6745},$$

(11)

where MAD denotes median absolute deviation

$$\mathrm{MAD}=\mathrm{Median}\left\{|y_n-\mathrm{Median}\left(y_n\right)|\right\},$$

(12)

and $y_n$ denotes the wavelet coefficients in the finest level [29]. In this paper, the estimation is carried out at the first level which which is often the finest level used for estimation in wavelet denoising problems and thus $y_n$ denotes the level 1 detail coefficients.

The single Laplace model will be discussed firstly and then extended to the mixture model. If the distribution of the wavelet coefficients $w$ follow single Laplacian pdf,

$$p\left(w\right)=\frac{1}{\sqrt{2}\sigma }{e}^{-\frac{\sqrt{2}}{\sigma }|w|}.$$

(13)

Plug the above equation into Equation (9), it can be found that

$$y=\widehat{w}+\frac{\sqrt{2}{\sigma }_n^2}{\sigma }\cdot \mathrm{sign}\left(\widehat{w}\right).$$

(14)

The relationship between $\widehat{w}$ and $y$ can be represented as

$$\widehat{w}=\{\begin{array}{cc}y+T,& y\le -T\\ 0,& |y|\le T\\ y-T,& y\ge T\end{array},$$

(15)

where $T=\sqrt{2}{\sigma }_n^2/\sigma$. The above function is the frequently used soft thresholding and the notation is often written as

$$\widehat{w}\left(y\right)=\mathrm{soft}\left(y,\frac{\sqrt{2}{\sigma }_n^2}{\sigma }\right).$$

(16)

Then, if the pdf of $w$ is a mixture of two Laplace pdfs as in Equation (3), the result is shown as

$$\widehat{w}\left(y\right)=p_{k_1}\left(y\right)\mathrm{soft}\left(y,\frac{\sqrt{2}{\sigma }_n^2}{{\sigma }_1}\right)+p_{k_2}\left(y\right)\mathrm{soft}\left(y,\frac{\sqrt{2}{\sigma }_n^2}{{\sigma }_2}\right),$$

(17)

where

$$p_{k_i}\left(y\right)=\frac{k_i{g}_i\left(y\right)}{k_1{g}_1\left(y\right)+k_2{g}_2\left(y\right)},i=1,2.$$

(18)

As to the model $y=w+n$, since $w$ and $n$ are independent of each other, the pdf of $y$ can be calculated as the convolution of $w$ and $n$ according to relevant knowledge of probability theory

$$g_i\left(y\right)=\frac{1}{\sqrt{2}{\sigma }_i}{e}^{-\frac{\sqrt{2}}{{\sigma }_i}|y|}\ast \frac{1}{\sqrt{2\pi }{\sigma }_n}{e}^{-\frac{y^2}{2{\sigma }_n^2}},i=1,2.$$

(19)

Through the aforementioned deductions, the pdf of $y$ is needed to improve the efficiency of calculations. That is, it is necessary to find the pdf of $y=w+n$, where $w$ is the noise-free value with a Laplace distribution and $n$ is the noise with a Gaussian distribution.

As $w$ is a Laplace random variable and assume that the standard deviation of $w$ is ${\sigma }_w$. The pdf of $w$ is:

$$p\left(w\right)=\frac{1}{\sqrt{2}{\sigma }_w}{e}^{-\frac{\sqrt{2}}{{\sigma }_w}|w|},$$

(20)

and $n$ denotes zero-mean Gaussian noise, assume the standard deviation of $n$ is ${\sigma }_n$. The pdf of $n$ is:

$$p\left(n\right)=\frac{1}{\sqrt{2\pi }{\sigma }_n}{e}^{-\frac{n^2}{2{\sigma }_n^2}}.$$

(21)

According to Hansen [30] and after verification, the pdf of $y$ is:

$$p\left(y\right)=\frac{1}{2\sqrt{2}{\sigma }_w}{e}^{-\frac{y^2}{2{\sigma }_n^2}}\cdot [\mathrm{erfcx}\left(\frac{{\sigma }_n}{{\sigma }_w}-\frac{|y|}{\sqrt{2}{\sigma }_n}\right)+\mathrm{erfcx}\left(\frac{{\sigma }_n}{{\sigma }_w}+\frac{|y|}{\sqrt{2}{\sigma }_n}\right)],$$

(22)

where

$$\mathrm{erfcx}\left(x\right)=e^{x^2}\cdot \mathrm{erfc}\left(x\right),$$

(23)

$$\mathrm{erfc}\left(x\right)=1-\mathrm{erf}\left(x\right),$$

(24)

$$\mathrm{erf}\left(x\right)=\frac{2}{\sqrt{\pi }}{\displaystyle {\int }_0^x{e}^{-t^2}dt}.$$

(25)

The notation $\mathrm{LapGauss}(y,{\sigma }_w,{\sigma }_n)$ is used to represent the pdf in Equation (22).

2.2. Parameter Estimation

After the model is established, parameters $k_1$, $k_2$, ${\sigma }_1$ and ${\sigma }_2$ need to be estimated according to the wavelet coefficients of the noisy bearing fault signals. At present, there are many mature iterative optimization algorithms to estimate the parameters of the model, such as gradient descent method, Gibbs sampling method and EM algorithm. The expectation-maximization (EM) algorithm has been successfully applied to the optimization of the sparse representation model and parameter estimation of mixture models in many fields [31]. Usually, the computational complexity of the EM algorithm for the one-dimensional signal is not very high. So, the EM algorithm is adopted in this paper to estimate the four parameters.

EM algorithm needs to set initial values for each parameter. In the proposed method, $k_1$, $k_2$, ${\sigma }_1$ and ${\sigma }_2$ is needed to be initialized. And then a series of expectation and maximization steps will alternate. The E-step updates the responsibility factors

$$q_i\left(n\right)\leftarrow \frac{k_i\mathrm{LapGauss}\left(y_n,{\sigma }_1,{\sigma }_n\right)}{k_1\mathrm{LapGauss}\left(y_n,{\sigma }_1,{\sigma }_n\right)+k_2\mathrm{LapGauss}\left(y_n,{\sigma }_2,{\sigma }_n\right)},i=1,2.$$

(26)

Next the updated parameters $k_1$, $k_2$, ${\sigma }_1$ and ${\sigma }_2$ are estimated in the M-step. $k_1$ and $k_2$ are estimated by

$$k_i\leftarrow \frac{1}{N}{\displaystyle \sum }_{n=1}^N{q}_i\left(n\right),i=1,2.$$

(27)

Then ${\sigma }_{y\left(i\right)}^2,i=1,2$ are estimated as weighted sum of the data values [32]. As to ${\sigma }_{y\left(i\right)}^2$, the weight for $y_n$ is the probability that the data point $y_n$ is generated by $\mathrm{LapGauss}(y,{\sigma }_w,{\sigma }_n)$

$${\sigma }_{y\left(i\right)}^2\leftarrow \frac{{\displaystyle {\displaystyle \sum }_{n=1}^N}{q}_i\left(n\right)y_n^2}{{\displaystyle {\displaystyle \sum }_{n=1}^N}{q}_i\left(n\right)},i=1,2.$$

(28)

Finally, as the wavelet coefficients $w$ and $n$ are independent, the standard deviation ${\sigma }_1$ and ${\sigma }_2$ can be calculated by

$${\sigma }_i^2={\sigma }_{y\left(i\right)}^2-{\sigma }_n^2,i=1,2.$$

(29)

And it is possible that the estimated variance ${\sigma }_1^2$ and ${\sigma }_2^2$ will be negative because these parameters are all estimated. So, in order to let the value under the radical sign not to be negative and the denominator not to be zero, ${\sigma }_1$ and ${\sigma }_2$ are estimated as [32]:

$${\sigma }_i=\sqrt{\max \left({\sigma }_{y\left(i\right)}^2-{\sigma }_n^2,\frac{{\sigma }_n^2}{250}\right)},i=1,2.$$

(30)

2.3. Method Summary

The flowchart of the proposed method is shown in Figure 1. Given an observed signal $g$, the corresponding wavelet coefficients $w$ are obtained by wavelet transform firstly. Next, the standard deviation of the independent Gaussian white noise ${\sigma }_n$ is estimated by Equation (11). Then the parameters $k_1$, $k_2$, ${\sigma }_1$, ${\sigma }_2$ are estimated by Equation (27) to Equation (30), and the new estimated wavelet coefficients $\widehat{w}$ are calculated by the Equation (17). Finally, the new wavelet coefficients $\widehat{w}$ are reconstructed to obtain the estimated signal $\widehat{x}$, and the fault feature can be obtained through envelope spectrum analysis.

3. Simulation Analysis

In this section, a simulated bearing fault signal is analyzed using the proposed method, to prove its effectiveness. Through the analysis of the failure mechanism, the bearing fault signal can be constructed as

$$y(t)={\displaystyle \sum _K\left[A{e}^{\frac{-\zeta }{\sqrt{1-{\zeta }^2}}2\pi f(t-\tau -KT)}\times 2\pi f(t-\tau -KT)\right]},$$

(31)

where $A$ denotes amplitude, $\xi$ denotes damping ratio, $f$ denotes frequency, $\tau$ denotes time index, $K$ denotes the number of impulses and $T$ denotes cyclic period [27]. The noisy signal is created by adding Gaussian noise to the simulated signal. We used the discrete Daubechies 5 wavelet, which is an orthogonal wavelet transform, with three decomposition levels, and detail coefficients at each level were used for the estimation of parameters $k_1$, $k_2$, ${\sigma }_1$ and ${\sigma }_2$ utilizing the EM algorithm. The corresponding results are shown in Figure 2.

The waveform of the simulated bearing fault signal in the time domain is shown in Figure 2a, where amplitude $A$ = 1.5, frequency $f$ = 1500 Hz, damping ratio $\xi$ = 0.080, cyclic period $T$ = 0.007 s and the corresponding theoretical fault feature frequency was $f_s$ = 135.1 Hz. The length of this signal was 51,200 and the sampling frequency was set at 51.2 kHz. The waveform of the noisy signal in the time domain is shown in Figure 2b. The standard deviation of the additive Gaussian white noise was ${\sigma }_n$ = 2.0.

The signal-to-noise ratio (SNR) is introduced to describe the influence of noisy components, which is defined as

$$\mathrm{SNR}=10{\log }_{10}\frac{{\Vert x\Vert }_2^2}{{\Vert n\Vert }_2^2}.$$

(32)

The SNR of the simulated signal was −16.082 dB. The envelope spectrum analysis of the noisy bearing fault signal is shown in Figure 2c. It is obvious that the fault feature frequency cannot be obtained directly. The time domain waveform of the estimated signal is shown in Figure 2d. The SNR of the estimated signal was −4.526 dB. We can see that the SNR improved significantly. The root mean square error (RMSE) was also used to evaluate our method. The RMSE is defined as

$$\mathrm{RMSE}=\sqrt{\frac{1}{N}{\Vert x-\widehat{x}\Vert }_2^2},$$

(33)

where $x$ denotes the noise-free simulated signal, $\widehat{x}$ denotes the estimated signal and $N$ denotes the length of the signal. The RMSE of the estimated signal is 0.5328. The noise level is reduced successfully.

The envelope spectrum analysis of the estimated bearing fault signal is shown in Figure 2e. Through the envelope analysis of the estimated signal, the fault feature frequency was $f_s$ = 135.1 Hz, which is equal to the simulated value (135.1 Hz). The comparison of Figure 2c,e shows that the proposed method can effectively extract the fault feature frequency from the noisy signal, even if the transient component is drowned by heavy noise.

The verification of the distribution of transient components and noisy bearing fault signals is shown in Figure 3. The distribution of transient components data is shown in Figure 3a. It is shown that the pdf of noise-free transient component was a mixture of two Laplace pdfs, where the estimated parameters were $k_1=0.416$, $k_2=0.584$, ${\sigma }_1=0.488$ and ${\sigma }_2=0.049$. The distribution of noisy data value is shown in Figure 3b. It can be found that the pdf of noisy value was a mixture of two LapGauss pdfs, where the estimated parameters $k_1=0.250$, $k_2=0.750$, ${\sigma }_1=0.447$ and ${\sigma }_2=0.186$.

4. Experimental Results

In this section, the cylindrical roller bearing installed on the test rig is taken as the test object. The test rig is shown in Figure 4. Three NJ208 (TMB) type cylindrical roller element bearing were used for testing. The vibration sensor model was a PCBM621B40 and the sensor had a reference sensitivity of 10 mv/g and the frequency range of measurable signals was 1.6 to 30,000 Hz. The data acquisition instrument model was an INV3018C. The test was carried out under the conditions of setting fault: a through crack fault with width and depth of 0.2 mm was set on the outer race, inner race and ball of the. rolling bearing respectively by using the wire-cut method. The load applied to the testing bearings was about 1 kN. The sampling frequency was 51.2 kHz and the Shaft rotating speed was 1496 r/min. The theoretical fault feature frequencies of the outer race, inner race and ball were 142.8 Hz, 206.2 Hz and 132.6 Hz respectively. And the corresponding periods were 7.00 ms, 4.85 ms, and 7.54 ms, respectively. As in the simulation, we used the discrete Daubechies 5 wavelet with three decomposition levels, and detail coefficients at each level were used for the estimation of parameters $k_1$, $k_2$, ${\sigma }_1$ and ${\sigma }_2$ utilizing the EM algorithm. The results for the outer race, inner race and ball are shown in Figure 5, Figure 6 and Figure 7 respectively.

4.1. Outer Race

In this section, the bearing outer race fault signal measured from the test rig is analyzed and the corresponding results are shown in Figure 5.

The time domain waveforms of the noisy outer race fault signal measured are shown in Figure 5a. The length of this signal was 76,800. The envelope spectrum analysis of the noisy bearing fault signal is shown in Figure 5b. Also, the fault feature frequency cannot be extracted directly. The results of the estimated transient components in the time domain is shown in Figure 5c. The envelope spectrum analysis of the estimated transient components is shown in Figure 5d. In Figure 5d, it can be seen that the fault feature frequency of the outer race was $f_1=142.8$ Hz, which is equal to the theoretical value (142.8 Hz).

4.2. Inner Race

In this section, the bearing inner race fault signal measured from the test rig is analyzed and the corresponding results are shown in Figure 6.

The time domain waveforms of the noisy inner race fault signal measured are shown in Figure 6a. The length of this signal was 76,800. The envelope spectrum analysis of the noisy bearing fault signal is shown in Figure 6b. The fault feature frequency also could not be obtained directly. The result of the estimated transient components in the time domain is shown in Figure 6c. The envelope spectrum analysis of the estimated transient components is shown in Figure 6d. From Figure 6d, it could be obtained that the fault feature frequency of the inner race was $f_2$ = 206 Hz, which is nearly equal to the theoretical value (206.2 Hz). Since there are always unavoidable errors in the measurement and calculation process, 206 Hz should be the correct required fault feature frequency.

4.3. Ball

In this section, the bearing ball fault signal measured from the test rig is analyzed and the corresponding results are shown in Figure 7.

The time domain waveforms of the noisy ball fault signal measured are shown in Figure 7a. The length of this signal was 76,800. The envelope spectrum analysis of the noisy bearing fault signal is shown in Figure 7b. The fault feature frequency also cannot be obtained directly. The result of the estimated transient components in the time domain is shown in Figure 7c. The envelope spectrum analysis of the estimated transient components is shown in Figure 7d. From Figure 7d, it can be obtained that the fault feature frequency of the ball was $f_3=132.7$ Hz, which is nearly equal to the theoretical value (132.6 Hz). Since there are always unavoidable errors in the measurement and calculation process, 132.7 Hz should be the correct required fault feature frequency.

Therefore, the proposed method can effectively extract the fault feature frequencies of three different types of faults from the corresponding noisy experimental signal. The aforementioned simulated and engineer tests were all implemented in MATLAB2016b on a computer whose operating system was Windows 10. The CPU parameter of the computer was Intel(R) Core(TM) i7-4790 at 3.60 GHz and the RAM was 12.0 GB, and the execution time of both simulated and engineering signals were less than 1 s. It can be seen that the time consumption of the proposed method is very low and stable.

4.4. Computational Complexity Comparison

Comparisons of the computational complexity between the proposed method and some other popular fault feature extraction methods are shown in

Table 1. Computational complexity comparison.

Methods	Spectral Kurtosis	Empirical Mode Decomposition	Time-Frequency Analysis	Sparse Representation	The Proposed Method
Computational complexity	$O(N\log N)$	$O(N\log N)$	$O(N\log N)$	$O(rN\log N)$	$O(N)$

.

Usually, the computational complexity of the spectral kurtosis method and empirical mode decomposition method are $O(N\log N)$. The frequently used time-frequency analysis methods include the short-time Fourier transform, continuous wavelet transform and Hilbert-Huang transform. The computational complexity of these methods is also $O(N\log N)$. The sparse representation method requires iterative optimization to solve many sub-problems of the whole model. At present, many mature iterative optimization algorithms have been applied to solve sparse representation models, such as the alternating direction method of multipliers (ADMM) [33], the split augmented Lagrangian shrinkage algorithm (SALSA) [34], majorization-minimization (MM) [35] and the forward-backward splitting (FBS) algorithm [36]. However, these iterative optimizations of sparse representation always involve a matrix-vector. This operation has $O(rN\log N)$ cost and $r$ is the redundancy factor [10]. Furthermore, calculating the inverse matrix may be involved sometimes. So, the computational complexity of the sparse representation method, if the construction of the sparse dictionary is not taken into account, is at least $O(rN\log N)$, which obviously limits the execution efficiency of the algorithm.

The EM algorithm used in the proposed method has computational complexity $O(N)$ when estimating model parameters iteratively. Besides, the computational complexity of the discrete wavelet transform is also $O(N)$. So, the computational complexity of the whole proposed method is $O(N)$ and the fault feature extraction results can be obtained quickly.

5. Discussion

The proposed method is expected to be applied to fault diagnosis of the rolling bearings used in modern mechanical equipment. In this paper, a simulated signal and three experimental signals with different fault types are analyzed, and the envelope spectrum of the estimated signal is compared with the envelope spectrum of the original noisy signal to determine whether the fault features can be effectively extracted. The results indicate that the proposed method can extract fault features accurately, although there is still some noise in the estimated signal. Once the fault feature is obtained, we can determine the type of failure in practical use. The signal data used in the experiments were collected from cylindrical roller bearings mounted on the test rig, which could effectively simulate the bearings working in mechanical equipment. Therefore, using these data for analysis is reliable and can be extended to practical bearing fault diagnosis after some improvements and testing. Furthermore, the proposed method has a great advantage, that is, the time required to extract the fault feature is very short (about 0.8 s), which is very important for quickly diagnosing bearing faults in practical applications. However, the proposed method has some inaccuracy in estimating the standard deviations of the mixture Laplace model, which may affect the noise level of the estimated bearing fault signal in the time domain. So, future research will focus on the application of the proposed method to actual mechanical equipment and finding a way to more accurately estimate the parameters of the mixture Laplace model while preserving the extraction time.

6. Conclusions

In this paper, an efficient method to quickly reduce noise level and extract fault feature frequencies of noisy bearing fault signals was proposed. This proposed method is based on the MAP estimator and pdfs of bearing fault source signal and noise, thereby having a reliable probabilistic interpretation. The proposed method first performed discrete wavelet transform on the noisy signal to obtain wavelet coefficients. Based on the wavelet coefficients obtained, the EM algorithm was adopted for parameter estimation of the derived model, and the newly estimated wavelet coefficients were then calculated using this model. Finally, the estimated transient components could be obtained after reconstruction. The computational complexity of this method was lower than some popular fault feature extraction methods and the fault feature extraction results could be obtained quickly. To validate the efficiency and accuracy of the proposed method, we analyzed a simulated bearing fault signal and three experimental signals with different types of faults using this method. The results indicate that this method can effectively and efficiently achieve the goal of noise level reduction and fault feature extraction in rolling element bearings.