# Applications of Fractional Lower Order Frequency Spectrum Technologies to Bearing Fault Analysis

Authors: Junbo Long; Haibin Wang; Peng Li
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

The traditional spectral analysis method is used to study the characteristics of bearing fault signals in frequency domain, which is reasonable and effective in general cases. However, it is proved that the fault signals have heavy tails in this paper, which areα stable distribution, and 1<α<2, and even the noises belong to α stable distribution. Then the conventional spectral analysis methods degenerate and even fail under α stable distribution environment. Several improved frequency spectral analysis methods are proposed employing fractional lower order covariation or fractional lower order covariance in this paper, including fractional lower order Blackman-Tukey covariation spectrum (FLOBTCS), fractional lower order periodogram covariation spectrum (FLOPCS), and fractional lower order welch covariation spectrum (FLOWCS). In order to suppress side lobe and improve resolution, we present novel fractional lower order autoregression (FLO-AR) and fractional lower order autoregressive moving average (FLO-ARMA) parameter model frequency spectrum methods, and the calculation steps are summarized. The proposed spectrum methods are compared with the existing methods based on second-order statistics under Gaussian and SαS distribution environments, and the results show that the new algorithms have better performance than the traditional methods. Finally, the improved methods are applied to estimate frequency spectrums of the normal and outer race fault signals, and it is demonstrated that they are effective for fault diagnosis.


---

## Body

## 1. Introduction

The mechanical fault vibration signals are a nonstationary and non-Gaussian process; their spectrum characteristic change with the time and their probability density function (PDF) have heavy tails. Hence, it is necessary to study their parameter characteristics and new methods suitable for the signals. Bearing fault diagnosis is a technology that can monitor the running state of mechanical equipment in real time. The frequency spectrum analysis is an important means for fault feature extraction in engineering application. The traditional spectrum estimation methods mainly include Blackman-Tukey (BT), Periodogram, Welch, multiple signal classification, wavelet spectrum, and parameter model spectrum estimation. The spectral analysis methods have been widely applied in radar signal processing [1], passive sonar signal processing [2], target recognition [3], seismic analysis [4], mechanical failure analysis, and other fields [5–7]. Recently, some improved methods were also applied to machine breakdown analysis, such that singular spectrum analysis was used for roller element bearing fault diagnosis [8, 9]. A new wavelet spectrum analysis was applied to incipient bearing fault diagnostics [10]. The frequency spectrum analysis method employing spectrum sparseness was proposed in [11], which was applied to the rolling bearing performance degradation assessment. Combining the segmentation point obtained from pretreatment and the spectrum analysis, SU Wei-Jun et al. proposed a new rolling bearing fault feature extraction method employing the local spectrum [12]. A new detection method of the change in amplitude was presented based on principal component analysis, which was used to detect the rolling element bearing fault [13]. Pandarakone and the others proposed a distinct motor bearing analysis method employing support vector machine for frequency spectrum determination [14], although the existing frequency spectrum methods have been well applied, which cannot work in peaked noise environment.The parameter model spectrum estimation method is a modern technology with practical significance, which mainly includes AR (autoregression), MA (moving average), and ARMA (autoregressive moving average) model methods. Recently, the methods are developed rapidly. A new method employing empirical mode decomposition (EMD) and AR model was proposed for fault feature extraction in [15]; firstly, EMD was applied to decompose the fault signals, and then each intrinsic mode function component was analyzed for feature extraction. The spectral kurtosis method based on AR model spectrum was applied to fault diagnosis and condition monitoring for rolling bearings [16]. A novel method based on fuzzy cluster analysis and AR model was presented for the bearing fault diagnosis; firstly, the parameters of the AR model were estimated based on the higher-order statistics, and then fuzzy cluster analysis was used for classification and pattern recognition [17]. The bispectrum analysis method based on ARMA model and fuzzy c-means was proposed in [18], which was applied to fault identification, and an algorithm employing AR neural networks and ARMA model has been applied to fault diagnosis in [19]. Although the traditional parameter model methods have been applied in many fields, which still have some defects, the methods degrade in performance under impulsive signal or noise environment and even fail. Therefore, it is of great significance to explore high performance and efficient ARMA model spectrum estimation algorithm.The conventional frequency spectrum estimation methods mentioned above are suitable for rotating machine condition monitoring and performance evaluation in general. However, probability density functions (PDFs) of the bearing fault signals have heavy trails in some special cases, even the same with the noise in the signals, which belong toα stable distribution [20–25]. Because α stable distribution has no finite second moment, the existing methods based on Gaussian hypothesis and second-order statistics degenerate under α stable distribution environment. We cannot reuse the frequency domain analysis and its extended analysis employing the subsistent power spectrum methods. Therefore, new spectral estimation methods applicable to α stable distribution are needed to seek. Recently, it is verified that the bearing fault signals belong to α stable distribution [26–28]. Furthermore, the method combining α stable distribution and support vector machine was proposed in [29], which was used in fault diagnosis of the gearbox. A new empirical mode decomposition method based on the parameters of α stable distribution was presented for fault diagnosis of the low-speed rolling bearings, which has better diagnosis accuracy and operational efficiency [30]. A short time Fourier transform local maxima (STFT-LM) spectrum based on FLOC was proposed for local damage detection [31]; the method combines the local maxima enhancement with dependency analysis in the frequency bands, which do not require second moment to be finite and can work in case when the sample distribution processes heavy tails.Aiming at the degradation of the traditional frequency spectrum methods underα stable distribution environment, the improved fractional lower order Blackman-Tukey covariation, fractional lower order Periodogram covariation, and fractional lower order Welch covariation methods are proposed based on fractional lower order covariation or fractional lower order covariance in this paper. The results show that the proposed frequency spectrum methods can better work in Gaussian noise and α stable distribution noise environment, which are robust, and their performance is better than the existing frequency spectrum methods based on second-order statistics. In order to suppress the side lobe and improve the resolution, fractional lower order autoregression (FLO-AR) and fractional lower order autoregressive moving average (FLO-ARMA) frequency spectrum methods are presented, which can effectively estimate out frequency spectrum of the signal from α stable distribution noise. Their mixed mean square error is significantly lower than that of the existing AR and ARMA methods under different characteristics index α and generalized signal noise ratio (GSNR); especially, when GSNR<19dB or α<1, the performance advantage of the FLO-AR and FLO-ARMA methods is more obvious. Finally, the improved methods are applied to analyze the outer race fault signal; the results show that the methods can better estimate out transient harmonic vibration components of the bearing fault signals under α stable noise environment, but the conventional methods degenerate.In this paper, several improved frequency spectrum estimation methods employing fractional lower order covariation or covariance are proposed for the bearing fault diagnosis under Gaussian orα stable distribution environment. This paper is structured in the following manner. α stable distribution and the bearing fault signals are introduced in Section 2. The improved fractional lower order frequency spectrum estimation methods are demonstrated, and simulation comparisons with the conventional spectrum analysis algorithms based on second-order statistics are performed to demonstrate justifiability of the improved methods in Section 3. The conjoint analysis of the outer race fault signals employing the proposed frequency spectrum and time-frequency representation methods in [28] are demonstrated in Section 4. Finally, conclusions and future research are given in Section 5.


## 2.α Stable Distribution and Bearing Fault Signals

### 2.1.α Stable Distribution

Probability density function (PDF) ofα stable distribution has heavy tail, and its characteristic function is defined as(1)φt=exp⁡jμt-γtα1+jβsign⁡tωτ,αwhere sign⁡(t)=1t>00t=0-1t<0,ω(τ,α)=tan⁡(απ/2)α≠12/πlog⁡τα=1. α stable distribution is a generalized Gaussian process, and α is its characteristic index, when α=2, which is Gaussian distribution, and when 0<α<2, which is lower order α stable distribution. μ is location parameter, and β and γ are symmetry and dispersion coefficient, respectively. When β=0, μ=0 and γ=1, which is called standard symmetric α stable distribution (SαS). The waveforms of SαS stable distribution are shown in Figure 1 under α = 0.5, 1.0, 1.5, and 2.0.Figure 1
The waveforms ofSαS stable distribution when α= 0.5, 1.0, 1.5, and 2.0.Figure2 is PDFs of SαS stable distribution under α=0.5, 1.0, 1.5, and 2.0, respectively. From Figure 2(a), we see that when x=0, the value of f(x) under α=2 (Gaussian distribution) is the smallest, and the value of f(x) under α=0.5 is the largest, which have pulse characteristics. Figure 2(b) shows that the value of f(x) under α=2 is equal to 0 when x is approximately 4.75, but the value of f(x) is not equal to 0 under α<2, which have tails. Hence, we can know that the smaller the characteristic index α of α stable distribution is, the thicker its tail is, and the more significant its pulse characteristic is. On the contrary, when the value of α increases, the tail becomes thinner and the pulse characteristic decreases.Figure 2
PDFs ofSαS stable distribution under different α.
(a) (b)


### 2.2. Bearing Fault Signals

The actual bearing fault signals data are obtained from the Case Western Reserve University (CWRU) bearing data center [32]. The test stand consists of a 2 hp motor (left), a torque transducer/encoder (center), a dynamometer (right), and control electronics (not shown). The test bearings support the motor shaft. Experiments are conducted using a 2 hp Reliance Electric motor, and acceleration data is measured at locations near to and remote from the motor bearings. The experimental equipment adopts 6205-2RS JEM SKF type bearing, the motor bearings are seeded with faults using electro-discharge machining, and the vibration data is collected using accelerometers, which are attached to the housing with magnetic bases. 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 out 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 centered at 6:00. The vibration signals are collected using a 16 channel DAT recorder, and the speed and horsepower data are collected using the torque transducer/encoder and were recorded by hand. The normal signals are given in Figure 3(a), and the fault signals of inner race, ball, and out race are shown in Figures 3(b), 3(c), and 3(d), respectively. We can know that the waveform of the fault signals has a certain impulse.Figure 3
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 have appliedα stable distribution statistical model to estimate the parameters of the normal signal and the fault signals of inner race, ball, and outer race, respectively, the results are given in Table 1. As it can be seen, the characteristic index of the normal signals is equal to 2, which are Gaussian distribution. However, the characteristic index of the bearing fault signals is greater than 1 but smaller than 2, which belongs to non-Gaussian α stable distribution (α<2).Table 1
Theα stable distribution parameters of the bearing fault signals.
Parameters  α β γ μ Normal DE 2.000 -0.2863 0.0532 0.0121 FE 2.000 1.000 0.0583 0.0236  Inner race fault BA 1.7682 0.0872 0.0590 0.0062 DE 1.4195 0.0155 0.2407 0.0175 FE 1.8350 0.0322 0.1495 0.0291  Ball fault BA 1.9790 0.0592 0.0293 0.0055 DE 1.8697 0.1215 0.0772 0.0193 FE 1.998 -0.0371 0.0674 0.0321  Outer race fault BA 1.6077 -0.1731 0.0530 0.0012 DE 1.1096 0.0433 0.1341 0.0367 FE 1.5435 -0.0169 0.0968 0.0296PDFs of the signals of inner race fault, ball fault, and outer race fault are shown in Figures4(a), 4(b), and 4(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 4 shows 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 4
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)


## 2.1.α Stable Distribution

Probability density function (PDF) ofα stable distribution has heavy tail, and its characteristic function is defined as(1)φt=exp⁡jμt-γtα1+jβsign⁡tωτ,αwhere sign⁡(t)=1t>00t=0-1t<0,ω(τ,α)=tan⁡(απ/2)α≠12/πlog⁡τα=1. α stable distribution is a generalized Gaussian process, and α is its characteristic index, when α=2, which is Gaussian distribution, and when 0<α<2, which is lower order α stable distribution. μ is location parameter, and β and γ are symmetry and dispersion coefficient, respectively. When β=0, μ=0 and γ=1, which is called standard symmetric α stable distribution (SαS). The waveforms of SαS stable distribution are shown in Figure 1 under α = 0.5, 1.0, 1.5, and 2.0.Figure 1
The waveforms ofSαS stable distribution when α= 0.5, 1.0, 1.5, and 2.0.Figure2 is PDFs of SαS stable distribution under α=0.5, 1.0, 1.5, and 2.0, respectively. From Figure 2(a), we see that when x=0, the value of f(x) under α=2 (Gaussian distribution) is the smallest, and the value of f(x) under α=0.5 is the largest, which have pulse characteristics. Figure 2(b) shows that the value of f(x) under α=2 is equal to 0 when x is approximately 4.75, but the value of f(x) is not equal to 0 under α<2, which have tails. Hence, we can know that the smaller the characteristic index α of α stable distribution is, the thicker its tail is, and the more significant its pulse characteristic is. On the contrary, when the value of α increases, the tail becomes thinner and the pulse characteristic decreases.Figure 2
PDFs ofSαS stable distribution under different α.
(a) (b)


## 2.2. Bearing Fault Signals

The actual bearing fault signals data are obtained from the Case Western Reserve University (CWRU) bearing data center [32]. The test stand consists of a 2 hp motor (left), a torque transducer/encoder (center), a dynamometer (right), and control electronics (not shown). The test bearings support the motor shaft. Experiments are conducted using a 2 hp Reliance Electric motor, and acceleration data is measured at locations near to and remote from the motor bearings. The experimental equipment adopts 6205-2RS JEM SKF type bearing, the motor bearings are seeded with faults using electro-discharge machining, and the vibration data is collected using accelerometers, which are attached to the housing with magnetic bases. 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 out 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 centered at 6:00. The vibration signals are collected using a 16 channel DAT recorder, and the speed and horsepower data are collected using the torque transducer/encoder and were recorded by hand. The normal signals are given in Figure 3(a), and the fault signals of inner race, ball, and out race are shown in Figures 3(b), 3(c), and 3(d), respectively. We can know that the waveform of the fault signals has a certain impulse.Figure 3
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 have appliedα stable distribution statistical model to estimate the parameters of the normal signal and the fault signals of inner race, ball, and outer race, respectively, the results are given in Table 1. As it can be seen, the characteristic index of the normal signals is equal to 2, which are Gaussian distribution. However, the characteristic index of the bearing fault signals is greater than 1 but smaller than 2, which belongs to non-Gaussian α stable distribution (α<2).Table 1
Theα stable distribution parameters of the bearing fault signals.
Parameters  α β γ μ Normal DE 2.000 -0.2863 0.0532 0.0121 FE 2.000 1.000 0.0583 0.0236  Inner race fault BA 1.7682 0.0872 0.0590 0.0062 DE 1.4195 0.0155 0.2407 0.0175 FE 1.8350 0.0322 0.1495 0.0291  Ball fault BA 1.9790 0.0592 0.0293 0.0055 DE 1.8697 0.1215 0.0772 0.0193 FE 1.998 -0.0371 0.0674 0.0321  Outer race fault BA 1.6077 -0.1731 0.0530 0.0012 DE 1.1096 0.0433 0.1341 0.0367 FE 1.5435 -0.0169 0.0968 0.0296PDFs of the signals of inner race fault, ball fault, and outer race fault are shown in Figures4(a), 4(b), and 4(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 4 shows 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 4
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 Frequency Spectrum Estimation Methods

### 3.1. Fractional Lower Order Covariation Spectrum Estimation Method

#### 3.1.1. Covariation Spectrum Principle

The autocorrelation function and power spectrum of a second-order random process can be expressed as(2)Rxxτ=Extx∗t-τ(3)Sxxω=∫-∞∞Rxxτe-jωτdτ(4)Rxxτ=12π∫-ππSxxωejωτdωFrom (2)-(4), we can see that power spectrum of x(t) is Fourier transform of its correlation function, which is called Wiener-Khinchin theorem. Rxx(τ) demonstrates second-order statistical properties of x(t) in time domain, and Sxx(ω) is the statistical properties of x(t) in frequency domain. The corresponding discrete correlation function and power spectrum are given by(5)R^xxm=1N∑n=0N-1xnx∗n-m(6)S^xxω=∑m=0N-1R^xxmejωm=1N∑m=0N-1∑n=0N-1xnx∗n-me-jωmWhenx(t) is α stable distribution process, the autocorrelation function and power spectrum based on second-order statistics are no longer applicable. Hence, fractional lower order covariation was proposed in [21], which was obtained by calculating fractional low order moments.(7)X,Yα=EXY<P-1>EYpγywhere p is a given real constant, 1≤p<α≤2 [19–22], and p denotes p order moment of x(t). When x(t) is a real signal, xp(t)=x(t)p·sign[x(t)], signx(t)=1x(t)>00x(t)=0-1x(t)<0, and when x(t) is a complex signal, xp(t)=x(t)p-1·x∗(t), ∗ denotes conjugate operation. γy is dispersion coefficient of Y. According to (7), the autocovariation function of x(t) can be expressed as (8)Cxxτ=xt,xt-τα=Extxt-τ<p-1>Ext-τpγxt-τ,1≤p<α≤2When x(t) is Gaussian process, then α=2 and p=2, and (8) degenerates into autocorrelation function(9)Rxxτ=xt,xt-τα=xt,xt-τ2=Extxt-τWe take Fourier transform on both sides of (8) and get(10)φxxω=∫-∞∞Cxxτe-jωτdτ(11)Cxxτ=12π∫-ππφxxωejωτdωEquation (10) is known as fractional low order covariation spectrum under α stable distribution process, which is corresponding to power spectrum under Gaussian distribution. When α=2 and p=2, (10) changes as power spectrum; hence, fractional low order covariation spectrum is a generalized power spectrum.


#### 3.1.2. Fractional Lower Order Blackman-Tukey Covariation Spectrum Estimation

We assume thatx(n),n=0,1,⋯,N-1 is N samples of a discrete α stable distribution random process, and the sample mean is 0, then(12)C^xxm=1N∑n=0N-1xnx<p-1>n-m(13)φ^xxω=∑m=-MMC^xxmejωm=1N∑m=-MM∑n=0N-1xnx<p-1>n-me-jωmEquation (12) can be applied to estimate the covariation function of x(n). By taking discrete Fourier transform of (12), Blackman-Tukey covariation spectrum estimation can be obtained, as shown in (13).


#### 3.1.3. Fractional Lower Order Periodogram Covariation Spectrum Estimation

LettingM=N-1 and k=n-m, substituting them into (13),(14)
ϕ
^
x
x
ω
=
1
N
∑
m
∑
k
x
n
x
<
p
-
1
>
k
e
-
j
ω
m
=
1
N
∑
n
=
0
N
-
1
x
n
e
-
j
ω
n
∑
k
=
0
N
-
1
x
<
p
-
1
>
k
e
j
ω
k
=
1
N
∑
n
=
0
N
-
1
x
n
e
-
j
ω
n
∑
k
=
0
N
-
1
x
<
p
-
1
>
k
e
-
j
ω
k
∗
=
1
N
∑
n
=
0
N
-
1
x
n
e
-
j
ω
n
∑
n
=
0
N
-
1
x
<
p
-
1
>
n
e
-
j
ω
n
∗
=
1
N
X
n
X
-
n
∗where X(ω)=∑n=0N-1x(n)e-jωn is discrete Fourier transform of x(n). The discrete Fourier transform of x<p-1>(n) is given by X-(n)=∑n=0N-1x<p-1>(n)e-jωn. From (14), we firstly compute discrete Fourier transform of x(n) and x<p-1>(n) and take their product; then fractional lower order Periodogram covariation spectrum estimation can be obtained.


#### 3.1.4. Fractional Lower Order Welch Covariation Spectrum Estimation

Covariation spectrums of the fractional lower order Blackman-Tukey and fractional lower order Welch methods can be gotten by Fourier transform, the principle is simple, and the calculation quantity is also small. However, FLOBTCS and FLOPCS have higher side lobe, and the variance is larger, which results in poor performance and affects their practical application. Therefore, the traditional Welch method needs to be improved.By dividingN discrete signal x(n),n=0,1,⋯,N-1 into partially overlapped K segments, the length of each segment is L, the overlap length between two adjacent segments is L/2, and then the ind segment can be expressed as(15)xin=xn+i-1L,i=1,2,⋯,K,n=0,1,⋯,L-1(16)K=N-L/2L/2By taking the windowed Fourier transform of x(n) and x<p-1>(n), the following equations are obtained:(17)Xin=∑n=0L-1xinhne-jωn,i=1,2,⋯,K(18)X-in=∑n=0N-1xi<p-1>nhne-jωn,i=1,2,⋯,KSubstituting (17)-(18) into (14), the ind covariation spectrum can be written as(19)φ^xiω=1NXinX-in∗,i=1,2,⋯,KCalculating the average of the covariation spectrum of all segments, then(20)φ^xxω=1K∑i=1Kφ^xiω=1K∑i=1K1NXinX-in∗,i=1,2,⋯,KEquation (20) is called fractional low order Welch covariation spectrum method. In order to reduce variance of the covariation spectrum estimation, the number of segments K can be increased appropriately, which is bound to reduce length of the segment L, so that the main lobe of the covariation spectrum is widened to reduce the resolution. In real application, K and L can be set according to the specific needs.


#### 3.1.5. Application Review

We will now apply the FLOBTCS and FLOPCS and FLO-Welch covariation spectrum methods to estimate frequency spectrum of the test signal. The test signaly(n) is defined as three sinusoidal signals with Gaussian or SαS stable distribution noise.(21)yn=Asin⁡2πf0n+sin⁡2πf1n+sin⁡2πf2n+vn=Axn+vnwhere f0=100, f1=300, f2=500, n=0,1,2,....,N-1, and N=2000. When v(n) is Gaussian noise, generalized signal to noise ratio (GSNR) is used to instead of SNR.(22)GSNR=10log10⁡Exn2γα==10log10⁡1Nγα∑n=0N-1xn2where γ is the dispersion coefficient of SαS stable distribution noise. According to the given GSNR, the amplitude of the signal can be written as(23)A=10GSNR/101/N∑n=0N-1xn2γα1/2LettingGSNR=20dB, the fractional lower order covariation spectrum methods and traditional power spectrum methods are used to estimate frequency spectrum of the signal x(n) under Gaussian distribution noise and SαS stable distribution noise; the simulation results of ten independent experiments are shown in Figures 5–8.Figure 5
Fractional lower order Blackman-Tukey covariation spectrums and power spectrums of the signalx(n) under Gaussian and SαS noise environment ((a) Gaussian noise environment, SNR=-5dB; (b) SαS noise environment, α=1.3, and GSNR=20dB).
(a) (b)Figure 6
Fractional lower order Periodogram covariation spectrums and power spectrums of the signalx(n) under Gaussian and SαS noise environment ((a) Gaussian noise environment, SNR=-5dB; (b) SαS noise environment, α=1.3, and GSNR=20dB).
(a) (b)Figure 7
Fractional lower order Welch covariation spectrums and power spectrums of the signalx(n) under Gaussian and SαS noise environment ((a) Gaussian noise environment, SNR=-5dB; (b) SαS noise environment, α=1.3, and GSNR=20dB).
(a) (b)Figure 8
Comparisons of fractional lower order Blackman-Tukey, Periodogram, and Welch covariation spectrums of the signalx(n) ((a) SαS (α=1.1) noise environment; (b) SαS (α=1.8) noise environment).
(a) (b)


#### 3.1.6. Remarks

Figure5(a) shows that the traditional Blackman-Tukey spectrum estimation method fails under SαS noise environment, but the improved FLOBTCS method in Figure 5(b) has good toughness. Figures 6(a) and 6(b) are frequency estimations of the signal x(n) employing the Periodogram and FLOPCS methods, respectively; the result reveals the advantage of the FLOPCS method. The Welch spectrum and FLO-Welch covariation spectrum of the signal x(n) are given in Figure 7; we can know that the Welch power spectrum method cannot work under SαS noise environment, but the proposed FLO-Welch covariation method demonstrates good performance. Figure 8 is frequency spectrum of the signal x(n) employing FLOBTCS, FLOPCS, and FLO-Welch covariation spectrum methods under SαS (α=1.1 and α=1.8) noise environment. The results show that the conventional power spectrum methods degenerate under SαS noise environment, but the improved fractional low order covariation spectrum methods can better estimate out the frequency spectrum of the signal x(n).The FLOBTCS and FLOPCS methods have a certain side lobe and little amount of calculation, which cause the weak signals to be submerged and lead to the distortion of the signals. FLO-Welch algorithm can suppress the side lobe effects by spectrum segmentation. However, the resolution and variance of the FLO-Welch method are controlled by the number of segments and length like the conventional Welch method. In applications, the shorter number of segments should be used when we would like to get higher frequency resolution, but variance gets larger, and if we wish to reduce estimation variance, a longer number of segments are preferred. Hence, we can compromise the number and length of the segments by the actual situation. As a result, the power spectrum methods are only suitable to analyze the signals under Gaussian environment, but the fractional low order covariation spectrum methods can be applied under Gaussian andSαS noise environment.


### 3.2. FLO-AR Model Spectrum Estimation Method

#### 3.2.1. Principle

An independent identically distributed (i.i.d)α process can be expressed by p order AR model(24)xn=a1xn-1+⋯+apxn-p+un=∑i=1paixn-i+unwhere ai(i=1,2,⋯,p) is the parameters of AR model. u(n) is an i.i.d process, and its characteristic index and dispersion coefficient are α and γu, respectively. When u(n) passes a finite impulse response (FIR) filter, then(25)xn=∑i=0qhiun-iTaking z transformation for (24)-(25),(26)Hz=∑i=0qhiz-i=11-∑i=1paiz-iwhere H(z) is transfer function of the filter. α spectrum of the signal x(n) was defined as [21](27)Sαz=Xn,Wznα=γuH1zα-1Hzα-1where(28)Wzn=∑i=-qi=qXn-iziWhen z=ejω, substituting z=1 and (26) into (27), AR model α spectrum of an α stable distribution process on the unit circle can be written as(29)Sαejω=γuH1ejωα-1Hejωα-1=γuHejωHejωα-1=γuHejωα=γu11-∑i=1paie-jωiαThen AR model α spectrum of a SαS (γu=1) process on the unit circle is given by(30)SSαSejω=Hejωα=11-∑i=1paie-jωiαAccording to the definition in [21], we can obtain fractional low order autocovariance of the signal x(n).(31)Dxxm=ΔExn<a>xn-m<b>,0≤a<α2,0≤b<α2The corresponding low order autocovariance estimation can be written as(32)D^xxm=∑n=L1+1L2xnaxn-mbsignxnxn-mL2-L1where L1=max⁡(0,-m), L2=min⁡(N-m,N), and N is number of the observation signal x(n).We defineA=[a^1,a^2,⋯,a^p]T, B=[D^xx(1),D^xx(2),⋯,D^xx(P)]T, and(33)Τ=D^xx0D^xx-1⋯D^xx1-pD^xx1D^xx0⋯D^xx2-p⋮⋮⋱⋮D^xxp-1D^xxp-2⋯D^xx0where Τ is call fractional low order autocovariance matrix (FLOACM). We apply FLOACM to fractional low order moment (FLOM) matrix in [21] and then get(34)ΤA=BAnd(35)D^xx0D^xx-1⋯D^xx1-pD^xx1D^xx0⋯D^xx2-p⋮⋮⋱⋮D^xxp-1D^xxp-2⋯D^xx0a^1a^2⋮a^p=D^xx1D^xx2⋮D^xxpWe call (35) fractional low order Euler Walker equation. The parameters a^i(i=1,2,⋯,p) can be gotten by solving (35). When a^i(i=1,2,⋯,p) are substituted into (26) or (29), fractional low order AR model spectrum can be obtained. 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. The calculation formula can be written as(36)FPEp=σa2N+pN-pwhere σa2 is the variance of residuals.We summarize the steps of fractional low order AR model spectrum method, as follows.Step 1.
Solve fractional low order AR model parametersai(i=1,2,⋯,p) of x(n) employing (35).Step 2.
Determine the orderp with the FPE criterion.Step 3.
Compute fractional low order AR model spectrum ofx(n) by substituting ai(i=1,2,⋯,p) into (29) or (30).


#### 3.2.2. Application Review

In this simulation, the test signal isx(n) in (21). We compare the performance of the traditional AR model method and the improved fractional low order AR model method under Gaussian distribution noise (SNR=-5dB) and SαS stable distribution noise (α=1.3, GSNR=20dB); the results of the ten independent experiments are shown in Figures 9-10. In order to further verify the effectiveness of the fractional low order AR method, we conduct comparative experiment on two methods under different α when GSNR=20, as show in Figure 11. When α=1.3 and GSNR changes from 14dB to 24dB, Figure 12 shows the change of Errors Power of the exiting AR model method and the improved FLO-AR method under α stable distribution noise environment.Figure 9
The traditional AR model spectrums and improved FLO-AR spectrums of the signalx(n) under Gaussian noise environment (SNR=-5dB) ((a) AR model spectrums of the signal x(n), where p=29; (b) FLO-AR spectrums of the signal x(n), where p=27 and a=b=0.8).
(a) (b)Figure 10
FLO-AR model and AR model spectrums of the signalx(n)under SαS (α=1.3) noise environment (GSNR=20dB) ((a) AR model spectrums of the signal x(n), where p=20; (b) FLO-AR model spectrums of the signal x(n), where p=21 and a=b=0.1).
(a) (b)Figure 11
Errors power of the AR model and FLO-AR model frequency spectrum algorithms under differentα (GSNR=20).Figure 12
Errors power of the AR model and FLO-AR model frequency spectrum algorithms under differentGSNR (α=1.3).


#### 3.2.3. Remarks

Figures9(a) and 9(b) show the traditional AR model spectrum and the improved FLO-AR model spectrum of the signal x(n), respectively. Both methods can estimate out the frequency spectrum of the signal x(n) well under Gaussian noise environment. Figure 10(a) reveals that the AR model spectrum method fails under SαS noise environment (GSNR=20dB), but the FLO-AR model spectrum method in Figure 10(b) has good toughness. As a result, the AR model spectrum method is only suitable to analyze the signals under Gaussian environment, but the FLO-AR model spectrum method can work under Gaussian and SαS noise environment, which is robust.Figure11 is Errors Power of the AR and FLO-AR model spectrum algorithms under different α. When the characteristic index α of SαS distribution noise changes from 0 to 2, Errors Power of the improved FLO-AR method is lower than the AR algorithm, which is very stable.WhenGSNR changes from 14dB to 24dB, the Errors Power of the AR model and FLO-AR model methods are given in Figure 12. The results show that Errors Power of the improved FLO-AR method is lower than that of the AR method, and they are very stable with the changing of GSNR. Especially, when GSNR<19dB, the performance superiority of the FLO-AR method is more obvious.


### 3.3. FLO-ARMA Model Spectrum Estimation Method

#### 3.3.1. Principle

An i.i.dα random process can be expressed by an ARMA (p,q).(37)xn=a1xn-1+⋯+apxn-p+b0un+b1un-1+⋯+bqun-q=∑i=1paixn-i+∑i=0qbiun-iwhere ai(i=1,2,⋯,p) and bi(i=0,1,2,⋯,q) are parameters of p order AR model and q order MA model, respectively. Taking Z transform on (37), obtaining(38)-∑i=0paixn-i=∑i=0qbiun-i(39)-∑i=0paiz-iXz=∑i=0qbiz-iUz(40)1-∑i=1paiz-iXz=1+∑i=1qbiz-iUz(41)Hz=XzUz=1+∑i=1qbiz-i1-∑i=1paiz-i=BzAz(42)hn↔Hz,xn=hn∗un=∑i=0qhiun-iwhere h(n) is impulse response function, H(z) is system transfer function, A(z) and B(z) are the FIR and IIR filters, respectively. From (41), FLO-ARMA model spectrum of α stable distribution process on the unit circle (z=ejω) is given by(43)Sαejω=γuHejωα==γu1+∑i=1qbie-jωi1-∑i=1paie-jωiαThen FLO-ARMA model spectrum of SαS (γu=1) process on the unit circle can be written as(44)SSαSejω=Hejωα=1+∑i=1qbie-jωi1-∑i=1pake-jωiαMultiplyingx(n-m) both sides of (38) and computing fractional low order covariance, then(45)Dxxm=∑i=1paiDxxm-i+∑i=0qbiDuxm-iwhere(46)Dxxm=Exn<a>xn-m<b>Dxxm-i=ΔExn<a>xn-m+i<b>Duxm-i=ΔEun<a>xn-m+i<b>Then we can get [20](47)Duxm=0m>0γuh-mα-1m<0When m>q, (45) changes as(48)Dxxm=∑i=1paiDxxm-i,m>qEquation (48) can be written as(49)DqDq-1⋯Dq-p+1Dq+1Dq⋯Dq-p+2⋮⋮⋱⋮Dq+p-1Dq+p-2⋯Dqa1a2⋮ap=Dq+1Dq+2⋮Dq+pwhere (49) is also a generalized Yule-Walker equation. Then we can obtain the parameters of the FLO-AR model by solving (49).When0≤m≤q, (45) changes as(50)Dxxm=∑i=0paiDxxm-i+∑i=0qbiDuxm-i=∑i=1paiDxxm-i+∑i=0qbihk-iα-1,0≤m≤qLetting p=0, then(51)Dxxm=γu∑i=0p-mbi+mbiα-1,0≤m≤qγu∑i=0p-mbibi-mα-1,-q≤m≤0By solving the nonlinear equation (51), the parameters bi(i=1,2,⋯,q) of the FLO-MA model can be obtained.An approximate infinitep′ order FLO-AR model can be used to be equivalent to a finite q′ order MA model (p′≫q′), then(52)am+∑i=1q′biam-i=emwhere e(m) is error. Letting e(m)=0, then(53)∑i=1q′biam-i=-amBy taking autocorrelation of both sides of (53) with a∗(n-m) one has(54)∑i=1q′biRm-i=-Rmwhere R(m)=Ea(m)a∗(n-m) and R(m-i)=Ea(m-i)a∗(n-m+i). Equation (38) can be written as(55)R0R-1⋯R1-q′R1R0⋯R2-q′⋮⋮⋱⋮Rq′-1Rq′-2⋯R0b1b2⋮bp=-R1R2⋮Rq′Equation (55) is a Toeplitz matrix. By solving the Yule-Walker equation (55), the parameters bi(i=1,2,⋯,q) of the FLO-MA model are gotten.We summarize the steps of the FLO-ARMA model spectrum method, as follows.Step 1.
Compute the parametersai(i=1,2,⋯,p) of the FLO-AR model of the signal x(n) employing (49).Step 2.
By solving the nonlinear equation (51) or the Yule-Walker equation (55), the parameters bi(i=1,2,⋯,q) of the FLO-MA model can be gotten.Step 3.
Determine the orderp and q employing FPE criterion.Step 4.
Substituteai(i=1,2,⋯,p) and bi(i=1,2,⋯,q) into (29) or (30); compute FLO-ARMA model spectrum of the signal x(n).


#### 3.3.2. Application Review

In this simulation,x(n) in (21) is used as the experimental test signal. The traditional ARMA method and modified FLO-ARMA method are applied to estimate the frequency spectrum of the signal x(n) under Gaussian noise (SNR=-5dB) and SαS noise environment (GSNR=20dB). The results of ten independent experiments are shown in Figures 13-14.Figure 13
The ARMA model and FLO-ARMA model spectrums of the signalx(n) under Gaussian noise environment (SNR=-5dB) ((a) ARMA model spectrums of the signal x(n), where p=29 and q=5; (b) FLO-ARMA model spectrums of the signal x(n), where p=27, q=6, and a=b=0.9).
(a) (b)Figure 14
FLO-ARMA model and ARMA model spectrums of the signalx(n) under SαS (α=1.3) noise environment (GSNR=20dB) ((a) ARMA model spectrums of the signal x(n), where p=20 and q=5; (b) FLO-ARMA model spectrums of the signal x(n), where p=21, q=5, and a=b=0.1).
(a) (b)


#### 3.3.3. Remarks

The ARMA frequency spectrum of the signalx(n) is shown in Figure 13(a), and Figure 13(b) is FLO-ARMA frequency spectrum of the signal x(n) under Gaussian(SNR=-5dB) noise environment. The results show that both methods can well estimate out the frequency spectrum of the signal x(n). Figure 14(a) shows that the traditional ARMA model spectrum method completely fails under SαS (α=1.3) noise environment (GSNR=20dB), but the improved FLO-ARMA model spectrum method in Figure 14(b) has very good toughness. Hence, the ARMA spectrum method is only suitable to estimate the frequency spectrum of the signals under Gaussian environment, but the FLO-ARMA model method can be applied in Gaussian and SαS noise environment, which can work in lower GSNR.The traditional parameter model spectrum estimation methods are equivalent to solving Yule-Walker equation and their extension equations based on the autocorrelation function. However, the fractional low order parameter model spectrum estimation methods construct the generalized Yule-Walker equation and their extension equations employing fractional low order covariance matrix, because fractional low order autocovariance is a generalized covariance, which changes into autocorrelation whena=b=2. As a result, the fractional low order parameter model spectrum estimation methods can be applied to both Gaussian and stable distribution noise environment.


### 3.4. Spectral Resolution of the FLO Spectrum Estimation Methods

In order to further verify the spectral resolution capability of the improved fractional low order spectrum methods, in this section, we apply the FLOBTCS, FLOPCS, FLO-Welch, FLOAR, and FLOARMA model spectrum methods to estimate the spectrum of the signal with closely spaced frequencies. The signaly(n) in (21) is used as the experimental test signal, where f0=100, f1=110, f2=500, n=0,1,2,…,N-1, and N=2000. v(n) is SαS stable distribution noise, α=1.4, and GSNR=22dB. The simulations of ten independent experiments are shown in Figure 15.Figure 15
Frequency spectrum estimations of the signal with closely spaced frequencies employing fractional lower order frequency spectrum methods ((a) fractional lower order Blackman-Tukey method; (b) fractional lower order Periodogram method; (c) fractional lower order Welch covariation method; (d) fractional lower order AR model method; (e) fractional lower order ARMA model method).
(a) (b) (c) (d) (e)Figures15(a) and 15(b) and Figure 15(c) are the FLOBTCS, FLOPCS, and FLO-Welch spectrums of the signal x(n), respectively. The FLOAR(90) model spectrum is shown in Figure 15(d), and Figure 15(e) is FLOAR(65, 3) model spectrum of the signal x(n). The results show that the improved methods can distinguish the closely spaced frequencies well. Hence, the proposed methods have better performance of spectral estimation and higher spectral resolution, and the FLOARMA model spectrum method performs best.


## 3.1. Fractional Lower Order Covariation Spectrum Estimation Method

### 3.1.1. Covariation Spectrum Principle

The autocorrelation function and power spectrum of a second-order random process can be expressed as(2)Rxxτ=Extx∗t-τ(3)Sxxω=∫-∞∞Rxxτe-jωτdτ(4)Rxxτ=12π∫-ππSxxωejωτdωFrom (2)-(4), we can see that power spectrum of x(t) is Fourier transform of its correlation function, which is called Wiener-Khinchin theorem. Rxx(τ) demonstrates second-order statistical properties of x(t) in time domain, and Sxx(ω) is the statistical properties of x(t) in frequency domain. The corresponding discrete correlation function and power spectrum are given by(5)R^xxm=1N∑n=0N-1xnx∗n-m(6)S^xxω=∑m=0N-1R^xxmejωm=1N∑m=0N-1∑n=0N-1xnx∗n-me-jωmWhenx(t) is α stable distribution process, the autocorrelation function and power spectrum based on second-order statistics are no longer applicable. Hence, fractional lower order covariation was proposed in [21], which was obtained by calculating fractional low order moments.(7)X,Yα=EXY<P-1>EYpγywhere p is a given real constant, 1≤p<α≤2 [19–22], and p denotes p order moment of x(t). When x(t) is a real signal, xp(t)=x(t)p·sign[x(t)], signx(t)=1x(t)>00x(t)=0-1x(t)<0, and when x(t) is a complex signal, xp(t)=x(t)p-1·x∗(t), ∗ denotes conjugate operation. γy is dispersion coefficient of Y. According to (7), the autocovariation function of x(t) can be expressed as (8)Cxxτ=xt,xt-τα=Extxt-τ<p-1>Ext-τpγxt-τ,1≤p<α≤2When x(t) is Gaussian process, then α=2 and p=2, and (8) degenerates into autocorrelation function(9)Rxxτ=xt,xt-τα=xt,xt-τ2=Extxt-τWe take Fourier transform on both sides of (8) and get(10)φxxω=∫-∞∞Cxxτe-jωτdτ(11)Cxxτ=12π∫-ππφxxωejωτdωEquation (10) is known as fractional low order covariation spectrum under α stable distribution process, which is corresponding to power spectrum under Gaussian distribution. When α=2 and p=2, (10) changes as power spectrum; hence, fractional low order covariation spectrum is a generalized power spectrum.


### 3.1.2. Fractional Lower Order Blackman-Tukey Covariation Spectrum Estimation

We assume thatx(n),n=0,1,⋯,N-1 is N samples of a discrete α stable distribution random process, and the sample mean is 0, then(12)C^xxm=1N∑n=0N-1xnx<p-1>n-m(13)φ^xxω=∑m=-MMC^xxmejωm=1N∑m=-MM∑n=0N-1xnx<p-1>n-me-jωmEquation (12) can be applied to estimate the covariation function of x(n). By taking discrete Fourier transform of (12), Blackman-Tukey covariation spectrum estimation can be obtained, as shown in (13).


### 3.1.3. Fractional Lower Order Periodogram Covariation Spectrum Estimation

LettingM=N-1 and k=n-m, substituting them into (13),(14)
ϕ
^
x
x
ω
=
1
N
∑
m
∑
k
x
n
x
<
p
-
1
>
k
e
-
j
ω
m
=
1
N
∑
n
=
0
N
-
1
x
n
e
-
j
ω
n
∑
k
=
0
N
-
1
x
<
p
-
1
>
k
e
j
ω
k
=
1
N
∑
n
=
0
N
-
1
x
n
e
-
j
ω
n
∑
k
=
0
N
-
1
x
<
p
-
1
>
k
e
-
j
ω
k
∗
=
1
N
∑
n
=
0
N
-
1
x
n
e
-
j
ω
n
∑
n
=
0
N
-
1
x
<
p
-
1
>
n
e
-
j
ω
n
∗
=
1
N
X
n
X
-
n
∗where X(ω)=∑n=0N-1x(n)e-jωn is discrete Fourier transform of x(n). The discrete Fourier transform of x<p-1>(n) is given by X-(n)=∑n=0N-1x<p-1>(n)e-jωn. From (14), we firstly compute discrete Fourier transform of x(n) and x<p-1>(n) and take their product; then fractional lower order Periodogram covariation spectrum estimation can be obtained.


### 3.1.4. Fractional Lower Order Welch Covariation Spectrum Estimation

Covariation spectrums of the fractional lower order Blackman-Tukey and fractional lower order Welch methods can be gotten by Fourier transform, the principle is simple, and the calculation quantity is also small. However, FLOBTCS and FLOPCS have higher side lobe, and the variance is larger, which results in poor performance and affects their practical application. Therefore, the traditional Welch method needs to be improved.By dividingN discrete signal x(n),n=0,1,⋯,N-1 into partially overlapped K segments, the length of each segment is L, the overlap length between two adjacent segments is L/2, and then the ind segment can be expressed as(15)xin=xn+i-1L,i=1,2,⋯,K,n=0,1,⋯,L-1(16)K=N-L/2L/2By taking the windowed Fourier transform of x(n) and x<p-1>(n), the following equations are obtained:(17)Xin=∑n=0L-1xinhne-jωn,i=1,2,⋯,K(18)X-in=∑n=0N-1xi<p-1>nhne-jωn,i=1,2,⋯,KSubstituting (17)-(18) into (14), the ind covariation spectrum can be written as(19)φ^xiω=1NXinX-in∗,i=1,2,⋯,KCalculating the average of the covariation spectrum of all segments, then(20)φ^xxω=1K∑i=1Kφ^xiω=1K∑i=1K1NXinX-in∗,i=1,2,⋯,KEquation (20) is called fractional low order Welch covariation spectrum method. In order to reduce variance of the covariation spectrum estimation, the number of segments K can be increased appropriately, which is bound to reduce length of the segment L, so that the main lobe of the covariation spectrum is widened to reduce the resolution. In real application, K and L can be set according to the specific needs.


### 3.1.5. Application Review

We will now apply the FLOBTCS and FLOPCS and FLO-Welch covariation spectrum methods to estimate frequency spectrum of the test signal. The test signaly(n) is defined as three sinusoidal signals with Gaussian or SαS stable distribution noise.(21)yn=Asin⁡2πf0n+sin⁡2πf1n+sin⁡2πf2n+vn=Axn+vnwhere f0=100, f1=300, f2=500, n=0,1,2,....,N-1, and N=2000. When v(n) is Gaussian noise, generalized signal to noise ratio (GSNR) is used to instead of SNR.(22)GSNR=10log10⁡Exn2γα==10log10⁡1Nγα∑n=0N-1xn2where γ is the dispersion coefficient of SαS stable distribution noise. According to the given GSNR, the amplitude of the signal can be written as(23)A=10GSNR/101/N∑n=0N-1xn2γα1/2LettingGSNR=20dB, the fractional lower order covariation spectrum methods and traditional power spectrum methods are used to estimate frequency spectrum of the signal x(n) under Gaussian distribution noise and SαS stable distribution noise; the simulation results of ten independent experiments are shown in Figures 5–8.Figure 5
Fractional lower order Blackman-Tukey covariation spectrums and power spectrums of the signalx(n) under Gaussian and SαS noise environment ((a) Gaussian noise environment, SNR=-5dB; (b) SαS noise environment, α=1.3, and GSNR=20dB).
(a) (b)Figure 6
Fractional lower order Periodogram covariation spectrums and power spectrums of the signalx(n) under Gaussian and SαS noise environment ((a) Gaussian noise environment, SNR=-5dB; (b) SαS noise environment, α=1.3, and GSNR=20dB).
(a) (b)Figure 7
Fractional lower order Welch covariation spectrums and power spectrums of the signalx(n) under Gaussian and SαS noise environment ((a) Gaussian noise environment, SNR=-5dB; (b) SαS noise environment, α=1.3, and GSNR=20dB).
(a) (b)Figure 8
Comparisons of fractional lower order Blackman-Tukey, Periodogram, and Welch covariation spectrums of the signalx(n) ((a) SαS (α=1.1) noise environment; (b) SαS (α=1.8) noise environment).
(a) (b)


### 3.1.6. Remarks

Figure5(a) shows that the traditional Blackman-Tukey spectrum estimation method fails under SαS noise environment, but the improved FLOBTCS method in Figure 5(b) has good toughness. Figures 6(a) and 6(b) are frequency estimations of the signal x(n) employing the Periodogram and FLOPCS methods, respectively; the result reveals the advantage of the FLOPCS method. The Welch spectrum and FLO-Welch covariation spectrum of the signal x(n) are given in Figure 7; we can know that the Welch power spectrum method cannot work under SαS noise environment, but the proposed FLO-Welch covariation method demonstrates good performance. Figure 8 is frequency spectrum of the signal x(n) employing FLOBTCS, FLOPCS, and FLO-Welch covariation spectrum methods under SαS (α=1.1 and α=1.8) noise environment. The results show that the conventional power spectrum methods degenerate under SαS noise environment, but the improved fractional low order covariation spectrum methods can better estimate out the frequency spectrum of the signal x(n).The FLOBTCS and FLOPCS methods have a certain side lobe and little amount of calculation, which cause the weak signals to be submerged and lead to the distortion of the signals. FLO-Welch algorithm can suppress the side lobe effects by spectrum segmentation. However, the resolution and variance of the FLO-Welch method are controlled by the number of segments and length like the conventional Welch method. In applications, the shorter number of segments should be used when we would like to get higher frequency resolution, but variance gets larger, and if we wish to reduce estimation variance, a longer number of segments are preferred. Hence, we can compromise the number and length of the segments by the actual situation. As a result, the power spectrum methods are only suitable to analyze the signals under Gaussian environment, but the fractional low order covariation spectrum methods can be applied under Gaussian andSαS noise environment.


## 3.1.1. Covariation Spectrum Principle

The autocorrelation function and power spectrum of a second-order random process can be expressed as(2)Rxxτ=Extx∗t-τ(3)Sxxω=∫-∞∞Rxxτe-jωτdτ(4)Rxxτ=12π∫-ππSxxωejωτdωFrom (2)-(4), we can see that power spectrum of x(t) is Fourier transform of its correlation function, which is called Wiener-Khinchin theorem. Rxx(τ) demonstrates second-order statistical properties of x(t) in time domain, and Sxx(ω) is the statistical properties of x(t) in frequency domain. The corresponding discrete correlation function and power spectrum are given by(5)R^xxm=1N∑n=0N-1xnx∗n-m(6)S^xxω=∑m=0N-1R^xxmejωm=1N∑m=0N-1∑n=0N-1xnx∗n-me-jωmWhenx(t) is α stable distribution process, the autocorrelation function and power spectrum based on second-order statistics are no longer applicable. Hence, fractional lower order covariation was proposed in [21], which was obtained by calculating fractional low order moments.(7)X,Yα=EXY<P-1>EYpγywhere p is a given real constant, 1≤p<α≤2 [19–22], and p denotes p order moment of x(t). When x(t) is a real signal, xp(t)=x(t)p·sign[x(t)], signx(t)=1x(t)>00x(t)=0-1x(t)<0, and when x(t) is a complex signal, xp(t)=x(t)p-1·x∗(t), ∗ denotes conjugate operation. γy is dispersion coefficient of Y. According to (7), the autocovariation function of x(t) can be expressed as (8)Cxxτ=xt,xt-τα=Extxt-τ<p-1>Ext-τpγxt-τ,1≤p<α≤2When x(t) is Gaussian process, then α=2 and p=2, and (8) degenerates into autocorrelation function(9)Rxxτ=xt,xt-τα=xt,xt-τ2=Extxt-τWe take Fourier transform on both sides of (8) and get(10)φxxω=∫-∞∞Cxxτe-jωτdτ(11)Cxxτ=12π∫-ππφxxωejωτdωEquation (10) is known as fractional low order covariation spectrum under α stable distribution process, which is corresponding to power spectrum under Gaussian distribution. When α=2 and p=2, (10) changes as power spectrum; hence, fractional low order covariation spectrum is a generalized power spectrum.


## 3.1.2. Fractional Lower Order Blackman-Tukey Covariation Spectrum Estimation

We assume thatx(n),n=0,1,⋯,N-1 is N samples of a discrete α stable distribution random process, and the sample mean is 0, then(12)C^xxm=1N∑n=0N-1xnx<p-1>n-m(13)φ^xxω=∑m=-MMC^xxmejωm=1N∑m=-MM∑n=0N-1xnx<p-1>n-me-jωmEquation (12) can be applied to estimate the covariation function of x(n). By taking discrete Fourier transform of (12), Blackman-Tukey covariation spectrum estimation can be obtained, as shown in (13).


## 3.1.3. Fractional Lower Order Periodogram Covariation Spectrum Estimation

LettingM=N-1 and k=n-m, substituting them into (13),(14)
ϕ
^
x
x
ω
=
1
N
∑
m
∑
k
x
n
x
<
p
-
1
>
k
e
-
j
ω
m
=
1
N
∑
n
=
0
N
-
1
x
n
e
-
j
ω
n
∑
k
=
0
N
-
1
x
<
p
-
1
>
k
e
j
ω
k
=
1
N
∑
n
=
0
N
-
1
x
n
e
-
j
ω
n
∑
k
=
0
N
-
1
x
<
p
-
1
>
k
e
-
j
ω
k
∗
=
1
N
∑
n
=
0
N
-
1
x
n
e
-
j
ω
n
∑
n
=
0
N
-
1
x
<
p
-
1
>
n
e
-
j
ω
n
∗
=
1
N
X
n
X
-
n
∗where X(ω)=∑n=0N-1x(n)e-jωn is discrete Fourier transform of x(n). The discrete Fourier transform of x<p-1>(n) is given by X-(n)=∑n=0N-1x<p-1>(n)e-jωn. From (14), we firstly compute discrete Fourier transform of x(n) and x<p-1>(n) and take their product; then fractional lower order Periodogram covariation spectrum estimation can be obtained.


## 3.1.4. Fractional Lower Order Welch Covariation Spectrum Estimation

Covariation spectrums of the fractional lower order Blackman-Tukey and fractional lower order Welch methods can be gotten by Fourier transform, the principle is simple, and the calculation quantity is also small. However, FLOBTCS and FLOPCS have higher side lobe, and the variance is larger, which results in poor performance and affects their practical application. Therefore, the traditional Welch method needs to be improved.By dividingN discrete signal x(n),n=0,1,⋯,N-1 into partially overlapped K segments, the length of each segment is L, the overlap length between two adjacent segments is L/2, and then the ind segment can be expressed as(15)xin=xn+i-1L,i=1,2,⋯,K,n=0,1,⋯,L-1(16)K=N-L/2L/2By taking the windowed Fourier transform of x(n) and x<p-1>(n), the following equations are obtained:(17)Xin=∑n=0L-1xinhne-jωn,i=1,2,⋯,K(18)X-in=∑n=0N-1xi<p-1>nhne-jωn,i=1,2,⋯,KSubstituting (17)-(18) into (14), the ind covariation spectrum can be written as(19)φ^xiω=1NXinX-in∗,i=1,2,⋯,KCalculating the average of the covariation spectrum of all segments, then(20)φ^xxω=1K∑i=1Kφ^xiω=1K∑i=1K1NXinX-in∗,i=1,2,⋯,KEquation (20) is called fractional low order Welch covariation spectrum method. In order to reduce variance of the covariation spectrum estimation, the number of segments K can be increased appropriately, which is bound to reduce length of the segment L, so that the main lobe of the covariation spectrum is widened to reduce the resolution. In real application, K and L can be set according to the specific needs.


## 3.1.5. Application Review

We will now apply the FLOBTCS and FLOPCS and FLO-Welch covariation spectrum methods to estimate frequency spectrum of the test signal. The test signaly(n) is defined as three sinusoidal signals with Gaussian or SαS stable distribution noise.(21)yn=Asin⁡2πf0n+sin⁡2πf1n+sin⁡2πf2n+vn=Axn+vnwhere f0=100, f1=300, f2=500, n=0,1,2,....,N-1, and N=2000. When v(n) is Gaussian noise, generalized signal to noise ratio (GSNR) is used to instead of SNR.(22)GSNR=10log10⁡Exn2γα==10log10⁡1Nγα∑n=0N-1xn2where γ is the dispersion coefficient of SαS stable distribution noise. According to the given GSNR, the amplitude of the signal can be written as(23)A=10GSNR/101/N∑n=0N-1xn2γα1/2LettingGSNR=20dB, the fractional lower order covariation spectrum methods and traditional power spectrum methods are used to estimate frequency spectrum of the signal x(n) under Gaussian distribution noise and SαS stable distribution noise; the simulation results of ten independent experiments are shown in Figures 5–8.Figure 5
Fractional lower order Blackman-Tukey covariation spectrums and power spectrums of the signalx(n) under Gaussian and SαS noise environment ((a) Gaussian noise environment, SNR=-5dB; (b) SαS noise environment, α=1.3, and GSNR=20dB).
(a) (b)Figure 6
Fractional lower order Periodogram covariation spectrums and power spectrums of the signalx(n) under Gaussian and SαS noise environment ((a) Gaussian noise environment, SNR=-5dB; (b) SαS noise environment, α=1.3, and GSNR=20dB).
(a) (b)Figure 7
Fractional lower order Welch covariation spectrums and power spectrums of the signalx(n) under Gaussian and SαS noise environment ((a) Gaussian noise environment, SNR=-5dB; (b) SαS noise environment, α=1.3, and GSNR=20dB).
(a) (b)Figure 8
Comparisons of fractional lower order Blackman-Tukey, Periodogram, and Welch covariation spectrums of the signalx(n) ((a) SαS (α=1.1) noise environment; (b) SαS (α=1.8) noise environment).
(a) (b)


## 3.1.6. Remarks

Figure5(a) shows that the traditional Blackman-Tukey spectrum estimation method fails under SαS noise environment, but the improved FLOBTCS method in Figure 5(b) has good toughness. Figures 6(a) and 6(b) are frequency estimations of the signal x(n) employing the Periodogram and FLOPCS methods, respectively; the result reveals the advantage of the FLOPCS method. The Welch spectrum and FLO-Welch covariation spectrum of the signal x(n) are given in Figure 7; we can know that the Welch power spectrum method cannot work under SαS noise environment, but the proposed FLO-Welch covariation method demonstrates good performance. Figure 8 is frequency spectrum of the signal x(n) employing FLOBTCS, FLOPCS, and FLO-Welch covariation spectrum methods under SαS (α=1.1 and α=1.8) noise environment. The results show that the conventional power spectrum methods degenerate under SαS noise environment, but the improved fractional low order covariation spectrum methods can better estimate out the frequency spectrum of the signal x(n).The FLOBTCS and FLOPCS methods have a certain side lobe and little amount of calculation, which cause the weak signals to be submerged and lead to the distortion of the signals. FLO-Welch algorithm can suppress the side lobe effects by spectrum segmentation. However, the resolution and variance of the FLO-Welch method are controlled by the number of segments and length like the conventional Welch method. In applications, the shorter number of segments should be used when we would like to get higher frequency resolution, but variance gets larger, and if we wish to reduce estimation variance, a longer number of segments are preferred. Hence, we can compromise the number and length of the segments by the actual situation. As a result, the power spectrum methods are only suitable to analyze the signals under Gaussian environment, but the fractional low order covariation spectrum methods can be applied under Gaussian andSαS noise environment.


## 3.2. FLO-AR Model Spectrum Estimation Method

### 3.2.1. Principle

An independent identically distributed (i.i.d)α process can be expressed by p order AR model(24)xn=a1xn-1+⋯+apxn-p+un=∑i=1paixn-i+unwhere ai(i=1,2,⋯,p) is the parameters of AR model. u(n) is an i.i.d process, and its characteristic index and dispersion coefficient are α and γu, respectively. When u(n) passes a finite impulse response (FIR) filter, then(25)xn=∑i=0qhiun-iTaking z transformation for (24)-(25),(26)Hz=∑i=0qhiz-i=11-∑i=1paiz-iwhere H(z) is transfer function of the filter. α spectrum of the signal x(n) was defined as [21](27)Sαz=Xn,Wznα=γuH1zα-1Hzα-1where(28)Wzn=∑i=-qi=qXn-iziWhen z=ejω, substituting z=1 and (26) into (27), AR model α spectrum of an α stable distribution process on the unit circle can be written as(29)Sαejω=γuH1ejωα-1Hejωα-1=γuHejωHejωα-1=γuHejωα=γu11-∑i=1paie-jωiαThen AR model α spectrum of a SαS (γu=1) process on the unit circle is given by(30)SSαSejω=Hejωα=11-∑i=1paie-jωiαAccording to the definition in [21], we can obtain fractional low order autocovariance of the signal x(n).(31)Dxxm=ΔExn<a>xn-m<b>,0≤a<α2,0≤b<α2The corresponding low order autocovariance estimation can be written as(32)D^xxm=∑n=L1+1L2xnaxn-mbsignxnxn-mL2-L1where L1=max⁡(0,-m), L2=min⁡(N-m,N), and N is number of the observation signal x(n).We defineA=[a^1,a^2,⋯,a^p]T, B=[D^xx(1),D^xx(2),⋯,D^xx(P)]T, and(33)Τ=D^xx0D^xx-1⋯D^xx1-pD^xx1D^xx0⋯D^xx2-p⋮⋮⋱⋮D^xxp-1D^xxp-2⋯D^xx0where Τ is call fractional low order autocovariance matrix (FLOACM). We apply FLOACM to fractional low order moment (FLOM) matrix in [21] and then get(34)ΤA=BAnd(35)D^xx0D^xx-1⋯D^xx1-pD^xx1D^xx0⋯D^xx2-p⋮⋮⋱⋮D^xxp-1D^xxp-2⋯D^xx0a^1a^2⋮a^p=D^xx1D^xx2⋮D^xxpWe call (35) fractional low order Euler Walker equation. The parameters a^i(i=1,2,⋯,p) can be gotten by solving (35). When a^i(i=1,2,⋯,p) are substituted into (26) or (29), fractional low order AR model spectrum can be obtained. 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. The calculation formula can be written as(36)FPEp=σa2N+pN-pwhere σa2 is the variance of residuals.We summarize the steps of fractional low order AR model spectrum method, as follows.Step 1.
Solve fractional low order AR model parametersai(i=1,2,⋯,p) of x(n) employing (35).Step 2.
Determine the orderp with the FPE criterion.Step 3.
Compute fractional low order AR model spectrum ofx(n) by substituting ai(i=1,2,⋯,p) into (29) or (30).


### 3.2.2. Application Review

In this simulation, the test signal isx(n) in (21). We compare the performance of the traditional AR model method and the improved fractional low order AR model method under Gaussian distribution noise (SNR=-5dB) and SαS stable distribution noise (α=1.3, GSNR=20dB); the results of the ten independent experiments are shown in Figures 9-10. In order to further verify the effectiveness of the fractional low order AR method, we conduct comparative experiment on two methods under different α when GSNR=20, as show in Figure 11. When α=1.3 and GSNR changes from 14dB to 24dB, Figure 12 shows the change of Errors Power of the exiting AR model method and the improved FLO-AR method under α stable distribution noise environment.Figure 9
The traditional AR model spectrums and improved FLO-AR spectrums of the signalx(n) under Gaussian noise environment (SNR=-5dB) ((a) AR model spectrums of the signal x(n), where p=29; (b) FLO-AR spectrums of the signal x(n), where p=27 and a=b=0.8).
(a) (b)Figure 10
FLO-AR model and AR model spectrums of the signalx(n)under SαS (α=1.3) noise environment (GSNR=20dB) ((a) AR model spectrums of the signal x(n), where p=20; (b) FLO-AR model spectrums of the signal x(n), where p=21 and a=b=0.1).
(a) (b)Figure 11
Errors power of the AR model and FLO-AR model frequency spectrum algorithms under differentα (GSNR=20).Figure 12
Errors power of the AR model and FLO-AR model frequency spectrum algorithms under differentGSNR (α=1.3).


### 3.2.3. Remarks

Figures9(a) and 9(b) show the traditional AR model spectrum and the improved FLO-AR model spectrum of the signal x(n), respectively. Both methods can estimate out the frequency spectrum of the signal x(n) well under Gaussian noise environment. Figure 10(a) reveals that the AR model spectrum method fails under SαS noise environment (GSNR=20dB), but the FLO-AR model spectrum method in Figure 10(b) has good toughness. As a result, the AR model spectrum method is only suitable to analyze the signals under Gaussian environment, but the FLO-AR model spectrum method can work under Gaussian and SαS noise environment, which is robust.Figure11 is Errors Power of the AR and FLO-AR model spectrum algorithms under different α. When the characteristic index α of SαS distribution noise changes from 0 to 2, Errors Power of the improved FLO-AR method is lower than the AR algorithm, which is very stable.WhenGSNR changes from 14dB to 24dB, the Errors Power of the AR model and FLO-AR model methods are given in Figure 12. The results show that Errors Power of the improved FLO-AR method is lower than that of the AR method, and they are very stable with the changing of GSNR. Especially, when GSNR<19dB, the performance superiority of the FLO-AR method is more obvious.


## 3.2.1. Principle

An independent identically distributed (i.i.d)α process can be expressed by p order AR model(24)xn=a1xn-1+⋯+apxn-p+un=∑i=1paixn-i+unwhere ai(i=1,2,⋯,p) is the parameters of AR model. u(n) is an i.i.d process, and its characteristic index and dispersion coefficient are α and γu, respectively. When u(n) passes a finite impulse response (FIR) filter, then(25)xn=∑i=0qhiun-iTaking z transformation for (24)-(25),(26)Hz=∑i=0qhiz-i=11-∑i=1paiz-iwhere H(z) is transfer function of the filter. α spectrum of the signal x(n) was defined as [21](27)Sαz=Xn,Wznα=γuH1zα-1Hzα-1where(28)Wzn=∑i=-qi=qXn-iziWhen z=ejω, substituting z=1 and (26) into (27), AR model α spectrum of an α stable distribution process on the unit circle can be written as(29)Sαejω=γuH1ejωα-1Hejωα-1=γuHejωHejωα-1=γuHejωα=γu11-∑i=1paie-jωiαThen AR model α spectrum of a SαS (γu=1) process on the unit circle is given by(30)SSαSejω=Hejωα=11-∑i=1paie-jωiαAccording to the definition in [21], we can obtain fractional low order autocovariance of the signal x(n).(31)Dxxm=ΔExn<a>xn-m<b>,0≤a<α2,0≤b<α2The corresponding low order autocovariance estimation can be written as(32)D^xxm=∑n=L1+1L2xnaxn-mbsignxnxn-mL2-L1where L1=max⁡(0,-m), L2=min⁡(N-m,N), and N is number of the observation signal x(n).We defineA=[a^1,a^2,⋯,a^p]T, B=[D^xx(1),D^xx(2),⋯,D^xx(P)]T, and(33)Τ=D^xx0D^xx-1⋯D^xx1-pD^xx1D^xx0⋯D^xx2-p⋮⋮⋱⋮D^xxp-1D^xxp-2⋯D^xx0where Τ is call fractional low order autocovariance matrix (FLOACM). We apply FLOACM to fractional low order moment (FLOM) matrix in [21] and then get(34)ΤA=BAnd(35)D^xx0D^xx-1⋯D^xx1-pD^xx1D^xx0⋯D^xx2-p⋮⋮⋱⋮D^xxp-1D^xxp-2⋯D^xx0a^1a^2⋮a^p=D^xx1D^xx2⋮D^xxpWe call (35) fractional low order Euler Walker equation. The parameters a^i(i=1,2,⋯,p) can be gotten by solving (35). When a^i(i=1,2,⋯,p) are substituted into (26) or (29), fractional low order AR model spectrum can be obtained. 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. The calculation formula can be written as(36)FPEp=σa2N+pN-pwhere σa2 is the variance of residuals.We summarize the steps of fractional low order AR model spectrum method, as follows.Step 1.
Solve fractional low order AR model parametersai(i=1,2,⋯,p) of x(n) employing (35).Step 2.
Determine the orderp with the FPE criterion.Step 3.
Compute fractional low order AR model spectrum ofx(n) by substituting ai(i=1,2,⋯,p) into (29) or (30).


## 3.2.2. Application Review

In this simulation, the test signal isx(n) in (21). We compare the performance of the traditional AR model method and the improved fractional low order AR model method under Gaussian distribution noise (SNR=-5dB) and SαS stable distribution noise (α=1.3, GSNR=20dB); the results of the ten independent experiments are shown in Figures 9-10. In order to further verify the effectiveness of the fractional low order AR method, we conduct comparative experiment on two methods under different α when GSNR=20, as show in Figure 11. When α=1.3 and GSNR changes from 14dB to 24dB, Figure 12 shows the change of Errors Power of the exiting AR model method and the improved FLO-AR method under α stable distribution noise environment.Figure 9
The traditional AR model spectrums and improved FLO-AR spectrums of the signalx(n) under Gaussian noise environment (SNR=-5dB) ((a) AR model spectrums of the signal x(n), where p=29; (b) FLO-AR spectrums of the signal x(n), where p=27 and a=b=0.8).
(a) (b)Figure 10
FLO-AR model and AR model spectrums of the signalx(n)under SαS (α=1.3) noise environment (GSNR=20dB) ((a) AR model spectrums of the signal x(n), where p=20; (b) FLO-AR model spectrums of the signal x(n), where p=21 and a=b=0.1).
(a) (b)Figure 11
Errors power of the AR model and FLO-AR model frequency spectrum algorithms under differentα (GSNR=20).Figure 12
Errors power of the AR model and FLO-AR model frequency spectrum algorithms under differentGSNR (α=1.3).


## 3.2.3. Remarks

Figures9(a) and 9(b) show the traditional AR model spectrum and the improved FLO-AR model spectrum of the signal x(n), respectively. Both methods can estimate out the frequency spectrum of the signal x(n) well under Gaussian noise environment. Figure 10(a) reveals that the AR model spectrum method fails under SαS noise environment (GSNR=20dB), but the FLO-AR model spectrum method in Figure 10(b) has good toughness. As a result, the AR model spectrum method is only suitable to analyze the signals under Gaussian environment, but the FLO-AR model spectrum method can work under Gaussian and SαS noise environment, which is robust.Figure11 is Errors Power of the AR and FLO-AR model spectrum algorithms under different α. When the characteristic index α of SαS distribution noise changes from 0 to 2, Errors Power of the improved FLO-AR method is lower than the AR algorithm, which is very stable.WhenGSNR changes from 14dB to 24dB, the Errors Power of the AR model and FLO-AR model methods are given in Figure 12. The results show that Errors Power of the improved FLO-AR method is lower than that of the AR method, and they are very stable with the changing of GSNR. Especially, when GSNR<19dB, the performance superiority of the FLO-AR method is more obvious.


## 3.3. FLO-ARMA Model Spectrum Estimation Method

### 3.3.1. Principle

An i.i.dα random process can be expressed by an ARMA (p,q).(37)xn=a1xn-1+⋯+apxn-p+b0un+b1un-1+⋯+bqun-q=∑i=1paixn-i+∑i=0qbiun-iwhere ai(i=1,2,⋯,p) and bi(i=0,1,2,⋯,q) are parameters of p order AR model and q order MA model, respectively. Taking Z transform on (37), obtaining(38)-∑i=0paixn-i=∑i=0qbiun-i(39)-∑i=0paiz-iXz=∑i=0qbiz-iUz(40)1-∑i=1paiz-iXz=1+∑i=1qbiz-iUz(41)Hz=XzUz=1+∑i=1qbiz-i1-∑i=1paiz-i=BzAz(42)hn↔Hz,xn=hn∗un=∑i=0qhiun-iwhere h(n) is impulse response function, H(z) is system transfer function, A(z) and B(z) are the FIR and IIR filters, respectively. From (41), FLO-ARMA model spectrum of α stable distribution process on the unit circle (z=ejω) is given by(43)Sαejω=γuHejωα==γu1+∑i=1qbie-jωi1-∑i=1paie-jωiαThen FLO-ARMA model spectrum of SαS (γu=1) process on the unit circle can be written as(44)SSαSejω=Hejωα=1+∑i=1qbie-jωi1-∑i=1pake-jωiαMultiplyingx(n-m) both sides of (38) and computing fractional low order covariance, then(45)Dxxm=∑i=1paiDxxm-i+∑i=0qbiDuxm-iwhere(46)Dxxm=Exn<a>xn-m<b>Dxxm-i=ΔExn<a>xn-m+i<b>Duxm-i=ΔEun<a>xn-m+i<b>Then we can get [20](47)Duxm=0m>0γuh-mα-1m<0When m>q, (45) changes as(48)Dxxm=∑i=1paiDxxm-i,m>qEquation (48) can be written as(49)DqDq-1⋯Dq-p+1Dq+1Dq⋯Dq-p+2⋮⋮⋱⋮Dq+p-1Dq+p-2⋯Dqa1a2⋮ap=Dq+1Dq+2⋮Dq+pwhere (49) is also a generalized Yule-Walker equation. Then we can obtain the parameters of the FLO-AR model by solving (49).When0≤m≤q, (45) changes as(50)Dxxm=∑i=0paiDxxm-i+∑i=0qbiDuxm-i=∑i=1paiDxxm-i+∑i=0qbihk-iα-1,0≤m≤qLetting p=0, then(51)Dxxm=γu∑i=0p-mbi+mbiα-1,0≤m≤qγu∑i=0p-mbibi-mα-1,-q≤m≤0By solving the nonlinear equation (51), the parameters bi(i=1,2,⋯,q) of the FLO-MA model can be obtained.An approximate infinitep′ order FLO-AR model can be used to be equivalent to a finite q′ order MA model (p′≫q′), then(52)am+∑i=1q′biam-i=emwhere e(m) is error. Letting e(m)=0, then(53)∑i=1q′biam-i=-amBy taking autocorrelation of both sides of (53) with a∗(n-m) one has(54)∑i=1q′biRm-i=-Rmwhere R(m)=Ea(m)a∗(n-m) and R(m-i)=Ea(m-i)a∗(n-m+i). Equation (38) can be written as(55)R0R-1⋯R1-q′R1R0⋯R2-q′⋮⋮⋱⋮Rq′-1Rq′-2⋯R0b1b2⋮bp=-R1R2⋮Rq′Equation (55) is a Toeplitz matrix. By solving the Yule-Walker equation (55), the parameters bi(i=1,2,⋯,q) of the FLO-MA model are gotten.We summarize the steps of the FLO-ARMA model spectrum method, as follows.Step 1.
Compute the parametersai(i=1,2,⋯,p) of the FLO-AR model of the signal x(n) employing (49).Step 2.
By solving the nonlinear equation (51) or the Yule-Walker equation (55), the parameters bi(i=1,2,⋯,q) of the FLO-MA model can be gotten.Step 3.
Determine the orderp and q employing FPE criterion.Step 4.
Substituteai(i=1,2,⋯,p) and bi(i=1,2,⋯,q) into (29) or (30); compute FLO-ARMA model spectrum of the signal x(n).


### 3.3.2. Application Review

In this simulation,x(n) in (21) is used as the experimental test signal. The traditional ARMA method and modified FLO-ARMA method are applied to estimate the frequency spectrum of the signal x(n) under Gaussian noise (SNR=-5dB) and SαS noise environment (GSNR=20dB). The results of ten independent experiments are shown in Figures 13-14.Figure 13
The ARMA model and FLO-ARMA model spectrums of the signalx(n) under Gaussian noise environment (SNR=-5dB) ((a) ARMA model spectrums of the signal x(n), where p=29 and q=5; (b) FLO-ARMA model spectrums of the signal x(n), where p=27, q=6, and a=b=0.9).
(a) (b)Figure 14
FLO-ARMA model and ARMA model spectrums of the signalx(n) under SαS (α=1.3) noise environment (GSNR=20dB) ((a) ARMA model spectrums of the signal x(n), where p=20 and q=5; (b) FLO-ARMA model spectrums of the signal x(n), where p=21, q=5, and a=b=0.1).
(a) (b)


### 3.3.3. Remarks

The ARMA frequency spectrum of the signalx(n) is shown in Figure 13(a), and Figure 13(b) is FLO-ARMA frequency spectrum of the signal x(n) under Gaussian(SNR=-5dB) noise environment. The results show that both methods can well estimate out the frequency spectrum of the signal x(n). Figure 14(a) shows that the traditional ARMA model spectrum method completely fails under SαS (α=1.3) noise environment (GSNR=20dB), but the improved FLO-ARMA model spectrum method in Figure 14(b) has very good toughness. Hence, the ARMA spectrum method is only suitable to estimate the frequency spectrum of the signals under Gaussian environment, but the FLO-ARMA model method can be applied in Gaussian and SαS noise environment, which can work in lower GSNR.The traditional parameter model spectrum estimation methods are equivalent to solving Yule-Walker equation and their extension equations based on the autocorrelation function. However, the fractional low order parameter model spectrum estimation methods construct the generalized Yule-Walker equation and their extension equations employing fractional low order covariance matrix, because fractional low order autocovariance is a generalized covariance, which changes into autocorrelation whena=b=2. As a result, the fractional low order parameter model spectrum estimation methods can be applied to both Gaussian and stable distribution noise environment.


## 3.3.1. Principle

An i.i.dα random process can be expressed by an ARMA (p,q).(37)xn=a1xn-1+⋯+apxn-p+b0un+b1un-1+⋯+bqun-q=∑i=1paixn-i+∑i=0qbiun-iwhere ai(i=1,2,⋯,p) and bi(i=0,1,2,⋯,q) are parameters of p order AR model and q order MA model, respectively. Taking Z transform on (37), obtaining(38)-∑i=0paixn-i=∑i=0qbiun-i(39)-∑i=0paiz-iXz=∑i=0qbiz-iUz(40)1-∑i=1paiz-iXz=1+∑i=1qbiz-iUz(41)Hz=XzUz=1+∑i=1qbiz-i1-∑i=1paiz-i=BzAz(42)hn↔Hz,xn=hn∗un=∑i=0qhiun-iwhere h(n) is impulse response function, H(z) is system transfer function, A(z) and B(z) are the FIR and IIR filters, respectively. From (41), FLO-ARMA model spectrum of α stable distribution process on the unit circle (z=ejω) is given by(43)Sαejω=γuHejωα==γu1+∑i=1qbie-jωi1-∑i=1paie-jωiαThen FLO-ARMA model spectrum of SαS (γu=1) process on the unit circle can be written as(44)SSαSejω=Hejωα=1+∑i=1qbie-jωi1-∑i=1pake-jωiαMultiplyingx(n-m) both sides of (38) and computing fractional low order covariance, then(45)Dxxm=∑i=1paiDxxm-i+∑i=0qbiDuxm-iwhere(46)Dxxm=Exn<a>xn-m<b>Dxxm-i=ΔExn<a>xn-m+i<b>Duxm-i=ΔEun<a>xn-m+i<b>Then we can get [20](47)Duxm=0m>0γuh-mα-1m<0When m>q, (45) changes as(48)Dxxm=∑i=1paiDxxm-i,m>qEquation (48) can be written as(49)DqDq-1⋯Dq-p+1Dq+1Dq⋯Dq-p+2⋮⋮⋱⋮Dq+p-1Dq+p-2⋯Dqa1a2⋮ap=Dq+1Dq+2⋮Dq+pwhere (49) is also a generalized Yule-Walker equation. Then we can obtain the parameters of the FLO-AR model by solving (49).When0≤m≤q, (45) changes as(50)Dxxm=∑i=0paiDxxm-i+∑i=0qbiDuxm-i=∑i=1paiDxxm-i+∑i=0qbihk-iα-1,0≤m≤qLetting p=0, then(51)Dxxm=γu∑i=0p-mbi+mbiα-1,0≤m≤qγu∑i=0p-mbibi-mα-1,-q≤m≤0By solving the nonlinear equation (51), the parameters bi(i=1,2,⋯,q) of the FLO-MA model can be obtained.An approximate infinitep′ order FLO-AR model can be used to be equivalent to a finite q′ order MA model (p′≫q′), then(52)am+∑i=1q′biam-i=emwhere e(m) is error. Letting e(m)=0, then(53)∑i=1q′biam-i=-amBy taking autocorrelation of both sides of (53) with a∗(n-m) one has(54)∑i=1q′biRm-i=-Rmwhere R(m)=Ea(m)a∗(n-m) and R(m-i)=Ea(m-i)a∗(n-m+i). Equation (38) can be written as(55)R0R-1⋯R1-q′R1R0⋯R2-q′⋮⋮⋱⋮Rq′-1Rq′-2⋯R0b1b2⋮bp=-R1R2⋮Rq′Equation (55) is a Toeplitz matrix. By solving the Yule-Walker equation (55), the parameters bi(i=1,2,⋯,q) of the FLO-MA model are gotten.We summarize the steps of the FLO-ARMA model spectrum method, as follows.Step 1.
Compute the parametersai(i=1,2,⋯,p) of the FLO-AR model of the signal x(n) employing (49).Step 2.
By solving the nonlinear equation (51) or the Yule-Walker equation (55), the parameters bi(i=1,2,⋯,q) of the FLO-MA model can be gotten.Step 3.
Determine the orderp and q employing FPE criterion.Step 4.
Substituteai(i=1,2,⋯,p) and bi(i=1,2,⋯,q) into (29) or (30); compute FLO-ARMA model spectrum of the signal x(n).


## 3.3.2. Application Review

In this simulation,x(n) in (21) is used as the experimental test signal. The traditional ARMA method and modified FLO-ARMA method are applied to estimate the frequency spectrum of the signal x(n) under Gaussian noise (SNR=-5dB) and SαS noise environment (GSNR=20dB). The results of ten independent experiments are shown in Figures 13-14.Figure 13
The ARMA model and FLO-ARMA model spectrums of the signalx(n) under Gaussian noise environment (SNR=-5dB) ((a) ARMA model spectrums of the signal x(n), where p=29 and q=5; (b) FLO-ARMA model spectrums of the signal x(n), where p=27, q=6, and a=b=0.9).
(a) (b)Figure 14
FLO-ARMA model and ARMA model spectrums of the signalx(n) under SαS (α=1.3) noise environment (GSNR=20dB) ((a) ARMA model spectrums of the signal x(n), where p=20 and q=5; (b) FLO-ARMA model spectrums of the signal x(n), where p=21, q=5, and a=b=0.1).
(a) (b)


## 3.3.3. Remarks

The ARMA frequency spectrum of the signalx(n) is shown in Figure 13(a), and Figure 13(b) is FLO-ARMA frequency spectrum of the signal x(n) under Gaussian(SNR=-5dB) noise environment. The results show that both methods can well estimate out the frequency spectrum of the signal x(n). Figure 14(a) shows that the traditional ARMA model spectrum method completely fails under SαS (α=1.3) noise environment (GSNR=20dB), but the improved FLO-ARMA model spectrum method in Figure 14(b) has very good toughness. Hence, the ARMA spectrum method is only suitable to estimate the frequency spectrum of the signals under Gaussian environment, but the FLO-ARMA model method can be applied in Gaussian and SαS noise environment, which can work in lower GSNR.The traditional parameter model spectrum estimation methods are equivalent to solving Yule-Walker equation and their extension equations based on the autocorrelation function. However, the fractional low order parameter model spectrum estimation methods construct the generalized Yule-Walker equation and their extension equations employing fractional low order covariance matrix, because fractional low order autocovariance is a generalized covariance, which changes into autocorrelation whena=b=2. As a result, the fractional low order parameter model spectrum estimation methods can be applied to both Gaussian and stable distribution noise environment.


## 3.4. Spectral Resolution of the FLO Spectrum Estimation Methods

In order to further verify the spectral resolution capability of the improved fractional low order spectrum methods, in this section, we apply the FLOBTCS, FLOPCS, FLO-Welch, FLOAR, and FLOARMA model spectrum methods to estimate the spectrum of the signal with closely spaced frequencies. The signaly(n) in (21) is used as the experimental test signal, where f0=100, f1=110, f2=500, n=0,1,2,…,N-1, and N=2000. v(n) is SαS stable distribution noise, α=1.4, and GSNR=22dB. The simulations of ten independent experiments are shown in Figure 15.Figure 15
Frequency spectrum estimations of the signal with closely spaced frequencies employing fractional lower order frequency spectrum methods ((a) fractional lower order Blackman-Tukey method; (b) fractional lower order Periodogram method; (c) fractional lower order Welch covariation method; (d) fractional lower order AR model method; (e) fractional lower order ARMA model method).
(a) (b) (c) (d) (e)Figures15(a) and 15(b) and Figure 15(c) are the FLOBTCS, FLOPCS, and FLO-Welch spectrums of the signal x(n), respectively. The FLOAR(90) model spectrum is shown in Figure 15(d), and Figure 15(e) is FLOAR(65, 3) model spectrum of the signal x(n). The results show that the improved methods can distinguish the closely spaced frequencies well. Hence, the proposed methods have better performance of spectral estimation and higher spectral resolution, and the FLOARMA model spectrum method performs best.


## 4. Application Simulations

In this simulation, the experiment signal adopts from the real normal signal and bearing outer race fault signal (DE) in Section2. 0.2 seconds data is selected as the test signal, which is collected at 12,000 samples per second, then N=2400. We apply the existing frequency spectrum methods and improved frequency spectrum methods to analyze the signals; the simulation results are shown in Figures 16–19.Figure 16
Frequency spectrum estimations of the normal bearing signal employing the conventional frequency spectrum methods ((a) Blackman-Tukey method; (b) Periodogram method; (c) Welch covariation method; (d) AR model method; (e) ARMA model method).
(a) (b) (c) (d) (e)Figure 17
Frequency spectrum estimations of the normal bearing signal employing fractional lower order frequency spectrum methods ((a) fractional lower order Blackman-Tukey method; (b) fractional lower order Periodogram method; (c) fractional lower order Welch covariation method; (d) fractional lower order AR model method; (e) fractional lower order ARMA model method).
(a) (b) (c) (d) (e)Figure 18
Frequency spectrum estimations of the bearing outer race fault signal in DE employing the conventional frequency spectrum ((a) Blackman-Tukey method; (b) Periodogram method; (c) Welch covariation method; (d) AR model method; (e) ARMA model method).
(a) (b) (c) (d) (e)Figure 19
Frequency spectrum estimations of the bearing outer race fault signal in DE employing fractional lower order frequency spectrums and time-frequency representation ((a) fractional lower order Blackman-Tukey method; (b) fractional lower order Periodogram method; (c) fractional lower order Welch covariation spectrum method; (d) fractional lower order AR model method; (e) fractional lower order ARMA model method; (f) fractional lower order Choi-Williams time-frequency representation).
(a) (b) (c) (d) (e) (f)Figures16 and 17 are frequency spectrum estimations of the normal bearing signal employing the conventional frequency spectrum methods and fractional lower order frequency spectrum methods, respectively. The figures show that the methods all can estimate out the frequency components of the normal signal; the transient harmonic vibration components is about 1060Hz. We can see that the BTCS, PCS, FLOBTCS, and FLOPCS methods have a certain side lobe in Figures 16(a), 16(b), 17(a), and 17(b), respectively, which lead to the distortion of the signal. The performance of the proposed methods is comparable to that of the traditional methods. The Welch and FLO-Welch methods in Figures 16(c) and 17(c) suppress the side lobe effects because of the spectrum segmentation. The FLO-AR and FLO-ARMA model methods in Figures 17(d) and 17(e) can better demonstrate the normal signal than the AR and ARMA method in Figures 16(d) and 16(e). Hence, the proposed methods have an important advantage in the normal bearing signal analysis.The frequency spectrums of the bearing outer race fault signal (DE) employing the traditional methods and the improved methods are shown in Figures18(a)–18(e) and Figures 19(a)–19(e), respectively, and Figure 19(f) is fractional lower order Choi-Williams time-frequency representation of the bearing outer race fault signal (DE). We can see that the transient harmonic vibration components about 600 Hz, 2800 Hz, and 3500 Hz dominate frequency domain from Figures 18(a)–18(e) and Figures 19(a)–19(e), and Figure 19(f) clearly reveals the gap between the impact. The interval between the pulse points A, B, C, D, E, and F is approximately 30ms, then the characteristic frequency of the bearing outer race fault is about 33.333Hz. However, the improved methods for the pulsive signals have better performance than the conventional methods in terms of variance and resolution; especially, the FLO-AR and FLO-ARMA model methods have higher resolution than the existing AR and ARMA methods, which can better demonstrate the frequency spectrums of the impulsive signals. Furthermore, the improved methods have their respective advantages and disadvantages, as follows. The FLOBTCS and FLOPCS methods in Figures 19(a) and 19(b) have some side lobes, the FLOBTCS method has lager variance, and their resolution is lower. The FLO-Welch method in Figure 19(c) has some improvement in suppressing side lobes and variance, but its resolution is not ideal. The parameter model frequency spectrum methods in Figures 19(d)–19(e) have better performance than the previous methods, which have higher resolution, less computation, and no side lobe interference.In order to further verify the advantages of the proposed methods,α distribution noise (α=1.1, GSNR=20dB) is added in the SαS stable distribution outer race fault signal as the actual working environment background noise. The improved methods and existing methods are applied to analyze frequency spectrum feature of the out race fault signal; ten independent experiments are shown in Figures 20 and 21, respectively. The simulation results show that the conventional methods in Figures 20(a)–20(e) fail, but the proposed methods in Figures 21(a)–21(e) have good performance, and we can know that the transient harmonic vibration components are about 600 Hz, 2800 Hz, and 3500 Hz.Figure 20
Frequency spectrum estimations of the bearing outer race fault signal in DE contaminated bySαS distribution noise employing the conventional frequency spectrum methods SαS distribution noise (α=1.1, GSNR=20dB) ((a) Blackman-Tukey method; (b) Periodogram method; (c) Welch covariation spectrum method; (d) AR model method; (e) ARMA model method).
(a) (b) (c) (d) (e)Figure 21
Frequency spectrum estimations of the bearing outer race fault signal in DE contaminated bySαS distribution noise employing the fractional lower order frequency spectrum methods (α=1.1, GSNR=20dB) ((a) fractional lower order Blackman-Tukey method; (b) fractional lower order Periodogram method; (c) fractional lower order Welch covariation spectrum method; (d) fractional lower order AR model method; (e) fractional lower order ARMA model method).
(a) (b) (c) (d) (e)


## 5. Conclusions

α stable distribution is a more appropriate statistical model for the bearing fault signals, and 1<α<2. We have proposed the improved FLOBTCS, FLOPCS, FLOWCS, FLOAR, and FLOARMA parameter model frequency spectrum analysis methods employing fractional low order statistics; the methods are suitable for Gaussian and α stable distribution noise environment and have good toughness; their performance is superior to the existing spectral estimation methods, which can effectively suppress SαS distribution noise and work in low GSNR. The methods are very stable in frequency estimation, which have better performance of spectral estimation and higher spectral resolution. The advantages of the FLOARMA model spectrum method are more obvious, which can effectively suppress side lobe and improve resolution. The methods have their respective advantages and disadvantages, the FLOBTCS and FLOPCS methods have a certain side lobe, and their resolutions are lower. The FLOBTCS method has lager variance, but it is very easy to compute it directly. The FLOWCS can effectively reduce the side lobe, which decreases the resolution of spectrum estimation. The FLOAR and FLOARMA parameter model frequency spectrum methods have advantages of specific sense of mathematics, higher resolution, and inexistent problem of side lobe interference, and the calculation is not very complicated. Hence, the parameter model frequency spectrum methods have high performance. In view of the α stable distribution bearing fault signal, even α stable distribution noise environment, the proposed methods can be applied to analyze their spectrum features. Combining the improved frequency spectrum methods with the fractional low order time-frequency methods, we can get the transient harmonic vibration components and fault characteristic frequency of the bearing fault signals and achieve their joint fault diagnosis. Therefore, the joint fault diagnosis method may melt the bright side of each, which is significance for the bearing fault diagnosis. In the future, we can also further study the methods, such as developing them into fractional low order multispectral, inverse spectrum, spectrum autocorrelation, and combined analysis based on EMD.



---

*Source: 7641383-2019-08-27.xml*