Sparse Fast Fourier Transform and its application in intelligent diagnosis system of train rolling bearing

2.1. Algorithm

SFFT output $\widehat{x}\mathrm{\text{'}}$ satisfies the following guarantee [4]:

where $ϵ$ is constant and $\delta =1/n^{O\left(1\right)}$ is an accuracy parameter, $y$ is the $k$-sparse vector minimizing ${‖\widehat{x}-y‖}_2$.

The main idea of SFFT [5] is shown below and a flow diagram is seen in Fig. 2:

1) Run location loops $L$ times and return $L$ sets of coordinates $I_1$, $I_2$,…, $I_L$.

2) Count every found coordinate $i$, of which the occurrences are added up to the number $s_i$: $s_i=\left|\right\{r|i\in I_r\}|$.

3) Retain the coordinates that appeared at least half of the location loops $I\mathrm{\text{'}}=\{i\in I_1\cup ...\cup I_L|s_i>L/2\}$.

4) On the set $I\mathrm{\text{'}}$, operate a number of $L$ estimation loops and return $L$ sets of frequency coefficients ${\widehat{x}}_{I^{\mathrm{\text{'}}}}^r$.

5) Estimate every frequency coefficient ${\widehat{x}}_i$ as${\widehat{x}}_i=\mathrm{m}\mathrm{e}\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{n}\left\{x_i^r|r\in \left\{1,...,L\right\}\right\}\text{,}$ where respectively take the median in real and imaginary part.

2.1.2. Estimation loops

The function of estimation loops, whose implementation is similar to location loops’, is to reconstruct the exact coefficient values by the given coordinates set $I$.

The first 3 steps share the same ones to location loops’.

Given a set of coordinates $I$, estimate $\widehat{x}$ as $\widehat{x}{\mathrm{\text{'}}}_i={\widehat{z}}_{h_{\sigma }\left(i\right)}{w}^{\tau i}/{\widehat{G}}_{o_{\sigma }\left(i\right)}$.

Fig. 2A flow diagram of SFFT

2.2. Implementation of SFFT

Randomly generate a group of signals $x$ with the length $n=\mathrm{}$512. Permuting $x$ in time domain with the permutation $P_{\sigma ,\tau }$ causes change in the frequency signal $\widehat{x}$ (only keep half of frequency domain), See Figs. 3, 4.

Set $k=$ 4, $B=\mathrm{}$16, compute the time domain signal $y=G{P}_{\sigma ,\tau }x$, of which the spectrum is large around the large coordinates of $P_{\sigma ,\tau }x$, such as Fig. 5(a). Next compute $\widehat{z}$, which is the rate $B$ subsampling of $\widehat{y}$ as plotted in Fig. 5(b), where red pots refer to data points which perform FFT after subsampling.

During estimation loops, the algorithm estimates ${\widehat{x}}_i$ based on the value of the nearest coordinate in $\widehat{z}$. During location loops, the algorithm chooses the largest $k$ coordinates of $\widehat{z}$, and selects the elements closest to those coordinates mapped the candidate interval $J$ as Fig. 5(c). It outputs the set $I$ of original signals of those elements. At this point, the inner loop ends. Finally combining the outer loop, location and value of large coefficient can be obtained in original signal spectrum as Fig. 5(d).

Fig. 3Permutation Pσ,τx of time-domain signal x (σ= 125, τ= 207)

a) Original signal

b) Permuted signal

Fig. 4Permutation P^σ,τx of frequency-domain signal x^ (σ=125, τ=207)

a) Original signal

b) Permuted signal

2.3. Compare SFFT with FFT

To evaluate the performance of SFFT, this paper compares SFFT with FFTW [6,7] from two aspects. Currently, FFTW is the fastest public implementation for computing the Discrete Fourier Transform and its computational complexity is ${\rm O}\left(n\mathrm{l}\mathrm{o}\mathrm{g}n\right)$.

Fix the sparsity parameter $k=\mathrm{}$50 and respectively perform SFFT as well as FFT several times. The runtime of the compared algorithms for different signal sizes $N$ is plotted in Fig. 6(a). The figure indicates that for signal sizes 100000$<N<\mathrm{}$67108864, SFFT is faster than FFT at $k=\mathrm{}$50, which proves the superiority of SFFT over FFT when processing enormous data.

Fix the data size $N=\mathrm{}$2²² and respectively perform SFFT as well as FFT several times. The runtime of the compared algorithms for different sparsity parameter $k$ is plotted in Fig. 6(b). The figure indicates that SFFT runs faster than FFT when $k$ ranges in certain sparse scope, that is to say, we can use less data points to realize better fast algorithm.

Fig. 5An example of signal permutation during inner loop and outer loop

a) Convolved signal

b) Subsampling and FFT

c) Estimate passing regions

d) Final result

Fig. 6Contrast result

a) Run time vs signal size ($k=\mathrm{}$50)

b) Run time vs sparsity ($N=$2²²))

3.1. Fault feature frequencies of rolling bearing

It’s very known that the bearing has four fault feature frequencies, corresponding to outer ring failure ($F_{ord}$), inner ring failure ($F_{\mathrm{i}rd}$), rolling body failure ($F_{bd}$) and retainer failure ($F_c$) [8]. These four kinds of fault feature frequencies are defined as Eqs. (2)-(5):

where $f_r$ is relative rotational frequency between inner ring and outer ring, $n$ is the number of rolling body, $d$ is the diameter of rolling body, $D$ is the pitch diameter of rolling bearing and $\varphi$ is contact angle (Fig. 7).

The datum used in this paper come from trains bearing testing table. The corresponding bearing parameters are $f_r=\mathrm{}$29.83 Hz, $n=\mathrm{}$9, $d=\mathrm{}$0.00794 m, $D=\mathrm{}$0.03904 m, $\varphi =\mathrm{}$0 rad. Then the four kinds of fault frequencies can be calculated as 106.93 Hz, 161.54 Hz, 128.68 Hz, 11.88 Hz respectively.

Fig. 7Bearing structure

3.2. Fault features extraction for rolling bearing

The validity of fault features has decisive influence on the fault diagnosis system. Another paper of authors Li et al. [9] has verified that amplitude ratios in frequency domain and kurtosis factor, crest factor, clearance factor, impulse factor, shape factor in time domain have good performance on fault diagnosis for train bearing operating in complex environment. In this paper, these features are also adopted to test the proposed method.

3.2.1. Feature in frequency domain

The feature extraction process is shown in Fig. 8. It contains vibration signal analysis with wavelet de-noising [10], envelope spectrum extraction [11] and SFFT.

In this paper, outer ring is taken as the example to test the proposed method. Fig. 8(d) is the spectrogram of the envelope signal of outer ring after SFFT. Fig. 8(e) is the first 600 Hz spectrogram in Fig. 8(d). It can be seen that the amplitude is the largest at 111.33 Hz and the energy is concentrated nearby, so this fault is estimated to be part of outer ring according to Eq. (2). Fig. 8(f) is spectrogram of envelope signal of outer ring performing the traditional FFT. Comparing Fig. 8(d) with Fig. 8(f), it can be seen that the spectrogram tendency of SFFT approximately keeps the same as FFT’s and verifies the validity of SFFT on fault feature extraction.

Feature amplitude ratios are defined as the ratios between the amplitude of outer or inner ring fault frequency, and the amplitude of rolling body fault frequency. In this paper, feature amplitude ratios are extracted on the basis of envelope signal after SFFT, showed in Fig. 8(d). Two feature amplitude ratios are listed as below:

6

$R_1=\frac{A\left(F_{ord}\right)}{A\left(F_{bd}\right)},$

7

$R_2=\frac{A\left(F_{ird}\right)}{A\left(F_{bd}\right)},$

where $A$ is the amplitude of fault feature frequency.

Respectively 20 groups of data samples under normal bearing, fault on outer ring, inner ring and rolling body are adopted. Then the fault feature amplitude ratios of every sample are computed by Eqs. (6), (7). Table 1 is one group of the feature amplitude ratios which is selected randomly.

Fig. 8Flow diagram for fault feature extraction of rolling bearing

It can be seen from Table 1, when bearing fault appears on outer ring, $R_1$ is a notable large value while $R_2$ is notably small. When bearing fault appears on inner ring, $R_1$ is a notable small value while $R_2$ is notably large. When fault appears on rolling body or bearing is normal, $R_1$ and $R_2$ are relatively steady but also can be distinguished obviously. Therefore, selecting feature amplitude ratio as the input vector of bearing fault classifier is effective.

Table 1Extract feature amplitude ratios

Type	Normal	Outer ring fault	Inner ring fault	Rolling body fault
$R_1$	1.2308	6.6330	0.5456	0.0262
$R_2$	0.6739	0.1286	10.9076	0.5529

3.3. Fault diagnosis of rolling bearing

An appropriate fault classifier is needed to be designed for fault diagnosis. There are several kinds of classifiers, of which Artificial Neural Network (ANN) is a highly interconnected network of a large number of processing elements called neurons. It quickly computes a large number of operations by parallel distributed processing [12]. Comparing with decision trees, Logistic Regression (LR) [13] model, ANN can better deal with ordered attributes in most cases [14] and approximate any complex nonlinear relationship and complicatedly partition feature space [15]. ANN is proved to be a powerful data-driven, self-adaptive computational tool which can accurately obtain nonlinear and complex underlying characteristics of physical process [16] and is widely used in fault diagnosing.

ANN attains knowledge about a problem by being trained with known examples. If the network is trained appropriately, it can be effectively used to solve some ‘unknown’ or ‘untrained’ instances of the problem [17]. ANN is constituted by multiple layers: an input layer, hidden layer (s) and an output layer. In each layer, neurons act as processing elements. Biases are added to all the input neurons. Network weights play an important role on the acceleration or deceleration of the input signals. In order to get the output of a neuron, the weighted sum of the inputs to each neuron should be passed through an activation function: threshold, piecewise-linear or sigmoid function.

There are dozens of neural networks models, of which ANN with back-propagation (BP) learning algorithm proposed by McClelland and Rumelhart [18] is widely used in variously classifying problems. This paper utilizes three-layered BP neural networks to realize the intelligent fault classification. BP neural networks must be firstly trained with training set, see its standard architecture in Fig. 9. Once the training pattern is completed, the network can be used for validation pattern. BP neural networks generally adopt error back propagation algorithm [19]. The main steps are as follows:

1) Initialize the synaptic weights and the biases of the input layer and output layer to random values.

2) A group of training examples are accepted by BP neural networks, signal is forward propagated through the network and finally computes the continuous output value for each output unit.

3) Compare the output value to the desired value and generate an error value, which can be evaluated by the mean square error.

4) Neural networks back propagate the error value to adjust the weights and biases of the neurons according to the principle of decreasing error value.

5) Step 2 to step 4 are iterated until the error value meets the minimum overall error.

Fig. 9Structure of BP neural network

To test the method proposed in this paper, the activation function of three-layered BP neural network is fixed as tangent-sigmoid function (tansig) at hidden layer, and a linear transfer function (purelin) at output layer is designed. The training function is selected as gradient descent algorithm (trainscg), which has a fast convergence speed. Assign the largest number of training to 1000, learning rate to 0.05 and target error to 0.01. Next, 80 groups of samples belonging to 4 kinds of healthy states (each kind has 20 groups) are obtained from the train bearing table. Then 2 feature vectors in frequency domain and 5 feature vectors in time domain mentioned above are calculated. The diagnosis process is described as following:

1) Constitute the input vector. Every 7 feature parameters constitute an input vector and all vectors in the same state merge into one vector set.

2) Fault type encoding. Normal bearing (1,0,0,0), outer ring fault (0,1,0,0), inner ring fault (0,0,1,0), rolling body fault (0,0,0,1).

3) Ensure parameters. Input 7 neurons referring to 7 feature parameters. Output 4 neurons referring to 4 kinds of healthy states. Choose 12 neurons of hidden layer according to empirical equations.

4) Train samples. For inspecting classification effect, choose 15 groups of feature vectors as training samples, and choose 5 groups as diagnosed samples in every healthy state. The error curve of BP neural network for training is showed in Fig. 10, where the best training performance that can be found is 0.005 at 1000 times, which takes 1 second and meets the accuracy requirement.

5) Fault classification. With the neural network having been trained, classify the diagnosed samples under 4 kinds of healthy states. Take the rolling bearing with fault on outer ring as the example, its accuracy of diagnosis result displayed in Table 3 reaches 100 %. After overall diagnosis, each accuracy of 4 kinds of operating states comes to be 100 %, 100 %, 100 %, 80 % in Table 4.

Fig. 10The training error curve of BP neural network. (Best training performance is 0.0052647 at epoch 1000)