An Optimal Parameter Selection Method for MOMEDA Based on EHNR and Its Spectral Entropy

As a vital component widely used in the industrial production field, rolling bearings work under complicated working conditions and are prone to failure, which will affect the normal operation of the whole mechanical system. Therefore, it is essential to conduct a health assessment of the rolling bearing. In recent years, Multipoint Optimal Minimum Entropy Deconvolution Adjusted (MOMEDA) is applied to the fault feature extraction for rolling bearings. However, the algorithm still has the following problems: (1) The selection of fault period T depends on prior knowledge. (2) The accuracy of signal denoising is affected by filter length L. To solve the limitations, an improved MOMEDA (IMOMEDA) method is proposed in this paper. Firstly, the envelope harmonic-to-noise ratio (EHNR) spectrum is adopted to estimate the fault period of MOMEDA. Then, the improved grid search method with EHNR spectral entropy as the objective function is constructed to calculate the optimal filter length used in the MOMEDA. Finally, a feature extraction method based on the improved MOMEDA (IMOMEDA) and Teager-Kaiser energy operator (TKEO) is applied in the field of rolling bearing fault diagnosis. The effectiveness and generalization performance of the proposed method is verified through comparison experiment with three data sets.


Introduction
Rolling bearings are vital components of mechanical systems and have been widely used in rotating machinery. However, rolling bearings often work in more demanding environments for a long time, such as corrosive, high temperature, large impact and high load conditions, which increases the uncertainty of the mechanical system, resulting in various fault defects of the bearings [1]. In order to avoid the environmental pollution, economic loss and even safety problems caused by rolling bearing shutdown, the fault diagnosis method based on signal analysis is widely used in industrial production, in which the bearing condition is recognized according to the effective information collected by the vibration signal [2,3], pressure signal [4], temperature signal [5] and acoustic signal [6,7].
The fault diagnosis technique based on vibration signal has been widely used in practice because of its intuitive operation and reliable diagnosis results, where the vibration signal of the bearing is monitored by the acceleration sensors installed in the appropriate direction of the bearing seat or the box, and the signal is analyzed and processed to determine the bearing condition. Furthermore, some advanced measurement systems are applied in rolling bearings. Wang [8] designed a data acquisition system of high speed railway bearing based on ZigBee technology. Gao et al. [9] designed a wireless sensor vibrating signal acquisition and processing system, and the researchers could collect and analyze the signal timely and conveniently. Li et al. [10] proposed a smart metasurface shaft (SMST) to achieve single-sensor identification of vibration sources, which showed potential applications in rotating machinery condition monitoring. Verstraete et al. [11] presented multi-sensor fusion systems for integrated remaining useful life prognostic capabilities. Meng et al. [12] developed a fully enclosed bearing-structured self-powered rotation sensor (SPRS) for multi-tasking motion measurement. Furthermore, the device displayed high precision, broad range, quick response, and excellent durability.
Fault feature extraction based on signal analysis is critical after the reliable vibration data source is available. In industrial production, when the bearing fails to work, the collected vibration signal will generate periodic transient pulses at a certain characteristic frequency [13]. However, in fact, in addition to the pulses caused by fault defects, there are also some modulation harmonics caused by external excitations, such as axial unbalance, misalignment, critical speed, structural resonance and uncertain factors on vibration signal transmission paths, which make the spectrum more complex and increases the difficulty of bearing fault diagnosis [14]. Therefore, recovering and enhancing effective fault period impulse components from vibration signals is of great significance and importance for fault feature extraction. Some signal processing methods have been successfully employed in capturing fault-related transient components in the non-stationary signal for rolling bearings [15], such as spectral kurtosis [16,17], wavelet transform theory [18] and so on. Unfortunately, most of the methods mainly concentrate on emphasize random impulse with large amplitude and thus ignores the periodicity of fault signal.
To tackle the limitation, the deconvolution principle has been introduced for bearing fault periodic composition extraction, which achieved satisfactory results. In this line, Wiggins [19] proposed the minimum entropy deconvolution (MED) in 1978, which took kurtosis as the objective function and the optimal inverse filter was iteratively found to enhance the transient ingredients in the original signals. In 2007, Randall successfully applied the MED to fault detection [20]. Since then, the MED has been widely used in the fault feature extraction of rotating machinery [21][22][23].
However, the MED method aims at maximizing the kurtosis value, and thus it is easily affected by strong noise. In the iterative process, only the pulse component in the vibration signal can be enhanced and this defect limits the development of MED method [24]. To address this issue, McDonald et al. [25] proposed the maximum correlation kurtosis deconvolution (MCKD), which took the correlation kurtosis as the objective function and iteratively searched for the optimal reverse filter to enhance the periodic transient components in the original signal. The MCKD algorithm has achieved remarkable results in the fault feature extraction of rotating machinery such as rolling bearings [26,27], gearboxes [25], and motor rotors [28]. Nonetheless, there are still some issues in the practical application of MCKD. In fact, its accuracy is affected by four parameters: filter length, fault period, the number of displacement and the number of iterations [29]. Moreover, the MCKD can only address the single period pulse signal.
To extract the continuous periodic pulse signal, McDonald et al. [30] proposed the multipoint optimal minimum entropy deconvolution adjusted (MOMEDA) method. This method is a non-iterative deconvolution process and utilizes infinite pulse sequence as the optimization object function to search for the optimal filter directly. With these advantages, the MOMEDA has been successfully applied to the rotating machinery fault characteristics extraction [31,32]. Ma et al. [33] introduced MOMEDA to detect planet bearings faults and the experiments show that MOMEDA presented specific applicability compared with the traditional kurtogram filtering method. Li et al. [34] suggested a method based on the MOMEDA and long and short time memory (LSTM) for fault diagnosis. Zhu et al. [35] introduced a method based on the MOMEDA and TKEO for the rolling bearing fault diagnosis with good experimental qualities. However, the following weaknesses of the MOMEDA method are also well noted in the community: (1) the performance of the MOMEDA mainly depends on the prior knowledge of fault period T. It is necessary to identify the periodic impulses signal through the multipoint kurtosis spectrum [30] in the default fault period T search interval. However, when the original signal is affected by complex background noise, the multipoint kurtosis spectrum may not provide correct fault period information according to the location of maximum multipoint kurtosis. (2) The effect of noise reduction for MOMEDA is also affected by the filter length L. Hence, unreasonable filter length L will weaken or enhance the energy of the original signal, and then cause misdiagnosis. To address the above-mentioned issues, scholars have introduced some optimization algorithms [36,37] to determine the fault period T and filter length L of MOMEDA, such as particle swarm optimization (PSO) algorithm [27,38], grasshopper optimization algorithm (GOA) [39] and so on. These parameter optimization algorithms mostly take kurtosis, envelope spectrum kurtosis (ESK) and other similar indexes as objective functions to improve MOMEDA. However, the methods also have some limitations. These algorithms are usually applied to multi-objective optimization, i.e., training parameters such as fault period T, filter length L and window function simultaneously, which make the implementation process complex and time-consuming. Therefore, it is necessary to explore other effective methods to quickly calculate the fault period T and filter length L of MOMEDA for improving the deconvolution performance.
An improved MOMEDA (IMOMEDA) method is investigated for fault feature extraction in this paper. The proposed method incorporates the envelope harmonic-to-noise ratio (EHNR) [40] into the MOMEDA scheme. By calculating the EHNR spectrum of the original signal, the fault period T can be quickly determined without using the prior information. Furthermore, kurtosis spectral entropy is proposed as a grid search objective function to search for the optimal filter length L of MOMEDA, for which the single objective optimization can greatly simplify the procedures [24]. However, kurtosis spectral entropy tends to emphasize random impulse with large amplitude (it inherits the essential characteristic of kurtosis) and thus ignores the periodicity of fault signals. Therefore, EHNR spectral entropy is presented as the target function to overcome this shortcoming and obtain the improved grid search method to select the optimal filter length L of MOMEDA. Compared with the existing methods, the overall approach has the following advantages: (1) A new index-based EHNR spectrum is proposed to calculate fault period T of MOMEDA and require less prior knowledge and calculation time.
(2) EHNR spectral entropy is introduced to construct the objective function of grid search method and select the optimal filter length L of MOMEDA.
(3) A novel feature extraction method-based IMOMEDA-TKEO is presented and applied in the field of rolling bearing fault diagnosis.
(4) The effectiveness and generalization performance of the presented method is verified through comparison experiment with three data sets.
The remaining part of the papers is organized as follows: Section 2 introduces the basic principles of MOMEDA and EHNR. Section 3 elaborates on the implementation process of IMOMEDA-TKEO. Section 4 presents the experimental results and analysis. Section 5 gives the discussion and summary.

MOMEDA Method
Suppose x is the vibration signal, y is the impulse sequence, h is the transfer function, e is noise, and "*" is the convolution product. The discrete-time response signal is as follows [41]: The MOMEDA algorithm aims at finding an optimal filter f in a non-iterative way to realize the reconstruction of the original impulse sequence y and try to reduce the impact of noise. The major process of MOMEDA algorithm is given as follows: Or in matrix form: Among them, k = 1, 2, . . . , N−L, according to the periodicity of the fault signal for rotating machinery, MOMEDA takes the multi D-norm of the filtered signal as the objective function, the concept of MOMEDA is presented in Equations (5) and (6): The target vector t is a constant vector that defines the location and weightings of the goal impulses to be deconvolved, which is suitable for the features extraction [19]. When t is completely consistent with original impulse signal y, the multi D-norm reaches the maximum value and the deconvolution effect works the best. Therefore, the non-iterative process of Equation (6) is formulated as follows: the Equation (8) can be simplified as: We then substitute y = X 0 T f into Equation (10), and then calculate: The f = (X 0 X 0 T −1 X 0 t can be denoted as a set of optimal filters, and then reconstruct the original impulse signal y: Finally, the multipoint kurtosis (MKurt) is defined as follows: However, the construction of target vector t needs to take T as a prior knowledge. In practical engineering, it is usually impossible to meet this condition due to the complex working environments, which in turn limits the application of the MOMEDA method. In addition, the periodic pulse energy is determined by the filter length L. If the parameter setting is unreasonable, it will weaken or enhance the energy of the periodic impulse signal which will lead to misdiagnosis. Therefore, how to calculate the fault period T and filter length L for the MOMEDA method quickly is one of the major issues, aiming at improving its performance.

EHNR Principle
In the field of condition-based maintenance, a large number of indicators are introduced to determine whether the measured signal is informative or not. These indicators are widely used, including kurtosis, entropy to L2/L1 norm, Gini index and other norm indexes [42]. Most of the indexes mainly concentrate on the peak value, non-gaussianity distribution and other general statistical distribution, thus ignoring the particularity of mechanical signals. Therefore, this paper introduces the EHNR method. According to the specific characteristics of rolling bearing fault signal, we can determine the parameter T of MOMEDA quickly.
The calculation process of EHNR is as follows [40]: (1) The envelope signal Env x (t) of vibration signal x(t) is obtained by Hilbert transform and the DC component is removed.
(2) The autocorrelation function of Env x (t) is calculated as: (3) The maximum position of the autocorrelation function of the signal is found in the lag domain, where the point is shown in Figure 1.  Then, the EHNR can be defined as Equation (18): where τmax is the time lag that makes the autocorrelation spectrum of Envx(t) reach its local maximum,  Then, the EHNR can be defined as Equation (18): where τ max is the time lag that makes the autocorrelation spectrum of Env x (t) reach its local maximum, r Env x (τ max ) is the energy of harmonic, r Env x (0) is the total energy of envelope signal.

Setting Filter Length L Based on Grid Search Method
The grid search [43] is an exhaustive search method for specifying the accurate parameter values. According to the definite step size, the grids are divided in a certain range and the constraint function values on each grid point are calculated. Then, the independent variables that make the objective function value optimal (maximum or minimum) are selected as the solution of the problem.
In this paper, the EHNR spectral entropy is proposed as the objective function and the expression is shown in Equation (19): Among them, where S R (ω) is the envelope spectrum of the original signal, N is the data length, EHNR(x) represents the EHNR value of the input signal, and EHNR_E(x) represents the EHNR spectral entropy of the signal. For the periodic pulse signal, the envelope spectrum is concentrated on the low-frequency region which results in a smaller E s value. EHNR(x) is a function to describe the periodic change of faults. The larger the EHNR(x) is, the stronger the periodicity. Therefore, EHNR_E(x) can reflect the uniformity of the signal period. In this paper, EHNR_E(x) is used as the objective function of grid search method, and then the optimal filter length L is determined.

TKEO Algorithm
The Teager-Kaiser energy operator is a nonlinear operator that can track the instantaneous energy of signal. The extraction of envelope by TKEO is to calculate the three adjacent sampling points of the measured waveform. TKEO can track the change of the measured signal waveform in real time. The algorithm has the advantages of excellent time resolution, simple and fast calculation. The specific mathematical expression and implementation process of the TKEO algorithm was proposed in [44].

The Proposed Method
Inspired by the above-introduced methods, in this section, an IMOMEDA-TKEO rolling bearing fault feature extraction method is proposed. The framework of the proposed method is shown in Figure 2. The general implementation procedure is summarized as follows: Step 1: The autocorrelation function r Env x (τ) = Env x (t)Env x (t + τ)dt of rolling bearing vibration signal x(t) is obtained by the autocorrelation analysis.
Step 2: The envelope harmonic-to-noise spectrum EHNR(x) of signal x(t) is obtained and the fault period T is calculated.
Step 3: According to [45], the filter length L meets the principle of L > 2f s /f c , where f s is the sampling frequency and f c is the resonance frequency. We initialize the input parameters of the grid search, set the search range of the filter length L to (250, 4000), the search step length as 50, and then the fault period T is obtained by Step 2.
Step 4: The parameters L are updated by EHNR spectral entropy EHNR_E(L). When EHNR_E(L) reaches maximum, the optimal filter length L is obtained. adjacent sampling points of the measured waveform. TKEO can track the change of the measured signal waveform in real time. The algorithm has the advantages of excellent time resolution, simple and fast calculation. The specific mathematical expression and implementation process of the TKEO algorithm was proposed in [44].

The Proposed Method
Inspired by the above-introduced methods, in this section, an IMOMEDA-TKEO rolling bearing fault feature extraction method is proposed. The framework of the proposed method is shown in Figure 2. The general implementation procedure is summarized as follows:

Experimental and Comparative Analysis
In order to verify the effectiveness of the IMOMEDA-TKEO method, the experimental database of rolling bearing from Case Western Reserve University (CWRU) [46], National Aeronautics and Space Administration (NASA) [47] which are authoritative in the field of bearing fault diagnosis and self-made experiment platforms (or KUST-SY for short) are used for experimental analysis. Comparative experiments among IMOMEDA-TKEO, MOMEDA-TKEO and fast kurtogram (FK) [48] are carried out to show the superiority of the proposed method. In each case, we utilize the same input data to complete each comparative experiment. The brief overview of the comparative experiment is shown in Table 1.
In each case, we utilize the same input data to complete each comparative experiment. Table 1. The brief overview of the comparative experiment.

Test 1 IMOMEDA-TKEO
The EHNR spectrum is utilized to calculate fault period T and the improved grid search method base on EHNR spectral entropy is proposed to determine optimal filter length L.
(a) CWRU data sets (b) NASA data sets (c) KUST-SY data sets Test 2 MOMEDA-TKEOclass 1 According to [38], set the default range of T is (0 Where f s is sampling frequency and f cα is fault frequency. The default filter length L [30] is put to 500.
The value range of T is consistent with Test 2 and the method to determine L follows Test 1.

Test 4 Fast kurtogram
Spectral kurtosis is taken as the function of Short Time Fourier Transform (STFT) window width to search the optimal filter parameters.
(a) Compared with the other experiments, the Test 1 demonstrates the effectiveness of the IMOMEDA-TKEO method.
(b) The comparison between Test 2 and Test 3 proves that it is necessary and effective to optimize the filter length of MOMEDA.
(c) Compared with fast kurtogram, the other tests verify that the MOMEDA method can enhance effective fault period impulse component from vibration signals.
(d) The demonstration of pseudo period in Test 3 further proves the importance to calculate fault period T exactly. Figure 3 shows the experiment platform and bearings. Some specific measured parameters are presented in Table 2. The test-bed consists of a two horsepower motor, a torque sensor/decoder, a power tester and an electronic controller. The vibration signals were collected by accelerometers, where the accelerometers are installed at 12 o'clock position at drive end of the motor and fan end of motor housing, respectively. The collected signals were stored in a 16-channel DAT recorder [46] and post-processed in MATLAB 2018a on an IntelI CITM) i5-7300HQ CPU @ 2.50 GHz Laptop. The motor speed is 1797rpm. The sampling frequency f s is 12 kHz, and the number of data points N is 4096. were collected by accelerometers, where the accelerometers are installed at 12 oʹclock position at drive end of the motor and fan end of motor housing, respectively. The collected signals were stored in a 16-channel DAT recorder [46] and post-processed in MATLAB 2018a on an IntelI CITM) i5-7300HQ CPU @ 2.50 GHz Laptop. The motor speed is 1797rpm. The sampling frequency fs is 12 kHz, and the number of data points N is 4096.

Case 1: CWRU Data Analysis
(a) (b) Outer race defect frequency:  Outer race defect frequency: Inner race defect frequency: According to the above parameters and formulas, BPFO = 107.36 Hz and BPFI = 162. 19 Hz are calculated, respectively. Next, we use the datasets in CWRU to complete the comparative experiments. Figure 4 shows the time-domain and frequency-domain analysis results of the original outer race signal, respectively. It is seen from Figure 4 that these waveforms cannot directly obtain the information related to the outer race fault. Further analysis is required to determine the fault characteristic frequency and other details. Therefore, IMOMEDA-TKEO, MOMEDA-TKEO and the FK methods are proposed to complete subsequent analysis.  (a) Using the EHNR spectrum to calculate fault period T. By carrying out the autocorrelation analysis of outer ring, the autocorrelation function (AC) spectrum of the original signal is obtained as shown in Figure 5. It can be seen from Figure 5 that the AC spectrum can reflect the periodic change of fault signal, which provides the basis for the estimation of fault period T. In order to estimate T, we further calculate the EHNR spectrum of signal as shown in Figure 6. It can be seen from Figure 6 that the fault signal shows obvious periodic change and gradually decays with time, which meets the trend of fault signal. Then, the fault period T = 112 is obtained.  (a) Using the EHNR spectrum to calculate fault period T. By carrying out the autocorrelation analysis of outer ring, the autocorrelation function (AC) spectrum of the original signal is obtained as shown in Figure 5. It can be seen from Figure 5 that the AC spectrum can reflect the periodic change of fault signal, which provides the basis for the estimation of fault period T. In order to estimate T, we further calculate the EHNR spectrum of signal as shown in Figure 6. It can be seen from Figure 6 that the fault signal shows obvious periodic change and gradually decays with time, which meets the trend of fault signal. Then, the fault period T = 112 is obtained.

Outer Race Fault Feature Extraction
(A) IMOMEDA-TKEO method; (a) Using the EHNR spectrum to calculate fault period T. By carrying out the autocorrelation analysis of outer ring, the autocorrelation function (AC) spectrum of the original signal is obtained as shown in Figure 5. It can be seen from Figure 5 that the AC spectrum can reflect the periodic change of fault signal, which provides the basis for the estimation of fault period T. In order to estimate T, we further calculate the EHNR spectrum of signal as shown in Figure 6. It can be seen from Figure 6 that the fault signal shows obvious periodic change and gradually decays with time, which meets the trend of fault signal. Then, the fault period T = 112 is obtained.  (b) Using the EHNR spectral entropy to get filter length L Given the fault period T = 112, the filter length L of the MOMEDA is quickly determined by the grid search method. Firstly, the filter length L satisfies L > 2fs /fc, where fs is the sampling frequency and fc is the resonance frequency [45]. To accurately cover the EHNR EHNR Figure 5. Autocorrelation function spectrum.
(A) IMOMEDA-TKEO method; (a) Using the EHNR spectrum to calculate fault period T. By carrying out the autocorrelation analysis of outer ring, the autocorrelation function (AC) spectrum of the original signal is obtained as shown in Figure 5. It can be seen from Figure 5 that the AC spectrum can reflect the periodic change of fault signal, which provides the basis for the estimation of fault period T. In order to estimate T, we further calculate the EHNR spectrum of signal as shown in Figure 6. It can be seen from Figure 6 that the fault signal shows obvious periodic change and gradually decays with time, which meets the trend of fault signal. Then, the fault period T = 112 is obtained.  (b) Using the EHNR spectral entropy to get filter length L Given the fault period T = 112, the filter length L of the MOMEDA is quickly determined by the grid search method. Firstly, the filter length L satisfies L > 2fs /fc, where fs is the sampling frequency and fc is the resonance frequency [45]. To accurately cover the EHNR EHNR Figure 6. EHNR spectrum.
(b) Using the EHNR spectral entropy to get filter length L Given the fault period T = 112, the filter length L of the MOMEDA is quickly determined by the grid search method. Firstly, the filter length L satisfies L > 2f s /f c , where f s is the sampling frequency and f c is the resonance frequency [45]. To accurately cover the whole fault frequency band, the initialization L search range is (250, 4000) and the search step is 50. Secondly, the EHNR spectral entropy is proposed as the target function. Finally, the trend of EHNR spectral entropy is obtained by constantly updating the parameter L as shown in Figure 7. When obtaining the maximum output of EHNR_E(L), the optimal filter length is 2050. whole fault frequency band, the initialization L search range is (250, 4000) and the search step is 50. Secondly, the EHNR spectral entropy is proposed as the target function. Finally, the trend of EHNR spectral entropy is obtained by constantly updating the parameter L as shown in Figure 7. When obtaining the maximum output of EHNR_E(L), the optimal filter length is 2050. (c) Outer race fault feature extraction According to (a) and (b), the fault period T and filter length L of MOMEDA are calculated as 112 and 2050, respectively. Then, the periodic impulse signal x cov (t) of the original signal x(t) is extracted by using the MOMEDA method, as shown in Figure 8. Then, we demodulate x cov (t) with TKEO and get the corresponding spectrum as shown in Figure 9. It can be seen from Figure 9 that the characteristic frequency (105.5 Hz), with its frequency doubling, approaches BPFO (2~9BPFO). Then, it can be inferred that the outer race fault has occurred.
(B) MOMEDA-TKEO method (c) Outer race fault feature extraction According to (a) and (b), the fault period T and filter length L of MOMEDA are calculated as 112 and 2050, respectively. Then, the periodic impulse signal xcov(t) of the original signal x(t) is extracted by using the MOMEDA method, as shown in Figure 8. Then, we demodulate xcov(t) with TKEO and get the corresponding spectrum as shown in Figure 9. It can be seen from Figure 9 that the characteristic frequency (105.5Hz), with its frequency doubling, approaches BPFO (2~9BPFO). Then, it can be inferred that the outer race fault has occurred.
(B) MOMEDA-TKEO method  According to [38] and [30], we analyze the outer race fault signal with the fault period T = (90,134) and filter length L = 500. The multi-point kurtosis diagram is shown in Figure 10. When the multi-point kurtosis reaches the maximum value, the fault period T is 111.5. The optimal deconvolution signal xcov(t) is presented in Figure 11.
The xcov(t) is demodulated by the TKEO and the corresponding spectrum is shown in Figure 12 where the result shows the characteristic frequency (108.4Hz), with its frequency doubling, approaches BPFO (2~9BPFO). Further, it is inferred from the original signal that there exists an outer ring fault.
According to [38] and [30], we analyze the outer race fault signal with the fault period T = (90,134) and filter length L = 500. The multi-point kurtosis diagram is shown in Figure 10. When the multi-point kurtosis reaches the maximum value, the fault period T is 111.5. The optimal deconvolution signal x cov (t) is presented in Figure 11.
The x cov (t) is demodulated by the TKEO and the corresponding spectrum is shown in Figure 12 where the result shows the characteristic frequency (108.4 Hz), with its frequency doubling, approaches BPFO (2~9BPFO). Further, it is inferred from the original signal that there exists an outer ring fault. T = (90,134) and filter length L = 500. The multi-point kurtosis diagram is shown in Figure 10. When the multi-point kurtosis reaches the maximum value, the fault period T is 111.5. The optimal deconvolution signal xcov(t) is presented in Figure 11.
The xcov(t) is demodulated by the TKEO and the corresponding spectrum is shown in Figure 12 where the result shows the characteristic frequency (108.4Hz), with its frequency doubling, approaches BPFO (2~9BPFO). Further, it is inferred from the original signal that there exists an outer ring fault.   According to [38] and [30], we analyze the outer race fault signal with the fault period T = (90,134) and filter length L = 500. The multi-point kurtosis diagram is shown in Figure 10. When the multi-point kurtosis reaches the maximum value, the fault period T is 111.5. The optimal deconvolution signal xcov(t) is presented in Figure 11.
The xcov(t) is demodulated by the TKEO and the corresponding spectrum is shown in Figure 12 where the result shows the characteristic frequency (108.4Hz), with its frequency doubling, approaches BPFO (2~9BPFO). Further, it is inferred from the original signal that there exists an outer ring fault.   (C) FK method Fast kurtogram (FK) [48] is a time-frequency analysis method for spectral kurtosis that shows sufficient efficiency in transient fault detection and has been widely used in fault diagnosis field. Therefore, FK is used as a reference for comparative analysis in this paper. The FK and fault characteristic spectrum are shown in Figure 13 and Figure 14, respectively. The result shows the characteristic frequency (107.2Hz) with its frequency doubling (2~6BPFO) can be extracted. (C) FK method Fast kurtogram (FK) [48] is a time-frequency analysis method for spectral kurtosis that shows sufficient efficiency in transient fault detection and has been widely used in fault diagnosis field. Therefore, FK is used as a reference for comparative analysis in this paper. The FK and fault characteristic spectrum are shown in Figures 13 and 14, respectively. The result shows the characteristic frequency (107.2 Hz) with its frequency doubling (2~6BPFO) can be extracted.
(C) FK method Fast kurtogram (FK) [48] is a time-frequency analysis method for spectral kurtosis that shows sufficient efficiency in transient fault detection and has been widely used in fault diagnosis field. Therefore, FK is used as a reference for comparative analysis in this paper. The FK and fault characteristic spectrum are shown in Figure 13 and Figure 14, respectively. The result shows the characteristic frequency (107.2Hz) with its frequency doubling (2~6BPFO) can be extracted.  (C) FK method Fast kurtogram (FK) [48] is a time-frequency analysis method for spectral kurtosis that shows sufficient efficiency in transient fault detection and has been widely used in fault diagnosis field. Therefore, FK is used as a reference for comparative analysis in this paper. The FK and fault characteristic spectrum are shown in Figure 13 and Figure 14, respectively. The result shows the characteristic frequency (107.2Hz) with its frequency doubling (2~6BPFO) can be extracted.   Compared with the MOMEDA-TKEO method and FK method, the amplitude obtained by the IMOMEDA-TKEO is clearer. Test 3 is carried out through the filter length optimized by the proposed method, without changing the fault period T search interval of traditional MOMEDA. Then, x cov (t) and its Teager-Kaiser energy spectrum are shown in Figures 15 and 16 Compared with the MOMEDA-TKEO method and FK method, the amplitude obtained by the IMOMEDA-TKEO is clearer. Test 3 is carried out through the filter length optimized by the proposed method, without changing the fault period T search interval of traditional MOMEDA. Then, xcov(t) and its Teager-Kaiser energy spectrum are shown in Figure 15 and Figure 16, respectively.    It can be seen from the Mkurt spectrum that the search interval of fault period T directly affects the denoising effect of MOMEDA. However, the interval of T depends on prior knowledge and has trial and error. When the original signal is affected by complex background noise, the Mkurt spectrum may not provide correct fault period information according to the location of maximum multipoint kurtosis. For example, we get the pseudo fault period as T = 97.2. Furthermore, the optimal periodic pulse signal xcov(t) cannot be obtained, as shown in Figure 17. The deconvolution effect is not ideal. The wrong fault characteristic frequency is 123Hz and its spectrum is shown in Figure 18. It can be seen from the Mkurt spectrum that the search interval of fault period T directly affects the denoising effect of MOMEDA. However, the interval of T depends on prior knowledge and has trial and error. When the original signal is affected by complex background noise, the Mkurt spectrum may not provide correct fault period information according to the location of maximum multipoint kurtosis. For example, we get the pseudo fault period as T = 97.2. Furthermore, the optimal periodic pulse signal x cov (t) cannot be obtained, as shown in Figure 17. The deconvolution effect is not ideal. The wrong fault characteristic frequency is 123Hz and its spectrum is shown in Figure 18.  In order to further confirm the effectiveness of the proposed method quantitatively, a new signal-to-noise ratio (SNR) [49] is used as an index for Teager-Kaiser energy spectrum to complete the quantitative comparative analysis among the IMOMEDA-TKEO, MOMEDA-TKEO and the FK. The new SNR is defined as follows.  In order to further confirm the effectiveness of the proposed method quantitatively, a new signal-to-noise ratio (SNR) [49] is used as an index for Teager-Kaiser energy spectrum to complete the quantitative comparative analysis among the IMOMEDA-TKEO, MOMEDA-TKEO and the FK. The new SNR is defined as follows. In order to further confirm the effectiveness of the proposed method quantitatively, a new signal-to-noise ratio (SNR) [49] is used as an index for Teager-Kaiser energy spectrum to complete the quantitative comparative analysis among the IMOMEDA-TKEO, MOMEDA-TKEO and the FK. The new SNR is defined as follows.
where f d is the fault characteristic frequency, M is the multiple of the fault characteristic frequency; p(f ) is the Teager-Kaiser energy spectrum [50]. It can be seen from Table 3 that the SNR fd index of the presented method is the highest, which indicates that the proposed method has better ability in weak feature representation.   Figure 19 shows the time-domain and frequency-domain analysis results of the original inner race signal, respectively. It is seen from Figure 19 that these waveforms cannot directly obtain the information related to the inner race fault. Further analysis is required to determine the fault characteristic frequency and other details. Therefore, IMOMEDA-TKEO and MOMEDA-TKEO are proposed to complete subsequent analysis.  The proposed method can achieve comparable effects as the MOMEDA-TKEO method as shown in Figure 20a,b. They both can extract the inner race fault features from the original vibration signals. It further proves the effectiveness of two methods.

Inner Race Fault Feature Extraction
Under the condition that the fault period T is invariant, filter length L is obtained by the proposed method and presented in Figure 20c. It can be seen from Figure 20c that the amplitude of TKEO spectrum apparently increases, which indicates filter length L directly affects the deconvolution performance of the MOMEDA method. In addition, it can be seen from Figure 20d that the selection of fault period T directly affects the performance of the method.
FK analysis results are shown in Figure 21. Compared with Figure 20, similar conclusions can be obtained, but the fault characteristic frequency with its frequency doubling characteristic obtained by IMOMEDA-TEO method are more prominent and clearer than those obtained by FK method.
From Table 4, the SNRfd for the IMOMEDA method is higher, which further implies the IMOMEDA has the qualitative and quantitative advantage over the original MOMEDA algorithm and FK method. The proposed method can achieve comparable effects as the MOMEDA-TKEO method as shown in Figure 20a,b. They both can extract the inner race fault features from the original vibration signals. It further proves the effectiveness of two methods.
Under the condition that the fault period T is invariant, filter length L is obtained by the proposed method and presented in Figure 20c. It can be seen from Figure 20c that the amplitude of TKEO spectrum apparently increases, which indicates filter length L directly affects the deconvolution performance of the MOMEDA method. In addition, it can be seen from Figure 20d that the selection of fault period T directly affects the performance of the method.
FK analysis results are shown in Figure 21. Compared with Figure 20, similar conclusions can be obtained, but the fault characteristic frequency with its frequency doubling characteristic obtained by IMOMEDA-TEO method are more prominent and clearer than those obtained by FK method.
From Table 4, the SNR fd for the IMOMEDA method is higher, which further implies the IMOMEDA has the qualitative and quantitative advantage over the original MOMEDA algorithm and FK method.

Case 2: NASA Data Analysis
In order to verify the good generalization ability of the IMOMEDA-TKEO method, the further experiments were carried out with NASA bearing datasets.
The test-bed shown in Figure 22 is composed of four bearings, rotating shaft, two acceleration sensors and electronic control unit. The test bearings are located in the motor shaft. Some accelerometers sensors are placed to the motor housing. The motor speed is 2000 rpm and the digital data is collected by sampling with the frequency of 20 kHz. The detailed measured parameters of the bearing are shown in Table 5. Experimental comparative analysis is complete by utilizing the vibration signals of outer race and inner race. Furthermore, the vibration signal time and frequency analysis is presented in Figure 23. It is seen from Figure 23 that these waveforms cannot directly determine the fault characteristic frequency and other details. In order to further analyze the fault information of the bearing, the same implementation process as Section 4.1 is taken to achieve the deconvolution results of the outer race and inner race, respectively, as shown in Figures 24 and 25.  From the bearing parameters shown in Table 5, according to Equation (19) and Equation (20), the fault characteristic frequency BPFO = 236.4 Hz for outer race and BPFI = 296.93 Hz for inner race can be calculated, respectively, and the data length N is 4096 points.
Experimental comparative analysis is complete by utilizing the vibration signals of outer race and inner race. Furthermore, the vibration signal time and frequency analysis is presented in Figure 23. It is seen from Figure 23 that these waveforms cannot directly determine the fault characteristic frequency and other details. In order to further analyze the fault information of the bearing, the same implementation process as Section 4.1 is taken to achieve the deconvolution results of the outer race and inner race, respectively, as shown in Figures 24 and 25. (a) (b) Figure 22. Bearing test rig of NASA data center (a) Test platform (b) Sensor layout.
Experimental comparative analysis is complete by utilizing the vibration signals of outer race and inner race. Furthermore, the vibration signal time and frequency analysis is presented in Figure 23. It is seen from Figure 23 that these waveforms cannot directly determine the fault characteristic frequency and other details. In order to further analyze the fault information of the bearing, the same implementation process as Section 4.1 is taken to achieve the deconvolution results of the outer race and inner race, respectively, as shown in Figures 24 and 25.    The experimental results of Figures 24 and 25 show that the proposed IMOMEDA-TKEO method effectively extracts fault features of outer race and inner race and achieves significantly better processing performance than the MOMEDA-TKEO method.
FK based method can extract fault features of outer race and inner race from the original vibration signals. Unfortunately, it can be seen from Figure 26 that the obtained fault characteristic frequency is not obvious (the amplitude is far less than the other two methods) and contains more noise components. Compared with Figures 24 and 25, the filtering performance is poorer than the IMOMEDA-TKEO method. The experimental results of Figures 24 and 25 show that the proposed IMOMEDA-TKEO method effectively extracts fault features of outer race and inner race and achieves significantly better processing performance than the MOMEDA-TKEO method.
FK based method can extract fault features of outer race and inner race from the original vibration signals. Unfortunately, it can be seen from Figure 26 that the obtained fault characteristic frequency is not obvious (the amplitude is far less than the other two methods) and contains more noise components. Compared with Figures 24 and 25, the filtering performance is poorer than the IMOMEDA-TKEO method. In order to intuitively compare the effectiveness of the method, the SNRfd index is calculated and shown in Table 6. By comparing and analyzing the data in the Table 6, the conclusions similar to Section 4.1 can be drawn.   In order to intuitively compare the effectiveness of the method, the SNR fd index is calculated and shown in Table 6. By comparing and analyzing the data in the Table 6, the conclusions similar to Section 4.1 can be drawn.

Application Testing
In order to further prove the validity of the proposed method, the practical test data from the self-made KUST-SY experiment platform is selected for subsequent analysis. The platform is able to simulate different failure conditions of rolling bearings. The structure of the test rig is shown in Figure 27, which consists of driving motor, rotating shaft, support bearing, hydraulic loading system, accelerometers and tested bearing. The device can realize radial hydraulic loading and the rotation speed is controlled by computer. Accelerometers are installed in the vertical and horizontal directions of the test bearing base to collect the bearing vibration signal. Because the load is applied in the horizontal direction, the accelerometer placed in this direction is capable of capturing more fault information of the tested bearings. Therefore, the horizontal vibration signals are selected to analyze bearing fault characteristics.  The failure behaviors of bearing system are diverse, including stripping, wear, fracture, crack, electrical erosion, indentation for the inner race and outer race. As displayed in Figure 28, in this experiment, inner ring crack and out ring crack of the bearing have been mainly simulated considering the existing experimental conditions. The raceway surface cracks (the fault diameter is 0.2 mm) of the inner ring and outer ring have an artificial linear cutting with fiber laser cutter. The sampling frequency is 51.2 kHz and the motor speed is 1797 rpm. The sampling interval is 20 s and the sampling time is 10 s. The tested bearing model is 6205-2RSH, and its relevant parameters are shown in Table 7. According to Table 7, the fault characteristic frequency BPFO = 107.22Hz for the outer race and BPFI = 162.33Hz for the inner race can be calculated, respectively. The failure behaviors of bearing system are diverse, including stripping, wear, fracture, crack, electrical erosion, indentation for the inner race and outer race. As displayed in Figure 28, in this experiment, inner ring crack and out ring crack of the bearing have been mainly simulated considering the existing experimental conditions. The raceway surface cracks (the fault diameter is 0.2 mm) of the inner ring and outer ring have an artificial linear cutting with fiber laser cutter. The sampling frequency is 51.2 kHz and the motor speed is 1797 rpm. The sampling interval is 20 s and the sampling time is 10 s. The tested bearing model is 6205-2RSH, and its relevant parameters are shown in Table 7. According to Table 7, the fault characteristic frequency BPFO = 107.22 Hz for the outer race and BPFI = 162.33 Hz for the inner race can be calculated, respectively. been mainly simulated considering the existing experimental conditions. The raceway surface cracks (the fault diameter is 0.2 mm) of the inner ring and outer ring have an artificial linear cutting with fiber laser cutter. The sampling frequency is 51.2 kHz and the motor speed is 1797 rpm. The sampling interval is 20 s and the sampling time is 10 s. The tested bearing model is 6205-2RSH, and its relevant parameters are shown in Table 7. According to Table 7, the fault characteristic frequency BPFO = 107.22Hz for the outer race and BPFI = 162.33Hz for the inner race can be calculated, respectively.  The time and frequency analysis is shown in Figure 29 for the vibration signal of the inner ring and outer ring (The length of data-set N is 32768). We cannot get valuable information related to fault characteristic frequency from Figure 29, which is the same as previous cases. Therefore, the fault information is determined by the same implementation process with Section 4.1 and plotted in Figure 30, Figure 31 and Figure 32, respectively.  The time and frequency analysis is shown in Figure 29 for the vibration signal of the inner ring and outer ring (The length of data-set N is 32768). We cannot get valuable information related to fault characteristic frequency from Figure 29, which is the same as previous cases. Therefore, the fault information is determined by the same implementation process with Section 4.1 and plotted in Figures 30-32, respectively.   It can be seen from Figures 30 and 31 that the fault features extracted by the proposed IMOMEDA-TKEO method are more obvious in the Teager energy spectrum. At the same time, it can be seen from Figures 30d and 31d that the improper period T will lead to erroneous results in the MOMEDA method.  The fault characteristic frequency for the inner race and outer race extracted by FK is shown in Figure 32. Comparing with the fault characteristic spectrum extracted by MOMEDA based method as shown in Figures 30 and 31, it can be seen that the fault characteristic frequencies can be accurately obtained by FK. However, there are many interference components in the spectrum, which result in the faint characteristic amplitude. In particular, the outer race seems to be more seriously disturbed by noise, and only three times of double frequency is extracted, while the other two methods can extract nine times of double frequency. In addition, through the comparison of performance indicator SNRfd presented in Table 8, it is also found that the IMOMEDA-TKEO method can achieve better performance than the MOMEDA-TKEO method and FK method in noise reduction. The fault characteristic frequency for the inner race and outer race extracted by FK is shown in Figure 32. Comparing with the fault characteristic spectrum extracted by MOMEDA based method as shown in Figures 30 and 31, it can be seen that the fault characteristic frequencies can be accurately obtained by FK. However, there are many interference components in the spectrum, which result in the faint characteristic amplitude. In particular, the outer race seems to be more seriously disturbed by noise, and only three times of double frequency is extracted, while the other two methods can extract nine times of double frequency. In addition, through the comparison of performance indicator SNR fd presented in Table 8, it is also found that the IMOMEDA-TKEO method can achieve better performance than the MOMEDA-TKEO method and FK method in noise reduction.

Comparative Analysis with Three Cases
Through the comparative experimental analysis of IMOMEDA-TKEO, MOMEDA-TKEO and FK, the following can be summarized: (a) Combining with three case studies, the experimental results show that the IMOMEDA-TKEO method is better than the other method described in Table 1, which enhances the application prospect of the proposed method.
(b) In this experiment, we select execution time and iterations to evaluate the complexity of the optimal parameter selection method and all the experiments are programmed in MATLAB 2018a on an Intel(R) Core(TM) i5-7300HQ CPU @ 2.50 GHz Laptop with 16.00 GB RAM. The parameter optimization time and iterations of IMOMEDA algorithm and the main program execution time of the comparative methods are shown in Tables 9 and 10, respectively. It is seen from Table 9 that we can calculate fault period T of MOMEDA quickly with little time. Compared with the traditional MOMEDA method in which the parameters are set to the default value, the grid search method based on EHNR spectral entropy for selecting optimal filter length L does consume some time. However, optimization algorithm itself is not the main reason of time-consuming. The main reason is that as the iterations increases, the MOMEDA filter time keeps accumulating, leading to more training time of filter length. The purpose of this paper is to accurately calculate the fault period T of MOMEDA and select the optimal filtering length L of MOMEDA. The rationality and necessity of MOMEDA parameter selection are demonstrated on three different datasets. Although the calculation needs a certain time, it is within the acceptable range.
(c) It can be seen from Table 10 that the FK method and the traditional MOMEDA method (default filter length and fault period) take less time. However, when dealing with abnormal or unknown vibration signals, the default parameters of MOMEDA are obviously not feasible. It is verified that the time-consuming by the trial and error to determine the MOMEDA parameters is almost the same as that of the proposed method. However, the result of unreasonable parameter setting is much worse than IMOMEDA, such as Figures 24d and 25d, which also proves the significance of optimizing MOMEDA parameters. The proposed IMOMEDA-TKEO method can effectively estimate the period T of vibration signal, and set the filter length L reasonably by using the improved grid search method. The execution time is within the acceptable range.

Discussion
(a) For the traditional MOMEDA method, the range of T has to be set by a trial and error method. Furthermore, the range of T selected in the above experiments is the best result obtained by multiple attempts. If the selected range changes, the pseudo period will cause misdiagnosis. The proposed method introduces the EHNR spectrum to calculate the signal period T quickly and optimizes the filter length L by the improved grid search method. It transforms the iterative optimization problem into a single dimension optimization problem. Hence, the proposed method increases the operation efficiency and improves the deconvolution performance of MOMEDA method.
(b) In the proposed method, the search range of the filter length L is set as [34]. Only the lower limit of the filter length L is given, while there is no upper limit, which makes the search range larger and the algorithm run for a long time. According to the specific characteristics of the signal, optimizing the search range can improve the efficiency of the algorithm. In addition, if we can obtain the corresponding relationship and change rule between the MOMEDA filter length L and period T on the basis of (a), it will further improve the operation efficiency and the practical application value of the proposed method.
(c) Compared with other types of signal characteristics, the practical test vibration signals for inner race from the rolling bearing fault simulation experiment platform are not only disturbed by noise, but also by harmonic signal of the outer race and rolling element. Under strong background noise, the IMOMEDA-TKEO method has been demonstrated to significantly outperform traditional signal processing methods in terms of noise reduction performance and fault feature representation robustness, which boosts the generalization ability.

Conclusions
In this paper, a rolling bearing fault feature extraction method based on the IMOMEDA-TKEO is proposed. Some conclusions are obtained as follows: (1) The EHNR method is introduced to obtain the EHNR spectrum of the original signal and solves the problem that the fault period T of MOMEDA depends on prior knowledge.
(2) The EHNR spectral entropy is proposed as the objective function, and the improved grid search method is introduced to solve the filter length L selection in MOMEDA.
(3) Comparative experimental analysis of the proposed method with the MOMEDA-TKEO method and the FK method are completed. The proposed method can achieve the best filtering performance, which verifies the effectiveness and feasibility of the IMOMEDA-TKEO.
In the proposed method, the upper limit of the filter length L is not given, which makes the execution time and iterations of the algorithm take a certain time (it is within the acceptable range). In future studies, we will focus on improving the calculation accuracy of fault period T for complex signals and optimizing the search range of filter length L. We will also explore the change rule between the filter length L and fault period T to improve the performance of MOMEDA. Institutional Review Board Statement: Not applicable.

Informed Consent Statement: Not applicable.
Data Availability Statement: Data is available on Reference [46] and Reference [47].