Fault Characteristic Extraction by Fractional Lower-Order Bispectrum Methods

.e generated signals generally contain a large amount of background noise when the mechanical bearing fails, and the fault signals present nonlinear and non-Gaussian feature, which have heavy tail and belong to α-stable distribution (1< α< 2); even the background noises are also α-stable distribution process. .en it is difficult to obtain reliable conclusion by using the traditional bispectral analysis method under α-stable distribution environment. Two improved bispectrum methods are proposed based on fractional lower-order covariation in this paper, including fractional low-order direct bispectrum (FLODB) method, fractional low-order indirect bispectrum (FLOIDB) method. In order to decrease the estimate variance and increase the bispectral flatness, the fractional lower-order autoregression (FLOAR) model bispectrum and fractional lower-order autoregressive moving average (FLOARMA) model bispectrum methods are presented, and their calculation steps are summarized. We compare the improved bispectrum methods with the conventional methods employing second-order statistics in Gaussian and SαS distribution environments; the simulation results show that the improved bispectrum methods have performance advantages compared to the traditional methods. Finally, we use the improved methods to estimate the bispectrum of the normal and outer race fault signal; the result indicates that they are feasible and effective for fault diagnosis.


Introduction
Bispectral analysis based on high-order statistics is an effective tool to solve nonlinear phase coupling and non-Gaussian fault diagnosis [1,2]. e traditional bispectrum methods include nonparametric bispectrum [3,4], parametric AR bispectrum [5,6], parametric ARMA bispectrum [7], and their improved bispectrum methods [8]. e bispectrum of the signal contains not only the amplitude information but also the phase information. Bispectrum can effectively suppress the influence of Gaussian background noise and extract the non-Gaussian features hidden in the signal, and the graphics are intuitive. e fault feature of different nonlinear coupling modes in the bispectrum can be applied to quickly identify the working state of the bearing. erefore, bispectrum can better extract the signal features than the traditional power spectrum, which has been widely used in mechanical fault diagnosis [9][10][11]. In recent years, the methods of mechanical fault signal analysis have been developed. A new family of model-based impulsive wavelets and their sparse representation method are presented for rolling bearing fault diagnosis in [12]. An SVD principle analysis method based on the correlation coefficient is proposed for the bearing fault diagnosis in [13]. Subsequently, Qin et al. proposed a K-SVD algorithm with adaptive transient dictionary and transient feature extraction by the improved orthogonal matching pursuit [14]. Guo et al. applied the resonance demodulation and vibration separation to tooth root crack detection of planet and sun gears [15].
Recently, some extension methods based on the traditional bispectrum have also been studied and applied to the rotating machinery signal analysis, such as the deterministic bispectrum [16], modulation signal bispectrum [17,18], principal component bispectrum [19,20], and cyclic bispectrum [21]. Cheng et al. proposed a mechanical fault location and diagnosis method based on two bispectra and fuzzy clustering, which can effectively diagnose the fault state and location [22]. A new fault diagnosis method based on wavelet packet decomposition and modulation signal bispectrum analysis was proposed in [23]. e method firstly reconstructs the wavelet packet energy signal in time-frequency domain and then carries out modulation bispectrum analysis on the reconstructed signal, which can realize early fault diagnosis. Wang et al. proposed a bispectrum image texture features manifold method based on the support vector machines and genetic optimization algorithms for the rolling bearing vibration signal analysis [24]. A new bispectrum analysis method based on the optimal scale shape slice was proposed in [25], and the method can extract the fault signal from the sensitive modal component; hence, it can better extract fault features. A new variant modulation signal bispectrum method was introduced in [26], which was used to measure and analyze the fault current signals of the different mechanical motors, and the results show that this method is better than the traditional bispectrum method. e traditional and improved bispectrum methods have been applied in the rotating machinery signal analysis fields, but the methods still have some defects, and their performances degrade in impulsive environment and even fail. erefore, it is of great significance to explore high efficiency and performance of bispectrum analysis methods.
In actual working conditions, the mechanical bearings work in poor environment, the generated signals generally contain a large amount of background noise when the fault occurs, and the fault signals have obvious nonlinear and non-Gaussian properties, which belong to α-stable distribution process, and even the same with the noise in the signals [27][28][29][30]. Hence, it is difficult to find a solid conclusion by using the traditional bispectrum analysis methods. erefore, improved bispectrum methods which can be applicable to α-stable distribution environment need to be explored. Recently, α-stable distribution model was used for statistical modeling of the ocean environmental noise [31]. e adaptive cumulative distribution detector and blind estimation of frequency hopping parameters methods were proposed based on α-stable distribution model in [32,33]. Several improved frequency spectrum analysis methods have been introduced for α-stable distribution environment in [34], and the improved time frequency representation algorithms are proposed in [35], which have been applied to mechanical fault signal analysis.
In view of the performance degradation of the conventional bispectrum methods in α-stable distribution environment, the improved fractional low-order direct bispectrum and fractional low-order indirect bispectrum methods have been proposed for α-stable distribution environment in this paper, and the improved fractional lowerorder autoregression model bispectrum and fractional lower-order autoregressive moving average model bispectrum methods are presented for decreasing the estimate variance and increasing the bispectral flatness. We also summarize their calculation steps. e improved bispectrum methods and the traditional bispectrum methods employing second-order statistics are compared under Gaussian and α-stable distribution environments; the simulation results show that the improved bispectrum methods have performance advantages compared to the traditional methods. Finally, we apply the improved methods to estimate the bispectrum of the normal and outer race fault signal; the result indicates that the proposed methods are feasible and effective for fault diagnosis.
In this paper, several improved bispectrum analysis methods based on fractional lower-order statistics are proposed for mechanical bearing fault diagnosis in Gaussian or α-stable distribution noise environment. e paper is structured in the following manner. α-stable distribution and the bearing fault signals are introduced in Section 2. e improved fractional lower-order bispectrum methods are demonstrated, and the simulation comparisons employing the traditional bispectrum methods and the improved bispectrum methods are performed to show the advantage of the proposed methods in Section 3. e conjoint application simulations of the actual bearing fault signals employing the proposed bispectrum and time-frequency distribution in [35] are demonstrated in Section 4. Finally, the conclusions and future research are given in Section 5.

Bearing Fault Signals
e actual bearing fault signals data are obtained from the Case Western Reserve University (CWRU) bearing data center [36]. e experimental equipment adopts 6205-2RS JEM SKF type bearing, the outer race diameter is 20.472 inches, and the inner race and the ball diameter are 0.9843 inches and 0.3126 inches, respectively. e bearing outer race thickness is 0.5906 inches, motor load is 0 HP, and motor speed is 1797 rpm. e bearing faults of inner race, ball, and outer race are set, and the fault diameters are all 0.021 inches. e fault data are collected at 12,000 samples per second, and the outer race position relative to load zone is centered at 6:00. e normal signals are given in Figure 1(a), and the fault signals of inner race, ball, and outer race are shown in Figures 1(b)-1(d), respectively. We can know that the waveform of the fault signals has a certain impulse.
In order to further verify the pulse characteristics of the bearing fault signals, we use α-stable distribution statistical model to estimate the parameters of the inner race fault, ball fault, and outer race fault signals, and the results are given in Table 1 [37,38]. As it can be seen, the characteristic index of the normal signals is equal to 2, which is Gaussian distribution. However, the characteristic index of the bearing fault signals is greater than 1 but smaller than 2, and it belongs to non-Gaussian α-stable distribution (α < 2).
PDFs of the signals of inner race fault, ball fault, and outer race fault are shown in Figures 2(a)-2(c), respectively. From the PDFs of normal and fault signals, we can see that PDFs of fault signals have heavy tails. Most of the parameters β are approximately equal to zero in Table 1, and Figure 2 shows that PDFs of the fault signals are near symmetric. Hence, SαS distribution is a more concise and accurate statistical model for the bearing fault signals.

Fractional Lower-Order Direct Bispectrum
Method.
x(n), n � 0, 1, . . . , N − 1 { } is N samples of the observation data; its discrete fractional lower-order Fourier transform is defined as where P is a given real constant and P < α ≤ 2. 〈P〉 denotes P order moment of x(t), when x(t) is a real signal, Dividing the data x(n) { } into L segments, each segment is M points, and two adjacent segments overlap (M/2) points; then L � (2(N − (M/2))/M) and x i (n), i � 1, 2, . . . , L; n � 1, 2, . . . , M . According to equation (1), the discrete fractional low-order Fourier transform of the ith segment can be written as    x where P 1 + P 2 + P 3 < α, and we let . e corresponding fractional low-order direct bispectrum of the ith segment is given by where ω 1 ′ � (Iω 1 /2πf s ) and ω 2 ′ � (Iω 2 /2πf s ). Averaging the bispectrum of those L segments, fractional low-order direct bispectrum of x(n) { } can be gotten.
Fractional low-order direct bispectral estimation process is to first calculate the discrete fractional low-order Fourier transform coefficient and then compute its fractional loworder triple correlation and finally average fractional loworder triple correlation of all segments.

Fractional Lower-Order Indirect Bispectrum Method.
By using the fractional lower-order moment, we define discrete fractional lower-order three-order cumulants where 0 ≤ P 1 < (α/2), 0 ≤ P 2 < (α/2), and 0 ≤ } are real, then the estimation of FLOTOC is given by Mathematical Problems in Engineering where From equations (7)-(9), e ith segment of FLOTOC of the signal x(n) { } can be given by where i � 1, 2, . . . , L, averaging fractional low-order thirdorder cumulant of the L segments; we can get Taking the windowed two-dimensional discrete Fourier transform of equation (11), fractional low-order indirect bispectrum estimation of the signal x(n) { } can be given by where w (m 1 , m 2 ) is a two-dimensional window function and K < M − 1. e bispectrum estimation process of fractional low-order indirect method is to first calculate the discrete fractional low-order third-order cumulant of each segment, followed by averaging the fractional low-order third-order cumulant of all segments, and finally compute the two-dimensional Fourier transform.

Application Review.
In this simulation, the test signal y(n) is defined as where x(n) is three cosinoidal signals and v(n) is additive Gaussian noise or SαS distribution noise.
When v(n) is additive Gaussian noise, signal to noise ratio (SNR) can be used. But v(n) is SαS distribution noise, SNR is inapplicable, and generalized signal to noise ratio (GSNR) is written as where c is the dispersion coefficient of SαS distribution noise. According to the given GSNR, the amplitude of the signal x(n) is written as Letting SNR � − 5 dB and GSNR � 20 dB, the traditional bispectrum direct and indirect methods and fractional lower-order bispectrum direct and indirect methods are applied to estimate the bispectrum of the signal x(n) under Gaussian distribution noise and SαS distribution noise; the simulation results are shown in Figures 3-6.

Remarks.
e direct bispectrum estimations of the signal x(n) under Gaussian noise environment (SNR � − 5 dB) are shown in Figure 3. e fractional lower-order direct and indirect bispectrum methods have larger variance; we can increase the number of segments and segment length by overlapping adjacent segments to reduce variance and add two-dimensional window function to improve the frequency resolution of bispectral estimation. e proposed fractional low-order bispectral (FLOB) estimation methods still have biperiodicity and symmetry, so they can be used to quickly calculate the bispectrum of the signal. Fractional lower-order bispectral biperiodicity properties are as follows: Fractional lower-order bispectral symmetry properties are as follows: Direct method bispectrum estimation

Fractional Lower-Order AR Model Bispectrum Method
Principle. e conventional real p order AR model process can be expressed as e traditional AR model bispectrum of the signal X(n) in (18) is given by where β ⌢ is the estimate of the third moment of the driving noise, H(ω) is the system transfer function, and A fractional lower-order AR model SαS process x(n) may be written as where p is order of the AR model and a m (i � 1, 2, . . . , p) are its parameters, which are real numbers. u(n) is an independent identically distributed (i.i.d) SαS random process, α is its characteristic index, and c u is its dispersion coefficient.
Indirect method bispectrum three-dimensional graph x(n) can be expressed by a finite impulse response (FIR) system [27].
where h(m) { } is system impulse response coefficient and Taking z transformation for equations (21) and (22) and obtaining the system transfer function, Taking α-order moment with equation (23), we have z � e jω , and |z| � 1 on the unit circle; then According to the definition of the AR model bispectrum in (19) and α-order moment in (24) and (25), we define fractional lower-order AR model bispectrum (FLOARB) as and FLOARB on the unit circle is written as where c u is the dispersion coefficient of the driving SαS random process u(n) and * is conjugate operation. According to the definition of the fractional lower-order three-order cumulants in where Φ is named as fractional low-order three-order cumulants matrix (FLOTOCM), which is Toeplitz. From solving the AR bispectrum coefficients equations and fractional low-order moment matrix equations in [30], we let en,
In this paper, we apply the final prediction error (FPE) criteria to determine the order p of fractional low-order AR model. When p increases gradually from 1, FPE will be the minimum at a certain p, which just is the most appropriate order. e calculation formula can be written as where σ 2 a is the variance of residuals. We summarize the steps of the FLOAR model bispectrum method as follows: Step 1: computing fractional lower-order three-order cumulants of the signal x(n) with equations (7)-(9) and constructing fractional low-order three-order cumulants matrix in (28).

Application Review.
In this simulation, the experimental signal is y(n) in equation (13). e performances of the traditional AR model bispectrum method and the improved FLOAR model bispectrum method are compared under Gaussian distribution noise (SNR � − 5 dB) and SαS distribution noise (α � 1.2; GSNR � 20 dB), and the simulation results are shown in Figures 7-10.
In order to further verify the advantages of FLOAR model bispectrum method, we conduct comparative experiment on two methods under different α when GSNR � 20 dB, and their parameter estimations are shown in Figure 11. When α � 1.3 and GSNR changes from 14 dB to 24 dB, we compare the changes of the errors power with the AR and FLO-AR model bispectrum methods under α-stable distribution noise environment, and the simulation is given in Figure 10. e AR and FLOAR model bispectrum estimations of the signal x(n) under SαS distribution noise (α � 1.2; GSNR � 20 dB) are given in Figure 8. e simulation result shows that the conventional AR bispectrum method fails under SαS noise environment in Figures 8(a) and 8(b), but the improved FLOAR bispectrum method shows good toughness in Figures 8(c) and 8(d). As a result, the AR bispectrum method is only suitable to analyze the signals in Gaussian environment, but the FLOAR bispectrum method can work in Gaussian and SαS noise environment, which is robust.

Fractional Lower-Order
where p and q are orders of the AR and MA model, respectively. a m (i � 1, 2, . . . , p) and b m (m � 0, 1, 2, . . . , q) are their parameters, which are real numbers. u(n) is an independent identically distributed (i.i.d) SαS random process, α is its characteristic index, and its dispersion coefficient is c u . Simplifying equation ( where a 0 � 1 and b 0 � 1. Taking z transform on both sides of equation (33), we obtain where h(m) { } is system impulse response coefficient,H(z) is system transfer function, and u(n − m) { } is i.i.d SαS random process.
e AR model parameters changes with a m (i � 1, 2, . . . , p); A(z) and B(z) are the FIR and IIR filters, respectively.
Taking α-order moment with equation (36), z � e jω , and |z| � 1 on the unit circle; then According to the definition of the ARMA model bispectrum and α-order moment of H(z) in (38) and (39), we define fractional lower-order ARMA model bispectrum (FLOARMAB) as and FLOARB on the unit circle is written as FLOARMAB e jω 1 , e jω 2 � c u H e jω 1 H e jω 2 H * e jω 1 + e jω 2 , where c u is the dispersion coefficient of the driving SαS random process u(n) and * is conjugate operation.
To obtain the coefficients of the fractional lower-order ARMA model a m (i � 0, 1, 2, . . . , p) and b m (m � 0, 1, 2, . . . , q), we should multiply x(n − i) on both sides of equation (32) taking fractional lower-order covariance to get where It can be shown [24] that, for the fractional lower-order covariance of the signal and noise in (45),   Equations (47) and (48) are a generalized Yule-Walker equation, and we can obtain the fractional lower-order AR model (FLOARM) coefficients a k (i � 0, 1, 2, . . . , p) by solving them.
If 0 ≤ m ≤ q, equation (42) changes as Letting p � 0, we have We can obtain the coefficients b k (k � 1, 2, . . . , q) of the fractional lower-order MA model by solving the nonlinear equation (35).
A finite q ′ -order FLOMA model can be equivalent by an approximate infinite p ′ -order FLOAR modelp ′ ≫ q ′ ; we have where e(m) is error. Letting e(m) � 0, we obtain We should multiply a 〈P〉 (n − i) on both sides of equation (52). We take fractional lower-order covariance to get where . We can also obtain the coefficients b i (i � 1, 2, . . . , q) of the FLOMA model by solving the Yule-Walker equation (54).
We summarize the steps of the FLOARMA model bispectrum method as follows: Step 1: solving fractional low-order bispectrum Yule-Walker equation in (48) and getting the coefficients a k (k � 0, 1, 2, . . . , p) of the FLOAR model.
Step 3: determining the order p and q employing FPE criterion.

Application Review.
In this simulation, y(n) in equation (14) is used as the experimental signal. e ARMA model and FLOARMA model bispectrum methods are compared to demonstrate their performance under Gaussian distribution noise (SNR � − 5 dB) and SαS distribution noise (α � 1.2; GSNR � 20 dB), and the simulations are shown in Figures 9-12.

Remarks.
e ARMA model and FLOARMA model bispectrum estimations of the signal x(n) are shown in Figure 9 under Gaussian noise environment (SNR � − 5 dB). Figures 9(a) and 9(b) are the existing ARMA model bispectral estimation and its three-dimensional graph estimation, respectively. FLOARMA model bispectral estimation and its three-dimensional graph are given in Figures 9(c) and 9(d), respectively. e results show that both methods can estimate out the bispectrum of the signal x(n) well under Gaussian noise. Figure 10 shows

Application Simulations
e impulse of the outer race fault signals in the vibration position of DE, FE, and BA is generated because of the local defects of rolling element bearings, and the waveforms are given in Figure 3(d) and Table 1. We can know that the fault signals are nonstationary and non-Gaussian α-stable distribution process. In this section, the experiment signal adopts from the normal signal and the bearing outer race fault signal in the vibration position of DE, and 0.2-second data is selected as the test signal, which is collected at 12,000 samples per second; then N � 2400. We apply the improved FLODB, FLOIDB, FLOAR, and FLOARMA bispectrum methods to analyze the normal and bearing outer race fault signals, and the simulations are shown in Figures 11 and 12. In this section, we have only extracted the first quadrant of the bispectral representation (f 1 > 0, f 2 > 0) to analyze the signals. Figures 11(a), 11(c), 11(e), and 11(g) are the bispectral contour maps of the bearing normal signal employing the FLODB, FLOIDB, FLOAR, and FLOARMA bispectrum methods, respectively. eir bispectral three-dimensional diagrams are shown in Figures 11(b), 11(d), 11(f ), and 11(h), respectively. It is observed that the improved methods can effectively suppress the noise and estimate out the frequency components of the normal signal, and the spectral peaks exist near the central frequency of 1060 Hz; hence, the transient harmonic vibration components of the normal signal are about 1060 Hz. Figures 12(a), 12(c), 12(e), and 12(g) are the bispectral contour maps of the bearing out race signal employing the FLODB, FLOIDB, FLOAR, and FLOARMA bispectrum methods, respectively, and their bispectral three-dimensional diagrams are shown in Figures 12(b), 12(d), 12(f), and 12(h), respectively. It is observed that the spectral peaks of Figures 12(a), 12(c), and 12(e) exist near the central frequencies of 600 Hz and 2800 Hz, and those in Figure 12(g) exist near the central frequencies of 600 Hz, 2800 Hz, and 3500 Hz. Hence, the outer race signal has fault harmonic vibration components of 600 Hz, 2800 Hz, and 3500 Hz, indicating that the bearing outer race is damaged. Figure 12(i) clearly shows the gap between the impacts regularly changing. e interval between impulses A, B, C, D, E, and F is approximately 30 ms; then the characteristic frequency of the bearing outer race fault is about 33.333 Hz. Figure 12(j) shows wave shape and spectrum of the envelope signal of the outer race fault signal using resonance-demodulation approach; we can clearly see the pulse separation. We can see that FLODB and FLOIDB methods have larger variance; the FLOARMA bispectrum method has lower variance and higher frequency resolution, and its performance is optimal.
To further verify the advantages of the proposed fractional low-order bispectrum methods, SαS distribution noise (α � 1.1; GSNR � 20 dB) is added in the α-stable distribution outer race fault signal as the actual working

Conclusions
e bearing fault signals are a non-Gaussian and nonstationary process, and α-stable distribution is a more appropriate statistical model for them. e improved FLODB, FLOIDB, FLOAR, and FLOARMA model bispectrum methods have been proposed for the fault signals employing fractional low-order statistics. e improved methods can be applicable to Gaussian and α-stable distribution noise environment, and their performances are superior to the existing direct bispectrum method, indirect bispectrum method, and AR and ARMA model bispectrum analysis methods. Fractional low-order nonparametric bispectrum estimation methods, FLODB and FLOIDB, require a large number of data samples and have a large estimation variance, but the fractional low-order parametric bispectrum estimation, FLOAR and FLOARMA model bispectrum, has small variance and produces fewer parameters that describe the characteristics of the target; hence, it can be directly used for target features. We can apply the improved methods to analyze the α-stable distribution bearing fault signal, even α-stable distribution noise environment, and the fault characteristic frequency, the dominant frequency, and the other fault frequency features of the fault signals can be clearly obtained. Combining the fractional low-order timefrequency methods, more fault characteristics can be obtained, and the joint diagnosis will be realized for the bearing fault signals. In the future, we can also apply the bispectrum diagonal slice to reflect the coupling information between

22
Mathematical Problems in Engineering fault signals so as to realize the recognition of the fault characteristics. e complete mechanical bearing fault state spectrum can be established based on fractional low-order bispectrum estimation for the fault signals, which can provide a new way for the fault diagnosis and online monitoring.
Data Availability e supplementary file in txt file format is the original experimental data of the paper.

Conflicts of Interest
e authors declare that they have no conflicts of interest.