Bearing Fault Diagnosis Based on Energy Spectrum Statistics and Modified Mayfly Optimization Algorithm

This study proposes a novel resonance demodulation frequency band selection method named the initial center frequency-guided filter (ICFGF) to diagnose the bearing fault. The proposed technology has a better performance on resisting the interference from the random impulses. More explicitly, the ICFGF can be summarized as two steps. In the first step, a variance statistic index is applied to evaluate the energy spectrum distribution, which can adaptively determine the center frequency of the fault impulse and suppress the interference from random impulse effectively. In the second step, a modified mayfly optimization algorithm (MMA) is applied to search the optimal resonance demodulation frequency band based on the center frequency from the first step, which has faster convergence. Finally, the filtered signal is processed by the squared envelope spectrum technology. Results of the proposed method for signals from an outer fault bearing and a ball fault bearing indicate that the ICFGF works well to extract bearing fault feature. Furthermore, compared with some other methods, including fast kurtogram, ensemble empirical mode decomposition, and conditional variance-based selector technology, the ICFGF can extract the fault characteristic more accurately.


Introduction
Rolling bearings have been widely used in modern industry, especially in rotation machinery. The operation pattern of the bearing determines it is a worn part. More explicitly, bearings in the engine or high-speed railway wear easily because of the poor working environment and unstable load. Consequently, bearing fault is a hot research topic in rotation machinery fault diagnosis [1][2][3]. The vibration signal can reflect the health status of rotation machinery effectively [4] and has been widely applied in fault diagnosis [5][6][7]. However, the response of the early localized damage is too weak to be detected directly and is usually submerged by noises from the environment and other running components. Consequently, a trustworthy signal process technology must be carried out.
For decades, many methods have been proposed for bearing fault diagnosis. Envelope analysis is one of the most popular approaches and can be designed as two steps [8][9][10]. The first step is to determine the optimal resonance frequency band (ORFB). The next step is to process the signal filtered by the demodulation technology to obtain the squared envelope (SE) or the squared envelope spectrum (SES), which can display the fault feature much more clearly than the raw feature. The key point in envelope analysis is to find the ORFB with the help of a suitable index [11,12].
For bearing fault diagnosis, the most popular approach is fast kurtogram (FK), whose index is kurtosis [13]. However, FK has two drawbacks. One is the frame inflexibility of the filter-banks, and the other is the vulnerable index. Many researchers have focused on solving these two problems. For the former, wavelet packet transform (WPT) [14] and a

The Proposed Method
This section briefly introduces some relative methods, including the conditional variance statistic index and the MMA optimization algorithm.

CV for Gaussian Data
The distribution for a set of Gaussian data can be described as where µ and σ 2 denote the mean and variance, respectively. From Hebda-Sobkowicz et al. [25] and Jaworski and Pitera [26], the conditional variance can be defined as Based on Jaworski and Pitera [26], Gaussian distribution can be divided into three parts, and each of them has the same variance. The ratio for each part is 20/60/20, which is displayed in Figure 1. where μ and 2 σ denote the mean and variance, respectively. From Hebda-Sobkowicz et al. [25] and Jaworski and Pitera [26], the conditional variance can be defined as Based on Jaworski and Pitera [26], Gaussian distribution can be divided into three parts, and each of them has the same variance. The ratio for each part is 20/60/20, which is displayed in Figure 1. Furthermore, Hebda-Sobkowicz et al. [25] and Jelito and Pitera [27] proposed to divide the distribution into seven parts, which is shown in Figure 2. The ratio for each is designed as 0.4/5.8/24.6/38.4/24.6/5.8/0.4. Similar to the 20/60/20 principle, the variance for each part is also equivalent. The power spectral density of Gaussian noise is a constant, which means the energy of Gaussian noise for the whole frequency domain should also be a constant. In practice, Gaussian noise will be filtered by the collection equipment, which means its energy spectrum will be presented as a Gaussian distribution. Based on the inheritance of data, a sub-part from the energy spectrum is also a Gaussian distribution. However, the energy for the impulse component will be concentrated on a series of resonance spectrum lines, which means its energy spectrum is not a Gaussian distribution. Consequently, we can distinguish it from the energy spectrum.
To identify the impulse component, a novel index is proposed in this study. Firstly, the signal is processed by FFT technology. Next, the square operation is applied to transform the FFT spectrum to the energy spectrum. Then, to suppress the effect from a harmonic component, areas A1 and A7 shown in Figure 2 are abandoned. Finally, by using the result of 20/60/20 for the remaining data, the index can be calculated as ( ) Furthermore, Hebda-Sobkowicz et al. [25] and Jelito and Pitera [27] proposed to divide the distribution into seven parts, which is shown in Figure 2. The ratio for each is designed as 0.4/5.8/24.6/38.4/24.6/5.8/0.4. Similar to the 20/60/20 principle, the variance for each part is also equivalent. where μ and 2 σ denote the mean and variance, respectively. From Hebda-Sobkowicz et al. [25] and Jaworski and Pitera [26], the conditional variance can be defined as Based on Jaworski and Pitera [26], Gaussian distribution can be divided into three parts, and each of them has the same variance. The ratio for each part is 20/60/20, which is displayed in Figure 1. Furthermore, Hebda-Sobkowicz et al. [25] and Jelito and Pitera [27] proposed to divide the distribution into seven parts, which is shown in Figure 2. The ratio for each is designed as 0.4/5.8/24.6/38.4/24.6/5.8/0.4. Similar to the 20/60/20 principle, the variance for each part is also equivalent. Sample for seven parts: A1, A2, A3, A4, A5, A6, and A7 corresponding to 0.4%, 5.8%, 24.6%, 38.4%, 24.6%, 5.8%, and 0.4%, respectively.
The power spectral density of Gaussian noise is a constant, which means the energy of Gaussian noise for the whole frequency domain should also be a constant. In practice, Gaussian noise will be filtered by the collection equipment, which means its energy spectrum will be presented as a Gaussian distribution. Based on the inheritance of data, a sub-part from the energy spectrum is also a Gaussian distribution. However, the energy for the impulse component will be concentrated on a series of resonance spectrum lines, which means its energy spectrum is not a Gaussian distribution. Consequently, we can distinguish it from the energy spectrum.
To identify the impulse component, a novel index is proposed in this study. Firstly, the signal is processed by FFT technology. Next, the square operation is applied to transform the FFT spectrum to the energy spectrum. Then, to suppress the effect from a harmonic component, areas A1 and A7 shown in Figure 2 are abandoned. Finally, by using the result of 20/60/20 for the remaining data, the index can be calculated as ( ) The power spectral density of Gaussian noise is a constant, which means the energy of Gaussian noise for the whole frequency domain should also be a constant. In practice, Gaussian noise will be filtered by the collection equipment, which means its energy spectrum will be presented as a Gaussian distribution. Based on the inheritance of data, a sub-part from the energy spectrum is also a Gaussian distribution. However, the energy for the impulse component will be concentrated on a series of resonance spectrum lines, which means its energy spectrum is not a Gaussian distribution. Consequently, we can distinguish it from the energy spectrum.
To identify the impulse component, a novel index is proposed in this study. Firstly, the signal is processed by FFT technology. Next, the square operation is applied to transform the FFT spectrum to the energy spectrum. Then, to suppress the effect from a harmonic component, areas A 1 and A 7 shown in Figure 2 are abandoned. Finally, by using the result of 20/60/20 for the remaining data, the index can be calculated as where A R and A M denote the areas shown in Figure 1, and N is the length of the data analyzed. From Equation (3), the value of CV should be small if the data to be analyzed are Gaussian. CV should be large for the impulse component. In this study, CV is combined with sliding technology. Consequently, the length of the window and its step size must be determined first. In this study, we suggest that the length of the window should meet the following conditions: (1) The length can divide the whole energy spectrum into 50-100 sub-parts, and the center frequency for each sub-part is mark as f c . (2) The length of the window should be wide enough to contain at least three fault feature spectrum lines. Considering the unknown fault type, the length can be set as two times larger than the inner fault feature frequency. (3) The overlap ratio between two windows should be more than 50%. In this study, the overlap ratio is set as 80%.
To verify the effectiveness of CV and the sliding window strategy, the simulation signal is designed as where the first item on the right of Equation (4) denotes the impulse component with a 20 Hz period. The second denotes the harmonic component, and the last item denotes Gaussian noise with 0.4 intensity. Three impulse interferences are added into the simulation signal and marked by red circles. The corresponding time waveform and FFT spectrum are shown in Figure 3. Based on Figure 3b, the impulse component is interfered by the noise seriously. The proposed energy spectrum statistic index is applied to process these data, and the result is shown in Figure 3c. From it, the frequency of the max CV is 1993 Hz, which is close to the design frequency. An important thing to keep in mind is that CV is normalized in the whole study.
where A R and A M denote the areas shown in Figure 1, and N is the length of the data analyzed. From Equation (3), the value of CV should be small if the data to be analyzed are Gaussian. CV should be large for the impulse component. In this study, CV is combined with sliding technology. Consequently, the length of the window and its step size must be determined first. In this study, we suggest that the length of the window should meet the following conditions: (1) The length can divide the whole energy spectrum into 50-100 sub-parts, and the center frequency for each sub-part is mark as c f .
(2) The length of the window should be wide enough to contain at least three fault feature spectrum lines. Considering the unknown fault type, the length can be set as two times larger than the inner fault feature frequency. (3) The overlap ratio between two windows should be more than 50%. In this study, the overlap ratio is set as 80%.
To verify the effectiveness of C V and the sliding window strategy, the simulation signal is designed as where the first item on the right of Equation (4) Figure 3. Based on Figure 3b, the impulse component is interfered by the noise seriously. The proposed energy spectrum statistic index is applied to process these data, and the result is shown in Figure 3c. From it, the frequency of the max CV is 1993 Hz, which is close to the design frequency. An important thing to keep in mind is that C V is normalized in the whole study. Amp. (a) Amp.
To verify the effectiveness of the index further, Equation (4) is repeated 1000 times to obtain 1000 simulation signals. The 1000 signals are analyzed by the CV index, and the results are shown in Figure 4. As shown in Figure 4a, most of the results are close to the designed resonant frequency, and the density curve of these results displayed in Figure   To verify the effectiveness of the index further, Equation (4) is repeated 1000 times to obtain 1000 simulation signals. The 1000 signals are analyzed by the CV index, and the results are shown in Figure 4. As shown in Figure 4a, most of the results are close to the designed resonant frequency, and the density curve of these results displayed in Figure   vs. c CV f .
To verify the effectiveness of the index further, Equation (4) is repeated 1000 times to obtain 1000 simulation signals. The 1000 signals are analyzed by the CV index, and the results are shown in Figure 4. As shown in Figure 4a, most of the results are close to the designed resonant frequency, and the density curve of these results displayed in Figure

The MMA Optimization Algorithm
The center frequency of the ORFB can be obtained by the index shown in the p vious section. In this current section, we determine the bandwidth for the ORFB throu an MMA optimization algorithm.
The mayfly algorithm (MA) was first proposed by Zervoudakis and Tsafarakis [ Its main steps can be summarized as follows.
(1) Initialization: The populations of males and females are first initiated

The MMA Optimization Algorithm
The center frequency of the ORFB can be obtained by the index shown in the previous section. In this current section, we determine the bandwidth for the ORFB through an MMA optimization algorithm.
The mayfly algorithm (MA) was first proposed by Zervoudakis and Tsafarakis [28]. Its main steps can be summarized as follows.
(1) Initialization: The populations of males and females are first initiated as x = [x 1 , · · · , x d ] and y = [y 1 , · · · , y d ], respectively. The corresponding velocity is (2) Males' movement: The global best position for the current iteration is determined by the object function f (x), which is marked as f (gbest). The position for the next iteration is updated based on f (gbest) and the Cartesian distance between the personal element and the global best agent gbest, which can be described as where x t ij corresponds to the agent i in dimension j at the current iteration t and v t ij denotes its velocity. a 1 and a 2 denote the global and personal learning coefficient, respectively. r g and r p denote the Cartesian distance for global and personal, respectively. The velocity of the best agent in the current iteration is updated as v t+1 = v t + d × r, where d denotes the nuptial dance and r is a random variable located in [−1, 1].
(3) Females' movement: Compared with the movement of males, females will fly to males for breeding. Thus, their velocities are updated based on the Cartesian distance between themselves and the males, which can be described as where y corresponds to the female agents and a 3 denotes the learning coefficient. β is the distance sight coefficient, which is a constant. r m f denotes the Cartesian distance between the male and female agents. The best female is arranged to match with the best male, and the second-best female is matched to the second-best male. Consequently, the positions of males are very important in MA. (4) Mating: Each couple produces two offspring in MA. One of them is added to the male population randomly, and the other is added to the female population. (5) Updating: The worst solution is replaced by the best solution, and the processes above are repeated until the stop criteria are met.
Based on the description above, understanding that the position of the male population is the key point in MA is easy. Their positions are initialized randomly in the original MA. However, in this study, we can design a specific initiation process based on the result of the energy spectrum statistic index, which can adjust the distribution of the males' positions corresponding to center frequencies reasonably. In this study, the number of agents is 10 for both the male population and female population. The range for the center frequency

The Proposed Method
Before introducing the proposed method, the optimization function must be determined. In this study, weight kurtosis is set as the objection function and can be calculated as where C denotes the correlation coefficient between the filtered signal and the raw signal, and K denotes the kurtosis of the filtered signal. Given that the goal of MA is to determine the global minimum, the objective function in MMA should be set as −CK.
The main steps of the proposed method can be summarized as follows: Step 1: The measured signal is first processed by FFT technology, and square operation is also applied to transform the FFT spectrum into an energy spectrum.
Step 2: The energy spectrum statistic index CV is applied to evaluate the center frequency for the ORFB.
Step 3: The width of the ORFB is determined by mining −CK with the MMA optimization algorithm.
Step 4: SES technology is applied to analyze the signal filtered with the optimal ORFB. To understand the proposed method comfortably, the flowchart of ICFGF is shown in Figure 5. ed as = × CK C K , (9) where C denotes the correlation coefficient between the filtered signal and the raw signal, and K denotes the kurtosis of the filtered signal. Given that the goal of MA is to determine the global minimum, the objective function in MMA should be set as CK − . The main steps of the proposed method can be summarized as follows: Step 1: The measured signal is first processed by FFT technology, and square operation is also applied to transform the FFT spectrum into an energy spectrum.
Step 2: The energy spectrum statistic index CV is applied to evaluate the center frequency for the ORFB.
Step 3: The width of the ORFB is determined by mining CK − with the MMA optimization algorithm.
Step 4: SES technology is applied to analyze the signal filtered with the optimal ORFB.
To understand the proposed method comfortably, the flowchart of ICFGF is shown in Figure 5.

Case Study
This section applies signals from an outer fault bearing and a ball fault bearing to verify the effectiveness of the proposed method. The highlight of the proposed method is shown by comparing with some existing methods, including FK, EEMD, and the CVB index technology.

Inner Fault Diagnosis
The data come from Curtin University and have been applied in Qin et al. [24]. The signal is produced by a machinery fault simulator, which is displayed in Figure 6.
The signal is collected by an accelerometer with the sample frequency of 51.2 kHz. The test bearing is MB ER-16K, and there is a local defect which exists in its outer race. The test shaft speed is 1740 rpm, and the corresponding ball pass frequency outer race (BPFO) is 103.6 Hz. The length of the signal applied in this part is 1 s.

.Inner Fault Diagnosis
The data come from Curtin University and have been applied in Qin et al. [24]. The signal is produced by a machinery fault simulator, which is displayed in Figure 6. The signal is collected by an accelerometer with the sample frequency of 51.2 kHz. The test bearing is MB ER-16K, and there is a local defect which exists in its outer race. The test shaft speed is 1740 rpm, and the corresponding ball pass frequency outer race (BPFO) is 103.6 Hz. The length of the signal applied in this part is 1 s.
The time waveform of the raw signal and its fast Fourier transform (FFT) spectrum are shown in Figure 7. As shown in Figure 7a, some random impulses (marked by red circles) are found. The FFT spectrum shown in Figure 7b tells us the fault feature is submerged by noise. Next, the proposed method is applied for the analysis.   The time waveform of the raw signal and its fast Fourier transform (FFT) spectrum are shown in Figure 7. As shown in Figure 7a, some random impulses (marked by red circles) are found. The FFT spectrum shown in Figure 7b tells us the fault feature is submerged by noise. Next, the proposed method is applied for the analysis.
signal is produced by a machinery fault simulator, which is displayed in Figure 6.    The center frequency is first evaluated by the CV index, and the result is shown in Figure 8. The center frequency is 2750 Hz. Then, the MMA is applied to find the ORFB. In this case, the range for searching center frequency is set as [2337, 3162] Hz, and the range for searching bandwidth is set as [124, 628] Hz. From the MMA, the center frequency of ORFB is 2806 Hz and the corresponding bandwidth is 628 Hz.
The center frequency is first evaluated by the CV index, and the result is shown in Figure 8. The center frequency is 2750 Hz. Then, the MMA is applied to find the ORFB. In this case, the range for searching center frequency is set as [2337, 3162] Hz, and the range for searching bandwidth is set as [124, 628] Hz. From the MMA, the center frequency of ORFB is 2806 Hz and the corresponding bandwidth is 628 Hz. The results are shown in Figure 9. From Figure 9a, the periodic impulses can be observed easily. Its SES depicted in Figure 9b shows the fault features (BPFO, 2BPFO, and 3BPFO) clearly. Both MMA and MA are re-run three times to show the superiority of MMA. Their average results are shown in Figure 10. From it, we can say the MMA has a faster convergence. The results are shown in Figure 9. From Figure 9a, the periodic impulses can be observed easily. Its SES depicted in Figure 9b shows the fault features (BPFO, 2BPFO, and 3BPFO) clearly.
The center frequency is first evaluated by the CV index, and the result is shown in Figure 8. The center frequency is 2750 Hz. Then, the MMA is applied to find the ORFB. In this case, the range for searching center frequency is set as [2337, 3162] Hz, and the range for searching bandwidth is set as [124, 628] Hz. From the MMA, the center frequency of ORFB is 2806 Hz and the corresponding bandwidth is 628 Hz. The results are shown in Figure 9. From Figure 9a, the periodic impulses can be observed easily. Its SES depicted in Figure 9b shows the fault features (BPFO, 2BPFO, and 3BPFO) clearly. Both MMA and MA are re-run three times to show the superiority of MMA. Their average results are shown in Figure 10. From it, we can say the MMA has a faster convergence.  Both MMA and MA are re-run three times to show the superiority of MMA. Their average results are shown in Figure 10. From it, we can say the MMA has a faster convergence.
The center frequency is first evaluated by the CV index, and the result is shown in Figure 8. The center frequency is 2750 Hz. Then, the MMA is applied to find the ORFB. In this case, the range for searching center frequency is set as [2337, 3162] Hz, and the range for searching bandwidth is set as [124,628] Hz. From the MMA, the center frequency of ORFB is 2806 Hz and the corresponding bandwidth is 628 Hz. The results are shown in Figure 9. From Figure 9a, the periodic impulses can be observed easily. Its SES depicted in Figure 9b shows the fault features (BPFO, 2BPFO, and 3BPFO) clearly. Both MMA and MA are re-run three times to show the superiority of MMA. Their average results are shown in Figure 10. From it, we can say the MMA has a faster convergence.  To highlight the superiority of the ICFGF, the outer fault signal is also processed by FK, EEMD, and CVB technology. The results by FK are displayed in Figure 11. From Figure 11a, the ORFB is [8533, 9599] Hz. The filtered signal and the corresponding SES are shown in Figure 11b, which illustrates that the SES is complex and the fault feature is interfered by the noise.
To highlight the superiority of the ICFGF, the outer fault signal is also processed by FK, EEMD, and CVB technology. The results by FK are displayed in Figure 11. From Figure 11a, the ORFB is [8533, 9599] Hz. The filtered signal and the corresponding SES are shown in Figure 11b, which illustrates that the SES is complex and the fault feature is interfered by the noise.  Finally, the outer fault signal is processed by the CVB technology, and the results are shown in Figure 13. From Hebda-Sobkowicz et al. [25], the ORFB is determined by the threshold value technology, which is a third smaller than the max value and is marked by the red line. Based on this, the ORFB is [9550, 10,700] Hz. The corresponding SES is shown in Figure 13c, illustrating that the fault feature is still interfered by the noise, and a large gap exists between it and Figure 9b.   Figure 12a, the second sub-signal has the max CK, and its SES is shown in Figure 12b. From Figure 12b, many interference components are found, and the feature is not as clear as the results shown in Figure 9b. To highlight the superiority of the ICFGF, the outer fault signal is also processed by FK, EEMD, and CVB technology. The results by FK are displayed in Figure 11. From Figure 11a, the ORFB is [8533, 9599] Hz. The filtered signal and the corresponding SES are shown in Figure 11b, which illustrates that the SES is complex and the fault feature is interfered by the noise.  Finally, the outer fault signal is processed by the CVB technology, and the results are shown in Figure 13. From Hebda-Sobkowicz et al. [25], the ORFB is determined by the threshold value technology, which is a third smaller than the max value and is marked by the red line. Based on this, the ORFB is [9550, 10,700] Hz. The corresponding SES is shown in Figure 13c, illustrating that the fault feature is still interfered by the noise, and a large gap exists between it and Figure 9b. Finally, the outer fault signal is processed by the CVB technology, and the results are shown in Figure 13. From Hebda-Sobkowicz et al. [25], the ORFB is determined by the threshold value technology, which is a third smaller than the max value and is marked by the red line. Based on this, the ORFB is [9550, 10,700] Hz. The corresponding SES is shown in Figure 13c, illustrating that the fault feature is still interfered by the noise, and a large gap exists between it and Figure 9b.

Ball Fault Diagnosis
The data come from the Case Western Reserve University Bearing Data Center Website [29] and had been applied in [1]. The test rig is displayed in Figure 14. The test ball fault bearing is 6205-2RS JEM SKF, which is installed at the drive end. The vibration data are collected by an accelerometer which is near the test bearing. The sample frequency is 48 kHz and the length of the signal used in this study is 1 s. The ball spin frequency (BSF) is 70.3 Hz because the shaft rotation speed is 1797 rpm. Importantly, the feature frequency of the ball fault should be 140.6 Hz because the local defect of the ball will touch with the inner race and outer race alternately.

Ball Fault Diagnosis
The data come from the Case Western Reserve University Bearing Data Center Website [29] and had been applied in [1]. The test rig is displayed in Figure 14. The test ball fault bearing is 6205-2RS JEM SKF, which is installed at the drive end. The vibration data are collected by an accelerometer which is near the test bearing. The sample frequency is 48 kHz and the length of the signal used in this study is 1 s. The ball spin frequency (BSF) is 70.3 Hz because the shaft rotation speed is 1797 rpm. Importantly, the feature frequency of the ball fault should be 140.6 Hz because the local defect of the ball will touch with the inner race and outer race alternately. The time waveform of the raw signal and its FFT spectrum are shown in Figure 15. The time waveform of the raw signal and its FFT spectrum are shown in Figure 15. The time waveform of the raw signal and its FFT spectrum are shown in Figure 15. The impulse component cannot be directly observed in both. The proposed metho is applied to process this signal. First, the center frequency of the ORFB is 3059 Hz base on the CV index, which is shown in Figure 16.  The impulse component cannot be directly observed in both. The proposed method is applied to process this signal. First, the center frequency of the ORFB is 3059 Hz based on the CV index, which is shown in Figure 16. The impulse component cannot be directly observed in both. The proposed method is applied to process this signal. First, the center frequency of the ORFB is 3059 Hz based on the CV index, which is shown in Figure 16.    Similar to the last part, the signal is also analyzed by FK, EEMD, and CVB technology. Figure 19 displays the results from FK. Based on Figure 19a, the ORFB is [13,124,13,499] Hz. The corresponding filtered signal and its SES are shown in Figure 19b. The SES is very clear. However, only the rotating frequency ( r f ) and its harmonics can be found in it. Consequently, FK fails to extract the fault feature.    Similar to the last part, the signal is also analyzed by FK, EEMD, and CVB technology. Figure 19 displays the results from FK. Based on Figure 19a, the ORFB is [13,124,13,499] Hz. The corresponding filtered signal and its SES are shown in Figure 19b. The SES is very clear. However, only the rotating frequency ( r f ) and its harmonics can be found in it. Consequently, FK fails to extract the fault feature.  Similar to the last part, the signal is also analyzed by FK, EEMD, and CVB technology. Figure 19 displays the results from FK. Based on Figure 19a, the ORFB is [13,124,13,499] Hz. The corresponding filtered signal and its SES are shown in Figure 19b. The SES is very clear. However, only the rotating frequency ( f r ) and its harmonics can be found in it. Consequently, FK fails to extract the fault feature.   Similar to the last part, the signal is also analyzed by FK, EEMD, and CVB technology. Figure 19 displays the results from FK. Based on Figure 19a, the ORFB is [13,124,13,499] Hz. The corresponding filtered signal and its SES are shown in Figure 19b. The SES is very clear. However, only the rotating frequency ( r f ) and its harmonics can be found in it. Consequently, FK fails to extract the fault feature.   The EEMD is applied to analyze the signal, and the results are shown in Figure 20. From Figure 20a, the first sub-signal has the max CK. Its SES is shown in Figure 20b. It is very difficult to observe the fault feature due to the existence of noise. The EEMD is applied to analyze the signal, and the results are shown in Figure 20. From Figure 20a, the first sub-signal has the max CK . Its SES is shown in Figure 20b. It is very difficult to observe the fault feature due to the existence of noise. Finally, CVB technology is used to extract the fault feature, and the results are shown in Figure 21. From Figure 21b, the ORFB should be [3234, 3421] Hz. The SES for the filtered signal is shown in Figure 21c, which illustrates the difficulty in observing the fault feature. Therefore, the three methods cannot catch up with the performance for extracting the bearing fault feature of the proposed method. Finally, CVB technology is used to extract the fault feature, and the results are shown in Figure 21. From Figure 21b, the ORFB should be [3234, 3421] Hz. The SES for the filtered signal is shown in Figure 21c, which illustrates the difficulty in observing the fault feature. Therefore, the three methods cannot catch up with the performance for extracting the bearing fault feature of the proposed method. The EEMD is applied to analyze the signal, and the results are shown in Figure 20. From Figure 20a, the first sub-signal has the max CK . Its SES is shown in Figure 20b. It is very difficult to observe the fault feature due to the existence of noise. Finally, CVB technology is used to extract the fault feature, and the results are shown in Figure 21. From Figure 21b, the ORFB should be [3234, 3421] Hz. The SES for the filtered signal is shown in Figure 21c, which illustrates the difficulty in observing the fault feature. Therefore, the three methods cannot catch up with the performance for extracting the bearing fault feature of the proposed method.

Conclusions
To extract the bearing fault feature, this study proposes a new method to determine the ORFB, which includes a novel energy spectrum statistic index and an MMA optimization algorithm. The novel spectrum index is built based on the differences between the

Conclusions
To extract the bearing fault feature, this study proposes a new method to determine the ORFB, which includes a novel energy spectrum statistic index and an MMA optimization algorithm. The novel spectrum index is built based on the differences between the energy spectrum of Gaussian noise and the energy spectrum of the periodic impulses. The energy spectrum statistic technology is applied to evaluate the distribution of the energy spectrum, which can suppress the interference from the harmonic component effectively. Thus, the center frequency of the ORFB can roughly be obtained by using the energy spectrum statistic index. Then, the MMA optimization algorithm, which has fast convergence, is proposed to determine the ORFB. Finally, the fault feature is obtained by SES technology. The effectiveness of the proposed method is verified by the signals from an outer fault bearing and a ball fault bearing. Furthermore, comparing it with FK, EEMD, and CVB technology, our method better extracts the bearing fault feature from a heavy noise signal. Consequently, this study offers a new method to extract the bearing fault feature, especially for a processing signal which includes non-cyclic impulse interferences (i.e., crushing machinery).