An Improved k-Value Symmetrical Difference Analytic Energy Operator With HTFIF and L-KCA for Bearing Fault Diagnosis

The opportune fault detection of the rolling element bearings can avoid serious equipment accidents or even casualties. However, the early fault features of the bearings are weak and often submerged in heavy background noise and interferences. This paper proposes a novel approach for bearing fault diagnosis based on the hard thresholding fast iterative filtering (HTFIF), IMF selection index integrating L- kurtosis, correlation coefficient and autocorrelation function impulse harmonic to noise ratio (L-KCA), and improved k-value symmetrical difference analytic energy operator (k-SDAEO). As an adaptive and fast time-frequency analysis method, HTFIF is first adopted to process the bearing vibration signal and obtain a series of IMFs. After that, an alternative fusion index L-KCA is developed to select the sensitive IMF. Finally, a novel k-SDAEO with strong noise robustness is presented to process the selected IMF. With this method, the weak bearing fault signatures can be identified from the energy spectrum. The performance of the proposed method comparing to the traditional methods are investigated by numerical simulation and experimental studies. The results show that the proposed method has better fault feature extraction capability and higher efficiency, which can be implemented in the on-line and real-time fault detection.

ABSTRACT The opportune fault detection of the rolling element bearings can avoid serious equipment accidents or even casualties. However, the early fault features of the bearings are weak and often submerged in heavy background noise and interferences. This paper proposes a novel approach for bearing fault diagnosis based on the hard thresholding fast iterative filtering (HTFIF), IMF selection index integrating L-kurtosis, correlation coefficient and autocorrelation function impulse harmonic to noise ratio (L-KCA), and improved k-value symmetrical difference analytic energy operator (k-SDAEO). As an adaptive and fast time-frequency analysis method, HTFIF is first adopted to process the bearing vibration signal and obtain a series of IMFs. After that, an alternative fusion index L-KCA is developed to select the sensitive IMF. Finally, a novel k-SDAEO with strong noise robustness is presented to process the selected IMF. With this method, the weak bearing fault signatures can be identified from the energy spectrum. The performance of the proposed method comparing to the traditional methods are investigated by numerical simulation and experimental studies. The results show that the proposed method has better fault feature extraction capability and higher efficiency, which can be implemented in the on-line and real-time fault detection.

FIF
fast iterative filtering HTFIF hard thresholding fast iterative filtering IMF intrinsic mode function ACFHNR autocorrelation function impulse harmonic to noise ratio L-KCA L-kurtosis combined with correlation coefficient and AFIHNR SDAEO symmetrical difference analytic energy operator The associate editor coordinating the review of this manuscript and approving it for publication was Dazhong Ma . The rolling element bearings are among the most fundamental elements in modern equipment, but they are prone to be damaged due to their frequent use and harsh working environments [1]. With the expansion of the damage, it will further lead to the failure of the equipment operation and even greater losses [2]. The function of a fault diagnosis system is to rapidly detect and determine the root causes of faults [3]. Hence, how to carry out on-line condition monitoring and real-time fault diagnosis of the rolling element bearings has become a key issue to ensure the good operation of the machine, and to reduce the maintenance cost and major accidents [3]- [5]. Fault diagnosis methods of the rolling element bearings can be categorized as model-based, signal-based, and datadriven method [6]. The model-based fault diagnostic method has a profound theoretical foundation but is only useful for the system with known mathematical models. It is difficult to obtain precise mathematical models in actual situations, which limits the use of this method [6]. Data-driven method, which does not require a prior mathematical model, is a hot spot in recent research [3], [6]- [9]. For instance, Cai et al. presented a dynamic Bayesian network (DBN)-based fault diagnosis methodology to identify the faulty component and distinguish the fault types, including transient fault (TF) and intermittent fault (IF) and Permanent fault (PF) [7]. Yang et al. proposed a distance metric named polynomial kernel induced maximum mean discrepancy (PK-MMD) to help construct a deep transfer residual network (TrResNet), which is expected to reuse the diagnosis knowledge from one machine to the other [8]. Cai et al. proposed a hybrid physicsmodel-based and data-driven remaining useful life (RUL) estimation methodology of structure systems considering the influence of multiple causes by using dynamic Bayesian networks (DBNs) [9]. Nevertheless, the data-driven methods generally have no clear physical meaning. Compared with above two methods, the signal-based method does not need to establish a mathematical model while has good physical significance. Moreover, the vibration signal of the bearing contains rich fault information as the periodic impulses will be generated when the rolling bearing has a localized defect [10]. Therefore, numerous vibration-based signal processing methods have been developed and applied to the bearing fault diagnosis, such as the short time Fourier transform (STFT), the wavelet transform (WT), the reassignment method (RM), the Wigner-Ville distribution (WVD), the sparsity assisted method and the empirical mode decomposition (EMD) [11], [12]. STFT and WT are effective and easy to implement but they both depend on the predetermined bases and cannot accomplish the high resolutions in time and frequency domains simultaneously [13], [14]. RM improves the time-frequency localization to some extent, but the original signal processed by RM cannot be recovered [15]. WVD has a higher time-frequency resolution, but it is restricted by the cross term interference [4], [15]. Sparse decomposition (e.g., matching pursuit (MP) and basis pursuit (BP)) can obtain the higher time-frequency resolution and self-adaptability, but expands the computation and storage costs [5].
As an adaptive decomposition method, EMD [16] is used to decompose a given signal into a set of intrinsic mode functions (IMFs) that have unique information [17]- [19]. It has been widely used in many fields, such as Medicine, Geophysics, Engineering, Information Technology and Economics [18]. As a powerful time-frequency analysis method, EMD also has exhibited some shortcomings, such as lack of mathematical foundation, un-stability of the sifting algorithm, noise susceptibility, intermittent and mode mixing [14], [17]. Over the past two decades, the widespread use and shortcomings of the technique have appealed to many researchers and a large number of improved and alternative approaches have been put forward. Deering and Kaiser [20] used a masking signal to improve the mode mixing and intermittent issue in the transient process, but its mode separation ability is affected by the signal magnitude [4], [21], [22]. Wu and Huang [23] and Yeh et al. [24] proposed the ensemble EMD (EEMD) and complementary EEMD (CEEMD) methods to improve the mode mixing and sensitivity of the standard EMD, but bring new issues of residue noise and spurious modes [25]. Lin et al. [14] used the convolution of the low-pass filter function and the signal itself as the ''signal moving average'' instead of the average between two envelopes in EMD and proposed the iterative filter (IF) algorithm. However, the compact support low-pass filters used in the IF are not smooth enough and may create artificial oscillations in the subsequent IMFs [26]- [28]. To solve this problem, Cicone et al. [13] constructed some smooth filters from the solutions of Fokker-Planck equations (FP filters) to replace the low-pass filters used in IF and obtained the modified IF (we call it MIF) [13], [26], [27]. Recently, Cicone [18] introduced the fast Fourier transform (FFT) and down-sampling techniques into MIF and proposed a hard thresholding fast iterative filtering (HTFIF). The HTFIF can simplify the process of MIF to one-step and make it closer to an on-line algorithm [28]. Therefore, the HTFIF decomposition method is adopted in this study to conduct on-line detection and fault diagnosis of the bearings.
After the vibration signals are decomposed by HTFIF, how to choose the sensitive IMF with rich fault information is still a tricky issue [17]. Recently, a lot of analysis methods for sensitive IMF selection have been reported. For example, Xiao et al. [29] selected the sensitive IMF by calculating the correlation coefficient between each IMF and the original signal. The correlation coefficient method has good adaptability but small discrimination and it is easily to be misjudged near the threshold boundary. Xu and Cai [30] selected sensitive IMF based on the maximum kurtosis method. Although the kurtosis has a high value for the impulse feature extracting, it is susceptible to the outliers, the noise and the sparsity of the impact [31]. In view of the above shortcomings of the traditional kurtosis method, Liu et al. [32] first applied the L-kurtosis to bearing fault diagnosis. Zheng et al. [31] proposed an autocorrelation function impulse harmonic to noise ratio (ACFHNR) index and its fusion method with Pearson correlation coefficient and kurtosis to measure the richness of the periodic pulse fault feature information [33]. In view of this, an alternative fusion index of L-kurtosis combined with correlation coefficient and ACFHNR (L-KCA) is developed in this paper to select the sensitive IMF.
However, when the working environment of the bearings is harsh, the noise and interference level is high. It is necessary to combine the mode decomposition method with other processing techniques. Energy operator demodulation technique, which has the characteristics of low computational complexity and high time resolution is typically a good option. Teager energy operator (TEO) is an effective demodulation method for the amplitude modulation-frequency modulation (AM-FM) signal and can adaptively detect the instantaneous changes of the signals [34]. However, it also has obvious disadvantages, such as low demodulation accuracy and noisesensitivity [35]. Choi and Kim [36], Choi et al. [37] proposed a multi-resolution Teager energy operator (MTEO) in the study of neural action potential detector, which has strong noise robustness compared with TEO. Feldman [38] found a new analytic energy operator (AEO) by integrating the concept of Hilbert transform and energy operator. Xu et al. [39] proposed a symmetrical difference analytic energy operator (SDAEO), which used the symmetric difference sequence with higher accuracy to replace the forward difference sequence in TEO and greatly reduced the end effect. Moreover, the symmetric difference sequence can also smooth the data and enhance the noise robustness of the energy operator. Therefore, to extract the fault features of the bearings from seriously contaminated signals quickly and accurately, an improved k-value symmetrical difference analysis energy operator (k-SDAEO), which not only retains the advantages of SDAEO but also has better robustness under low signal-to-noise ratio (SNR) is presented in this work.
The fault signal of the rolling element bearings is often accompanied by the heavy background noise. Meanwhile, the traditional fault diagnosis methods are usually timeconsuming. To accurately and fast extract the fault features, a new bearing fault diagnosis approach based on HTFIF and k-SDAEO is proposed. The main contributions of this work can be summarized as follows: The HTFIF is first adopted for bearing fault diagnosis in this paper, and it is able to decompose a given signal into a set of IMFs. Moreover, an alternative fusion indicator L-KCA is developed to select the sensitive IMF with rich fault information. Finally, a k-SDAEO is presented to apply the spectrum analysis on the selected IMF. As a result, the fault characteristics can be easily and quickly identified.
The remainder of this paper is organized as follows: Section 2 introduces the fundamental theory of the proposed method. A numerical simulation is carried out on a complex signal with heavy noise and interference harmonics in Section 3. In Section 4, two experiments are conducted to validate the effectiveness of the proposed method. Finally, the conclusions are drawn in Section 5.

II. THEORETICAL ANALYSIS A. HARD THRESHOLDING FAST ITERATIVE FILTERING
The structure of the Fast Iterative Filtering (FIF) resembles EMD and contains two nested loops: the outer loop to generate all the IMFs and the inner loop to compute each single IMF. The key step of the algorithm is the iteration equation of the inner loop [18]. For a discrete signal s(x) supported on [0, 1], the mth iteration of the inner loop is given by: and its matrix form is expressed as: in which and w m is a Fokker-Planck filter described in the article [13] with compact support on [l m , −l m ], where l m is the filter length, whose value can be tuned by parameter ζ and once set will be kept constant throughout the entire inner loop iteration process. Hence, for every m ≥ 1, l m = l 1 = l and W m = W . Since W represents the discrete convolution operator, in the case of periodicity of the signal at the boundaries, W is a circulant matrix. In this case, the eigenvalues of the circulant matrix W = w pq p,q=0,1,...,n−1 are given by: w 1q e −2πip q n , p = 0, 1, . . . , n − 1 (4) VOLUME 9, 2021 and the corresponding eigenvectors are: Given a matrix U , whose columns are the eigenvectors u p , then U is a unitary matrix. In addition, W is diagonalizable by U . Hence W = UDU T , where D is a diagonal matrix whose diagonal elements are the eigenvalues of W and U T . Assume that N 0 is the number of iterations needed by FIF method to derive the first IMF computation. Based on a predefined stopping criterion, the first IMF is given by: As we know, U T s is the discrete Fourier transform (DFT) of s, which can be computed using the FFT algorithm. Multiplying the matrix U on the left is equivalent to calculating the inverse DFT (IDFT), which can be computed using the inverse FFT. Hence, assuming we can evaluate a priori the value of N 0 needed in equation (6), then the IMF can be derived in one step as: Let M be the number of IMFs already computed. The subsequent IMFs are generated by iteratively applying the same procedure described above to the remainder signal r = s − M k=1 IMF k . The computation of s m+1 is the critical step of the algorithm, which is now made fast significantly by exploiting the FFT. When r becomes the trend signal, i.e. it does not contain any more oscillations, the algorithm stops.
The FIF method, which is originally an iterative procedure, can be made to be a direct method, based on equation (7). In order to do so, it needs to compute a priori the number of iterations N 0 required in the inner loop of the algorithm to extract each IMF. However this task is not easy at all and currently the only possibility proposed in the literature [18] is a formula which allows only to upper bound the exact value of N 0 . For this reason in [18] an alternative approach is proposed. In particular in [18] it is pointed out that the purpose of the inner loop iterations in the FIF method is to send, after N 0 iterations, all eigenvalues in I − D that are close to zero to be zero, up to machine precision, and keep the eigenvalues bigger enough than zero unchanged [18], [40]. Therefore, the idea is to send the eigenvalues of matrix I − D that are smaller than the threshold τ to 0, while preserving the other eigenvalues unchanged. This is a direct method, and it is called the hard thresholding fast iterative filtering (HTFIF).

B. COMPARISON OF DIFFERENT MODE DECOMPOSITION METHODS
To validate the effectiveness and superiority of the innovative HTFIF algorithm, the excellent signal decomposition methods newly proposed in recent years, such as time-varying filter EMD (TVF-EMD) and variational mode decomposition (VMD) are employed to compare with HTFIF. Besides, HTFIF is developed on the basis of MIF. The comparison between MIF and HTFIF is also conducted.
Given an artificial signal f 1 (t), 0 ≤ t ≤ 5, which is composed of an intermittent signal and two harmonics. Perturbed by the white noise n(t) ∼ N (0, t), the signals are considered and expressed as follows: The mixed-signal f 1 (t) and its real components are displayed in fig. 1. The decomposition methods of the VMD with the number of the mode K to be 4, TVF-EMD with the bandwidth threshold to be 0.4 and the B-spline order to be 20, MIF with the filter length adjustment parameter ζ to be 4 and HTFIF with the filter length adjustment parameter ζ to be 4 and the threshold τ to be 0.8 are performed on f 1 (t) respectively. The decomposition results are displayed in Fig. 2. It can be clearly seen that severe mode mixing occurs in the signals decomposed by VMD. TVF-EMD, MIF and HTFIF are all able to separate the real components from the mixed-signal with the first and second IMF revealing the noise and intermittent signal respectively, and the last two IMFs corresponding to the harmonics s 2 and s 1 respectively. The modes decomposed by the three methods are almost equal that is hard to judge which one is better from the waveform.
To further explore the decomposition performance of TVF-EMD, MIF and HTFIF, a relative error introduced in [41] is used here. And for the harmonic component s 2 , the relative error can be expressed as follows: The decomposition results are shown in Table. 1. It indicates that the HTFIF decomposition method has higher accuracy.   What's more, the operating time of the three methods above is computed. The experimental hardware conditions are as follows: Intel(R) Core(TM) i5-8265 CPU @1.6 GHz and 8 GB RAM. The soft-ware conditions are as follows: Windows 10, 64-bit operating system, X64-based processor, and MATLAB R2017b. The calculation results are shown in Table. 2. Obviously, the running time of HTFIF decomposition is 0.0128s, far less than that of TVF-EMD and MIF, which indicates that HTFIF has excellent efficiency in signal decomposition.

C. L-KCA INDICTOR FOR SENSITIVE IMF SELECTION
To choose the sensitive feature component from a series of IMFs that decomposed by HTFIF, a fusion indicator L-KCA on the basis of L-kurtosis, autocorrelation function impulse harmonic to noise ratio index (ACFHNR) and correlation coefficient (Corre. Coe.) is developed in this research.

1) L-KURTOSIS
Given a real valued random variable X , let its associated cumulative distribution function be F(x) = P(X ≤ x) and the inverse function of F(x) be x(F). Let X 1:n ≤ X 2:n ≤ ···X n:n be the statistics of the random sample X with size n in ascending order. Hosking [42] defined the rth L-moment of X as: where E (X r−k:r ) is described as: Then the first four order L-moments are calculated as follows: According to the above results, L-kurtosis can be obtained as: The energy of the autocorrelation function of the general periodic signal is periodically distributed, while that of the white noise signal is mainly concentrated in the zero point and becomes zero at the non-zero points. Therefore, Zheng et al. [33] proposed a new autocorrelation function periodic impulse harmonic to noise ratio index (ACFHNR) to detect the periodic pulse fault characteristics. The calculation process is as follows: Let x be the Hilbert transform of the original signal x(t) with the length of N , and y be the absolute value of x. And they can be calculated as follows: Then remove the direct current component of y. It should be note that for the purpose of nondimensionalize, during the process, y with its average removed need to divide by its standard deviation, which can be expressed as: Next, compute the autocorrelation function R(µ) of s(t): where E{·} is the expectation operator and s(n) * is the conjugation of s(n). Finally, the ACFHNR index is described as follows: where N T is the number of the sampling points in a pulse period, R en (kN T ) is the energy value of the autocorrelation function R(µ) at kth period impulse, R en (0) is the energy value of the autocorrelation function R(µ) at the zero point. And k = 3 is suggested [31].

3) L-KCA INDICATOR
Compared with kurtosis, L-kurtosis is robust to outliers and more sensitive to the fault pulse sequence itself, but it still mainly focuses on the general statistical distribution of a signal and may ignore the specific characteristics of the vibration signals. The ACFHNR index can indicate the ratio of period impulse harmonics and noise components in the vibration signals, but it is sensitive to interference harmonics. According to the correlation coefficients between each IMF and the original vibration signal, the IMFs can be divided into clusters of sensitive and insensitive IMFs [33]. Make the above three indicators have complementary advantages to obtain a fusion index L-KCA for the sensitive IMF selection, which is calculated as follows: First, the L-kurtosis and ACFHNR values of each IMF and the correlation coefficient between each IMF and the original vibration signal are calculated. Then the following Min-Max scaling formula (22) is used to normalize them in the same scale: Finally, the L-KCA indictor of each IMF is calculated as: where λ L−kurtosis (k), λ Corre.Coe. (k) and λ ACFHNR (k) denote the normalization value of L-kurtosis, correlation coefficient and ACFHNR of the kth IMF, respectively. Obviously, the limit [Min, Max] of L-KCA indicator is set to [0, 1], and the larger the L-KCA indicator of the IMF is, the more periodic impact features is contained in it.
In order to illustrate the effectiveness and superiority of the proposed sensitive IMF selection indicator, the comparisons of L-KCA with kurtosis, L-kurtosis, ACFHNR and correlation coefficient are carried out. The tested signal S(t) can be described as follow: where x(t) is the fault transient impulse sequence and n(t) is the noise signal. The SNR of the mixed-signal S(t) is set to −1dB to −15dB. The indicators of the sensitive IMF decomposed by HTFIF at different SNR are shown in Fig. 3. It can be seen that when the input SNR of the mixed signal changed from −1dB to −15 dB, all the indicators are gradually decreased, which indicates that all the indicators can distinguish the sensitive IMF. The L-KCA index has the best resolution among the indicators, indicating that the proposed L-KCA index outperforms other indicators in identifying the periodic impulse feature of the rolling bearings. In addition, it is worth noting that the kurtosis indicator changes suddenly at higher noise level, which is due to the fact that the kurtosis metric is susceptible to outliers caused by random events. This could mislead the detection decision. L-kurtosis is more robust and reliable than the traditional kurtosis, which is the reason why we choose L-kurtosis instead of kurtosis in the fusion index L-KCA.

D. IMPROVED k-VALUE SYMMETRICAL DIFFERENCE ANALYTIC ENERGY OPERATOR
The concept of analytic energy operator is developed from the analytic signal. Suppose a modulated signal: whose analytic form is: where A (t) is the instantaneous amplitude. Let x be the Hilbert transform of x(t), then the instantaneous amplitude can be obtained by: ϕ(t) is the instantaneous phase, which can be expressed as: The instantaneous frequency can be expressed by the first derivative of the instantaneous phase, that is: Hence, Obviously, it includes amplitude demodulation and frequency demodulation, so it is called analytical energy operator (AEO) and represented by {x(t)}, that is: where x (t) is the first derivative of the signal x(t); x (t) is the first derivative of the Hilbert transform of x(t). The discrete form of equation (30) can be expressed as: At the same time, for the discrete time signal x(n), Teager energy operator adopts the forward difference • x(n) = x(n + 1) − x(n) to approximate the time derivative. However, the maximum relative error can reach 11% after taking the derivative by the forward difference, resulting in low demodulation accuracy [34]. Therefore, to smooth the original discrete signal and improve the accuracy of demodulation, instead of forward difference the central finite difference sequence is adopted in the symmetrical difference analytic energy operator (SDAEO), which is expressed as: In order to enhance the noise robustness of the energy operator, inspired by the multi-resolution Teager energy operator (MTEO) [36], [37]. This paper uses three sampling points that spaced k apart in the signal to improve the symmetrical difference sequence and obtains: Then the improved k-value Symmetric Difference Analytic Energy Operator (k-SDAEO) is constructed by combining Equation (31) with Equation (33), and expressed as: where h [n] = H [x (n)], represents the discrete Hilbert transformation.

E. THE SNR IMPROVEMENT OF THE k-SDAEO
The demodulation performance of the optimized k-SDAEO is hard to evaluate due to its nonlinearity. It is clear that higher SNR improvement (SNRI) brings about better demodulation results. Therefore, we use SNRI indicator to evaluate the demodulation performance of k-SDAEO indirectly. Let s(n) with a sample size N be the original signal and ν(n) be the signal denoised by k-SDAEO, and we define the SNRI of K-SDAEO as follows:

F. THE PROCEDURE OF THE PROPOSED HTFIF-k-SDAEO METHOD
Motivated by the advantages of HTFIF, L-KCA indictor and the improved k-SDAEO demodulation technology, a new HTFIF-k-SDAEO method for rolling element bearing fault diagnosis is proposed. The detailed procedure of the HTFIFk-SDAEO is illustrated in Fig.4 and the explanation is as follows: Step 1: Decompose the collected vibration signal into a series of IMFs by HTFIF; Step 2: Calculate the L-KCA indictor of each IMF and select the one with the maximum L-KCA as the sensitive IMF; Step 3: Calculate the SNRI of each k-SDAEO of the sensitive IMF at different resolution parameter k and select the one with the maximum SNRI as the optimal k-SDAEO; Step4: Perform the optimal k-SDAEO on the sensitive IMF component to extract the fault characteristics of the bearing.

III. NUMERICAL SIMULATION
Due to the vibration signal acquired from a defective bearing is usually a mixture of the bearing fault-induced signal, the interferences are transmitted from other components and background noise. To verify the proposed method, the fault signal with heavy background noise and harmonic interferences was created by Equations (37), (38) and (39).
where I (t) is the interference signal; A h is the amplitude of the interference harmonics; H is the total number of the each interference component; K is the number of harmonics of the ith interference; B ij is the amplitude of the jth harmonic of the ith interference;f hi is the fundamental frequency of the ith interference, here three vibration interference frequencies (73Hz, 135 and 183Hz) are added. x(t) is the transient impulse fault signal; A m is the amplitude of the fault impulse and is set to be 3; α is the damping coefficient and is set to be 0.1; f n with a value of 2,000Hz is the resonance frequency. The fault characteristic frequency is f c = 105Hz. n(t) is the additive white Gaussian noise; here the SNR of the fault signal is set to −9.5dB. The number of the sampling points of each simulated signal is set to be 12,000 and the sampling frequency is set to be 12kHz. The heavily contaminated signal s(t) and its corresponding envelope spectrum are displayed in Fig. 5. It is found that the time-domain signal is dominated by noise, and the envelope spectrum cannot reveal any information about the impulse fault features.
To increase the SNR of the spectrum and enhance the spectral peaks of the fault characteristic frequency, the TEO, SDAEO and k-SDAEO are performed on the simulated signal respectively, and the results are shown in Fig. 6. Since the fault characteristic frequency can hardly be extracted when k is greater than 5 and k-SDAEO degenerates to SDAEO when k=1, this research only takes the cases of k = 2, 3, 4 and 5 into consideration. The calculated SNRIs of SDAEO and k-SDAEO at different k are listed in Table. 3. As shown in Fig. 6, it is difficult for these energy operators to extract the fault impulse characteristics from strong background noise. But it can be seen that 2-SDAEO is able to identify several harmonics of the fault characteristic frequency and is slightly better than the other energy operators, which is consistent with the case that the value of SNRI of k-SDAEO at k=2 is slightly higher.
Hence, the proposed method is performed on the simulated signal to identify the fault impulse features. Firstly, HTFIF is employed to decompose the simulated signal s(t), and eleven components are obtained as shown in Fig. 7.   Secondly, in order to select the sensitive component which contains rich fault feature information, the value of the kurtosis, correlation coefficient, ACFHNR, and the developed L-KCA indicator for each IMF component are calculated respectively, as shown in Fig. 8.
Obviously, the kurtosis, ACFHNR and correlation coefficient values of IMF1 are respectively larger than the ones of other components, while IMF2 has maximum L-KCA value. To confirm the effective component, the envelop analyses of IMF1 and IMF2 are carried out, as shown in Fig. 9. It can be found that the periodic impulse features in the envelop spectrum of IMF1 are totally submerged in background noise. However, several harmonics of fault characteristic frequency can be identified from the envelop spectrum of IMF2, which shows that IMF2 is the sensitive component and the developed L-KCA indicator is superior to the kurtosis, ACFHNR and correlation coefficient indexes in selecting the sensitive IMF.   Compared the envelope spectrum of the mixed signal ( Fig. 5(b)) with that of the optimal IMF2 ( Fig. 9(b)), it can be found that after HTFIF processing, the noise has been reduced and the impact characteristics have been enhanced. Several harmonics of the fault characteristic frequency can be identified, but the fault characteristic frequency itself is still submerged in noise. To select the optimal resolution parameter k, the SNRI of each k-SDAEO is calculated as shown in Table. 4. It is evident that the maximum SNRI value appears at k = 2, so the optimal resolution parameter k is set to be 2 here. Therefore, TEO, SDAEO and 2-SDAEO are performed on the sensitive IMF2 respectively and the results are shown in Fig. 10. It can be seen that compared with TEO and SDAEO, 2-SDAEO is able to extract more obvious fault characteristic frequencies and more harmonics, which gives another proof that k-SDAEO outperforms TEO, SDAEO.
In order to further prove the superiority of the proposed method, the excellent signal decomposition methods TVF-EMD and VMD newly proposed in recent years are employed for comparison and analysis. Six IMFs are obtained by VMD with the number of the mode K to be 6 and fourteen IMFs are obtained by TVF-EMD with the bandwidth threshold to be 0.25 and the B-spline order to be 27. The third IMF in VMD and the second one in TVFEMD are respectively selected as the sensitive IMF according to the maximum L-KCA indicators, as shown in Fig. 11. Then the SNRIs of k-SDAEO of the sensitive IMFs for the two methods are calculated as shown in Table. 5. 4-SDAEO and 2-SDAEO are the optimal k-SDAEO for VMD and TVFEMD respectively according to the maximum SNRI.
Perform the optimal k-SDAEO on the sensitive IMFs derived from VMD and TVF-EMD and the results are illustrated in Fig.12. It can be seen that although the fault  characteristic frequency and its some harmonics could be identified by the methods based on VMD and TVFEMD. However, compared with the analysis results of the proposed method ( Fig. 10(c)), the observed periodic pulse characteristics are not so clear due to the strong noise contamination.
For the purpose of quantitative evaluation, this paper selects the proportion of the spectrum energy as an evaluation function, which is expressed as: where f c is fault characteristic frequency; i represents the harmonic order of the fault characteristic frequency; m is the total number of the harmonics of the fault characteristic frequency that can be found from the spectrum within a frequency band (0∼1000Hz); ε is the error caused by bearing mounting, E if c +ε is the frequency band energy sum of observable fault characteristic frequency and its harmonics; E is the total sum of the frequency band (0∼1000Hz); η is the proportion of the spectrum energy, and the greater the value, the better the detection results. Then the proportion of the spectrum energy of the HTFIFk-SDAEO, VMD-k-SDAEO and TVFEMD-k-SDAEO method are calculated to be 3.02%, 2.13% and 2.68% respectively, which verifies that the proposed HTFIF-k-SDAEO method has the best detection result.
Finally, MATLAB 2017b was used as the operating software to compare the computing efficiency of HTFIF-k-SDAEO, VMD-k-SDAEO and TVFEMD-k-SDAEO on a laptop configured with Intel Core i5-8265U CPU and 8.0GB RAM under windows10 and 64-bit operating system. The computation time of the three methods is shown in Table. 6. It can be seen that the calculation time of HTFIF-k-SDAEO is 0.0989s, which is significantly smaller than that of VMDk-SDAEO (38.9011s) and TVFEMD-k-SDAEO (40.3258s).
Therefore, it can be seen from the above analysis that although TEO, SDAEO and k-SDAEO have   high computation efficiency, they cannot identify periodic impulses from the composite signals containing heavy noise and harmonic interferences. VMD-k-SDAEO and TVFEMDk-SDAEO can extract the periodic impulse features, but the effects are not obvious and are time-consuming. However, the proposed HTFIF-k-SDAEO method not only has strong feature extraction performance, but also has high computational efficiency.

IV. EXPERIMENTAL EVALUATIONS
To further validate the effectiveness and superiority of the proposed HTFIF-k-SDAEO method in practical fault diagnosis, the experiments of the vibrating screen bearings are conducted. In this part, two types of the bearing faults are considered: the outer race fault and the inner race fault.

A. EXPERIMENTAL SETUP
The experiments are performed on the multi-function vibrating screen (SDM00) with two shafts and two motors, in mechanical vibration laboratory of Xi'an University of Architecture and Technology. The DYTRAN piezoelectric accelerometers with the sensitivity of 97.3 MV/g and sampling frequency of 20 kHz are adopted. The experimental setup is shown in Fig.13. Type 1308 aligning ball bearings are applied to these real tests and their specifications are displayed in Table.7. The sampling frequency of the signals is 20 KHz and 20000 points are collected for the outer and inner race analysis. The bearing speed is n = 1000 r/min and the rotating frequency is f r = 16.7 Hz. According to the formulas of bearing theoretical fault frequency, the fault characteristic frequencies of the outer race and inner race are calculated as f o = 104.18 Hz and f i = 145.87 Hz, respectively.

B. FAULT DETECTION FOR THE OUTER RACE
The time-domain waveform of the outer race fault signal and its corresponding envelope spectrum are presented in Fig.14 (a) and Fig.14 (b), respectively.
From the time domain diagram we cannot see any fault impulse characteristics with uniform interval. The envelope spectrum is abound with interference spectral lines and noise, which makes the fault features difficult to be extracted. Then, VOLUME 9, 2021    perform the energy operators TEO, SDAEO and k-SDAEO on the original vibration signal for demodulate analysis and their corresponding energy spectra are obtained as shown in Fig.15. According to the maximum SNRI (Table.8), k is set to be 4.
As shown in the TEO spectrum, the fault features are totally submerged in the heavy background noise. SDAEO cannot extract the outer race fault characteristic frequency f o although its second and fourth harmonics can be revealed. 4-SDAEO can both extract the outer race fault characteristic frequency and its several harmonics. Therefore, among these energy operators, 4-SDAEO has the best fault extraction performance, which verifies the superiority of the proposed k-SDAEO method. However, there are also some interference frequencies and background noise in the 4-SDAEO spectrum and the order of the extracted harmonics are discontinuous.  Therefore, the proposed method is performed on the outer race fault signal. Firstly, the HTFIF is employed to decompose the vibration signal adaptively, and ten components are obtained as shown in Fig.16.
To select the sensitive IMF, the values of kurtosis, correlation coefficient and L-KCA indictor for each IMF component are calculated respectively as shown in Fig.17. It can be found that the components with the maximum values of kurtosis, correlation coefficient and L-KCA indictor are IMF7, IMF1 and IMF2, respectively. As can be seen from fig.16, IMF7 does not contain any periodic impact characteristics. The envelop spectra of IMF1 and IMF2 are illustrated in Fig.18. It is easy to discover that IMF2, which contain rich According to the maximum SNRI of k-SDAEO for IMF2 (Table.9), k is set to be 3. Then the sensitive IMF2 is processed by TEO, SDAEO, and 3-SDAEO and the results are illustrated in Fig.19. It can be seen that the fault characteristic frequency and the associated harmonics can be detected from the IMF2 by TEO, SDAEO and 3-SDAEO all. But compared with 3-SDAEO spectrum, the amplitudes are small in SDAEO spectrum and the harmonic orders are discontinuous in TEO spectrum. Accordingly, the k-SDAEO outperforms TEO and SDAEO under the same conditions.
To further verify the advantages of the HTFIF-based method, the vibration signal is also decomposed by VMD and TVF-EMD. Seven IMFs are obtained by VMD with the number of mode k to be 7 and eight IMFs are obtained by TVFEMD with the bandwidth threshold to be 0.28 and the B-spline order to be 26. The third IMF in VMD and the fourth IMF in TVFEMD are respectively selected as the sensitive VOLUME 9, 2021    IMFs according to the maximum L-KCA indicators, as shown in Fig.20. The SNRI of k-SDAEO for the selected IMFs obtained by the two methods are shown in Table.10.
It is can be seen that both VMD and TVFEMD have maximum SNRI at k=5. The 5-SDAEO spectra of the selected IMFs for the two methods are shown in Fig.21. It is clear to be seen that the results are inferior to the 3-SDAEO spectrum of the sensitive IMF obtained by HTFIF (Fig.19(c)).
For the purpose of quantitative evaluation, the proportion of the spectrum energy of the HTFIF-k-SDAEO, VMD-k-SDAEO and TVFEMD-k-SDAEO method are calculated to be 16.6%, 4.58% and 7.06% respectively, which verifies that the proposed HTFIF-k-SDAEO method has the best detection result.
Finally, MATLAB2017b was used as the operating software to compare the computing efficiency of HTFIF-k-SDAEO, VMD-k-SDAEO and TVFEMD-k-SDAEO on a laptop configured with Intel Core i5-8265U CPU and 8.0GB RAM under windows10, 64-bit operating system. And the computational time of the three methods is shown in Table.11. It can be seen that the calculated time of HTFIF-k-SDAEO is only 0.0609s which is far less than that of VMD-k-SDAEO (17.1528s) and TVFEMD-KSDAEO (34.3194). Therefore, the proposed HTFIF-k-SDAEO method can be used for real-time monitoring and on-line fault diagnosis.

C. FAULT DETECTION FOR THE INNER RACE
The time domain waveform of the inner race fault signal and its corresponding envelope spectrum are presented in Fig. 22(a) and Fig. 22(b), respectively. Similarly, the fault impulses cannot be seen from the time domain signal, and the envelope spectrum of the measured signal is dominated by noise and interferences.    As before, according to the SNRI of each k-SDAEO (Table. 12) the optimal k-SDAEO presents at k=5.
Perform TEO, SDAEO and 5-SDAEO on the raw signal and the results are illustrated in fig. 23. It clear to be seen that TEO spectrum is dominated by noise. In the SDAEO spectrum only the fourth and sixth harmonics of the fault characteristic frequency can be found. 5-SDAEO can both identify the inner race fault characteristic frequency f i and its 5-fold and 6-fold frequency, which is better than the results of the above two energy operators. Nevertheless, there are also too many interference spectral lines in 5-SDAEO spectrum and the orders of the extracted harmonics are discontinuous.
Then, the HTFIF, VMD and TVFEMD are applied to decompose this vibration signal and the L-KCA indictors of each IMF component for the three methods are calculated as shown in Fig. 24. It can be seen that the components with the maximum value of L-KCA are IMF2, IMF3, and IMF2, respectively for the three methods. According to the maximum SNRI of k-SDAEO for the sensitive IMF (Table.13), the optimal resolution parameter k is set to be 2, 1, and 2 VOLUME 9, 2021  respectively for the three methods. It should be noted that k-SDAEO degenerates to SDAEO when k = 1.
The sensitive IMFs for the three methods and their corresponding optimal k-SDAEO spectra are illustrated in Fig. 25. It can be seen that both VMD-and TVFEMD-based methods can extract several harmonics of the fault characteristic frequency with small amplitudes and accompanied by strong noise and harmonic interferences. However, the fault characteristic frequency and its harmonics can be identified evidently from the 2-SDAEO spectrum of IMF2 in HTFIF. In addition, the execution time has to be taken into account. Thus, the proposed method is more likely to be further appreciated in practice.

V. CONCLUSION
A novel rolling element bearing fault diagnosis approach HTFIF-k-SDAEO is proposed in this study. It consists of three steps: hard thresholding fast iterative filtering (HTFIF) adopted for signal decomposition, L-KCA indicator developed for sensitive IMF selection and the improved k-value symmetrical difference analytic energy operator (k-SDAEO) presented for demodulation analysis. The results of numerical simulation and experiments demonstrated the superiority of the proposed method. The main conclusions are summarized as follows: (1) HTFIF was first applied to the bearing fault diagnosis. The comparison of HTFIF with VMD, TVFEMD and MIF based on simulation signal highlighted that the HTFIF decomposition method has higher decomposition accuracy and more favorable decomposition efficiency. Furthermore, to select the sensitive IMF, the L-KCA indicator combining the L-kurtosis, correlation coefficient and ACFHNR was developed and verified to be superior to the indexes of kurtosis, ACFHNR and correlation coefficient. Finally, a k-SDAEO with an optimal k value selected according to SNRI for different signals was presented to extract the bearing fault features from the sensitive IMF. (2) Through numerical simulation and experimental studies, the performance of the proposed method is compared with that of the non-preprocessing energy operator demodulation, HTFIF-TEO, HTFIF-SDAEO, VMD-and TVFEMD-based methods. The results show that the proposed method achieves a better fault feature extraction capability and higher efficiency and it is more likely to be an effective tool for bearing fault real-time detection and on-line diagnosis.
In addition, some parameters of the algorithms in this paper are the choices when the best decomposition results are obtained after repeated debugging according to experience.
In the succeeding research, we will study some parameter optimization algorithms to realize the automatic adjustment of the parameters and improve the engineering applicability of the algorithms.