# Application of FIFD-OMGD-1.5D Teager Energy to Extract Microfault Features of Rolling Bearing

Authors: Fangming Liu; Jie Ma
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

Deconvolution-related methods are the mainstream choice when it comes to enhancing the pulse impact of bearing fault and reducing noise interference. Kurtogram algorithm is used to optimize the minimum generalized Lp/Lq deconvolution to improve the nonconvexity of other optimization criteria. However, it has low computational efficiency and poor diagnostic accuracy under strong background noise. The paper proposes an optimized method using protrugram algorithm that combines fast iterative filter decomposition (FIFD) with minimum generalized Lp/Lq deconvolution (OMGD) for the 1.5-dimension Teager energy spectrum demodulation. Here is the specific process of the application: Fast iterative filtering (FIF) was used to reduce noise interference before using the maximum kurtosis to obtain the center frequency and frequency band and optimize the filter design, which was for the MGD initialization operation to prevent the result from falling into the local optimal solution and check the interference of impulse noise to a certain extent. The 1.5-dimension Teager energy spectrum was then used for demodulation analysis to extract small fault features of rolling bearings. The verification of simulation signals and actual data showed that this method was better in terms of extraction effect and efficiency than the use of fast kurtogram algorithm to optimize minimum generalized Lp/Lq deconvolution when it comes to extracting microfault features with high interference of background noise.


---

## Body

## 1. Introduction

The bearing is an indispensable component in the mechanical system. It is prone to all types of failure because of its fast motion speed, large load, and complex structure [1]. All failures start from small faults (early faults or initial faults) according to the final analysis. However, when the bearing signals are being collected from sensors, the signals will inevitably be mixed with other vibration interference components. Besides, one of the characteristics of microfault is small amplitude, which can make the fault impulse response very weak. Thus, it is easy to drown in the strong background noise, and it is too hard to find and solve timely. Over time, the equipment may be damaged or shut down, even bringing significant personal casualty and huge economic loss. For example, from September 2003 to October 2004, major accidents of freight train derailment caused by bearing fatigue fracture occurred many times in some important sections of China, resulting in economic losses of more than 2 billion yuan [2]; in 2010, a Russian airliner had an accident during landing due to motor shaft failure, resulting in 132 deaths [3]; in 2015, a serious leakage accident occurred in a petrochemical company in China due to the rupture of the ball caused by the distortion and fracture of the bearing inner ring [4]. In addition, the operational safety of modern engineering system which has large scale and complexity has gradually reduced [5]. Therefore, the study on microfaults of rolling bearings is helpful to prevent equipment accidents and ensure the operation safety of mechanical systems [6–9]. However, how to detect and enhance fault impulse response and extract fault features effectively has always been an important concern in bearing microfault diagnosis, and it is also a popular subject in the study of nonlinear vibration signals.Due to the background noise and interference, noise reduction filtering is a necessary step before extracting the features of microfaults of rolling bearings. Wavelet can be used to denoise nonstationary signals [10]. However, its transform results depend on wavelet basis and lack adaptability. Furthermore, energy leakage occurs easily during the signal transformation [11], which leads to the nonideal results of denoising. Therefore, some decomposition methods are proposed to denoise the signal, such as empirical mode decomposition (EMD) [12], ensemble empirical mode decomposition (EEMD) [13], complementary ensemble empirical mode decomposition (CEEMD) [14], complete ensemble empirical mode decomposition with adaptive noise (CEEMDAN) [15], and variational mode decomposition (VMD) [16]. Most decomposition technologies break down the signal into some IMFs. In order to achieve noise reduction, the components with more useful information are retained and the those with more noise are discarded. However, such decomposition methods only work better for the signals that have very short data points. When long data points are involved, they will have various problems, such as poor decomposition effect and long decomposition time [17]. Considering that the iterative filtering (IF) method is an iterative decomposition method inspired by EMD and that it can ensure convergence and stability [18], paper [19] has proposed an iterative decomposition method based on fast Fourier transform (FFT), namely, fast iterative filtering decomposition (FIFD) method. This method has excellent rapidity and accuracy in signal decomposition.Furthermore, because the generalized Lp/Lq norm applied in blind deconvolution performs well in extracting sparse features from noise signals [20], the generalized Lp/Lq norm sparse filtering method for pulse feature enhancement is applied to rolling bearing fault diagnosis successfully [21]. Based on the contributions of the above papers, paper [22] has proposed the minimum correlation generalized Lp/Lq deconvolution method, which can successfully extract the fault features and complete the detection of composite faults. Because the objective function of the minimum generalized Lp/Lq deconvolution (MGD) is nonconvex, the result may fall into the local optimal solution if improper initial values are selected. Therefore, paper [9] has used the filter designed by fast kurtogram to provide appropriate initialization and achieve adaptive adjustment. It is helpful to find the frequency band of the best filter and optimize MGD. Therefore, it is very important to adopt an appropriate method based on deconvolution to enhance the impact of bearing fault pulse.Moreover, since the nonlinear signal can modulate amplitude and modulation frequency, the nonlinear broadband component can be converted into narrowband component through demodulation technology; that is, the concentrated spectrum can be obtained [23]. Therefore, the traditional demodulation technology, such as Hilbert transform (HT) [24] and Teager energy operator (TEO) demodulation [25], is usually applied to the bearing fault signal that can modulate amplitude and modulation frequency after noise reduction and brings out the fault features. However, considering the poor noise resistance of these demodulation methods, the paper [26] has proposed that the 1.5-dimension energy spectrum be applied to analyze the bearing fault, since it can reduce noise well and recover nonlinear features. Thus, it works well in extracting the frequency doubling information of fault feature frequency.To sum up, in order to solve the computational efficiency loss in the above optimization methods, this paper has proposed MGD optimization based on the protrugram algorithm (referred to as OMGD herein). It is used to obtain the center frequency and the corresponding frequency band at the maximal kurtosis [27, 28]. Moreover, we can use the filter parameters obtained from this algorithm to design the filter and then to complete MGD optimization. This method does not achieve calculation accuracy at the expense of calculation efficiency. On the contrary, it ensures not only the accuracy of the diagnosis results, but also the calculation speed. However, the optimization method alone is not enough to process the signal under strong background noise. Therefore, in order to maintain the performance of deconvolution when the signal is affected by strong background noise, this paper has proposed a combination of OMGD and FIFD. First, we used FIFD to eliminate the noise and other signal interference, improve the signal-to-noise ratio, and reduce the number of objective functions converging to the local optimal solution. Then, the deconvolution method was used to make it easier to converge to the global optimal solution and then to enhance the pulse characteristics of small faults. Finally, we used the 1.5-dimension Teager energy spectrum for demodulation analysis to complete the feature extraction of microfaults. With the proposed FIFD-OMGD-1.5-dimension Teager energy spectrum, we could extract the microfault features of rolling bearing efficiently and accurately even with serious interference of noise.The rest of the content is arranged in the following order: first, an introduction of the formula derivation of FIFD method is presented in Section2. Second, an introduction of the derivation process of OMGD is provided in Section 3. Finally, an introduction of the 1.5-dimension energy spectrum is given in Section 4. In Sections 3 and 4, we use the abovementioned methods to analyze and verify the analog signal and actual signal, respectively. Finally, we summarize the study in Section 5.


## 2. Methodology

### 2.1. Fast Iterative Filtering Decomposition Method

Fast iterative filter decomposition (FIFD) is essentially the discrete version of iterative filter decomposition (IFD) based on fast Fourier transform (FFT).Assuming that the signal issx,x∈R, we sample it at n points xj=j/n−1 to ensure that the sampling rate captures all its details. The main objective of FIFD is that vector IMFs are obtained by vector sxjj=0n−1. In addition, it is assumed that sxj2=1 to make the method more general.For convenience,sxjj=0n−1 is replaced by symbol s to simplify the formula.There are two “while” loops in the discrete version of the iterative filter decomposition (DIFD) algorithm. The first “while” loop is known as an outer loop, and the second one is an inner loop.In the inner loop, the firstIMF can be obtained through IMF1=limm⟶∞I−Cmsm, where the matrix Cm=ωmxi−xji,j=0n−1 represents the discrete convolution operator, and its structure is determined by the semi-support length lm at each step m. A nonnegative number is used for ωm. The closely supported symmetric filter ωm within area 1 is generated by the convolution of the symmetric filter hm and itself. We can use the following equation to calculate ωmu:(1)ωmu=ωg−1u∫−11ωtgm′tdt,where gm is the scaling function, which is considered monotonically reversible. By linear scaling, g−1u=u/lm,gm′t=lm≥0. Thus, we can rewrite (1) as(2)ωmu=ωu/lmlm.Suppose that for any numberm and m≥1, there is C=Cm. Given the semi-simple eigenvalue λpp=0,…,n−1 of C and its eigenvector upp=0,…,n−1, we can define the column of matrix U as the eigenvector up. In addition, suppose there are k zero eigenvalues in C, where k∈0,1,…,n−1.Therefore, in the external circulation, the firstIMF is obtained by the following formula:(3)IMF1=limm⟶0I−Cms=UQUTs,where U is known as unitary matrix. Besides, the k elements in the diagonal of the diagonal matrix Q are equal to 1.To make DIFD reach the accuracyδ in calculating an IMF, we need to approximate the number of iterations through the following equation:(4)N0N0N0+1N0+1<δs˜∞n−1−k,N0∈N,where s˜=UTs.Therefore, we can calculateIMF from the following formula:(5)IMF=UP01−λ1N0⋱1−λn−1−kN01⋱1×PTUTs,where the permutation matrix obtained by reordering the columns of U is represented by P. Since the column of U corresponds to the eigenvector of C, λ represents the corresponding eigenvalues in descending order.UTs is obtained by DFT of s using FFT algorithm. In addition, its computational complexity is nlogn. Multiplying the matrix U on the left side of UTs is equivalent to calculating the IDFT, which can be completed by using the IFFT. Therefore, IMF can use the following formula for rapid calculation:(6)IMF=∑k=0n−1uk1−λkσkN0=IDFTI−DN0DFTs,where D is a diagonal matrix in which its diagonal includes the eigenvalue of C, and σk is the k−th element obtained after DFT of s.In order to shorten the execution time of the algorithm, a stop criterion is introduced here, that is,(7)MIsx−MI+1sxL2<δ,∀I≥N0,where Msm=sm−Lmsm=sm+1 and Lmsmx=∫−lmlmsmx+tωmtdt.In order to avoidm unlimited growth, we use the above stopping criterion to intervene, shown as follows:(8)SD≔sm+1−sm2sm2.When the value ofSD is less than or equal to δ, the process can be stopped, so as to greatly shorten the calculation time.


### 2.2. OMGD Method

#### 2.2.1. Generalized Lp/Lq Norm

In some current literature, the concise sparse representation of signal is expressed by the following formula:(9)minθθl0,x−Aθ≤ε,where A represents the dictionary matrix or base matrix, θ represents the expected sparse representation, x represents the input signal, and ε represents the estimation error which is an infinitesimal positive number.However, because it is difficult to optimizel0, some literature has proposed to use lp norm 0<p≤1 to represent the sparsity of signal. Because signal sparsity regularization was not applicable in the study, the study in [29] has proposed to make the signal sparse due to several advantages of l1/l2 in calculation. Its advantages are as follows: First, the scale remains unchanged. Second, it has all the input gradient information and will give an easy minimization algorithm. According to that, l1/l2 norm is a modified sparse index. We further extend l1/l2 to the generalized lp/lq norm, which is a generalized sparse measure. Its definition can be expressed by the following formula:(10)Jp,qx=logqp⋅xlpxlqp=logqp⋅∑iNxip∑iNxiqp/q,where logq/p is a constant with a sign. According to the above explanations, the minimization method can be used to replace the optimization method.


#### 2.2.2. OMGD

The method based on the generalizedlp/lq norm can optimize the sparsity of the input characteristic matrix. Each increasing diagonal element in the characteristic matrix of the input signal obtained by this method is a constant. This matrix is called Hankel matrix, shown as follows:(11)H=x1x2⋯xlx2x3⋯xl+1⋮⋮⋮⋮xN−l+1xN−l+2⋯xN,where x∈RN represents the input signal vector and H∈RN−L×L represents the Hankel matrix derived from x. We use f=H⋅W to describe the signal features, where W∈RL is the weighted vector and f∈RN−L+1 represents the feature vector after sparse learning. The sparse constraint applied to the eigenvector f is expressed by the following equation:(12)minfJp,qf=minf∑iNfip∑iNfiqp/qsubject to:f=H⋅W,WL2=1,p<q.Because sparse filtering and deconvolution can be equated mathematically, the sparse filtering of generalizedlp/lq norm is renamed as minimum generalized lp/lq deconvolution (hereinafter referred to as MGD) in [8, 9]. We can use the following equation to represent deconvolution:(13)y=x∗ginv,where x indicates the input signal, ginv is the inverse filter that we desire, y is the filtered signal, and ∗ is the convolution operator.In MGD, according to the relationship between sparse filtering and deconvolution, we can rewrite the convolution process of the above formula as the following formula by using Hankel matrix:(14)y=H⋅W,where W∈RL is equivalent to ginv in the formula. At this time, the filtered signal y to which the sparse constraint is applied has sparse features, which can be indicated as follows:(15)miny,wJp,qy,subject toy=H⋅W,W2=1,p<q.The reason why we cannot minimizeJp,qy is that the objective function in (15) is not smooth and the function also has nonconvexity. Therefore, MGD needs to be optimized:(1)
The nonsmooth problem can be solved by replacingyn with the soft absolute function ci=yn2+ξ, where ξ is a small number, which is greater than 0 and is generally specified as ξ=10−8.(2)
We choose to use the gradient descent algorithm to optimize the nonconvexity of the function. The gradient term of the optimization algorithm is as follows:(16)∂Jp,q∂wj=∑i∂Jp,q∂ci⋅∂ci∂yi⋅∂yi∂wj=p∑icip−1∑iciqp/q−∑ici⋅ciq−1∑iciqp/q+1⋅fifi2+ε⋅xi+j−1.However, it is indispensable to confirm the initial value before gradient descent. A good initial value is very important to solve the optimal solution. The initial value optimization method adopted in [9] is to realize the design of the filter through the fast kurtogram algorithm and then take the filter coefficient as the initial value. However, when the environmental noise has a great impact and the signal-to-noise ratio is low, the fast kurtogram can only provide a large bandwidth. Therefore, it may result in the lack of accuracy of the center frequency and filter bandwidth and then make it unable to find the parameters correctly.Therefore, in this paper, we used the protrugram algorithm to find the optimal initial value of MGD. The algorithm selects the optimal frequency band based on the kurtosis of the signal narrowband envelope spectrum, which solves the problem of low parameter accuracy obtained by fast kurtogram algorithm. The general steps are as follows:(1)
Thep and q parameter values of the generalized Lp/Lq norm, the number of decomposition layers k, and the filter length L are set; the coefficients of the stage 0 all-pass filter are initialized to 0,1, 0, …, 0(2)
Taking three times of the fault characteristic frequency of the rolling bearing as the fixed bandwidth value of the protrugram method, the kurtosis value corresponding to each center frequency is calculated according to the protrugram algorithm(3)
The bandwidth and center frequency with the largest kurtosis are selected as filter parameters to design the corresponding filter(4)
The designed filter coefficient is used as the initial value of MGD to complete the optimization work


### 2.3. 1.5-Dimension Teager Energy Spectrum

#### 2.3.1. Review of the Teager Energy Operator

We useyt to represent the signal and define the energy operator ψ by(17)Ψyt=dytdt2−ytd2ytdt2,where Ψyt is the instantaneous energy signal, and dyt/dt and d2yt/dt2 are the first derivative and second derivative of signal yt, respectively. Teager energy operator has good time resolution. Its output can quickly and accurately track the change of total energy required by the signal and enhance the transient features.


#### 2.3.2. Definitions of 1.5-Dimension Spectrum

For the stationary random signalyt, the third-order cumulant diagonal slice R3yT,TT1=T2=T is(18)R3yT,T=Eynyn+Tyn+T.By performing one-dimension Fourier transform onR3yT,T, the 1.5-dimension spectrum Bω of signal yt is obtained as(19)Bω=∫−∞∞R3yT,Te−jωTdT.In addition, the 1.5-dimension spectrum has the properties of strengthening the fundamental frequency component, suppressing Gaussian white noise, and detecting harmonic components [30].


#### 2.3.3. 1.5-Dimension Teager Energy Spectrum

1.5-dimension energy spectrum combines Teager energy operator and 1.5-dimension spectrum. Its principle is as follows: Firstly, the energy operator is used to demodulate the signalyt to obtain the instantaneous energy signal so that the transient impact of the signal becomes obvious. Then, a 1.5-dimensional spectrum of Ψyt is made to suppress the noise component. Therefore, the 1.5-dimension Teager energy spectrum of yt is(20)Eω=∫−∞∞R3Ψτ,τe−jωτdτ.Because the 1.5-dimension energy spectrum can not only effectively highlight the transient impact features, but also suppress the interference of noise, it is very appropriate for extracting the impact features of bearing fault signal.


### 2.4. Proposed Method for Bearing Microfault Diagnosis under Strong Background Noise

In this paper, a bearing microfault feature extraction method under strong background noise has been proposed. The flowchart of the proposed method is shown in Figure1. It consists of three main steps: First, the original vibration signal polluted by interference and noise is decomposed by FIFD method, and the useful components are selected by the combination of correlation coefficient and kurtosis criterion to complete the signal reconstruction. In the second step, OMGD is performed on the reconstructed signal to enhance the pulse characteristics. During this step, we use the parameters obtained by the protrugram algorithm to optimize the filter and then take the filter coefficient as the initial value of MGD to optimize MGD. In the third step, the 1.5-dimension Teager energy spectrum is used for demodulation analysis to extract the characteristics of bearing microfaults under strong background noise.Figure 1
Flowchart of the proposed bearing microfault diagnosis.


## 2.1. Fast Iterative Filtering Decomposition Method

Fast iterative filter decomposition (FIFD) is essentially the discrete version of iterative filter decomposition (IFD) based on fast Fourier transform (FFT).Assuming that the signal issx,x∈R, we sample it at n points xj=j/n−1 to ensure that the sampling rate captures all its details. The main objective of FIFD is that vector IMFs are obtained by vector sxjj=0n−1. In addition, it is assumed that sxj2=1 to make the method more general.For convenience,sxjj=0n−1 is replaced by symbol s to simplify the formula.There are two “while” loops in the discrete version of the iterative filter decomposition (DIFD) algorithm. The first “while” loop is known as an outer loop, and the second one is an inner loop.In the inner loop, the firstIMF can be obtained through IMF1=limm⟶∞I−Cmsm, where the matrix Cm=ωmxi−xji,j=0n−1 represents the discrete convolution operator, and its structure is determined by the semi-support length lm at each step m. A nonnegative number is used for ωm. The closely supported symmetric filter ωm within area 1 is generated by the convolution of the symmetric filter hm and itself. We can use the following equation to calculate ωmu:(1)ωmu=ωg−1u∫−11ωtgm′tdt,where gm is the scaling function, which is considered monotonically reversible. By linear scaling, g−1u=u/lm,gm′t=lm≥0. Thus, we can rewrite (1) as(2)ωmu=ωu/lmlm.Suppose that for any numberm and m≥1, there is C=Cm. Given the semi-simple eigenvalue λpp=0,…,n−1 of C and its eigenvector upp=0,…,n−1, we can define the column of matrix U as the eigenvector up. In addition, suppose there are k zero eigenvalues in C, where k∈0,1,…,n−1.Therefore, in the external circulation, the firstIMF is obtained by the following formula:(3)IMF1=limm⟶0I−Cms=UQUTs,where U is known as unitary matrix. Besides, the k elements in the diagonal of the diagonal matrix Q are equal to 1.To make DIFD reach the accuracyδ in calculating an IMF, we need to approximate the number of iterations through the following equation:(4)N0N0N0+1N0+1<δs˜∞n−1−k,N0∈N,where s˜=UTs.Therefore, we can calculateIMF from the following formula:(5)IMF=UP01−λ1N0⋱1−λn−1−kN01⋱1×PTUTs,where the permutation matrix obtained by reordering the columns of U is represented by P. Since the column of U corresponds to the eigenvector of C, λ represents the corresponding eigenvalues in descending order.UTs is obtained by DFT of s using FFT algorithm. In addition, its computational complexity is nlogn. Multiplying the matrix U on the left side of UTs is equivalent to calculating the IDFT, which can be completed by using the IFFT. Therefore, IMF can use the following formula for rapid calculation:(6)IMF=∑k=0n−1uk1−λkσkN0=IDFTI−DN0DFTs,where D is a diagonal matrix in which its diagonal includes the eigenvalue of C, and σk is the k−th element obtained after DFT of s.In order to shorten the execution time of the algorithm, a stop criterion is introduced here, that is,(7)MIsx−MI+1sxL2<δ,∀I≥N0,where Msm=sm−Lmsm=sm+1 and Lmsmx=∫−lmlmsmx+tωmtdt.In order to avoidm unlimited growth, we use the above stopping criterion to intervene, shown as follows:(8)SD≔sm+1−sm2sm2.When the value ofSD is less than or equal to δ, the process can be stopped, so as to greatly shorten the calculation time.


## 2.2. OMGD Method

### 2.2.1. Generalized Lp/Lq Norm

In some current literature, the concise sparse representation of signal is expressed by the following formula:(9)minθθl0,x−Aθ≤ε,where A represents the dictionary matrix or base matrix, θ represents the expected sparse representation, x represents the input signal, and ε represents the estimation error which is an infinitesimal positive number.However, because it is difficult to optimizel0, some literature has proposed to use lp norm 0<p≤1 to represent the sparsity of signal. Because signal sparsity regularization was not applicable in the study, the study in [29] has proposed to make the signal sparse due to several advantages of l1/l2 in calculation. Its advantages are as follows: First, the scale remains unchanged. Second, it has all the input gradient information and will give an easy minimization algorithm. According to that, l1/l2 norm is a modified sparse index. We further extend l1/l2 to the generalized lp/lq norm, which is a generalized sparse measure. Its definition can be expressed by the following formula:(10)Jp,qx=logqp⋅xlpxlqp=logqp⋅∑iNxip∑iNxiqp/q,where logq/p is a constant with a sign. According to the above explanations, the minimization method can be used to replace the optimization method.


### 2.2.2. OMGD

The method based on the generalizedlp/lq norm can optimize the sparsity of the input characteristic matrix. Each increasing diagonal element in the characteristic matrix of the input signal obtained by this method is a constant. This matrix is called Hankel matrix, shown as follows:(11)H=x1x2⋯xlx2x3⋯xl+1⋮⋮⋮⋮xN−l+1xN−l+2⋯xN,where x∈RN represents the input signal vector and H∈RN−L×L represents the Hankel matrix derived from x. We use f=H⋅W to describe the signal features, where W∈RL is the weighted vector and f∈RN−L+1 represents the feature vector after sparse learning. The sparse constraint applied to the eigenvector f is expressed by the following equation:(12)minfJp,qf=minf∑iNfip∑iNfiqp/qsubject to:f=H⋅W,WL2=1,p<q.Because sparse filtering and deconvolution can be equated mathematically, the sparse filtering of generalizedlp/lq norm is renamed as minimum generalized lp/lq deconvolution (hereinafter referred to as MGD) in [8, 9]. We can use the following equation to represent deconvolution:(13)y=x∗ginv,where x indicates the input signal, ginv is the inverse filter that we desire, y is the filtered signal, and ∗ is the convolution operator.In MGD, according to the relationship between sparse filtering and deconvolution, we can rewrite the convolution process of the above formula as the following formula by using Hankel matrix:(14)y=H⋅W,where W∈RL is equivalent to ginv in the formula. At this time, the filtered signal y to which the sparse constraint is applied has sparse features, which can be indicated as follows:(15)miny,wJp,qy,subject toy=H⋅W,W2=1,p<q.The reason why we cannot minimizeJp,qy is that the objective function in (15) is not smooth and the function also has nonconvexity. Therefore, MGD needs to be optimized:(1)
The nonsmooth problem can be solved by replacingyn with the soft absolute function ci=yn2+ξ, where ξ is a small number, which is greater than 0 and is generally specified as ξ=10−8.(2)
We choose to use the gradient descent algorithm to optimize the nonconvexity of the function. The gradient term of the optimization algorithm is as follows:(16)∂Jp,q∂wj=∑i∂Jp,q∂ci⋅∂ci∂yi⋅∂yi∂wj=p∑icip−1∑iciqp/q−∑ici⋅ciq−1∑iciqp/q+1⋅fifi2+ε⋅xi+j−1.However, it is indispensable to confirm the initial value before gradient descent. A good initial value is very important to solve the optimal solution. The initial value optimization method adopted in [9] is to realize the design of the filter through the fast kurtogram algorithm and then take the filter coefficient as the initial value. However, when the environmental noise has a great impact and the signal-to-noise ratio is low, the fast kurtogram can only provide a large bandwidth. Therefore, it may result in the lack of accuracy of the center frequency and filter bandwidth and then make it unable to find the parameters correctly.Therefore, in this paper, we used the protrugram algorithm to find the optimal initial value of MGD. The algorithm selects the optimal frequency band based on the kurtosis of the signal narrowband envelope spectrum, which solves the problem of low parameter accuracy obtained by fast kurtogram algorithm. The general steps are as follows:(1)
Thep and q parameter values of the generalized Lp/Lq norm, the number of decomposition layers k, and the filter length L are set; the coefficients of the stage 0 all-pass filter are initialized to 0,1, 0, …, 0(2)
Taking three times of the fault characteristic frequency of the rolling bearing as the fixed bandwidth value of the protrugram method, the kurtosis value corresponding to each center frequency is calculated according to the protrugram algorithm(3)
The bandwidth and center frequency with the largest kurtosis are selected as filter parameters to design the corresponding filter(4)
The designed filter coefficient is used as the initial value of MGD to complete the optimization work


## 2.2.1. Generalized Lp/Lq Norm

In some current literature, the concise sparse representation of signal is expressed by the following formula:(9)minθθl0,x−Aθ≤ε,where A represents the dictionary matrix or base matrix, θ represents the expected sparse representation, x represents the input signal, and ε represents the estimation error which is an infinitesimal positive number.However, because it is difficult to optimizel0, some literature has proposed to use lp norm 0<p≤1 to represent the sparsity of signal. Because signal sparsity regularization was not applicable in the study, the study in [29] has proposed to make the signal sparse due to several advantages of l1/l2 in calculation. Its advantages are as follows: First, the scale remains unchanged. Second, it has all the input gradient information and will give an easy minimization algorithm. According to that, l1/l2 norm is a modified sparse index. We further extend l1/l2 to the generalized lp/lq norm, which is a generalized sparse measure. Its definition can be expressed by the following formula:(10)Jp,qx=logqp⋅xlpxlqp=logqp⋅∑iNxip∑iNxiqp/q,where logq/p is a constant with a sign. According to the above explanations, the minimization method can be used to replace the optimization method.


## 2.2.2. OMGD

The method based on the generalizedlp/lq norm can optimize the sparsity of the input characteristic matrix. Each increasing diagonal element in the characteristic matrix of the input signal obtained by this method is a constant. This matrix is called Hankel matrix, shown as follows:(11)H=x1x2⋯xlx2x3⋯xl+1⋮⋮⋮⋮xN−l+1xN−l+2⋯xN,where x∈RN represents the input signal vector and H∈RN−L×L represents the Hankel matrix derived from x. We use f=H⋅W to describe the signal features, where W∈RL is the weighted vector and f∈RN−L+1 represents the feature vector after sparse learning. The sparse constraint applied to the eigenvector f is expressed by the following equation:(12)minfJp,qf=minf∑iNfip∑iNfiqp/qsubject to:f=H⋅W,WL2=1,p<q.Because sparse filtering and deconvolution can be equated mathematically, the sparse filtering of generalizedlp/lq norm is renamed as minimum generalized lp/lq deconvolution (hereinafter referred to as MGD) in [8, 9]. We can use the following equation to represent deconvolution:(13)y=x∗ginv,where x indicates the input signal, ginv is the inverse filter that we desire, y is the filtered signal, and ∗ is the convolution operator.In MGD, according to the relationship between sparse filtering and deconvolution, we can rewrite the convolution process of the above formula as the following formula by using Hankel matrix:(14)y=H⋅W,where W∈RL is equivalent to ginv in the formula. At this time, the filtered signal y to which the sparse constraint is applied has sparse features, which can be indicated as follows:(15)miny,wJp,qy,subject toy=H⋅W,W2=1,p<q.The reason why we cannot minimizeJp,qy is that the objective function in (15) is not smooth and the function also has nonconvexity. Therefore, MGD needs to be optimized:(1)
The nonsmooth problem can be solved by replacingyn with the soft absolute function ci=yn2+ξ, where ξ is a small number, which is greater than 0 and is generally specified as ξ=10−8.(2)
We choose to use the gradient descent algorithm to optimize the nonconvexity of the function. The gradient term of the optimization algorithm is as follows:(16)∂Jp,q∂wj=∑i∂Jp,q∂ci⋅∂ci∂yi⋅∂yi∂wj=p∑icip−1∑iciqp/q−∑ici⋅ciq−1∑iciqp/q+1⋅fifi2+ε⋅xi+j−1.However, it is indispensable to confirm the initial value before gradient descent. A good initial value is very important to solve the optimal solution. The initial value optimization method adopted in [9] is to realize the design of the filter through the fast kurtogram algorithm and then take the filter coefficient as the initial value. However, when the environmental noise has a great impact and the signal-to-noise ratio is low, the fast kurtogram can only provide a large bandwidth. Therefore, it may result in the lack of accuracy of the center frequency and filter bandwidth and then make it unable to find the parameters correctly.Therefore, in this paper, we used the protrugram algorithm to find the optimal initial value of MGD. The algorithm selects the optimal frequency band based on the kurtosis of the signal narrowband envelope spectrum, which solves the problem of low parameter accuracy obtained by fast kurtogram algorithm. The general steps are as follows:(1)
Thep and q parameter values of the generalized Lp/Lq norm, the number of decomposition layers k, and the filter length L are set; the coefficients of the stage 0 all-pass filter are initialized to 0,1, 0, …, 0(2)
Taking three times of the fault characteristic frequency of the rolling bearing as the fixed bandwidth value of the protrugram method, the kurtosis value corresponding to each center frequency is calculated according to the protrugram algorithm(3)
The bandwidth and center frequency with the largest kurtosis are selected as filter parameters to design the corresponding filter(4)
The designed filter coefficient is used as the initial value of MGD to complete the optimization work


## 2.3. 1.5-Dimension Teager Energy Spectrum

### 2.3.1. Review of the Teager Energy Operator

We useyt to represent the signal and define the energy operator ψ by(17)Ψyt=dytdt2−ytd2ytdt2,where Ψyt is the instantaneous energy signal, and dyt/dt and d2yt/dt2 are the first derivative and second derivative of signal yt, respectively. Teager energy operator has good time resolution. Its output can quickly and accurately track the change of total energy required by the signal and enhance the transient features.


### 2.3.2. Definitions of 1.5-Dimension Spectrum

For the stationary random signalyt, the third-order cumulant diagonal slice R3yT,TT1=T2=T is(18)R3yT,T=Eynyn+Tyn+T.By performing one-dimension Fourier transform onR3yT,T, the 1.5-dimension spectrum Bω of signal yt is obtained as(19)Bω=∫−∞∞R3yT,Te−jωTdT.In addition, the 1.5-dimension spectrum has the properties of strengthening the fundamental frequency component, suppressing Gaussian white noise, and detecting harmonic components [30].


### 2.3.3. 1.5-Dimension Teager Energy Spectrum

1.5-dimension energy spectrum combines Teager energy operator and 1.5-dimension spectrum. Its principle is as follows: Firstly, the energy operator is used to demodulate the signalyt to obtain the instantaneous energy signal so that the transient impact of the signal becomes obvious. Then, a 1.5-dimensional spectrum of Ψyt is made to suppress the noise component. Therefore, the 1.5-dimension Teager energy spectrum of yt is(20)Eω=∫−∞∞R3Ψτ,τe−jωτdτ.Because the 1.5-dimension energy spectrum can not only effectively highlight the transient impact features, but also suppress the interference of noise, it is very appropriate for extracting the impact features of bearing fault signal.


## 2.3.1. Review of the Teager Energy Operator

We useyt to represent the signal and define the energy operator ψ by(17)Ψyt=dytdt2−ytd2ytdt2,where Ψyt is the instantaneous energy signal, and dyt/dt and d2yt/dt2 are the first derivative and second derivative of signal yt, respectively. Teager energy operator has good time resolution. Its output can quickly and accurately track the change of total energy required by the signal and enhance the transient features.


## 2.3.2. Definitions of 1.5-Dimension Spectrum

For the stationary random signalyt, the third-order cumulant diagonal slice R3yT,TT1=T2=T is(18)R3yT,T=Eynyn+Tyn+T.By performing one-dimension Fourier transform onR3yT,T, the 1.5-dimension spectrum Bω of signal yt is obtained as(19)Bω=∫−∞∞R3yT,Te−jωTdT.In addition, the 1.5-dimension spectrum has the properties of strengthening the fundamental frequency component, suppressing Gaussian white noise, and detecting harmonic components [30].


## 2.3.3. 1.5-Dimension Teager Energy Spectrum

1.5-dimension energy spectrum combines Teager energy operator and 1.5-dimension spectrum. Its principle is as follows: Firstly, the energy operator is used to demodulate the signalyt to obtain the instantaneous energy signal so that the transient impact of the signal becomes obvious. Then, a 1.5-dimensional spectrum of Ψyt is made to suppress the noise component. Therefore, the 1.5-dimension Teager energy spectrum of yt is(20)Eω=∫−∞∞R3Ψτ,τe−jωτdτ.Because the 1.5-dimension energy spectrum can not only effectively highlight the transient impact features, but also suppress the interference of noise, it is very appropriate for extracting the impact features of bearing fault signal.


## 2.4. Proposed Method for Bearing Microfault Diagnosis under Strong Background Noise

In this paper, a bearing microfault feature extraction method under strong background noise has been proposed. The flowchart of the proposed method is shown in Figure1. It consists of three main steps: First, the original vibration signal polluted by interference and noise is decomposed by FIFD method, and the useful components are selected by the combination of correlation coefficient and kurtosis criterion to complete the signal reconstruction. In the second step, OMGD is performed on the reconstructed signal to enhance the pulse characteristics. During this step, we use the parameters obtained by the protrugram algorithm to optimize the filter and then take the filter coefficient as the initial value of MGD to optimize MGD. In the third step, the 1.5-dimension Teager energy spectrum is used for demodulation analysis to extract the characteristics of bearing microfaults under strong background noise.Figure 1
Flowchart of the proposed bearing microfault diagnosis.


## 3. Simulation Validation

### 3.1. Construction of Bearing Fault Simulation Signal

In this part, we use the same method to diagnose the analog signal to further verify its effectiveness. In this part, we construct the simulation model of defective bearing as follows:(21)Y=y0e−2πgfnt0sin2πfs1−g2t0+nt,where y0 is the displacement constant, the natural frequency is shown by fn, g is the damping coefficient, fs is the sampling frequency, and t0 is the single cycle sampling time. Besides, nt is Gaussian white noise and its amplitude is 1.5 [31].Sety0=2,g=0.1,fn=2kHz,andfs=20kHz; the repetition period is 0.1 s; the number of sampling points is 10000; and the fault feature frequency of the simulation signal is 100 Hz.The waveform diagram and envelope spectrum from the above simulation model are shown in Figure2. Due to the interference of noise, we cannot identify the fault impact component of the signal. Direct envelope analysis of simulation signals cannot directly find the correct fault frequency.Figure 2
Simulation signal: (a) its waveform and (b) its corresponding envelope spectrum.
(a)(b)


### 3.2. Verification of the Proposed Method

Now, the proposed method is used to extract the fault characteristic frequency from the severely damaged signal. Fast iterative filter decomposition (FIFD) is implemented on the original simulation signal to obtain 8 IMFs. Besides, the correlation degree between each component and bearing signal is evaluated by calculating the cross-correlation coefficient. The effective component with more information has more correlation with the signal, and the calculated correlation value is large [32]. Of course, kurtosis is often used to characterize bearing faults. The larger the kurtosis value is, the more obvious the fault impact component is, and the easier it is to extract fault information [33]. However, because kurtosis is particularly susceptible to abnormal components, it is difficult to distinguish random pulses from cycle pulses [9, 34, 35]. Therefore, we use the screening method on cross-correlation coefficient and kurtosis value to select the sensitive IMF and then to complete the filtered signal reconstruction.The cross-correlation coefficients and kurtosis values are shown in Tables1 and 2. According to the screening criteria, we select IMF1 and IMF2 components for signal reconstruction after filtering. The waveform of the reconstructed signal is shown in Figure 3. We can find out that compared with the original simulation model, the reconstructed signal after preliminary filtering shows obvious impact, which proves the effect of the fast iterative filtering decomposition method in noise reduction filtering.Table 1
Kurtosis values of components obtained by decomposing the simulation model.
IMF1IMF2IMF3IMF4IMF5IMF6IMF7IMF82.87793.36223.06153.04913.31313.04184.62112.6637Table 2
Correlation values of components obtained by decomposing the simulation model.
IMF1IMF2IMF3IMF4IMF5IMF6IMF7IMF80.83290.59760.33520.21690.14810.08470.07320.0506Figure 3
Reconstructed signal after FIFD filtering.However, we can also see that after fast iterative filtering, there are still some noise and interference in the reconstructed signal. Therefore, its periodicity is not shown, and the amplitude of the filtered signal is significantly reduced. To solve these problems, we will use OMGD based on protrugram algorithm to filter the reconstructed signal again and enhance the fault pulse features, so as to reduce the difficulty of extracting small faults in the filtering process. Firstly, we set the parameters of the protrugram algorithm responsible for MGD initialization. We set the iteration step to 100; select the triple fault characteristic frequency as the fixed bandwidth, that is,bw=300Hz; and determine the center frequency as CF=bw/2:step:Fs/2−bw/2. According to the protrugram algorithm, the center frequency corresponding to the maximum kurtosis value is 1050 Hz. According to the above values, the corresponding band-pass filter is designed, and the coefficient of the band-pass filter is taken as the initial value of MGD iteration. Then, we set the filter length to 100 in MGD, p = 1, q = 2. Analyzing the reconstructed signal by using OMGD based on protrugram algorithm, the processing result is shown in Figure 4(a). We can see from the figure that after OMGD processing the signal shows obvious periodicity and impact, the filtering effect is very obvious and the signal amplitude has been significantly improved. In addition, the signal processed by OMGD is analyzed by 1.5-dimension Teager energy spectrum. We can see that the fault feature frequency and frequency doubling feature of the simulation signal are effectively extracted. From the final envelope spectrum, we can know that using the FIFD-OMGD combined method for noise reduction filtering in addition to 1.5-dimension Teager energy spectrum analysis has a better effect in extracting small fault features of rolling bearing.Figure 4
The reconstructed signal based on protrugram algorithm processing after OMGD: (a) its waveform and (b) its 1.5D Teager energy spectrum.
(a)(b)


### 3.3. The Comparison of the Proposed Method with FIFD-OMGD (Based on Fast Kurtogram Algorithm), CEEMDAN-OMGD, and FIFD-MCKD

We compare the proposed method with FIFD-OMGD (based on fast kurtogram algorithm), CEEMDAN-OMGD, and FIFD-MCKD and analyze the results to evaluate the effectiveness of the proposed method.Referring to the idea of control variable method, firstly, we apply OMGD based on fast kurtogram algorithm to the reconstructed signal after FIFD decomposition. The processing result is shown in Figure5. In this method, we keep the parameter settings in MGD consistent with the parameters of OMGD based on protrugram algorithm; that is, the filter length in MGD is set to 100, p=1, and q = 2. The only difference is that the two methods provide two different initialization techniques for the initial value of MGD. Among them, OMGD based on fast kurtogram algorithm uses fast kurtogram to search the resonance band of the bearing, designs the filter by laying the frequency plane in the case of 1/3 binary tree, and finally completes the initialization operation via taking the filter coefficient as the initial value.Figure 5
The signal processed by OMGD using fast kurtogram algorithm: (a) its waveform and (b) envelope spectrum.
(a)(b)By comparing Figures4 and 5, we can see that the OMGD method using the fast kurtogram algorithm has exposed many problems in feature extraction of signals seriously polluted by noise. For example, when the measured signal is tampered with by some strong interference signals, the fault impulse enhancement effect has been reduced greatly. In contrast, it achieves more accuracy and robustness in extracting microfault features under strong background interference.Through the comparison of the above contents, we can see that OMGD method based on protrugram algorithm has more advantages in enhancing the fault pulse impact of contaminated signal. Therefore, in the following comparison, we choose OMGD method based on protrugram algorithm. In order to highlight the advantages of FIFD algorithm in the proposed method, this paper selects the method of combining CEEMDAN and OMGD to process the contaminated signal. The processing result is shown in Figure6.Figure 6
Filtered signal by CEEMDAN-OMGD: (a) its waveform and (b) its 1.5D Teager energy spectrum.
(a)(b)Before using the above decomposition tool, we need to set some predefined parameters. The key parameters are as follows: for CEEMDAN, the noise standard deviation is set to 0.2, the number of times noise is added is 100, and the maximum number of iterations is 100. The number of IMFs is set to 8. For FIFD, the mask length is 2 in this paper and the number of IMFs is 8. Unless otherwise specified, these parameter values no longer change.In order to reflect the decomposition accuracy of the two decomposition methods better, we calculate the RMSE (root mean square error) between the eight IMF components of the two methods and the ground truth. It can be clearly seen from Figure7 that the RMSE of IMFs decomposed by FIFD is significantly smaller than that of IMFs decomposed by CEEMDAN algorithm. Therefore, FIFD is obviously better than CEEMDAN algorithm in decomposition accuracy.Figure 7
RMSE between 8 IMFs and the ground truth.In addition, by comparing Figures4 and 6, we see that the signal after being denoised and filtered using the CEEMDAN-OMGD method presents obvious impact and periodicity. Compared with the signal filtered by the other method, the signal is sparser, which makes it sensible to suspect that using this method may result in losing some useful information in the process of denoising and filtering. Finally, the filtered signal was analyzed using the 1.5-dimension Teager energy spectrum. It can be observed from the energy spectrum that only the triple frequency was extracted using the method. In addition, by comparison with the effect shown from the quintuple fault feature frequency extracted in Figure 4(b), we have more ground to suspect that the above method has lost some useful information in filtering. This also reveals the robustness of the FIFD-OMGD method.Next, we choose to use FIFD algorithm to denoise the simulation signal, and then compare the denoised signal with MCKD method and OMGD method based on protrugram algorithm to verify the performance of OMGD. Similarly, we need to set the parameters before using MCKD. For MCKD, the number of iterations is set to 30, and the filter length is set to 100. Then, we use MCKD to enhance the fault pulse impact of the reconstructed signal processed by FIFD algorithm. From Figure8, we see that the signal impact and periodicity after the signal was filtered by FIFD-MCKD are not obvious, which indicates that filtering using FIFD-OMGD combined with noise reduction is better than that using FIFD-MCKD. To further expand on that, in Figure 8(b), the feature frequency extraction effect of this method and the performance of frequency doubling are not ideal. In addition, the amplitudes of the first and second frequency doubling of the extracted fault feature frequency are extremely small compared with the frequency amplitude extracted by the FIFD-OMGD method. Thus, we conclude that the performance of OMGD is better in fault pulse features enhancement.Figure 8
Filtered signal by FIFD-MCKD: (a) its waveform and (b) its 1.5D Teager energy spectrum.
(a)(b)In order to illustrate the effectiveness of 1.5-dimension Teager energy spectrum, we performed envelope demodulation analysis on the signal after FIFD-OMGD noise reduction directly. The analysis result is shown in Figure9. According to the comparison between Figures 4(b) and 9, we see that the effect of enveloping the joint filtered signal directly is not very satisfactory. However, using the 1.5-dimension Teager energy spectrum can well suppress the interference of noise and highlight the transient impact characteristics of the signal. The superiority of demodulation analysis using 1.5-dimension Teager energy spectrum is confirmed. Due to spatial constraints, we use 1.5-dimension Teager energy spectrum for demodulation analysis directly in the next experimental data verification part.Figure 9
Envelope spectrum of the filtered signal after FIFD-OMGD processing.


### 3.4. Computational Complexity

This section will compare the three methods mentioned with the proposed method in terms of calculation time. Figure10 shows the relationship between data length and execution time, where the sampling data interval is 10000. It can be clearly seen from the figure that the calculation efficiency of the proposed method is the highest among the four methods. Moreover, when the data length increases from 10000 to 50000, the calculation time curve of the proposed method is always below the other three curves, which also shows that the calculation efficiency of the proposed method is still in the leading position among the four methods no matter how long the data length is.Figure 10
The relationship between data length and execution time.


### 3.5. Conservativeness of the Proposed Method

Because the proposed method in this paper is aimed at the signal polluted by strong background noise, we add different degrees of noise to the simulation signal in this section. The conservativeness of the proposed method is verified by comparing the fault feature frequency extracted by different methods with the ideal fault feature frequency under different SNR and calculating the average relative error between them. The results are obtained as shown in Table3. It can be seen from Table 3 that with different SNR, the accuracy of fault feature frequency extraction by the proposed method is higher than that of other algorithms. In addition, as more and more noise is added, the ability of the proposed method to extract fault feature frequency is not affected, which also verifies that the method has good conservativeness.Table 3
Fault feature frequency extracted by different methods and average relative error with different SNRs.
Fault feature Frequency (Hz)SNRAlgorithm−3−5−7−9−11−13−15Average relative errorTrue100100100100100100100FIFD-OMGD (based on protrugram algorithm)99.899.899.8100.2100.299.8100.20.000200FIFD-OMGD (based on fast kurtogram algorithm)10198.9898.9810198.9898.9898.980.010143CEEMDAN-OMGD101101102.2101101.698.981020.014029FIFD-MCKD100.1100.110198.6107.3111.889.40.046142


## 3.1. Construction of Bearing Fault Simulation Signal

In this part, we use the same method to diagnose the analog signal to further verify its effectiveness. In this part, we construct the simulation model of defective bearing as follows:(21)Y=y0e−2πgfnt0sin2πfs1−g2t0+nt,where y0 is the displacement constant, the natural frequency is shown by fn, g is the damping coefficient, fs is the sampling frequency, and t0 is the single cycle sampling time. Besides, nt is Gaussian white noise and its amplitude is 1.5 [31].Sety0=2,g=0.1,fn=2kHz,andfs=20kHz; the repetition period is 0.1 s; the number of sampling points is 10000; and the fault feature frequency of the simulation signal is 100 Hz.The waveform diagram and envelope spectrum from the above simulation model are shown in Figure2. Due to the interference of noise, we cannot identify the fault impact component of the signal. Direct envelope analysis of simulation signals cannot directly find the correct fault frequency.Figure 2
Simulation signal: (a) its waveform and (b) its corresponding envelope spectrum.
(a)(b)


## 3.2. Verification of the Proposed Method

Now, the proposed method is used to extract the fault characteristic frequency from the severely damaged signal. Fast iterative filter decomposition (FIFD) is implemented on the original simulation signal to obtain 8 IMFs. Besides, the correlation degree between each component and bearing signal is evaluated by calculating the cross-correlation coefficient. The effective component with more information has more correlation with the signal, and the calculated correlation value is large [32]. Of course, kurtosis is often used to characterize bearing faults. The larger the kurtosis value is, the more obvious the fault impact component is, and the easier it is to extract fault information [33]. However, because kurtosis is particularly susceptible to abnormal components, it is difficult to distinguish random pulses from cycle pulses [9, 34, 35]. Therefore, we use the screening method on cross-correlation coefficient and kurtosis value to select the sensitive IMF and then to complete the filtered signal reconstruction.The cross-correlation coefficients and kurtosis values are shown in Tables1 and 2. According to the screening criteria, we select IMF1 and IMF2 components for signal reconstruction after filtering. The waveform of the reconstructed signal is shown in Figure 3. We can find out that compared with the original simulation model, the reconstructed signal after preliminary filtering shows obvious impact, which proves the effect of the fast iterative filtering decomposition method in noise reduction filtering.Table 1
Kurtosis values of components obtained by decomposing the simulation model.
IMF1IMF2IMF3IMF4IMF5IMF6IMF7IMF82.87793.36223.06153.04913.31313.04184.62112.6637Table 2
Correlation values of components obtained by decomposing the simulation model.
IMF1IMF2IMF3IMF4IMF5IMF6IMF7IMF80.83290.59760.33520.21690.14810.08470.07320.0506Figure 3
Reconstructed signal after FIFD filtering.However, we can also see that after fast iterative filtering, there are still some noise and interference in the reconstructed signal. Therefore, its periodicity is not shown, and the amplitude of the filtered signal is significantly reduced. To solve these problems, we will use OMGD based on protrugram algorithm to filter the reconstructed signal again and enhance the fault pulse features, so as to reduce the difficulty of extracting small faults in the filtering process. Firstly, we set the parameters of the protrugram algorithm responsible for MGD initialization. We set the iteration step to 100; select the triple fault characteristic frequency as the fixed bandwidth, that is,bw=300Hz; and determine the center frequency as CF=bw/2:step:Fs/2−bw/2. According to the protrugram algorithm, the center frequency corresponding to the maximum kurtosis value is 1050 Hz. According to the above values, the corresponding band-pass filter is designed, and the coefficient of the band-pass filter is taken as the initial value of MGD iteration. Then, we set the filter length to 100 in MGD, p = 1, q = 2. Analyzing the reconstructed signal by using OMGD based on protrugram algorithm, the processing result is shown in Figure 4(a). We can see from the figure that after OMGD processing the signal shows obvious periodicity and impact, the filtering effect is very obvious and the signal amplitude has been significantly improved. In addition, the signal processed by OMGD is analyzed by 1.5-dimension Teager energy spectrum. We can see that the fault feature frequency and frequency doubling feature of the simulation signal are effectively extracted. From the final envelope spectrum, we can know that using the FIFD-OMGD combined method for noise reduction filtering in addition to 1.5-dimension Teager energy spectrum analysis has a better effect in extracting small fault features of rolling bearing.Figure 4
The reconstructed signal based on protrugram algorithm processing after OMGD: (a) its waveform and (b) its 1.5D Teager energy spectrum.
(a)(b)


## 3.3. The Comparison of the Proposed Method with FIFD-OMGD (Based on Fast Kurtogram Algorithm), CEEMDAN-OMGD, and FIFD-MCKD

We compare the proposed method with FIFD-OMGD (based on fast kurtogram algorithm), CEEMDAN-OMGD, and FIFD-MCKD and analyze the results to evaluate the effectiveness of the proposed method.Referring to the idea of control variable method, firstly, we apply OMGD based on fast kurtogram algorithm to the reconstructed signal after FIFD decomposition. The processing result is shown in Figure5. In this method, we keep the parameter settings in MGD consistent with the parameters of OMGD based on protrugram algorithm; that is, the filter length in MGD is set to 100, p=1, and q = 2. The only difference is that the two methods provide two different initialization techniques for the initial value of MGD. Among them, OMGD based on fast kurtogram algorithm uses fast kurtogram to search the resonance band of the bearing, designs the filter by laying the frequency plane in the case of 1/3 binary tree, and finally completes the initialization operation via taking the filter coefficient as the initial value.Figure 5
The signal processed by OMGD using fast kurtogram algorithm: (a) its waveform and (b) envelope spectrum.
(a)(b)By comparing Figures4 and 5, we can see that the OMGD method using the fast kurtogram algorithm has exposed many problems in feature extraction of signals seriously polluted by noise. For example, when the measured signal is tampered with by some strong interference signals, the fault impulse enhancement effect has been reduced greatly. In contrast, it achieves more accuracy and robustness in extracting microfault features under strong background interference.Through the comparison of the above contents, we can see that OMGD method based on protrugram algorithm has more advantages in enhancing the fault pulse impact of contaminated signal. Therefore, in the following comparison, we choose OMGD method based on protrugram algorithm. In order to highlight the advantages of FIFD algorithm in the proposed method, this paper selects the method of combining CEEMDAN and OMGD to process the contaminated signal. The processing result is shown in Figure6.Figure 6
Filtered signal by CEEMDAN-OMGD: (a) its waveform and (b) its 1.5D Teager energy spectrum.
(a)(b)Before using the above decomposition tool, we need to set some predefined parameters. The key parameters are as follows: for CEEMDAN, the noise standard deviation is set to 0.2, the number of times noise is added is 100, and the maximum number of iterations is 100. The number of IMFs is set to 8. For FIFD, the mask length is 2 in this paper and the number of IMFs is 8. Unless otherwise specified, these parameter values no longer change.In order to reflect the decomposition accuracy of the two decomposition methods better, we calculate the RMSE (root mean square error) between the eight IMF components of the two methods and the ground truth. It can be clearly seen from Figure7 that the RMSE of IMFs decomposed by FIFD is significantly smaller than that of IMFs decomposed by CEEMDAN algorithm. Therefore, FIFD is obviously better than CEEMDAN algorithm in decomposition accuracy.Figure 7
RMSE between 8 IMFs and the ground truth.In addition, by comparing Figures4 and 6, we see that the signal after being denoised and filtered using the CEEMDAN-OMGD method presents obvious impact and periodicity. Compared with the signal filtered by the other method, the signal is sparser, which makes it sensible to suspect that using this method may result in losing some useful information in the process of denoising and filtering. Finally, the filtered signal was analyzed using the 1.5-dimension Teager energy spectrum. It can be observed from the energy spectrum that only the triple frequency was extracted using the method. In addition, by comparison with the effect shown from the quintuple fault feature frequency extracted in Figure 4(b), we have more ground to suspect that the above method has lost some useful information in filtering. This also reveals the robustness of the FIFD-OMGD method.Next, we choose to use FIFD algorithm to denoise the simulation signal, and then compare the denoised signal with MCKD method and OMGD method based on protrugram algorithm to verify the performance of OMGD. Similarly, we need to set the parameters before using MCKD. For MCKD, the number of iterations is set to 30, and the filter length is set to 100. Then, we use MCKD to enhance the fault pulse impact of the reconstructed signal processed by FIFD algorithm. From Figure8, we see that the signal impact and periodicity after the signal was filtered by FIFD-MCKD are not obvious, which indicates that filtering using FIFD-OMGD combined with noise reduction is better than that using FIFD-MCKD. To further expand on that, in Figure 8(b), the feature frequency extraction effect of this method and the performance of frequency doubling are not ideal. In addition, the amplitudes of the first and second frequency doubling of the extracted fault feature frequency are extremely small compared with the frequency amplitude extracted by the FIFD-OMGD method. Thus, we conclude that the performance of OMGD is better in fault pulse features enhancement.Figure 8
Filtered signal by FIFD-MCKD: (a) its waveform and (b) its 1.5D Teager energy spectrum.
(a)(b)In order to illustrate the effectiveness of 1.5-dimension Teager energy spectrum, we performed envelope demodulation analysis on the signal after FIFD-OMGD noise reduction directly. The analysis result is shown in Figure9. According to the comparison between Figures 4(b) and 9, we see that the effect of enveloping the joint filtered signal directly is not very satisfactory. However, using the 1.5-dimension Teager energy spectrum can well suppress the interference of noise and highlight the transient impact characteristics of the signal. The superiority of demodulation analysis using 1.5-dimension Teager energy spectrum is confirmed. Due to spatial constraints, we use 1.5-dimension Teager energy spectrum for demodulation analysis directly in the next experimental data verification part.Figure 9
Envelope spectrum of the filtered signal after FIFD-OMGD processing.


## 3.4. Computational Complexity

This section will compare the three methods mentioned with the proposed method in terms of calculation time. Figure10 shows the relationship between data length and execution time, where the sampling data interval is 10000. It can be clearly seen from the figure that the calculation efficiency of the proposed method is the highest among the four methods. Moreover, when the data length increases from 10000 to 50000, the calculation time curve of the proposed method is always below the other three curves, which also shows that the calculation efficiency of the proposed method is still in the leading position among the four methods no matter how long the data length is.Figure 10
The relationship between data length and execution time.


## 3.5. Conservativeness of the Proposed Method

Because the proposed method in this paper is aimed at the signal polluted by strong background noise, we add different degrees of noise to the simulation signal in this section. The conservativeness of the proposed method is verified by comparing the fault feature frequency extracted by different methods with the ideal fault feature frequency under different SNR and calculating the average relative error between them. The results are obtained as shown in Table3. It can be seen from Table 3 that with different SNR, the accuracy of fault feature frequency extraction by the proposed method is higher than that of other algorithms. In addition, as more and more noise is added, the ability of the proposed method to extract fault feature frequency is not affected, which also verifies that the method has good conservativeness.Table 3
Fault feature frequency extracted by different methods and average relative error with different SNRs.
Fault feature Frequency (Hz)SNRAlgorithm−3−5−7−9−11−13−15Average relative errorTrue100100100100100100100FIFD-OMGD (based on protrugram algorithm)99.899.899.8100.2100.299.8100.20.000200FIFD-OMGD (based on fast kurtogram algorithm)10198.9898.9810198.9898.9898.980.010143CEEMDAN-OMGD101101102.2101101.698.981020.014029FIFD-MCKD100.1100.110198.6107.3111.889.40.046142


## 4. Experimental Data Validation

The experimental data in this paper come from the accelerated life test data set of XJTU-SY rolling bearing conducted by Xi'an Jiao tong University. The data set collects the vibration signals of bearings switching from normal state to failure state, which can be used to diagnose the microfaults of rolling bearings [36]. The bearing tested is LDK UER204 rolling bearing, and the geometric parameters are shown in Table 4. During data acquisition, the sampling frequency was set to 25.6 kHz, the sampling interval was 1 min, and the sampling duration of each time was 1.28 s [36].Table 4
Geometric parameters of LDK UER204 rolling bearing.
Number of ballsm (piece)Ball diameterd (mm)Bearing pitch diameterD (mm)Contact angleα (°)87.9234.550In this section, the fault vibration signal of the bearing outer ring was selected under the working conditions of 2400 r/min of rotation speed and 10 kN radial force. The fault feature frequency of the used bearing outer ring was 123.32 Hz, which is calculated by the calculation formula in Table5.Table 5
Calculation formula of fault characteristic frequency of rolling bearing.
Fault typeCalculation formula of fault characteristic frequency of rolling bearingOuter ring faultf0=n/601−d/DCOSαm/2The waveform and envelope of the outer ring vibration signal are shown intuitively in Figure11. It can be seen that the original vibration signal contains a lot of random noise and interference and that the signal features are weak. The bearing fault feature frequency extracted by envelope analysis is not obvious.Figure 11
Original signal: (a) its waveform and (b) its envelope spectrum.
(a)(b)As in the previous section, the proposed method was used for fast iterative filter decomposition (FIFD) to obtain 8 IMFs on the original outer ring signal. The cross-correlation coefficients and kurtosis values of the IMFs from decomposing the outer ring signal are shown in Tables6 and 7. According to the screening criteria, we select IMF1 and IMF3 components for signal reconstruction after filtering. The waveform of the reconstructed signal is shown in Figure 12. It is obvious that most random interference in the reconstructed signal after fast iterative filtering has been removed. It also shows obvious periodicity and impact, which also shows the effectiveness of fast iterative filtering (FIF) in noise reduction filtering.Table 6
Kurtosis values of components obtained by decomposing the outer ring signal.
IMF1IMF2IMF3IMF4IMF5IMF6IMF7IMF83.20363.08183.25932.30821.72197.75113.05941.8473Table 7
Correlation values of the components obtained by decomposing the outer ring signal.
IMF1IMF2IMF3IMF4IMF5IMF6IMF7IMF80.81160.38010.47400.30510.11940.01740.00720.0056Figure 12
Reconstructed signal.However, in Figure12, we see that the amplitude of the reconstructed signal after fast iterative filtering is small. To solve this problem, we use OMGD based on protrugram algorithm to increase the fault feature of reconstructed signal. In Figure 13, the signal amplitude was significantly improved, and the 1.5-dimension envelope of the signal shows that the fault feature frequency and its doubled frequency are effectively and obviously extracted. The extracted fault feature frequency is 124.1 Hz.Figure 13
Filtered signal after processing reconstructed signal by OMGD: (a) its waveform and (b) its 1.5D Teager energy spectrum.
(a)(b)Similarly, we apply the proposed method, FIFD-OMGD (based on fast kurtogram algorithm), CEEMDAN-OMGD, and FIFD-MCKD to the actual fault data of bearing outer ring and compare their processing results (the parameter settings of these algorithms have been described earlier and will not be described here). The effectiveness of the proposed method is further verified by analyzing the results.Firstly, we apply FIFD-OMGD (based on fast kurtogram algorithm) to the bearing outer ring signal, and the processing results are shown in Figure14. By comparing Figures 13 and 14, we can find that FIFD-OMGD (based on fast kurtogram algorithm) is not outstanding in feature extraction of signals seriously polluted by noise. This also verifies the superiority of the proposed method indirectly.Figure 14
Filtered signal by FIFD-OMGD based on fast kurtogram algorithm: (a) its waveform and (b) its 1.5D Teager energy spectrum.
(a)(b)Secondly, in order to highlight the advantages of FIFD algorithm in the proposed method, this paper selects the method of combining CEEMDAN with OMGD based on protrugram algorithm to process the bearing outer ring signal. The processing results are shown in Figure15. It can be seen from Figure 15(b) that the fault feature frequency extracted by CEEMDAN-OMGD method is not complete, and the frequency doubling of the fault feature frequency is submerged in the noise. Comparing this with Figure 13, we can see that FIFD-OMGD method is significantly better than CEEMDAN-OMGD method.Figure 15
Filtered signal by CEEMDAN-OMGD: (a) its waveform and (b) its 1.5D Teager energy spectrum.
(a)(b)Finally, we choose to use the FIFD algorithm to denoise the bearing outer ring signal, and then use the MCKD method to process the denoised signal. The processing results are shown in Figure16. It can be seen from Figure 16 that the signal impact and periodicity after filtering with FIFD-MCKD are not obvious, which shows that the effect of FIFD-OMGD combined with noise reduction is better than that of FIFD-MCKD. Moreover, in Figure 16(b), the extraction effect of this method on the fault feature frequency and its frequency doubling are not very ideal. Furthermore, the fault feature frequency and its frequency doubling of the outer ring signal are not clearly extracted. It can be seen that OMGD based on protrugram algorithm is obviously better in enhancing fault pulse feature.Figure 16
Filtered signal by FIFD-MCKD: (a) its waveform and (b) its 1.5D Teager energy spectrum.
(a)(b)In addition, we can clearly see from Table8 that the computational efficiency of FIFD-OMGD (based on protrugram algorithm) is still better than that of the other three methods. Although the relative error of FIFD-MCKD in extracting fault feature frequency is slightly less than that of the proposed method, the extraction effect of the proposed method is still the best in combination with the energy spectrum of the above method. Therefore, taking everything into consideration, the proposed method in this paper has good accuracy, high efficiency, and good conservativeness in extracting the microfault feature of seriously polluted signals.Table 8
Computational time and error analysis of fault feature frequency for different algorithms.
AlgorithmComputational time (s)FCF (Hz)Relative errorTrue—123.32FIFD-OMGD (based on protrugram algorithm)42.653022124.10.006325FIFD-OMGD (based on fast kurtogram algorithm)133.73336038.780.685534CEEMDAN-OMGD174.489326121.50.014758FIFD-MCKD107.101329122.90.003406


## 5. Conclusion

In this paper, the combination of FIFD-OMGD and 1.5-dimension energy spectrum was used to extract the features of microfault pulse under massive noise. In order to solve the problem of low precision and low efficiency of OMGD method using the fast kurtogram algorithm to extract bearing fault features under strong background noise, the protrugram algorithm was proposed to determine the filter parameters for the filter design, and the designed filter coefficient was adopted as the initial value of MGD to achieve more efficient pulse feature enhancement. After using different methods to analyze and compare the signal of the outer ring of the simulation bearing and the actual bearing, the following conclusions are drawn:(1)
By comparing the FIFD-OMGD with the OMGD method that uses fast kurtogram algorithm, we can see that the proposed method has better accuracy, noise robustness, and efficiency in fault feature extraction.(2)
By comparing CEEMDAN-OMGD and FIFD-OMGD, it can be seen that the decomposition speed of FIFD is significantly faster than CEEMDAN. This efficiency is very important for processing large-scale data. In addition, for strong background noise, FIFD method has better noise reduction effect and stronger robustness.(3)
By comparing FIFD-OMGD with FIFD-MCKD, we can see that the effect of noise reduction of FIFD-OMGD combined with noise reduction filter is better. Moreover, FIFD-OMGD is obviously better in enhancing fault pulse features.(4)
By comparing the analysis results of direct envelope demodulation and 1.5-dimension Teager energy spectrum demodulation, we can see that the performance of the latter method is better in that it highlights the transient impact features of signal and restrains the interference of noise, which further affirms the superiority of 1.5-dimension Teager energy spectrum in demodulation analysis.To sum up, to combine FIFD-OMGD with noise reduction filter and 1.5-dimension Teager energy spectrum is more effective as a practical tool in the microfault feature detection of rolling bearings. In addition, because it is outstanding in microfault feature extraction with strong background interference, this method can be applied to a more general range of targets, such as planetary gear vibration signal seriously polluted by noise. Moreover, signal processing is combined with optimized deconvolution to enhance the performance of the latter, which provides a new idea for optimizing the deconvolution during the processing of strong interference signals in the future.



---

*Source: 7648288-2022-02-27.xml*