# Fault Characteristic Extraction by Fractional Lower-Order Bispectrum Methods

Authors: Haibin Wang; Junbo Long; Zeliang Liu; Fang You
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

The 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. Then 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.


---

## Body

## 1. Introduction

Bispectral analysis based on high-order statistics is an effective tool to solve nonlinear phase coupling and non-Gaussian fault diagnosis [1, 2]. The traditional bispectrum methods include nonparametric bispectrum [3, 4], parametric AR bispectrum [5, 6], parametric ARMA bispectrum [7], and their improved bispectrum methods [8]. The 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. The fault feature of different nonlinear coupling modes in the bispectrum can be applied to quickly identify the working state of the bearing. Therefore, bispectrum can better extract the signal features than the traditional power spectrum, which has been widely used in mechanical fault diagnosis [9–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]. The 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. The 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. Therefore, 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–30]. Hence, it is difficult to find a solid conclusion by using the traditional bispectrum analysis methods. Therefore, 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].The 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 lower-order 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. The 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. The paper is structured in the following manner. α-stable distribution and the bearing fault signals are introduced in Section 2. The 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. The 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.


## 2. Bearing Fault Signals

The actual bearing fault signals data are obtained from the Case Western Reserve University (CWRU) bearing data center [36]. The 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. The bearing outer race thickness is 0.5906 inches, motor load is 0 HP, and motor speed is 1797 rpm. The bearing faults of inner race, ball, and outer race are set, and the fault diameters are all 0.021 inches. The fault data are collected at 12,000 samples per second, and the outer race position relative to load zone is centered at 6:00. The 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.Figure 1
Waveforms of the bearing fault signals. (a) Normal signals in DE and FE. (b) Inner race fault signals in BA, DE, and FE. (c) Ball fault signals in BA, DE, and FE. (d) Outer race fault signals in BA, DE, and FE.
(a)(b)(c)(d)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).Table 1
The -stable distribution parameters of the bearing fault signals.
ParametersαβγμNormalDE2.000−0.28630.05320.0121FE2.0001.0000.05830.0236Inner race faultBA1.76820.08720.05900.0062DE1.41950.01550.24070.0175FE1.83500.03220.14950.0291Ball faultBA1.97900.05920.02930.0055DE1.86970.12150.07720.0193FE1.998−0.03710.06740.0321Outer race faultBA1.6077−0.17310.05300.0012DE1.10960.04330.13410.0367FE1.5435−0.01690.09680.0296PDFs of the signals of inner race fault, ball fault, and outer race fault are shown in Figures2(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.Figure 2
PDFs of the normal and bearing fault signals. (a) PDFs of the normal and inner race fault signals in DE, FE, and BA. (b) PDFs of the normal and ball fault signals in DE, FE, and BA. (c) PDFs of the normal and outer race fault signals in DE, FE, and BA.
(a)(b)(c)


## 3. Fractional Lower-Order Nonparametric Bispectrum Methods

### 3.1. Fractional Lower-Order Direct Bispectrum Method

xn,n=0,1,…,N−1 is N samples of the observation data; its discrete fractional lower-order Fourier transform is defined as (1)X˜ω=∑n=0N−1xPne−jωn,where P is a given real constant and P<α≤2. P denotes P order moment of xt, when xt is a real signal, xPt=xtP⋅signxt, signxt=1,xt>0,0,xt=0,−1,xt<0,, and when xt is a complex signal, xPt=xtP−1⋅x∗t, where ∗ denotes conjugate operation.Dividing the dataxn into L segments, each segment is M points, and two adjacent segments overlap M/2 points; then L=2N−M/2/M and xin,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(2)X˜iω=∑n=1MxiPne−jωn.The corresponding discrete fractional low-order Fourier transform coefficient is defined as(3)X˜iφ=1M∑n=1MxiPne−j2πnφ/M,φ=0,1,…,M2.Let the sampling frequency of the signalxn be fs and the total frequency sampling points be I; then K=M/I−1/2. Combining equation (3), fractional low-order triple correlation of xn can be expressed as(4)B^iφ1,φ2=I2fs2∑k1=−KK∑k2=−KKX˜i′φ1+k1X˜i″φ2+k2X˜i″′−φ1−φ2−k1−k2,X˜i′φ1+k1=1M∑n=1MxiP1ne−j2πnφ1+k1/MX˜i″φ2+k2=1M∑n=1MxiP2ne−j2πnφ2+k2/M,X˜i″′−φ1−φ2−k1−k2=1M∑n=1MxiP3ne−j2πn−φ1−φ2−k1−k2/M,where P1+P2+P3<α, and we let P1=P2=P3<α/3, 0≤φ2≤φ1, and φ2+φ1≤fs/2. The corresponding fractional low-order direct bispectrum of the ith segment is given by(5)FLODBiω1,ω2=I2fs2∑k1=−KK∑k2=−KKX˜i′ω1′+k1X˜i″ω2′+k2X˜i″′−ω1′−ω2′−k1−k2,where ω1′=Iω1/2πfs and ω2′=Iω2/2πfs. Averaging the bispectrum of those L segments, fractional low-order direct bispectrum of xn can be gotten.(6)FLODBω1,ω2=1L∑i=1LFLODBiω1,ω2.Fractional low-order direct bispectral estimation process is to first calculate the discrete fractional low-order Fourier transform coefficient and then compute its fractional low-order triple correlation and finally average fractional low-order triple correlation of all segments.


### 3.2. Fractional Lower-Order Indirect Bispectrum Method

By using the fractional lower-order moment, we define discrete fractional lower-order three-order cumulants (FLOTOC) ofxn as(7)Cm1,m2=ΔExmP1xm+m1P2xm+m2P3,where 0≤P1<α/2, 0≤P2<α/2, and 0≤P3<α/2. If xn are real, then the estimation of FLOTOC is given by(8)C^m1,m2=1M∑m=M1M2xmP1xm+m1P2xm+m2P3signxmxm+m1xm+m2,and if xn are complex, then(9)C^m1,m2=1M∑m=M1M2xmP1−1xm+m1P2−1xm+m2P3−1x∗mx∗m+m1x∗m+m2,where M1=max0,−m1,−m2 and M2=minM−1,M−1−m1,M−1−m2.From equations (7)–(9), The ith segment of FLOTOC of the signal xn can be given by(10)C^im1,m2=1M∑m=M1M2ximP1xim+m1P2xim+m2P3,where i=1,2,…,L, averaging fractional low-order third-order cumulant of the L segments; we can get(11)C^m1,m2=1L∑i=1LC^im1,m2.Taking the windowed two-dimensional discrete Fourier transform of equation (11), fractional low-order indirect bispectrum estimation of the signal xn can be given by(12)FLOIDBω1,ω2=∑m1=−KK∑m2=−KKC^m1,m2wm1,m2e−jm1ω1+m2ω2,where wm1,m2 is a two-dimensional window function and K<M−1. The 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.


### 3.3. Application Review

In this simulation, the test signalyn is defined as(13)yn=Acos2πf1n+φ1+cos2πf2n+φ2+cos2πf3n+φ3+vn,=Axn+vn,where xn is three cosinoidal signals and vn is additive Gaussian noise or SαS distribution noise. f1=0.1, f2=0.15, f3=0.2, φ1=π/3, φ2=π/4, φ3=π/5, n=0,1,2,…,N−1, and N=512. When vn is additive Gaussian noise, signal to noise ratio (SNR) can be used. But vn is SαS distribution noise, SNR is inapplicable, and generalized signal to noise ratio (GSNR) is written as(14)GSNR=10log10Exn2γα=10log101Nγα∑n=0N−1xn2,where γ is the dispersion coefficient of SαS distribution noise. According to the given GSNR, the amplitude of the signal xn is written as(15)A=10GSNR/101/N∑n=0N−1xn2γα1/2.LettingSNR=−5dB and GSNR=20dB, 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 xn under Gaussian distribution noise and SαS distribution noise; the simulation results are shown in Figures 3–6.Figure 3
The direct bispectrum estimations of the signalx(n) under Gaussian noise environment. (a) The direct method bispectrum estimation of the signal x(n). (b) The direct method bispectrum three-dimensional graph estimation of the signal x(n). (c) The FLO direct method bispectrum estimation of the signal x(n). (d) The FLO direct method bispectrum three-dimensional graph estimation of the signal x(n)).
(a)(b)(c)(d)Figure 4
The direct bispectrum estimations of the signalx(n) under SαS noise environment. (a) The direct method bispectrum estimation of the signal x(n). (b) The direct method bispectrum three-dimensional graph estimation of the signal x(n). (c) The FLO direct method bispectrum estimation of the signal x(n). (d) The FLO direct method bispectrum three-dimensional graph estimation of the signal x(n).
(a)(b)(c)(d)Figure 5
The indirect bispectrum estimations of the signalx(n) under Gaussian noise environment. (a). The indirect method bispectrum estimation of the signal x(n). (b) The indirect method bispectrum three-dimensional graph estimation of the signal x(n). (c) The FLO indirect method bispectrum estimation of the signal x(n). (d) The FLO indirect method bispectrum three-dimensional graph estimation of the signal x(n)).
(a)(b)(c)(d)Figure 6
The indirect bispectrum estimations of the signalx(n) under SαS noise environment. (a) The indirect method bispectrum estimation of the signal x(n). (b) The indirect method bispectrum three-dimensional graph estimation of the signal x(n). (c) The FLO indirect method bispectrum estimation of the signal x(n). (d) The FLO indirect method bispectrum three-dimensional graph estimation of the signal x(n)).
(a)(b)(c)(d)


### 3.4. Remarks

The direct bispectrum estimations of the signalx(n) under Gaussian noise environment (SNR=−5dB) are shown in Figure 3. Figures 3(a) and 3(b) are the traditional direct bispectral estimation and its three-dimensional graph estimation, respectively. Fractional lower-order direct bispectral estimation and its three-dimensional graph are given in Figures 3(c) and 3(d), respectively. The result shows that both methods can estimate the bispectrum of the signal x(n) well. Figure 4 shows the direct bispectral estimations of the signal x(n) under SαS distribution noise environment (GSNR=20dB; α=1.3). We can know that the traditional direct bispectral estimation fails in Figures 4(a) and 4(b), but the proposed fractional lower-order direct bispectral estimation method in Figures 4(c) and 4(d) shows good performance. Figures 5 and 6 are the traditional indirect bispectrum and fractional lower-order indirect bispectrum of the signal x(n) under Gaussian noise environment (SNR=−5dB) and SαS distribution noise environment (GSNR=20dB; α=1.3), respectively. The results show that both methods have better performance under Gaussian noise environment, but the conventional indirect bispectrum method degenerates under SαS noise environment, and fractional lower-order indirect bispectrum method can better estimate out the bispectrum of the signal x(n). Hence, the fractional lower-order direct and indirect bispectrum methods are robust.The 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. The 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:(16)FLOBω1,ω2=FLOBω1+2π,ω2+2π.Fractional lower-order bispectral symmetry properties are as follows:(17)FLOBω1,ω2=FLOBω2,ω1=FLOB−ω1−ω2,ω1=FLOB−ω1−ω2,ω2=FLOBω1,−ω1−ω2=FLOBω2,−ω1−ω2=FLOB∗−ω1,−ω2=FLOB∗−ω2,−ω1=FLOB∗−ω1,ω1+ω2=FLOB∗−ω2,ω1+ω2=FLOB∗ω1+ω2,−ω1=FLOB∗ω1+ω2,−ω2.


## 3.1. Fractional Lower-Order Direct Bispectrum Method

xn,n=0,1,…,N−1 is N samples of the observation data; its discrete fractional lower-order Fourier transform is defined as (1)X˜ω=∑n=0N−1xPne−jωn,where P is a given real constant and P<α≤2. P denotes P order moment of xt, when xt is a real signal, xPt=xtP⋅signxt, signxt=1,xt>0,0,xt=0,−1,xt<0,, and when xt is a complex signal, xPt=xtP−1⋅x∗t, where ∗ denotes conjugate operation.Dividing the dataxn into L segments, each segment is M points, and two adjacent segments overlap M/2 points; then L=2N−M/2/M and xin,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(2)X˜iω=∑n=1MxiPne−jωn.The corresponding discrete fractional low-order Fourier transform coefficient is defined as(3)X˜iφ=1M∑n=1MxiPne−j2πnφ/M,φ=0,1,…,M2.Let the sampling frequency of the signalxn be fs and the total frequency sampling points be I; then K=M/I−1/2. Combining equation (3), fractional low-order triple correlation of xn can be expressed as(4)B^iφ1,φ2=I2fs2∑k1=−KK∑k2=−KKX˜i′φ1+k1X˜i″φ2+k2X˜i″′−φ1−φ2−k1−k2,X˜i′φ1+k1=1M∑n=1MxiP1ne−j2πnφ1+k1/MX˜i″φ2+k2=1M∑n=1MxiP2ne−j2πnφ2+k2/M,X˜i″′−φ1−φ2−k1−k2=1M∑n=1MxiP3ne−j2πn−φ1−φ2−k1−k2/M,where P1+P2+P3<α, and we let P1=P2=P3<α/3, 0≤φ2≤φ1, and φ2+φ1≤fs/2. The corresponding fractional low-order direct bispectrum of the ith segment is given by(5)FLODBiω1,ω2=I2fs2∑k1=−KK∑k2=−KKX˜i′ω1′+k1X˜i″ω2′+k2X˜i″′−ω1′−ω2′−k1−k2,where ω1′=Iω1/2πfs and ω2′=Iω2/2πfs. Averaging the bispectrum of those L segments, fractional low-order direct bispectrum of xn can be gotten.(6)FLODBω1,ω2=1L∑i=1LFLODBiω1,ω2.Fractional low-order direct bispectral estimation process is to first calculate the discrete fractional low-order Fourier transform coefficient and then compute its fractional low-order triple correlation and finally average fractional low-order triple correlation of all segments.


## 3.2. Fractional Lower-Order Indirect Bispectrum Method

By using the fractional lower-order moment, we define discrete fractional lower-order three-order cumulants (FLOTOC) ofxn as(7)Cm1,m2=ΔExmP1xm+m1P2xm+m2P3,where 0≤P1<α/2, 0≤P2<α/2, and 0≤P3<α/2. If xn are real, then the estimation of FLOTOC is given by(8)C^m1,m2=1M∑m=M1M2xmP1xm+m1P2xm+m2P3signxmxm+m1xm+m2,and if xn are complex, then(9)C^m1,m2=1M∑m=M1M2xmP1−1xm+m1P2−1xm+m2P3−1x∗mx∗m+m1x∗m+m2,where M1=max0,−m1,−m2 and M2=minM−1,M−1−m1,M−1−m2.From equations (7)–(9), The ith segment of FLOTOC of the signal xn can be given by(10)C^im1,m2=1M∑m=M1M2ximP1xim+m1P2xim+m2P3,where i=1,2,…,L, averaging fractional low-order third-order cumulant of the L segments; we can get(11)C^m1,m2=1L∑i=1LC^im1,m2.Taking the windowed two-dimensional discrete Fourier transform of equation (11), fractional low-order indirect bispectrum estimation of the signal xn can be given by(12)FLOIDBω1,ω2=∑m1=−KK∑m2=−KKC^m1,m2wm1,m2e−jm1ω1+m2ω2,where wm1,m2 is a two-dimensional window function and K<M−1. The 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.


## 3.3. Application Review

In this simulation, the test signalyn is defined as(13)yn=Acos2πf1n+φ1+cos2πf2n+φ2+cos2πf3n+φ3+vn,=Axn+vn,where xn is three cosinoidal signals and vn is additive Gaussian noise or SαS distribution noise. f1=0.1, f2=0.15, f3=0.2, φ1=π/3, φ2=π/4, φ3=π/5, n=0,1,2,…,N−1, and N=512. When vn is additive Gaussian noise, signal to noise ratio (SNR) can be used. But vn is SαS distribution noise, SNR is inapplicable, and generalized signal to noise ratio (GSNR) is written as(14)GSNR=10log10Exn2γα=10log101Nγα∑n=0N−1xn2,where γ is the dispersion coefficient of SαS distribution noise. According to the given GSNR, the amplitude of the signal xn is written as(15)A=10GSNR/101/N∑n=0N−1xn2γα1/2.LettingSNR=−5dB and GSNR=20dB, 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 xn under Gaussian distribution noise and SαS distribution noise; the simulation results are shown in Figures 3–6.Figure 3
The direct bispectrum estimations of the signalx(n) under Gaussian noise environment. (a) The direct method bispectrum estimation of the signal x(n). (b) The direct method bispectrum three-dimensional graph estimation of the signal x(n). (c) The FLO direct method bispectrum estimation of the signal x(n). (d) The FLO direct method bispectrum three-dimensional graph estimation of the signal x(n)).
(a)(b)(c)(d)Figure 4
The direct bispectrum estimations of the signalx(n) under SαS noise environment. (a) The direct method bispectrum estimation of the signal x(n). (b) The direct method bispectrum three-dimensional graph estimation of the signal x(n). (c) The FLO direct method bispectrum estimation of the signal x(n). (d) The FLO direct method bispectrum three-dimensional graph estimation of the signal x(n).
(a)(b)(c)(d)Figure 5
The indirect bispectrum estimations of the signalx(n) under Gaussian noise environment. (a). The indirect method bispectrum estimation of the signal x(n). (b) The indirect method bispectrum three-dimensional graph estimation of the signal x(n). (c) The FLO indirect method bispectrum estimation of the signal x(n). (d) The FLO indirect method bispectrum three-dimensional graph estimation of the signal x(n)).
(a)(b)(c)(d)Figure 6
The indirect bispectrum estimations of the signalx(n) under SαS noise environment. (a) The indirect method bispectrum estimation of the signal x(n). (b) The indirect method bispectrum three-dimensional graph estimation of the signal x(n). (c) The FLO indirect method bispectrum estimation of the signal x(n). (d) The FLO indirect method bispectrum three-dimensional graph estimation of the signal x(n)).
(a)(b)(c)(d)


## 3.4. Remarks

The direct bispectrum estimations of the signalx(n) under Gaussian noise environment (SNR=−5dB) are shown in Figure 3. Figures 3(a) and 3(b) are the traditional direct bispectral estimation and its three-dimensional graph estimation, respectively. Fractional lower-order direct bispectral estimation and its three-dimensional graph are given in Figures 3(c) and 3(d), respectively. The result shows that both methods can estimate the bispectrum of the signal x(n) well. Figure 4 shows the direct bispectral estimations of the signal x(n) under SαS distribution noise environment (GSNR=20dB; α=1.3). We can know that the traditional direct bispectral estimation fails in Figures 4(a) and 4(b), but the proposed fractional lower-order direct bispectral estimation method in Figures 4(c) and 4(d) shows good performance. Figures 5 and 6 are the traditional indirect bispectrum and fractional lower-order indirect bispectrum of the signal x(n) under Gaussian noise environment (SNR=−5dB) and SαS distribution noise environment (GSNR=20dB; α=1.3), respectively. The results show that both methods have better performance under Gaussian noise environment, but the conventional indirect bispectrum method degenerates under SαS noise environment, and fractional lower-order indirect bispectrum method can better estimate out the bispectrum of the signal x(n). Hence, the fractional lower-order direct and indirect bispectrum methods are robust.The 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. The 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:(16)FLOBω1,ω2=FLOBω1+2π,ω2+2π.Fractional lower-order bispectral symmetry properties are as follows:(17)FLOBω1,ω2=FLOBω2,ω1=FLOB−ω1−ω2,ω1=FLOB−ω1−ω2,ω2=FLOBω1,−ω1−ω2=FLOBω2,−ω1−ω2=FLOB∗−ω1,−ω2=FLOB∗−ω2,−ω1=FLOB∗−ω1,ω1+ω2=FLOB∗−ω2,ω1+ω2=FLOB∗ω1+ω2,−ω1=FLOB∗ω1+ω2,−ω2.


## 4. Fractional Lower-Order Parametric Bispectrum Methods

### 4.1. Fractional Lower-Order AR Model Bispectrum Method

#### 4.1.1. Principle

The conventional realp order AR model process can be expressed as(18)Xn=∑i=1paiXn−i+Wn.The traditional AR model bispectrum of the signalXn in (18) is given by(19)Bω1,ω2=β⌢Hω1Hω2H∗ω1+ω2,where β⌢ is the estimate of the third moment of the driving noise, Hω is the system transfer function, and(20)Hejω=11−∑i=1paie−jωi.A fractional lower-order AR modelSαS process xn may be written as(21)xn=∑m=1pamxn−m+un,where p is order of the AR model and ami=1,2,…,p are its parameters, which are real numbers. un is an independent identically distributed (i.i.d) SαS random process, α is its characteristic index, and γu is its dispersion coefficient. xn can be expressed by a finite impulse response (FIR) system [27].(22)xn=∑m=0qhmun−m,where hm is system impulse response coefficient and un−m is i.i.d SαS random process. Taking z transformation for equations (21) and (22) and obtaining the system transfer function,(23)Hz=∑m=0qhmz−m=11−∑m=1pamz−m.Takingα-order moment with equation (23), we have(24)H¯z=Hzα=∑m=0qhmz−mα=11−∑m=1pamz−mα,z=ejω, and z=1 on the unit circle; then(25)H¯ejω=11−∑m=1pame−jωmα.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(26)FLOARBω1,ω2=γuH¯ω1H¯ω2H¯∗ω1+ω2,and FLOARB on the unit circle is written as(27)FLOARBejω1,ejω2=γuH¯ejω1H¯ejω2H¯∗ejω1+ejω2,where γu is the dispersion coefficient of the driving SαS random process un and ∗ is conjugate operation.According to the definition of the fractional lower-order three-order cumulants in (8), we define A=a1,a2,…,apT, C˜=C^1,1,C^2,2,…,C^p,pT, and(28)Φ=ΔC^0,0C^−1,−1…C^1−p,1−pC^2,2C^0,0…C^2−p,2−p⋮⋮⋱⋮C^p−1,p−1C^p−2,p−2…C^0,0,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(29)ΦA=C˜.Then,(30)C^0,0C^−1,−1…C^1−p,1−pC^2,2C^0,0…C^2−p,2−p⋮⋮⋱⋮C^p−1,p−1C^p−2,p−2…C^0,0a1a2⋮ap=C^1,1C^2,2⋮C^p,p.Equation (30) is a fractional low-order bispectrum Yule-Walker equation. The AR model parameters a^ii=1,2,…,p can be gotten by solving equation (30). When a^ii=1,2,…,p are substituted into equation (27), we can get fractional low-order AR model bispectrum estimation.In this paper, we apply the final prediction error (FPE) criteria to determine the orderp 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. The calculation formula can be written as(31)FPEp=σa2N+pN−p,where σa2 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 signalxn with equations (7)–(9) and constructing fractional low-order three-order cumulants matrix in (28).Step 2: solving fractional low-order bispectrum Yule-Walker equation in equation (30) and getting the coefficients a^ii=1,2,…,p.Step 3: computing fractional low-order AR model bispectrum ofxn by substituting a^ii=1,2,…,p into equation (19).


#### 4.1.2. Application Review

In this simulation, the experimental signal isyn in equation (13). The performances of the traditional AR model bispectrum method and the improved FLOAR model bispectrum method are compared under Gaussian distribution noise (SNR=−5dB) and SαS distribution noise (α=1.2; GSNR=20dB), and the simulation results are shown in Figures 7–10.Figure 7
The AR model and FLOAR bispectrum estimations of the signalx(n) under Gaussian noise environment. (a) The AR model bispectrum estimation of the signal x(n). (b) The AR model bispectrum three-dimensional graph estimation of the signal x(n). (c) The FLOAR model bispectrum estimation of the signal x(n). (d) The FLOAR model bispectrum three-dimensional graph estimation of the signal x(n).
(a)(b)(c)(d)Figure 8
The AR model and FLOAR bispectrum estimations of the signalx(n) under SαS noise environment. (a) The AR model bispectrum estimation of the signal x(n). (b) The AR model bispectrum three-dimensional graph estimation of the signal x(n). (c) The FLOAR model bispectrum estimation of the signal x(n). (d) The FLOAR model bispectrum three-dimensional graph estimation of the signal x(n).
(a)(b)(c)(d)Figure 9
The ARMA model and FLOARMA bispectrum estimations of the signalx(n) under Gaussian noise environment. (a) The ARMA(2, 6) model bispectrum estimation of the signal x(n). (b) The ARMA(2, 6) model bispectrum three-dimensional graph estimation of the signal x(n). (c) The FLOARMA(2, 6) model bispectrum estimation of the signal x(n). (d) The FLOARMA(2, 6) model bispectrum three-dimensional graph estimation of the signal x(n).
(a)(b)(c)(d)Figure 10
The ARMA model and FLOARMA bispectrum estimations of the signalx(n) under SαS noise environment. (a) The ARMA(2, 6) model bispectrum estimation of the signal x(n). (b) The ARMA(2, 6) model bispectrum three-dimensional graph estimation of the signal x(n). (c) The FLOARMA(2, 6) model bispectrum estimation of the signal x(n). (d) The FLOARMA(2, 6) model bispectrum three-dimensional graph estimation of the signal x(n).
(a)(b)(c)(d)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.Figure 11
The fractional lower-order bispectrum estimations of the outer race normal signal. (a) FLODB estimation. (b) FLODB three-dimensional graph. (c) FLOIDB estimation. (d) FLOIDB three-dimensional graph. (e) FLOARB estimation. (f) FLOARB three-dimensional graph. (g) FLOARMAB estimation. (h) FLOARMAB three-dimensional graph.
(a)(b)(c)(d)(e)(f)(g)(h)


#### 4.1.3. Remarks

Figure7 is the AR model and FLOAR bispectrum estimations of the signal x(n) under Gaussian noise environment (SNR=−5dB). The traditional AR model bispectral estimation and its three-dimensional graph estimation are shown in Figures 7(a) and 7(b), respectively. Fractional lower-order AR model bispectral estimation and its three-dimensional graph are given in Figures 7(c) and 7(d), respectively. The result shows that both methods can estimate the bispectrum of the signal x(n) well under Gaussian noise environment.The AR and FLOAR model bispectrum estimations of the signalx(n) under SαS distribution noise (α=1.2; GSNR=20dB) are given in Figure 8. The 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.


### 4.2. Fractional Lower-Order ARMA Model Bispectrum Method

#### 4.2.1. Principle

A fractional lower-order autoregressive moving average (FLOARMA) model i.i.dSαS process xn can be given by(32)xn=−∑m=1paixn−m+∑m=0qbiun−m,where p and q are orders of the AR and MA model, respectively. ami=1,2,…,p and bmm=0,1,2,…,q are their parameters, which are real numbers. un is an independent identically distributed (i.i.d) SαS random process, α is its characteristic index, and its dispersion coefficient is γu.Simplifying equation (32), we have(33)∑m=0paixn−m=∑m=0qbiun−m,where a0=1 and b0=1. Taking z transform on both sides of equation (33), we obtain(34)∑m=0paiz−mXz=∑m=0qbiz−mUz,(35)1+∑m=1paiz−mXz=1+∑m=1qbiz−mUz,(36)Hz=XzUz=1+∑m=1qbiz−m1+∑m=1paiz−m=BzAz,(37)hm↔Hz,xn=∑m=0qhmun−m,where hm is system impulse response coefficient,Hz is system transfer function, and un−m is i.i.d SαS random process. The AR model parameters changes with ami=1,2,…,p; Az and Bz are the FIR and IIR filters, respectively.Takingα-order moment with equation (36),(38)H˜z=Hzα=∑m=0qhmz−mα=1+∑m=1pbmz−m1+∑m=1pamz−mα,z=ejω, and z=1 on the unit circle; then(39)H˜ejω=1+∑m=1pbme−jωm1+∑m=1pame−jωmα.According to the definition of the ARMA model bispectrum andα-order moment of Hz in (38) and (39), we define fractional lower-order ARMA model bispectrum (FLOARMAB) as(40)FLOARMABω1,ω2=γuH˜ω1H˜ω2H˜∗ω1+ω2,and FLOARB on the unit circle is written as(41)FLOARMABejω1,ejω2=γuH˜ejω1H˜ejω2H˜∗ejω1+ejω2,where γu is the dispersion coefficient of the driving SαS random process un and ∗ is conjugate operation.To obtain the coefficients of the fractional lower-order ARMA modelami=0,1,2,…,p and bmm=0,1,2,…,q, we should multiply xn−i on both sides of equation (32) taking fractional lower-order covariance to get(42)Dxxi=−∑k=1pakDxxi−k+∑k=0qbkDuxi−k,where(43)Dxxi=ΔExnP1xn−iP2,(44)Dxxi−k=ΔExnP1xn−i+kP2,(45)Duxi−k=ΔEunP1xn−i+kP2.It can be shown [24] that, for the fractional lower-order covariance of the signal and noise in (45),(46)Duxi=0,i>0,γuh−iα−1,i<0.Then, ifi>q, we have(47)Dxxi=∑k=1pakDxxi−k,i>q,(48)DxxqDxxq−1…Dxxq−p+1Dxxq+1Dxxq…Dxxq−p+2⋮⋮⋱⋮Dxxq+p−1Dxxq+p−2…Dxxqa1a2⋮ap=Dxxq+1Dxxq+2⋮Dxxq+p.Equations (47) and (48) are a generalized Yule-Walker equation, and we can obtain the fractional lower-order AR model (FLOARM) coefficients aki=0,1,2,…,p by solving them.If0≤m≤q, equation (42) changes as(49)Dxxi=∑k=0pakDxxi−k+∑k=0qbkDuxi−k=∑k=1pakDxxi−k+∑k=0qbkhi−kα−1,0≤i≤q.Lettingp=0, we have(50)Dxxm=γu∑k=0p−ibk+ibkα−1,0≤i≤q,γu∑k=0p−ibkbk−iα−1,−q≤i≤0.We can obtain the coefficientsbkk=1,2,…,q of the fractional lower-order MA model by solving the nonlinear equation (35).A finiteq′-order FLOMA model can be equivalent by an approximate infinite p′-order FLOAR modelp′≫q′; we have(51)ai+∑k=1q′bkai−k=ei,where em is error. Letting em=0, we obtain(52)∑k=1q′bkai−k=−ai.We should multiplyaPn−i on both sides of equation (52). We take fractional lower-order covariance to get(53)∑k=1q′bkR˜i−k=−R˜i,(54)R˜0R˜−1…R˜1−q′R˜1R˜0…R˜2−q′⋮⋮⋱⋮R˜q′−1R˜q′−2…R˜0b1b2⋮bp=−R˜1R˜2⋮R˜q′,where R˜i=EaiaPn−i and R˜i−k=Eai−kaPn−i+k. We can also obtain the coefficients bii=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 akk=0,1,2,…,p of the FLOAR model.Step 2: obtaining the coefficientsbkk=1,2,…,q of the FLOMA mode by solving the nonlinear equation (50) or the Yule-Walker equation in (54).Step 3: determining the orderp and q employing FPE criterion.Step 4: computing fractional low-order ARMA model bispectrum ofxn by substituting akk=0,1,2,…,p and bkk=1,2,…,q into equation (40) or equation (41).


#### 4.2.2. Application Review

In this simulation,yn in equation (14) is used as the experimental signal. The ARMA model and FLOARMA model bispectrum methods are compared to demonstrate their performance under Gaussian distribution noise (SNR=−5dB) and SαS distribution noise (α=1.2; GSNR=20dB), and the simulations are shown in Figures 9–12.Figure 12
The fractional lower-order bispectrum estimations of the outer race fault signal. (a) FLODB estimation. (b) FLODB three-dimensional graph. (c) FLOIDB estimation. (d) FLOIDB three-dimensional graph. (e) FLOARB estimation. (f) FLOARB three-dimensional graph. (g) FLOARMAB estimation. (h) FLOARMAB three-dimensional graph. (i) Fractional lower-order Choi-Williams distribution time-frequency representation method. (j) Wave shape and spectrum of the envelope signal.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)


#### 4.2.3. Remarks

The ARMA model and FLOARMA model bispectrum estimations of the signalx(n) are shown in Figure 9 under Gaussian noise environment (SNR=−5dB). 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. The results show that both methods can estimate out the bispectrum of the signal x(n) well under Gaussian noise.Figure10 shows the ARMA and FLOARMA model bispectrum estimations of the signal x(n) under SαS distribution noise (α=1.2; GSNR=20dB). The simulation results show that the existing ARMA bispectrum method in Figures 10(a) and 10(b) fails under SαS noise environment, but the improved FLOARMA model bispectrum method in Figures 10(c) and 10(d) has good toughness. Hence, the ARMA bispectrum method only works for Gaussian noise environment, but the FLOARMA bispectrum method can be applied in Gaussian and SαS noise environment, which is robust.


## 4.1. Fractional Lower-Order AR Model Bispectrum Method

### 4.1.1. Principle

The conventional realp order AR model process can be expressed as(18)Xn=∑i=1paiXn−i+Wn.The traditional AR model bispectrum of the signalXn in (18) is given by(19)Bω1,ω2=β⌢Hω1Hω2H∗ω1+ω2,where β⌢ is the estimate of the third moment of the driving noise, Hω is the system transfer function, and(20)Hejω=11−∑i=1paie−jωi.A fractional lower-order AR modelSαS process xn may be written as(21)xn=∑m=1pamxn−m+un,where p is order of the AR model and ami=1,2,…,p are its parameters, which are real numbers. un is an independent identically distributed (i.i.d) SαS random process, α is its characteristic index, and γu is its dispersion coefficient. xn can be expressed by a finite impulse response (FIR) system [27].(22)xn=∑m=0qhmun−m,where hm is system impulse response coefficient and un−m is i.i.d SαS random process. Taking z transformation for equations (21) and (22) and obtaining the system transfer function,(23)Hz=∑m=0qhmz−m=11−∑m=1pamz−m.Takingα-order moment with equation (23), we have(24)H¯z=Hzα=∑m=0qhmz−mα=11−∑m=1pamz−mα,z=ejω, and z=1 on the unit circle; then(25)H¯ejω=11−∑m=1pame−jωmα.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(26)FLOARBω1,ω2=γuH¯ω1H¯ω2H¯∗ω1+ω2,and FLOARB on the unit circle is written as(27)FLOARBejω1,ejω2=γuH¯ejω1H¯ejω2H¯∗ejω1+ejω2,where γu is the dispersion coefficient of the driving SαS random process un and ∗ is conjugate operation.According to the definition of the fractional lower-order three-order cumulants in (8), we define A=a1,a2,…,apT, C˜=C^1,1,C^2,2,…,C^p,pT, and(28)Φ=ΔC^0,0C^−1,−1…C^1−p,1−pC^2,2C^0,0…C^2−p,2−p⋮⋮⋱⋮C^p−1,p−1C^p−2,p−2…C^0,0,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(29)ΦA=C˜.Then,(30)C^0,0C^−1,−1…C^1−p,1−pC^2,2C^0,0…C^2−p,2−p⋮⋮⋱⋮C^p−1,p−1C^p−2,p−2…C^0,0a1a2⋮ap=C^1,1C^2,2⋮C^p,p.Equation (30) is a fractional low-order bispectrum Yule-Walker equation. The AR model parameters a^ii=1,2,…,p can be gotten by solving equation (30). When a^ii=1,2,…,p are substituted into equation (27), we can get fractional low-order AR model bispectrum estimation.In this paper, we apply the final prediction error (FPE) criteria to determine the orderp 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. The calculation formula can be written as(31)FPEp=σa2N+pN−p,where σa2 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 signalxn with equations (7)–(9) and constructing fractional low-order three-order cumulants matrix in (28).Step 2: solving fractional low-order bispectrum Yule-Walker equation in equation (30) and getting the coefficients a^ii=1,2,…,p.Step 3: computing fractional low-order AR model bispectrum ofxn by substituting a^ii=1,2,…,p into equation (19).


### 4.1.2. Application Review

In this simulation, the experimental signal isyn in equation (13). The performances of the traditional AR model bispectrum method and the improved FLOAR model bispectrum method are compared under Gaussian distribution noise (SNR=−5dB) and SαS distribution noise (α=1.2; GSNR=20dB), and the simulation results are shown in Figures 7–10.Figure 7
The AR model and FLOAR bispectrum estimations of the signalx(n) under Gaussian noise environment. (a) The AR model bispectrum estimation of the signal x(n). (b) The AR model bispectrum three-dimensional graph estimation of the signal x(n). (c) The FLOAR model bispectrum estimation of the signal x(n). (d) The FLOAR model bispectrum three-dimensional graph estimation of the signal x(n).
(a)(b)(c)(d)Figure 8
The AR model and FLOAR bispectrum estimations of the signalx(n) under SαS noise environment. (a) The AR model bispectrum estimation of the signal x(n). (b) The AR model bispectrum three-dimensional graph estimation of the signal x(n). (c) The FLOAR model bispectrum estimation of the signal x(n). (d) The FLOAR model bispectrum three-dimensional graph estimation of the signal x(n).
(a)(b)(c)(d)Figure 9
The ARMA model and FLOARMA bispectrum estimations of the signalx(n) under Gaussian noise environment. (a) The ARMA(2, 6) model bispectrum estimation of the signal x(n). (b) The ARMA(2, 6) model bispectrum three-dimensional graph estimation of the signal x(n). (c) The FLOARMA(2, 6) model bispectrum estimation of the signal x(n). (d) The FLOARMA(2, 6) model bispectrum three-dimensional graph estimation of the signal x(n).
(a)(b)(c)(d)Figure 10
The ARMA model and FLOARMA bispectrum estimations of the signalx(n) under SαS noise environment. (a) The ARMA(2, 6) model bispectrum estimation of the signal x(n). (b) The ARMA(2, 6) model bispectrum three-dimensional graph estimation of the signal x(n). (c) The FLOARMA(2, 6) model bispectrum estimation of the signal x(n). (d) The FLOARMA(2, 6) model bispectrum three-dimensional graph estimation of the signal x(n).
(a)(b)(c)(d)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.Figure 11
The fractional lower-order bispectrum estimations of the outer race normal signal. (a) FLODB estimation. (b) FLODB three-dimensional graph. (c) FLOIDB estimation. (d) FLOIDB three-dimensional graph. (e) FLOARB estimation. (f) FLOARB three-dimensional graph. (g) FLOARMAB estimation. (h) FLOARMAB three-dimensional graph.
(a)(b)(c)(d)(e)(f)(g)(h)


### 4.1.3. Remarks

Figure7 is the AR model and FLOAR bispectrum estimations of the signal x(n) under Gaussian noise environment (SNR=−5dB). The traditional AR model bispectral estimation and its three-dimensional graph estimation are shown in Figures 7(a) and 7(b), respectively. Fractional lower-order AR model bispectral estimation and its three-dimensional graph are given in Figures 7(c) and 7(d), respectively. The result shows that both methods can estimate the bispectrum of the signal x(n) well under Gaussian noise environment.The AR and FLOAR model bispectrum estimations of the signalx(n) under SαS distribution noise (α=1.2; GSNR=20dB) are given in Figure 8. The 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.


## 4.1.1. Principle

The conventional realp order AR model process can be expressed as(18)Xn=∑i=1paiXn−i+Wn.The traditional AR model bispectrum of the signalXn in (18) is given by(19)Bω1,ω2=β⌢Hω1Hω2H∗ω1+ω2,where β⌢ is the estimate of the third moment of the driving noise, Hω is the system transfer function, and(20)Hejω=11−∑i=1paie−jωi.A fractional lower-order AR modelSαS process xn may be written as(21)xn=∑m=1pamxn−m+un,where p is order of the AR model and ami=1,2,…,p are its parameters, which are real numbers. un is an independent identically distributed (i.i.d) SαS random process, α is its characteristic index, and γu is its dispersion coefficient. xn can be expressed by a finite impulse response (FIR) system [27].(22)xn=∑m=0qhmun−m,where hm is system impulse response coefficient and un−m is i.i.d SαS random process. Taking z transformation for equations (21) and (22) and obtaining the system transfer function,(23)Hz=∑m=0qhmz−m=11−∑m=1pamz−m.Takingα-order moment with equation (23), we have(24)H¯z=Hzα=∑m=0qhmz−mα=11−∑m=1pamz−mα,z=ejω, and z=1 on the unit circle; then(25)H¯ejω=11−∑m=1pame−jωmα.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(26)FLOARBω1,ω2=γuH¯ω1H¯ω2H¯∗ω1+ω2,and FLOARB on the unit circle is written as(27)FLOARBejω1,ejω2=γuH¯ejω1H¯ejω2H¯∗ejω1+ejω2,where γu is the dispersion coefficient of the driving SαS random process un and ∗ is conjugate operation.According to the definition of the fractional lower-order three-order cumulants in (8), we define A=a1,a2,…,apT, C˜=C^1,1,C^2,2,…,C^p,pT, and(28)Φ=ΔC^0,0C^−1,−1…C^1−p,1−pC^2,2C^0,0…C^2−p,2−p⋮⋮⋱⋮C^p−1,p−1C^p−2,p−2…C^0,0,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(29)ΦA=C˜.Then,(30)C^0,0C^−1,−1…C^1−p,1−pC^2,2C^0,0…C^2−p,2−p⋮⋮⋱⋮C^p−1,p−1C^p−2,p−2…C^0,0a1a2⋮ap=C^1,1C^2,2⋮C^p,p.Equation (30) is a fractional low-order bispectrum Yule-Walker equation. The AR model parameters a^ii=1,2,…,p can be gotten by solving equation (30). When a^ii=1,2,…,p are substituted into equation (27), we can get fractional low-order AR model bispectrum estimation.In this paper, we apply the final prediction error (FPE) criteria to determine the orderp 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. The calculation formula can be written as(31)FPEp=σa2N+pN−p,where σa2 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 signalxn with equations (7)–(9) and constructing fractional low-order three-order cumulants matrix in (28).Step 2: solving fractional low-order bispectrum Yule-Walker equation in equation (30) and getting the coefficients a^ii=1,2,…,p.Step 3: computing fractional low-order AR model bispectrum ofxn by substituting a^ii=1,2,…,p into equation (19).


## 4.1.2. Application Review

In this simulation, the experimental signal isyn in equation (13). The performances of the traditional AR model bispectrum method and the improved FLOAR model bispectrum method are compared under Gaussian distribution noise (SNR=−5dB) and SαS distribution noise (α=1.2; GSNR=20dB), and the simulation results are shown in Figures 7–10.Figure 7
The AR model and FLOAR bispectrum estimations of the signalx(n) under Gaussian noise environment. (a) The AR model bispectrum estimation of the signal x(n). (b) The AR model bispectrum three-dimensional graph estimation of the signal x(n). (c) The FLOAR model bispectrum estimation of the signal x(n). (d) The FLOAR model bispectrum three-dimensional graph estimation of the signal x(n).
(a)(b)(c)(d)Figure 8
The AR model and FLOAR bispectrum estimations of the signalx(n) under SαS noise environment. (a) The AR model bispectrum estimation of the signal x(n). (b) The AR model bispectrum three-dimensional graph estimation of the signal x(n). (c) The FLOAR model bispectrum estimation of the signal x(n). (d) The FLOAR model bispectrum three-dimensional graph estimation of the signal x(n).
(a)(b)(c)(d)Figure 9
The ARMA model and FLOARMA bispectrum estimations of the signalx(n) under Gaussian noise environment. (a) The ARMA(2, 6) model bispectrum estimation of the signal x(n). (b) The ARMA(2, 6) model bispectrum three-dimensional graph estimation of the signal x(n). (c) The FLOARMA(2, 6) model bispectrum estimation of the signal x(n). (d) The FLOARMA(2, 6) model bispectrum three-dimensional graph estimation of the signal x(n).
(a)(b)(c)(d)Figure 10
The ARMA model and FLOARMA bispectrum estimations of the signalx(n) under SαS noise environment. (a) The ARMA(2, 6) model bispectrum estimation of the signal x(n). (b) The ARMA(2, 6) model bispectrum three-dimensional graph estimation of the signal x(n). (c) The FLOARMA(2, 6) model bispectrum estimation of the signal x(n). (d) The FLOARMA(2, 6) model bispectrum three-dimensional graph estimation of the signal x(n).
(a)(b)(c)(d)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.Figure 11
The fractional lower-order bispectrum estimations of the outer race normal signal. (a) FLODB estimation. (b) FLODB three-dimensional graph. (c) FLOIDB estimation. (d) FLOIDB three-dimensional graph. (e) FLOARB estimation. (f) FLOARB three-dimensional graph. (g) FLOARMAB estimation. (h) FLOARMAB three-dimensional graph.
(a)(b)(c)(d)(e)(f)(g)(h)


## 4.1.3. Remarks

Figure7 is the AR model and FLOAR bispectrum estimations of the signal x(n) under Gaussian noise environment (SNR=−5dB). The traditional AR model bispectral estimation and its three-dimensional graph estimation are shown in Figures 7(a) and 7(b), respectively. Fractional lower-order AR model bispectral estimation and its three-dimensional graph are given in Figures 7(c) and 7(d), respectively. The result shows that both methods can estimate the bispectrum of the signal x(n) well under Gaussian noise environment.The AR and FLOAR model bispectrum estimations of the signalx(n) under SαS distribution noise (α=1.2; GSNR=20dB) are given in Figure 8. The 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.


## 4.2. Fractional Lower-Order ARMA Model Bispectrum Method

### 4.2.1. Principle

A fractional lower-order autoregressive moving average (FLOARMA) model i.i.dSαS process xn can be given by(32)xn=−∑m=1paixn−m+∑m=0qbiun−m,where p and q are orders of the AR and MA model, respectively. ami=1,2,…,p and bmm=0,1,2,…,q are their parameters, which are real numbers. un is an independent identically distributed (i.i.d) SαS random process, α is its characteristic index, and its dispersion coefficient is γu.Simplifying equation (32), we have(33)∑m=0paixn−m=∑m=0qbiun−m,where a0=1 and b0=1. Taking z transform on both sides of equation (33), we obtain(34)∑m=0paiz−mXz=∑m=0qbiz−mUz,(35)1+∑m=1paiz−mXz=1+∑m=1qbiz−mUz,(36)Hz=XzUz=1+∑m=1qbiz−m1+∑m=1paiz−m=BzAz,(37)hm↔Hz,xn=∑m=0qhmun−m,where hm is system impulse response coefficient,Hz is system transfer function, and un−m is i.i.d SαS random process. The AR model parameters changes with ami=1,2,…,p; Az and Bz are the FIR and IIR filters, respectively.Takingα-order moment with equation (36),(38)H˜z=Hzα=∑m=0qhmz−mα=1+∑m=1pbmz−m1+∑m=1pamz−mα,z=ejω, and z=1 on the unit circle; then(39)H˜ejω=1+∑m=1pbme−jωm1+∑m=1pame−jωmα.According to the definition of the ARMA model bispectrum andα-order moment of Hz in (38) and (39), we define fractional lower-order ARMA model bispectrum (FLOARMAB) as(40)FLOARMABω1,ω2=γuH˜ω1H˜ω2H˜∗ω1+ω2,and FLOARB on the unit circle is written as(41)FLOARMABejω1,ejω2=γuH˜ejω1H˜ejω2H˜∗ejω1+ejω2,where γu is the dispersion coefficient of the driving SαS random process un and ∗ is conjugate operation.To obtain the coefficients of the fractional lower-order ARMA modelami=0,1,2,…,p and bmm=0,1,2,…,q, we should multiply xn−i on both sides of equation (32) taking fractional lower-order covariance to get(42)Dxxi=−∑k=1pakDxxi−k+∑k=0qbkDuxi−k,where(43)Dxxi=ΔExnP1xn−iP2,(44)Dxxi−k=ΔExnP1xn−i+kP2,(45)Duxi−k=ΔEunP1xn−i+kP2.It can be shown [24] that, for the fractional lower-order covariance of the signal and noise in (45),(46)Duxi=0,i>0,γuh−iα−1,i<0.Then, ifi>q, we have(47)Dxxi=∑k=1pakDxxi−k,i>q,(48)DxxqDxxq−1…Dxxq−p+1Dxxq+1Dxxq…Dxxq−p+2⋮⋮⋱⋮Dxxq+p−1Dxxq+p−2…Dxxqa1a2⋮ap=Dxxq+1Dxxq+2⋮Dxxq+p.Equations (47) and (48) are a generalized Yule-Walker equation, and we can obtain the fractional lower-order AR model (FLOARM) coefficients aki=0,1,2,…,p by solving them.If0≤m≤q, equation (42) changes as(49)Dxxi=∑k=0pakDxxi−k+∑k=0qbkDuxi−k=∑k=1pakDxxi−k+∑k=0qbkhi−kα−1,0≤i≤q.Lettingp=0, we have(50)Dxxm=γu∑k=0p−ibk+ibkα−1,0≤i≤q,γu∑k=0p−ibkbk−iα−1,−q≤i≤0.We can obtain the coefficientsbkk=1,2,…,q of the fractional lower-order MA model by solving the nonlinear equation (35).A finiteq′-order FLOMA model can be equivalent by an approximate infinite p′-order FLOAR modelp′≫q′; we have(51)ai+∑k=1q′bkai−k=ei,where em is error. Letting em=0, we obtain(52)∑k=1q′bkai−k=−ai.We should multiplyaPn−i on both sides of equation (52). We take fractional lower-order covariance to get(53)∑k=1q′bkR˜i−k=−R˜i,(54)R˜0R˜−1…R˜1−q′R˜1R˜0…R˜2−q′⋮⋮⋱⋮R˜q′−1R˜q′−2…R˜0b1b2⋮bp=−R˜1R˜2⋮R˜q′,where R˜i=EaiaPn−i and R˜i−k=Eai−kaPn−i+k. We can also obtain the coefficients bii=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 akk=0,1,2,…,p of the FLOAR model.Step 2: obtaining the coefficientsbkk=1,2,…,q of the FLOMA mode by solving the nonlinear equation (50) or the Yule-Walker equation in (54).Step 3: determining the orderp and q employing FPE criterion.Step 4: computing fractional low-order ARMA model bispectrum ofxn by substituting akk=0,1,2,…,p and bkk=1,2,…,q into equation (40) or equation (41).


### 4.2.2. Application Review

In this simulation,yn in equation (14) is used as the experimental signal. The ARMA model and FLOARMA model bispectrum methods are compared to demonstrate their performance under Gaussian distribution noise (SNR=−5dB) and SαS distribution noise (α=1.2; GSNR=20dB), and the simulations are shown in Figures 9–12.Figure 12
The fractional lower-order bispectrum estimations of the outer race fault signal. (a) FLODB estimation. (b) FLODB three-dimensional graph. (c) FLOIDB estimation. (d) FLOIDB three-dimensional graph. (e) FLOARB estimation. (f) FLOARB three-dimensional graph. (g) FLOARMAB estimation. (h) FLOARMAB three-dimensional graph. (i) Fractional lower-order Choi-Williams distribution time-frequency representation method. (j) Wave shape and spectrum of the envelope signal.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)


### 4.2.3. Remarks

The ARMA model and FLOARMA model bispectrum estimations of the signalx(n) are shown in Figure 9 under Gaussian noise environment (SNR=−5dB). 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. The results show that both methods can estimate out the bispectrum of the signal x(n) well under Gaussian noise.Figure10 shows the ARMA and FLOARMA model bispectrum estimations of the signal x(n) under SαS distribution noise (α=1.2; GSNR=20dB). The simulation results show that the existing ARMA bispectrum method in Figures 10(a) and 10(b) fails under SαS noise environment, but the improved FLOARMA model bispectrum method in Figures 10(c) and 10(d) has good toughness. Hence, the ARMA bispectrum method only works for Gaussian noise environment, but the FLOARMA bispectrum method can be applied in Gaussian and SαS noise environment, which is robust.


## 4.2.1. Principle

A fractional lower-order autoregressive moving average (FLOARMA) model i.i.dSαS process xn can be given by(32)xn=−∑m=1paixn−m+∑m=0qbiun−m,where p and q are orders of the AR and MA model, respectively. ami=1,2,…,p and bmm=0,1,2,…,q are their parameters, which are real numbers. un is an independent identically distributed (i.i.d) SαS random process, α is its characteristic index, and its dispersion coefficient is γu.Simplifying equation (32), we have(33)∑m=0paixn−m=∑m=0qbiun−m,where a0=1 and b0=1. Taking z transform on both sides of equation (33), we obtain(34)∑m=0paiz−mXz=∑m=0qbiz−mUz,(35)1+∑m=1paiz−mXz=1+∑m=1qbiz−mUz,(36)Hz=XzUz=1+∑m=1qbiz−m1+∑m=1paiz−m=BzAz,(37)hm↔Hz,xn=∑m=0qhmun−m,where hm is system impulse response coefficient,Hz is system transfer function, and un−m is i.i.d SαS random process. The AR model parameters changes with ami=1,2,…,p; Az and Bz are the FIR and IIR filters, respectively.Takingα-order moment with equation (36),(38)H˜z=Hzα=∑m=0qhmz−mα=1+∑m=1pbmz−m1+∑m=1pamz−mα,z=ejω, and z=1 on the unit circle; then(39)H˜ejω=1+∑m=1pbme−jωm1+∑m=1pame−jωmα.According to the definition of the ARMA model bispectrum andα-order moment of Hz in (38) and (39), we define fractional lower-order ARMA model bispectrum (FLOARMAB) as(40)FLOARMABω1,ω2=γuH˜ω1H˜ω2H˜∗ω1+ω2,and FLOARB on the unit circle is written as(41)FLOARMABejω1,ejω2=γuH˜ejω1H˜ejω2H˜∗ejω1+ejω2,where γu is the dispersion coefficient of the driving SαS random process un and ∗ is conjugate operation.To obtain the coefficients of the fractional lower-order ARMA modelami=0,1,2,…,p and bmm=0,1,2,…,q, we should multiply xn−i on both sides of equation (32) taking fractional lower-order covariance to get(42)Dxxi=−∑k=1pakDxxi−k+∑k=0qbkDuxi−k,where(43)Dxxi=ΔExnP1xn−iP2,(44)Dxxi−k=ΔExnP1xn−i+kP2,(45)Duxi−k=ΔEunP1xn−i+kP2.It can be shown [24] that, for the fractional lower-order covariance of the signal and noise in (45),(46)Duxi=0,i>0,γuh−iα−1,i<0.Then, ifi>q, we have(47)Dxxi=∑k=1pakDxxi−k,i>q,(48)DxxqDxxq−1…Dxxq−p+1Dxxq+1Dxxq…Dxxq−p+2⋮⋮⋱⋮Dxxq+p−1Dxxq+p−2…Dxxqa1a2⋮ap=Dxxq+1Dxxq+2⋮Dxxq+p.Equations (47) and (48) are a generalized Yule-Walker equation, and we can obtain the fractional lower-order AR model (FLOARM) coefficients aki=0,1,2,…,p by solving them.If0≤m≤q, equation (42) changes as(49)Dxxi=∑k=0pakDxxi−k+∑k=0qbkDuxi−k=∑k=1pakDxxi−k+∑k=0qbkhi−kα−1,0≤i≤q.Lettingp=0, we have(50)Dxxm=γu∑k=0p−ibk+ibkα−1,0≤i≤q,γu∑k=0p−ibkbk−iα−1,−q≤i≤0.We can obtain the coefficientsbkk=1,2,…,q of the fractional lower-order MA model by solving the nonlinear equation (35).A finiteq′-order FLOMA model can be equivalent by an approximate infinite p′-order FLOAR modelp′≫q′; we have(51)ai+∑k=1q′bkai−k=ei,where em is error. Letting em=0, we obtain(52)∑k=1q′bkai−k=−ai.We should multiplyaPn−i on both sides of equation (52). We take fractional lower-order covariance to get(53)∑k=1q′bkR˜i−k=−R˜i,(54)R˜0R˜−1…R˜1−q′R˜1R˜0…R˜2−q′⋮⋮⋱⋮R˜q′−1R˜q′−2…R˜0b1b2⋮bp=−R˜1R˜2⋮R˜q′,where R˜i=EaiaPn−i and R˜i−k=Eai−kaPn−i+k. We can also obtain the coefficients bii=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 akk=0,1,2,…,p of the FLOAR model.Step 2: obtaining the coefficientsbkk=1,2,…,q of the FLOMA mode by solving the nonlinear equation (50) or the Yule-Walker equation in (54).Step 3: determining the orderp and q employing FPE criterion.Step 4: computing fractional low-order ARMA model bispectrum ofxn by substituting akk=0,1,2,…,p and bkk=1,2,…,q into equation (40) or equation (41).


## 4.2.2. Application Review

In this simulation,yn in equation (14) is used as the experimental signal. The ARMA model and FLOARMA model bispectrum methods are compared to demonstrate their performance under Gaussian distribution noise (SNR=−5dB) and SαS distribution noise (α=1.2; GSNR=20dB), and the simulations are shown in Figures 9–12.Figure 12
The fractional lower-order bispectrum estimations of the outer race fault signal. (a) FLODB estimation. (b) FLODB three-dimensional graph. (c) FLOIDB estimation. (d) FLOIDB three-dimensional graph. (e) FLOARB estimation. (f) FLOARB three-dimensional graph. (g) FLOARMAB estimation. (h) FLOARMAB three-dimensional graph. (i) Fractional lower-order Choi-Williams distribution time-frequency representation method. (j) Wave shape and spectrum of the envelope signal.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)


## 4.2.3. Remarks

The ARMA model and FLOARMA model bispectrum estimations of the signalx(n) are shown in Figure 9 under Gaussian noise environment (SNR=−5dB). 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. The results show that both methods can estimate out the bispectrum of the signal x(n) well under Gaussian noise.Figure10 shows the ARMA and FLOARMA model bispectrum estimations of the signal x(n) under SαS distribution noise (α=1.2; GSNR=20dB). The simulation results show that the existing ARMA bispectrum method in Figures 10(a) and 10(b) fails under SαS noise environment, but the improved FLOARMA model bispectrum method in Figures 10(c) and 10(d) has good toughness. Hence, the ARMA bispectrum method only works for Gaussian noise environment, but the FLOARMA bispectrum method can be applied in Gaussian and SαS noise environment, which is robust.


## 5. Application Simulations

The 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 Figure3(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 (f1>0, f2>0) to analyze the signals.Figures11(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. Their 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.Figures12(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. The 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=20dB) is added in the α-stable distribution outer race fault signal as the actual working environment background noise. The conventional bispectrum methods including the direct bispectrum, indirect bispectrum method, AR model, ARMA bispectrum methods, and the improved fractional lower-order bispectrum methods including FLODB, FLOIDB, FLOARB, and FLOARMAB methods are used to analyze the outer race bearing fault signal under SαS distribution background noise environment. The simulations are given in Figures 13 and 14. The simulation results show that the existing bispectrum methods in Figures 13(a)–14(h) fail, but the proposed fractional low-order bispectrum methods in Figures 14(a)–14(h) have good performance, and it is observed that the transient harmonic vibration components are near the central frequencies of 600 Hz and 2800 Hz.Figure 13
The conventional bispectrum estimations of the outer race fault signal contaminated bySαS distribution noise. (a) The conventional direct bispectrum estimation. (b) The conventional direct bispectrum three-dimensional graph. (c) The conventional indirect bispectrum estimation. (d) The conventional indirect bispectrum three-dimensional graph. (e) The conventional AR model bispectrum estimation. (f) The conventional AR model bispectrum three-dimensional graph. (g) The conventional ARMA model bispectrum estimation. (h) The conventional ARMA model bispectrum three-dimensional graph.
(a)(b)(c)(d)(e)(f)(g)(h)Figure 14
The improved fractional lower-order bispectrum estimations of the outer race fault signal contaminated bySαS distribution noise. (a) The improved FLODB estimation. (b) The improved FLODB three-dimensional graph. (c) The improved FLOIDB estimation. (d) The improved FLOIDB three-dimensional graph. (e) The improved FLOAR model bispectrum estimation. (f) The improved FLOAR model bispectrum three-dimensional graph. (g) The improved FLOARMA model bispectrum estimation. (h) The improved FLOARMA model bispectrum three-dimensional graph.
(a)(b)(c)(d)(e)(f)(g)(h)


## 6. Conclusions

The bearing fault signals are a non-Gaussian and nonstationary process, andα-stable distribution is a more appropriate statistical model for them. The improved FLODB, FLOIDB, FLOAR, and FLOARMA model bispectrum methods have been proposed for the fault signals employing fractional low-order statistics. The 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 time-frequency 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 fault signals so as to realize the recognition of the fault characteristics. The 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.



---

*Source: 8823389-2020-12-31.xml*