Sparse Envelope Spectra for Feature Extraction of Bearing Faults Based on NMF

Periodic impulses and the oscillation response signal are the vital feature indicators of rolling bearing faults. However, finding the suitable feature frequency band is usually difficult due to the interferences of other components and multiple resonance regions. According to the characteristics of non-negative matrix factorization (NMF) on a spectrogram, the feature extraction method from a sparse envelope spectrum for rolling bearing faults is proposed in this paper. On the basis of the time–frequency distribution (TFD) of the periodic transient oscillations, the basic matrix can be interpreted as the spectral bases, and the time weight matrix corresponding to spectral bases can be extracted by NMF. Because the bases and the weights have a one-to-one correspondence, the frequency band filtering with the basic component and the time domain envelope of the weight vector are calculated respectively. Then, the sparse envelope spectrum can be derived by the inner product of the above results. The effectiveness of the proposed method is verified by simulations and experiments. Compared with band-pass filtering and spectral kurtosis methods, and considering the time weights and corresponding the spectral bases for the periodic transient oscillations, the weak fault-rated feature can be enhanced in the sparse spectrum, while other components and noise are weakened. Therefore, the proposed method can reduce the requirement of selecting frequency band filtering.


Introduction
As an important component in rotating machinery, the rolling bearing is a frequent area in machinery failure which can bring about economic losses and casualties [1].Therefore, the condition monitoring and fault diagnosis of the rolling bearings are important means of health assessment, and the feature extraction based on vibration signals [2,3] is of significant interest to researchers.
Theoretically, once the rolling bearing has a localized fault (such as pitting, spalling, etc.), which creates the periodic impulses during the rotation of rotor, the periodic resonance of structures can be excited in the vibration signal.Thus, the frequency of the transient oscillation is decided by the natural frequency of the mechanical system, and the repetition frequency of the transient oscillation corresponds to the fault-related impulses [4].The periodic vibration response signal is clearly non-stationary because the magnitude of oscillations decreases exponentially, and the signal having a discretized spectrum and wide frequency band of the high rotational speed of bearing is shown in Figure 1a,b.The discrete interval of frequency is the repetition period of the impact excitation.However, in incipient stages, due to the lower amplitude of the transient oscillation, the large noise and other components often interfere with the identification of discrete frequencies [5], as shown in Figure 1c.Therefore, the feature extraction of periodic impulses is proposed with various detection methods, such as spectrum analysis, cepstral analysis [6], envelope analysis, and time-frequency analysis.
Appl.Sci.2019, 9, x FOR PEER REVIEW 2 of 23 natural frequency of the mechanical system, and the repetition frequency of the transient oscillation corresponds to the fault-related impulses [4].The periodic vibration response signal is clearly nonstationary because the magnitude of oscillations decreases exponentially, and the signal having a discretized spectrum and wide frequency band of the high rotational speed of bearing is shown in Figure 1a,b.The discrete interval of frequency is the repetition period of the impact excitation.However, in incipient stages, due to the lower amplitude of the transient oscillation, the large noise and other components often interfere with the identification of discrete frequencies [5], as shown in Figure 1c.Therefore, the feature extraction of periodic impulses is proposed with various detection methods, such as spectrum analysis, cepstral analysis [6], envelope analysis, and time-frequency analysis.
(a) (b) (c) According to the frequency representation of the vibration response caused by the fault-related impulses, the envelope demodulation can be used to extract the repetitive frequency of the transient oscillations with band-pass filtering and envelope transformation [7].However, when the impact excitation has low energy, the magnitude of the exponential decay signal is also small and leads to difficulty in the selection of a resonance region.Moreover, it often interacts with the frequency band of other components and the excitation response of random impulses.It is obvious that the optimal filtering frequency band cannot be obtained easily.Thus, to capture the signature of repetitive transients, the spectral kurtosis (SK) is proposed [8,9].SK selects the best band-pass filter based on spectral kurtosis variation.However, the kurtosis index is susceptible to various outliers that are independent of the fault [10], so the selected frequency band could be unable to extract valid fault information.
Unlike the response of random impact excitation, the periodic transient oscillation corresponds to the discrete components of the resonance frequency band.However, due to the influence of noise Amplitude Amplitude According to the frequency representation of the vibration response caused by the fault-related impulses, the envelope demodulation can be used to extract the repetitive frequency of the transient oscillations with band-pass filtering and envelope transformation [7].However, when the impact excitation has low energy, the magnitude of the exponential decay signal is also small and leads to difficulty in the selection of a resonance region.Moreover, it often interacts with the frequency band of other components and the excitation response of random impulses.It is obvious that the optimal filtering frequency band cannot be obtained easily.Thus, to capture the signature of repetitive transients, the spectral kurtosis (SK) is proposed [8,9].SK selects the best band-pass filter based on spectral kurtosis variation.However, the kurtosis index is susceptible to various outliers that are independent of the fault [10], so the selected frequency band could be unable to extract valid fault information.
Unlike the response of random impact excitation, the periodic transient oscillation corresponds to the discrete components of the resonance frequency band.However, due to the influence of noise and other components, it is difficult to observe the discretization spectra.In this case, the time-frequency analysis provides a localized feature representation for the repetitive transients, and the time-frequency distribution (TFD) of the signal has been adopted to identify the response of the system, such as the natural frequency and damping ratio [11][12][13][14].Moreover, due to splitting the signal into different frequency bands and central frequencies, the wavelet transform is also used for frequency band optimization.In related research, Qiu [5] constructed a suitable filter for envelope analysis by selecting the optimal parameters of the wavelet filter, and Su [15] proposed the optimal Morlet wavelet filter.Yiakopoulos [16] also adopted the wavelet transform of the original signal to envelope the low-frequency feature.The results denote that the performances are very dependent on the wavelet basis and the practical applicability is limited.
In addition, the short time Fourier transform (STFT) has also been adopted to represent the non-stationary information of the vibration signal [17][18][19].In view of the non-negative property of the TFD, the non-negative matrix factorization (NMF) is suitable for the analysis of vibration signals due to the ability of the parts-based representation.Delong [19] proposed the constrained NMF of the spectral matrix and used the decomposed basis matrix to extract features.After utilizing the semi-binary NMF of the spectrum, Wodecki [18] selected the appropriate decomposition frequency band for slicing analysis based on the principle of maximum kurtosis.Meanwhile, Wodecki [17] also adopted standard NMF on spectrograms and utilized the basis vector to filter frequency bands.In addition to vibration signals, the speech signal has also been analyzed while performing NMF [20][21][22].Unlike traditional band-pass filtering, the NMF-based methods can decompose frequency bands according to time-frequency distribution adaptively, and the interference component in the bands cannot be eliminated.
By applying NMF to the TFD, the frequency distribution of the vibration response can be decomposed into the base matrix, which can be interpreted as the spectral bases and the weight matrix which is the time envelope corresponding to the spectral bases.For the periodic transient oscillations, the spectral bases can reflect the natural frequency band, while the magnitude change of the periodic transient oscillations is also represented by the weight matrix.Due to the influence of noise and the decomposed frequency band, neither the band-pass filtering based on the frequency base matrix nor the envelope of the time weight matrix can clearly extract the impact features.However, since the bases and the weights have a one-to-one correspondence, it means that the same feature component with different amplitudes can exist in the band-pass filtering and the envelope of time weights.Thus, the inner product of the above results can be used to obtain a sparse envelope spectrum, which can enhance the components of fault features and weaken interference components.
The rest of this paper is organized as follows.Section 2 gives a brief description of NMF and provides the characteristics of NMF of time-frequency distribution, and then introduces the principles and processes of the proposed method.Section 3 verifies the proposed method with simulated signals and fault bearing data.Conclusions are drawn in Section 4.

NMF
With the non-negative constraints, NMF is also widely used in image recognition, data mining etc. [23][24][25].Given a non-negative matrix where m and n are the numbers of samples and features respectively, NMF is designed to decompose the matrix V into the product of two non-negative matrices W = , such that Equation (1) holds: where W and H are the basis matrix and the weight matrix, respectively, and r is the low rank of the embedding space.Normally, the choice of r is to satisfy (m + n) < mn.W can be regarded as a set of bases that linearly approximate the observed array V, and H is the non-negative projection weight of the observed array V on the basis W.
Due to the purely additive expression, NMF provides a good explanation of the local characteristics of data.In addition, from the perspective of clustering, the decomposition result of NMF is clustering for the non-negative data [26].Ding [27] proved the consistency among NMF, K-means and spectral clustering methods.We can interpret the columns of W as the posterior probabilities for row clustering and the rows of H as the posterior probabilities for column clustering.
For the solution of Equation ( 1), the objective function based on Euclidean distance is written as: where ||•|| 2 F denotes the Frobenius norm.Using the gradient descent method for this optimization problem, the classical multiplicative iterative rules are given by: The multiplicative iterative algorithm is simple and easy to use, but the poor performance and convergence are also pointed out in paper [28,29].Therefore, improved algorithms are proposed, such as the alternating least squares algorithm [30], the projected quasi-Newton algorithm [31], the projected gradient algorithm [28], and the block principal pivoting algorithm [32].The block principal pivoting algorithm is developed with the non-negative least squares (NNLS) framework.
To obtain good convergence, the main iterative rules are as follows: Meanwhile, block principal rotation [33] is used to improve the convergence of the algorithm.Hence, the block principal pivoting algorithm is chosen for analyzing complex vibration signals in this paper.

Time Frequency Separation of Time-Frequency Matrix
For the original signal x(t), its STFT coefficient is defined as: where w(t − k) is a window function, t k is the time corresponding to each slip window and N is the overall number of the computed frames.Considering the frequency band selection and the strong noise in vibration signals, the window function with low side-lobes is applied, such as the hanning window, hamming window, etc. From the coefficient STFT(t, f ), the TFD is set to V = |STFT(t, f )|, and then the NMF algorithm is adopted to decompose V.According to the clustering, the sub-matrix W = w 1 w 2 • • • w r will be interpreted as spectral bases; that is, the cluster of the frequency domain.Each w i includes some frequency regions.Meanwhile, the magnitude envelope of signals belonging to the frequency region is represented by weight vector h i .This correspondence can be explained by the bilinear NMF model [34].In this model, the non-negative data matrix is approximately represented by a sum of the linear combination of rank-one non-negative matrices: where E represents error or noise.Note that the frequency base vector w i and the time weight vector h i have a one-to-one correspondence.It can be considered that the feature extraction of signals is realized either by spectral basis or by the time weight matrix.
Figure 2 depicts the correspondence between w i and h i .The original signal shown in Figure 2 is mixed with the two non-continuous sinusoidal components of f 1 = 1000 Hz and f 2 = 2300 Hz, where the component of f 1 contains three non-zero segments and the component of f 2 contains four segments.Applying NMF with r = 2, two spectral bases and two time weight vectors are obtained from the TFD of the original signal.It can be clearly seen that w 1 and h 1 are the descriptions of the sinusoidal signal of frequency f 1 in the frequency domain and the time envelope respectively.Likewise, the sinusoidal signal of the frequency f 2 is described in both the time domain and the frequency domain by the basis w 2 and the weight vector h 2 .Therefore, NMF on a spectrogram is a powerful method to describe each transient oscillation and its position in time with a set of spectral bases and weights.
Appl.Sci.2019, 9, x FOR PEER REVIEW 5 of 23 is represented by weight vector hi.This correspondence can be explained by the bilinear NMF model [34].In this model, the non-negative data matrix is approximately represented by a sum of the linear combination of rank-one non-negative matrices: where E represents error or noise.Note that the frequency base vector wi and the time weight vector hi have a one-to-one correspondence.It can be considered that the feature extraction of signals is realized either by spectral basis or by the time weight matrix.
Figure 2 depicts the correspondence between wi and hi.The original signal shown in Figure 2 is mixed with the two non-continuous sinusoidal components of f1 = 1000 Hz and f2 = 2300 Hz, where the component of f1 contains three non-zero segments and the component of f2 contains four segments.Applying NMF with r = 2, two spectral bases and two time weight vectors are obtained from the TFD of the original signal.It can be clearly seen that w1 and h1 are the descriptions of the sinusoidal signal of frequency f1 in the frequency domain and the time envelope respectively.Likewise, the sinusoidal signal of the frequency f2 is described in both the time domain and the frequency domain by the basis w2 and the weight vector h2.Therefore, NMF on a spectrogram is a powerful method to describe each transient oscillation and its position in time with a set of spectral bases and weights.

Feature Extraction Based on Sparse Envelope Spectrum
As a spectral basis, each column vector of W represents different frequency regions.Thus, the original signal x(t) can be filtered with the base wi containing the fault characteristic information, and the filtered spectrum is given by: This can be followed by performing the inverse Fourier transform (IFT) on the filtered spectrum;

Feature Extraction Based on Sparse Envelope Spectrum
As a spectral basis, each column vector of W represents different frequency regions.Thus, the original signal x(t) can be filtered with the base w i containing the fault characteristic information, and the filtered spectrum is given by: This can be followed by performing the inverse Fourier transform (IFT) on the filtered spectrum; the filtered waveform is calculated as follows: Next, the envelope signal is obtained from the filtered signal x f (t), which is defined as: where x f (t) is the Hilbert transform of the signal x f (t), and t−τ dτ.Finally, removing the zero-frequency component, the spectrum of the envelope signal is computed by Fourier transform (FT): After computing the envelope spectrum, the feature information of the fault frequency can be observed.However, the acquired vibration signal includes not only the impact excitation caused by bearing defects, but also the external impacts.The envelope spectrum usually contains a lot of noise and other interference components.Therefore, to improve performance, it is necessary to further analyze the selected filter band.
Generally speaking, each spectral basis w i describes the frequency region of the transient oscillation, while the weight vector h i indicates its time envelope.Thus, performing Fourier transform on h i can display the frequency components of the time envelope: Obviously, the band-pass filtering is realized by the spectral base w i for the transform of the whole signal.The weight vector h i is obtained from the TFD directly, where the frequency component of each slip moment is calculated by the local signal.This difference leads to the different amplitude distributions of the two envelope spectra.However, the periodic impact component, as the main characteristic, has a larger amplitude in both envelope spectra.Therefore, to enhance the fault-related frequency, the product of the two squared spectra is defined as: The squares of the components represent the energy of the frequency component, which can highlight a larger component than that of the amplitude spectrum [3].The inner product, denoted by "•", can strengthen the common frequency of the envelope spectra and reduce other frequencies.As a result, the envelope spectrum becomes sparse, and the feature frequency of bearing faults can be clearly observed naturally.The algorithm of the sparse envelope spectrum extracted by combining W and H is described as follows.
Step 1: Calculating the STFT coefficients of the vibration signal with a low side-lobes window function and obtaining its TFD matrix.To better extract the impact response waveform, it is recommended to use a smaller window length; Step 2: Applying NMF on the TFD matrix with low rank r.Considering that the frequency spectrum often contains multiple frequency regions, r is set to k ≤ r ≤ k + 3, where k is the estimated number of resonance regions.By properly increasing the value of r, the similar frequency distribution in the decomposition results can be avoided; Step 3: Removing spectral basis with the low frequency region in W; Step 4: Computing the filtered spectrum according to Equation ( 9) by the remaining spectral basis; Step 5: Calculating the filtered signal and its envelope spectrum H W ( f ); Step 6: Obtaining the envelop spectrum H H ( f ) of each weight vector's corresponding frequency basis; Step 7: Getting the sparse envelope spectra with Equation ( 12); Step 8: Selecting the sparse spectrum with the maximum amplitude at the fault frequency in H( f ) as the final result.

Numerical Experiment
To evaluate the performance of the proposed method, according to the proposed vibration model of a rolling bearing [35,36], the synthetic signal x(t) is generated by: The exponentially decaying component usually indicates a localized fault in rolling bearings or external excitation, such as electromagnetic interference.A k is the amplitude of impulse.T is the period of impulse, f h = 1/T is the fault feature frequency, f k is the natural frequency and B k and ϕ k are the damping ratio and phase respectively.It is obvious that the vibration response excited by a bearing fault is periodic.The second item of Equation ( 16) corresponds to the rotor fault, which is represented by the frequency harmonic component f o .P m is the amplitude of the harmonic, and the last item n(t) is Gaussian white noise.
Considering the low frequency position for the rotor harmonic fault, the impact excitations are only computed in an experiment.The outer ring fault-related frequency f h is set to 80 Hz, and the excited natural frequency is set to 2200 Hz.Meanwhile, three non-periodic transient oscillations are added, where the carrier frequencies are 1000 Hz, 2600 Hz and 3000 Hz, and the fractional bandwidths are set to 0.1, 0.3 and 0.1.The phases are all zero.Finally, the Gaussian white noise with signal-to-noise ratio (SNR) = 0 is added, and the vibration response signal is sampled at 10,240 Hz.The components of the simulated signal are described in Figure 3, and its waveform and spectrum are shown in Figure 4, where the red line denotes the envelope of the valid fault component.In Figure 4, because of the strong background noise and the random transient interference, it is difficult to directly observe the impact response caused by the local fault and its frequency component.
The Hanning window with 33 samples and 97% overlap rate is used for the STFT of the signals.Then, NMF is performed on the TFD of the simulated signal, where the low rank dimension r = 5. Figure 5 illustrates the decomposed spectral basis and the time weight vectors.For the distribution of spectral basis W, the basis w 5 has a wide frequency band of about 2000 Hz, which can contain the periodic transient oscillations of 2200 Hz and the transient impact excitation component of 2600 Hz.The frequency basis w 2 focuses on the frequency band around 3000 Hz and corresponds to the transient oscillation components of natural frequency 3000 Hz.The frequency basis w 3 is concentrated in the frequency band of 1000 Hz, which is the frequency feature of the other transient impact excitation components.In addition, the bases w 1 and w 4 correspond to the noise components.As for the distribution of the time weight matrix H, the two peaks h 2 at approximately 0.14 s and 0.22 s correspond to the two impulses whose natural frequency is 3000 Hz; the peak h 3 at approximately 0.35 s corresponds to the impulse whose natural frequency is 1000 Hz.Useful feature components of other time weight vectors are difficult to observe.
bandwidths are set to 0.1, 0.3 and 0.1.The phases are all zero.Finally, the Gaussian white noise with signal-to-noise ratio (SNR) = 0 is added, and the vibration response signal is sampled at 10,240 Hz.The components of the simulated signal are described in Figure 3, and its waveform and spectrum are shown in Figure 4, where the red line denotes the envelope of the valid fault component.In Figure 4, because of the strong background noise and the random transient interference, it is difficult to directly observe the impact response caused by the local fault and its frequency component.The Hanning window with 33 samples and 97% overlap rate is used for the STFT of the signals.Then, NMF is performed on the TFD of the simulated signal, where the low rank dimension r = 5. Figure 5   The spectral bases w 2 , w 3 , and w 5 , which contain the impact excitations, are selected to extract impulse responses.The sparse envelope spectra constructed with three frequency bases and weight vectors are given in Figure 6, in which most of the feature components are significantly enhanced and noise is almost eliminated.In the spectrum of Figure 6c, the frequency component of 80 Hz (denoted by the arrow) is effectively extracted in the sparse envelope spectrum constructed with the bases w 5 and h 5 .The envelope spectra extracted by other bases such as w 2 and w 3 do not have the certain feature components because both point to aperiodic impact.It is noted that the impulsive features are extracted, which confirms that the sparse envelope spectrum offers huge advantages for the representation of feature information.To compare the effects of the proposed method, the results of the envelope demodulation with band-pass filtering [17] are also presented in Figure 7.The feature component of 80 Hz can be observed by the frequency band filtering with basis w 5 in Figure .Compared with Figure 6c, it can be seen that no clear fault features can be observed due to the influence of in-band noise although the frequency basis w 5 corresponds to the frequency band of 2200 Hz.This illustrates that the position and bandwidth of band-pass filtering can seriously affect the envelope effect.Considering the time domain envelope of the impulse series, as shown in Figure 8, the proposed method can obtain a sparse envelope spectrum.Therefore, the frequency feature of band-pass filtering is greatly enhanced, and is eliminated for in-band noise.
Appl.Sci.2019, 9, x FOR PEER REVIEW 9 of 23 vectors are given in Figure 6, in which most of the feature components are significantly enhanced and noise is almost eliminated.In the spectrum of Figure 6c, the frequency component of 80 Hz (denoted by the arrow) is effectively extracted in the sparse envelope spectrum constructed with the bases w5 and h5.The envelope spectra extracted by other bases such as w2 and w3 do not have the certain feature components because both point to aperiodic impact.It is noted that the impulsive features are extracted, which confirms that the sparse envelope spectrum offers huge advantages for the representation of feature information.To compare the effects of the proposed method, the results of the envelope demodulation with band-pass filtering [17] are also presented in Figure 7.The feature component of 80 Hz can be observed by the frequency band filtering with basis w5 in Figure.
Compared with Figure 6c, it can be seen that no clear fault features can be observed due to the influence of in-band noise although the frequency basis w5 corresponds to the frequency band of 2200 Hz.This illustrates that the position and bandwidth of band-pass filtering can seriously affect the envelope effect.Considering the time domain envelope of the impulse series, as shown in Figure 8, the proposed method can obtain a sparse envelope spectrum.Therefore, the frequency feature of band-pass filtering is greatly enhanced, and is eliminated for in-band noise.For comparison, Figure 9a presents the kurtogram [9] of the same signal, from which the optimal demodulation band is 960-1280 Hz and the waveform and spectrum of the filtered signal are presented in Figure 9b,c.The extracted result indicates that the filtered signal with maximum kurtosis still has noise interference and the transient feature of the waveform is not very clear.When envelope analysis is applied to this signal, little fault frequency can be obtained, as shown in Figure 9d.In fact, we know that kurtosis is a measure of the peakedness, so it easily tends to highlight the outliers of the signal caused by noise.For comparison, Figure 9a presents the kurtogram [9] of the same signal, from which the optimal demodulation band is 960-1280 Hz and the waveform and spectrum of the filtered signal are presented in Figure 9b,c.The extracted result indicates that the filtered signal with maximum kurtosis still has noise interference and the transient feature of the waveform is not very clear.When envelope analysis is applied to this signal, little fault frequency can be obtained, as shown in Figure 9d.In fact, we know that kurtosis is a measure of the peakedness, so it easily tends to highlight the outliers of the signal caused by noise.

Case 1
To further explore the effectiveness of the proposed method, the data from the Case Western Reserve University Bearing Data Center [37] is used for the verification.The test stand consists of a two-horsepower motor with a torque transducer and a dynamometer to apply different loads.Vibration data was collected using accelerometers, which were attached to the housing with magnetic

Case 1
To further explore the effectiveness of the proposed method, the data from the Case Western Reserve University Bearing Data Center [37] is used for the verification.The test stand consists of a two-horsepower motor with a torque transducer and a dynamometer to apply different loads.Vibration data was collected using accelerometers, which were attached to the housing with magnetic

Case 1
To further explore the effectiveness of the proposed method, the data from the Case Western Reserve University Bearing Data Center [37] is used for the verification.The test stand consists of a two-horsepower motor with a torque transducer and a dynamometer to apply different loads.Vibration data was collected using accelerometers, which were attached to the housing with magnetic bases.Smith [38] conducted a comprehensive analysis of each type of fault in the data set and summarized all the data into three types of faults, namely, obvious (Y), weak (P), and cannot be diagnosed (N).Faults marked as "P" or "N" are recommended for the performance analysis of the newly proposed algorithm.
The fault of the inner ring bearing numbered 277DE, marked as "P", is adopted for identification.The sensor for 277DE data is located at the motor driving end, and far from the fan end of the fault bearing, which means the test signal is coupled with a large amount of interference components.The analyzed data set consists of 4096 sampling data with a sampling frequency of 12 kHz, and its waveform and frequency spectrum are shown in Figure 10.According to the structural parameters of the bearing and the rotation speed, the characteristic frequency of the fault bearing is f i = 142.9Hz.In the frequency spectrum of Figure 10b, distinct frequency bands can be observed except for the harmonic frequency of the rotation speed, and the high frequency region is located in the range of 1900-4000 Hz.
Appl.Sci.2019, 9, x FOR PEER REVIEW 12 of 23 bases.Smith [38] conducted a comprehensive analysis of each type of fault in the data set and summarized all the data into three types of faults, namely, obvious (Y), weak (P), and cannot be diagnosed (N).Faults marked as "P" or "N" are recommended for the performance analysis of the newly proposed algorithm.The fault of the inner ring bearing numbered 277DE, marked as "P", is adopted for identification.The sensor for 277DE data is located at the motor driving end, and far from the fan end of the fault bearing, which means the test signal is coupled with a large amount of interference components.The analyzed data set consists of 4096 sampling data with a sampling frequency of 12 kHz, and its waveform and frequency spectrum are shown in Figure 10.According to the structural parameters of the bearing and the rotation speed, the characteristic frequency of the fault bearing is fi = 142.9Hz.In the frequency spectrum of Figure 10b, distinct frequency bands can be observed except for the harmonic frequency of the rotation speed, and the high frequency region is located in the range of 1900-4000 Hz.After STFT is computed with the same parameters in the simulation experiment, the TFD matrix is decomposed by NMF with low rank r = 5; the basis and weight vectors are shown in Figure 11.Because the high frequency bands usually have small noise interference, for the five decomposed frequency regions in Figure 11, the spectral bases w3, w4, w5 are located at the frequency band of 1900-4000 Hz and selected to extract the impulse response.The three selected frequency regions and the sub-matrix of corresponding weight are adopted to calculate the inner product of the squared envelope, and the sparse envelope spectra are displayed in Figure 12.
Amplitude Amplitude After STFT is computed with the same parameters in the simulation experiment, the TFD matrix is decomposed by NMF with low rank r = 5; the basis and weight vectors are shown in Figure 11.Because the high frequency bands usually have small noise interference, for the five decomposed frequency regions in Figure 11, the spectral bases w 3, w 4, w 5 are located at the frequency band of 1900-4000 Hz and selected to extract the impulse response.The three selected frequency regions and the sub-matrix of corresponding weight are adopted to calculate the inner product of the squared envelope, and the sparse envelope spectra are displayed in Figure 12.
In Figure 12a, constructed by frequency basis w 3 and weight vector h 3 , the fault characteristic frequency f i can be observed clearly.The sparse envelope spectrum constructed by the spectral bases w 4 and w 5 are shown in Figure 12b,c.Thereinto, an obvious rotational frequency of 29.2 Hz and double harmonic components are found, while the frequency of the inner ring fault is smaller.Although the frequency band regions of the three frequency bases are all located in the 1900-4000 Hz region, the three sparse envelope spectra are different.It can be supposed that w 3 corresponds to the vibration response excited by fault, while w 5 can be generated by other components, such as motor vibration excitation or bearing friction.Due to the overlapping of multi-frequency bands, the bases w 3 and w 5 are very close, as shown in Figure 13.Therefore, except for w 3 , the sparse envelope spectrum constructed by the basis w 5 and the weight h 5 also contains part of the features of the periodic impact response.In Figure 12a, constructed by frequency basis w3 and weight vector h3, the fault characteristic frequency fi can be observed clearly.The sparse envelope spectrum constructed by the spectral bases  In Figure 12a, constructed by frequency basis w3 and weight vector h3, the fault characteristic frequency fi can be observed clearly.The sparse envelope spectrum constructed by the spectral bases region, the three sparse envelope spectra are different.It can be supposed that w3 corresponds to the vibration response excited by fault, while w5 can be generated by other components, such as motor vibration excitation or bearing friction.Due to the overlapping of multi-frequency bands, the bases w3 and w5 are very close, as shown in Figure 13.Therefore, except for w3, the sparse envelope spectrum constructed by the basis w5 and the weight h5 also contains part of the features of the periodic impact response.To compare the effects, band-pass filtering is also applied.Figure 14 shows the envelope demodulation results with the three frequency bands.Because of the large in-band noise, the fault characteristic frequency is influenced by other frequency components.Compared with the result of Figure 12, due to the combination of the time domain envelope of the impulse response shown in Figure 15, the proposed method can obtain the sparse envelope spectrum constructed with the frequency basis w3 and the weight h3, thereby effectively extracting the fault characteristic frequency.To compare the effects, band-pass filtering is also applied.Figure 14 shows the envelope demodulation results with the three frequency bands.Because of the large in-band noise, the fault characteristic frequency is influenced by other frequency components.Compared with the result of Figure 12, due to the combination of the time domain envelope of the impulse response shown in Figure 15, the proposed method can obtain the sparse envelope spectrum constructed with the frequency basis w 3 and the weight h 3 , thereby effectively extracting the fault characteristic frequency.
vibration response excited by fault, while w5 can be generated by other components, such as motor vibration excitation or bearing friction.Due to the overlapping of multi-frequency bands, the bases w3 and w5 are very close, as shown in Figure 13.Therefore, except for w3, the sparse envelope spectrum constructed by the basis w5 and the weight h5 also contains part of the features of the periodic impact response.To compare the effects, band-pass filtering is also applied.Figure 14 shows the envelope demodulation results with the three frequency bands.Because of the large in-band noise, the fault characteristic frequency is influenced by other frequency components.Compared with the result of Figure 12, due to the combination of the time domain envelope of the impulse response shown in Figure 15, the proposed method can obtain the sparse envelope spectrum constructed with the frequency basis w3 and the weight h3, thereby effectively extracting the fault characteristic frequency.Additionally, the fast spectral kurtosis method is also used to analyze the signal.The results are shown in Figure 16.The optimal demodulation band is located in the 2625-3125 Hz region, and the waveform and spectrum of the filtered signal are presented in Figure 16b,c.Compared with the proposed method, although the envelope spectrum obtained by the fast spectral kurtosis method contains fault characteristic frequencies, it also contains many interference components, which increases the difficulty of identification and is not as simple and intuitive as the sparse spectrum.Additionally, the fast spectral kurtosis method is also used to analyze the signal.The results are shown in Figure 16.The optimal demodulation band is located in the 2625-3125 Hz region, and the waveform and spectrum of the filtered signal are presented in Figure 16b,c.Compared with the proposed method, although the envelope spectrum obtained by the fast spectral kurtosis method contains fault characteristic frequencies, it also contains many interference components, which increases the difficulty of identification and is not as simple and intuitive as the sparse spectrum.

Case 2
Another experiment is carried out on the test bench of mechanical equipment failure simulation, and its structure is shown in Figure 17.The rotating shaft is driven by an alternating current (AC) motor, and fault simulation can be conducted by replacing the fault motors.The motor bearing is 6203.The inner ring peeling fault is simulated on the inner ring of the motor drive end bearing.When the rotation frequency f r is 20 Hz, the fault characteristic frequency f i is 99.5 Hz.The accelerometer is installed at the vertical direction of the motor end bearing.The vibration signal is acquired with the sampling frequency of 20,480 Hz and the sample size of 4096. Figure 18 shows the waveform and frequency spectrum of the acquired signal.In Figure 18a, weak impulse features are buried in interference and noise, while it is shown that from the spectrum shown in Figure 18b, there are multiple resonance frequency bands, and the rotating characteristic frequencies are submerged.

Case 2
Another experiment is carried out on the test bench of mechanical equipment failure simulation, and its structure is shown in Figure 17.The rotating shaft is driven by an alternating current (AC) motor, and fault simulation can be conducted by replacing the fault motors.The motor bearing is 6203.The inner ring peeling fault is simulated on the inner ring of the motor drive end bearing.When the rotation frequency fr is 20 Hz, the fault characteristic frequency fi is 99.5 Hz.The accelerometer is installed at the vertical direction of the motor end bearing.The vibration signal is acquired with the sampling frequency of 20,480 Hz and the sample size of 4096. Figure 18 shows the waveform and frequency spectrum of the acquired signal.In Figure 18a, weak impulse features are buried in interference and noise, while it is shown that from the spectrum shown in Figure 18b, there are multiple resonance frequency bands, and the rotating characteristic frequencies are submerged.First, the proposed method is applied to the sampled signal.After computing STFT with the same parameters, the TFD matrix from the STFT coefficients is decomposed by NMF with low rank r = 5, and the decomposed bases are shown in Figure 19.The five spectral bases correspond to different frequency bands of 895-7400 Hz.Afterwards, the corresponding five band-pass filtering with frequency bases from w 1 to w 5 and the time envelope based on the weight matrix H are computed.The resulting sparse envelope spectra are shown in Figure 20.In Figure 20b, the fault characteristic frequency f i is effectively extracted, in which the frequency basis w 2 corresponds to the resonance frequency band from 1400-3800 Hz.In addition to the frequency basis w 3 , the other filtering bands cannot extract the corresponding fault feature frequency.First, the proposed method is applied to the sampled signal.After computing STFT with the same parameters, the TFD matrix from the STFT coefficients is decomposed by NMF with low rank r = 5, and the decomposed bases are shown in Figure 19.The five spectral bases correspond to different frequency bands of 895-7400 Hz.Afterwards, the corresponding five band-pass filtering with frequency bases from w1 to w5 and the time envelope based on the weight matrix H are computed.The resulting sparse envelope spectra are shown in Figure 20.In Figure 20b, the fault characteristic frequency fi is effectively extracted, in which the frequency basis w2 corresponds to the resonance frequency band from 1400-3800 Hz.In addition to the frequency basis w3, the other filtering bands cannot extract the corresponding fault feature frequency.Furthermore, the band-pass filtering for the five bands and the fast spectral kurtosis method are also analyzed.The results are shown in Figures 21 and 22, respectively.In the ranges of the five decomposed frequency bands, the fault component is not displayed in w 1, w 3, w 4, and w 5 due to the noise component, and the feature frequency in w 2 is also inconspicuous.From the result extracted by the fast spectral kurtosis method, it is difficult to detect the existence of fault characteristic frequencies because of its large bandwidth and interference components.Compared with Figure 20d, the advantage of the proposed sparse envelope spectrum is shown again.Besides the characteristic frequency, the spectrum of band-pass filtering also contains a large number of other components, which are easy to interfere with the judgment of fault.Additionally, the inaccurate frequency band extracts the small amplitude of the characteristic frequency.Obviously, the sparse envelope spectrum has an excellent visualization effect, so the product of time weight and filtering envelopes also can reduce the requirement of parameter selection for band-pass filtering.advantage of the proposed sparse envelope spectrum is shown again.Besides the characteristic frequency, the spectrum of band-pass filtering also contains a large number of other components, which are easy to interfere with the judgment of fault.Additionally, the inaccurate frequency band extracts the small amplitude of the characteristic frequency.Obviously, the sparse envelope spectrum has an excellent visualization effect, so the product of time weight and filtering envelopes also can reduce the requirement of parameter selection for band-pass filtering.(e)

Conclusions
In this paper, based on the non-negative matrix factorization of the time-frequency distribution of periodic transient oscillations, we propose feature extraction for the characteristic frequencies of rolling bearing faults, which takes into account the non-stationary characteristic of impact excitation for localized faults.The application of simulations and experiments shows that, based on timefrequency distribution, non-negative matrix factorization not only decomposes the frequency band

Conclusions
In this paper, based on the non-negative matrix factorization of the time-frequency distribution of periodic transient oscillations, we propose feature extraction for the characteristic frequencies of rolling bearing faults, which takes into account the non-stationary characteristic of impact excitation for localized faults.The application of simulations and experiments shows that, based on time-frequency distribution, non-negative matrix factorization not only decomposes the frequency band adaptively, but also extracts the time domain envelope of the signal of the corresponding frequency band.Using the inner product of the band-pass filtering and time weight envelopes, the sparse envelope spectrum is realized to enhance the fault-related frequency component.Thus, the strict requirements of band-pass filtering on the position and bandwidth of the filter band are reduced.Compared with band-pass filtering and fast spectral kurtosis methods, the weak fault-rated features are easier to observe in the sparse spectrum.In practice, the running speed of many pieces of equipment, such as a wind-powered gearbox, is not stable.Therefore, fault feature extraction under variable speed is the work of further research.

Figure 1 .
Figure 1.Simulation signal of bearing fault.The carrier frequency is 2200 Hz and the fault feature frequency (fh) is 80 Hz: (a) Simulation signal waveform; (b) Discretized frequency spectrum of simulation signal; (c) Spectrum of simulation signal with added noise.

Figure 1 .
Figure 1.Simulation signal of bearing fault.The carrier frequency is 2200 Hz and the fault feature frequency (f h ) is 80 Hz: (a) Simulation signal waveform; (b) Discretized frequency spectrum of simulation signal; (c) Spectrum of simulation signal with added noise.

Figure 2 .
Figure 2. Non-negative matrix factorization (NMF) on time-frequency distribution (TFD).The lower left part is the basis matrix W, and the upper right part is the weight matrix H.

Figure 2 .
Figure 2. Non-negative matrix factorization (NMF) on time-frequency distribution (TFD).The lower left part is the basis matrix W, and the upper right part is the weight matrix H.

Figure 4 .Figure 4 .
Figure 4. (a) Waveform of the simulated signal; (b) Spectrum of the simulated signal.

Figure 5 .Figure 5 .
Figure 5. NMF on the TFD of the simulated signal.The spectral bases w2, w3, and w5, which contain the impact excitations, are selected to extract impulse responses.The sparse envelope spectra constructed with three frequency bases and weight

Figure 6 . 5 Figure 6 .
Figure 6.Sparse envelop spectrums is obtained with decomposed sub-matrices: (a) Basis w 2 and weight vector h 2 ; (b) Basis w 3 and weight vector h 3 ; (c) Basis w 5 and weight vector h 5 .

Figure 9 .
Figure 9. Fault characteristics extracted by the fast spectral kurtosis method: (a) Kurtogram of the simulated signal, the global maximum is achieved for f c = 960 Hz and Level = 6; (b) Waveform of the filtered signal; (c) Spectrum of the filtered signal; (d) Envelope spectrum of the filtered signal.

Figure 10 .
Figure 10.Inner ring fault bearing: (a) Waveform of the fault signal; (b) Frequency spectrum of the fault signal.

Figure 10 .
Figure 10.Inner ring fault bearing: (a) Waveform of the fault signal; (b) Frequency spectrum of the fault signal.

Figure 11 .Figure 12 .
Figure 11.NMF on the TFD of the inner ring fault bearing.

3 w and weight vector 3 h 5 w and weight vector 5 h
; (b) Basis 4 w and weight vector 4 h ; (c) Basis .fi = fault characteristic frequency.

3 w and weight vector 3 h 5 w and weight vector 5 h
; (b) Basis 4 w and weight vector 4 h ; (c) Basis .fi = fault characteristic frequency.

Figure 12 .
Figure 12.Spare envelope spectra are obtained with decomposed sub-matrices: (a) Basis w 3 and weight vector h 3 ; (b) Basis w 4 and weight vector h 4 ; (c) Basis w 5 and weight vector h 5 .f i = fault characteristic frequency.

Figure 13 .
Figure 13.Position of three frequency band filtering.

3 Figure 13 .
Figure 13.Position of three frequency band filtering.

Figure 13 .
Figure 13.Position of three frequency band filtering.

Figure 16 .f
Figure 16.Fault characteristics extracted by the fast spectral kurtosis method: (a) Kurtogram of the simulated signal, the global maximum is achieved for = 2625 c f Hz and Level = 4.6; (b) Waveform of the filtered signal; (c) Spectrum of the filtered signal; (d) Envelope spectrum of the filtered signal.

Figure 16 .
Figure 16.Fault characteristics extracted by the fast spectral kurtosis method: (a) Kurtogram of the simulated signal, the global maximum is achieved for f c = 2625 Hz and Level = 4.6; (b) Waveform of the filtered signal; (c) Spectrum of the filtered signal; (d) Envelope spectrum of the filtered signal.

Figure 20 .
Figure 20.Obtained sparse spectrum with decomposed sub-matrix: (a) Basis w 1 and weight vector h 1 ; (b) Basis w 2 and weight vector h 2 ; (c) Basis w 3 and weight vector h 3 ; (d) Basis w 4 and weight vector h 4 ; (e) Basis w 5 and weight vector h 5 .

Figure 22 .
Figure 22.Fault characteristics extracted by fast spectral kurtosis method: (a) Kurtogram of the simulated signal, the global maximum is achieved for f c = 6400 Hz and Level = 3.6; (b) Waveform of the filtered signal; (c) Spectrum of the filtered signal; (d) Envelope spectrum of the filtered signal.