A Novel Method for Adaptive Multiresonance Bands Detection Based on VMD and Using MTEO to Enhance Rolling Element Bearing Fault Diagnosis

Vibration signals of the defect rolling element bearings are usually immersed in strong background noise, which make it difficult to detect the incipient bearing defect. In our paper, the adaptive detection of the multiresonance bands in vibration signal is firstly considered based on variational mode decomposition (VMD). As a consequence, the novel method for enhancing rolling element bearing fault diagnosis is proposed. Specifically, the method is conducted by the following three steps. First, the VMD is introduced to decompose the raw vibration signal. Second, the one ormoremodeswith the information of fault-related impulses are selected through the kurtosis index.Third,MultiresolutionTeager EnergyOperator (MTEO) is employed to extract the fault-related impulses hidden in the vibration signal and avoid the negative value phenomenon of Teager Energy Operator (TEO). Meanwhile, the physical meaning of MTEO is also discovered in this paper. In addition, an idea of combining the multiresonance bands is constructed to further enhance the fault-related impulses. The simulation studies and experimental verifications confirm that the proposed method is effective for identifying the multiresonance bands and enhancing rolling element bearing fault diagnosis by comparing with Hilbert transform, EMD-based demodulation, and fast Kurtogram analysis.


Introduction
Rolling element bearings are widely used in rotating machinery to support rotating shafts, and the major cause of machinery breakdown is the bearing failure.Hence, it is necessary to detect bearing faults at an early stage.However, the rolling element bearing early incipient defect feature is very weak for reasons of being buried in the strong background noises and the interference of the rotating frequency.Besides, there exist the severe signal attenuation phenomenon between the fault source and the sensor collecting the fault signal if the sensor is placed far from the fault-related location.Currently, the fault diagnosis of rolling bearing early weak fault is not only a hot area but also a difficult area [1].
Rolling element bearings usually consist of an inner race, an outer race, several rollers, and a cage.When the surface of one or more of these components develops a localized fault, the impacts generated excite the resonant frequencies of the bearing and adjacent components and induce a modulating phenomenon [2,3], which is the basis of bearing fault diagnosis.Vibration signals collected from bearings carry the rich information on machine health conditions.Therefore, the vibration-based methods have received intensive study during the past decades.It is possible to obtain vital characteristic information from the vibration signals through the use of signal processing techniques [4,5].In order to extract the transient features from the vibration signals, different signal processing techniques have been developed in the area of rotary machine fault diagnosis, such as the wavelet analysis [6], empirical mode decomposition (EMD) [7], time-frequency analysis (TFA) [8], sparse decomposition [9][10][11], manifold learning [12], spectral kurtosis (SK) [13,14], cyclic spectral analysis [15,16], and envelope analysis [17][18][19][20].In fact, the essence of most methods for the weak fault diagnosis is to detect the resonance bands excited by bearing defect.However, by the optimal frequency band selection methods, only the noises outside the selected frequency band are removed from the original signal, while those inside the 2 Shock and Vibration selected frequency band cannot be wiped off effectively.As a result, the performances of these methods may be poor in the presence of low signal-to-noise ratio and each of these methods leaves something to be desired and some even perform badly in analyzing weak vibration signal.For example, the same class decomposition models of EMD are mostly limited by their algorithmic ad hoc nature lacking mathematical theory, the recursive sifting which does not allow for backward error correction, and the inability to properly cope with noise [21].As for the wavelet transform, if the wavelet function is selected properly, the defect-related features may be well extracted [5].Therefore, the constructing wavelet function that adaptively matches the defect-related characteristics of vibration signal is one of the key issues.However, the way of constructing wavelet function needs the prior knowledge and the hard band-limits of wavelet transform is also another inevitable shortcoming of it.In addition, the resonance bands excited by bearing defect are often more than one, while the conventional methods may fail to uncover the accurate multiresonance bands.
Recently, Dragomiretskiy and Zosso [22] proposed a new decomposition model called VMD that determined the relevant bands adaptively and estimated the corresponding modes concurrently, thus properly balancing errors between them.VMD is a fully intrinsic and adaptive variational method which is motivated by the narrow-band properties corresponding to the current common intrinsic mode function (IMF) definition and looks for an ensemble of modes that reconstruct the given input signal optimally (either exactly, or in a least-squares sense).Each mode being band-limited about a center frequency is estimated on-line.Specifically, the method can address the presence of noise in the input signal in which the tight relations of it to the Wiener filter actually suggest that it has some optimality in dealing with noise.Subsequently, Mohanty et al. [23] employed the VMD to decompose the vibration signal in several modes, extracted the energy feature of it to diagnose whether the bearing doped with sand or not, and compared the method with EMD to verify its performance.Wang et al. [24] used the VMD to detect the rub-impact fault of the rotor system and verify that its performance outperforms the conventional decomposition methods such as empirical wavelet transform (EWT) [25], EMD [21], and EEMD [26].To our knowledge, there is no report in the literature so far on its applications to weak bearing vibration signal analysis and its property of detecting multiresonance bands is also not discovered up till now.Thus, the method of VMD is firstly applied to detect the resonance bands for taking full advantage of the resonance information in vibration signal of defective bearing in our paper.
After the resonance bands are determined in the bearing vibration signal where the fault-related periodic impulse is a modulator to the high natural frequencies of the machine, the demodulation techniques should be used to demodulate the impact impulses from the resonance modes obtained by VMD.Many demodulation methods have been studied, such as FFT-based Hilbert transform [21,27], wavelet-based [28,29], and TEO [17,18].Among them, the TEO method is an attractive demodulation method proposed and developed by Maragos et al. [30][31][32], is a sort of time-frequency analyzer, and has been used for many applications such as speech processing, image processing, and AM/FM demodulation.Compared with HT method, TEO method is totally based on the local differential operation without involving integral transform, so it has a better localization property and lower computational complexity.The TEO is also known to be sensitive to spikes, where a spike means that a signal is concentrated in a short time interval and at a high frequency band.With these advantages, the TEO method has been also introduced into machinery fault diagnosis.Lin et al. [33] utilized TEO for resonance demodulation analysis to extract fault characteristic of roller bearings.Junsheng et al. [34] presented a TEO demodulation approach based on EMD to diagnose machinery fault.Liang and Bozchalooi [17] introduced a repetitive application of TEO on detecting the fault characteristic frequency in the spectrum of the energytransformed signal.However, under a low signal-to-noise ratio (SNR) or in background noise at high frequencies, the TEO gets more sensitive to high noisy peaks than to the true fault-related impulses, and the performance of TEO as an impulse detector degrades rapidly.Besides, the extracted envelope waveform by these demodulation methods is at a single scale [35].Thus, Choi et al. [35,36] proposed the MTEO method to detect the action potential of neural signals which could make up for the weakness by tuning the TEO to the frequency of action potentials with the resolution parameter.The analysis results of experimental data in [35,36] have shown that the MTEO method outperforms TEO and the other conventional demodulation methods in handling both noise and interferences.However, the MTEO physical meaning is not given in paper [35,36].We will further develop and uncover the MTEO physical meaning and then firstly employ MTEO to enhance the fault-related impulses in bearing vibration signal.And yet, the enhanced vibration signal still confronts the contamination of in-band noise.The optimal smoothing window, hamming window, together with MTEO is used to further enhance the faultrelated impulses.
Furthermore, we can consider that the noise components are less correlated with the different resonance bands while the impulses of envelope waveform are still correlated.Therefore, the combination of more possible resonance bands could be further beneficial to in-band noise removal.Eventually, we propose the novel method for enhancing faultrelated impulses in bearing vibration signal by combining the VMD and MTEO.The simulation studies and experiment verification on an experimental rolling element bearing will also be conducted to test the improved performance of the proposed method and compare with the conventional FFTbased Hilbert transform, EMD-based demodulation [27], and fast Kurtogram analysis [2,13].
This paper is organized as follows.Section 2 describes the details on the theoretical background of the proposed method for the enhanced fault diagnosis of rotating machines.The MTEO physical meaning is also given and uncovered in this section.In Section 3, the simulation studies of the proposed method are studied where there are two cases, single resonance band and double resonance bands.Then in Section 4, the practical applications to bearing defect identification are conducted to verify the effectiveness of the proposed method.A discussion is given in Section 5.The conclusions are finally drawn in Section 6.

Variational Mode Decomposition (VMD).
The VMD [22] is a recently developed methodology for adaptive signal decomposition which decomposes an input signal into  discrete number of subsignals (modes)   .These modes have specific sparsity properties while reproducing the input, where each mode has limited bandwidth in the spectral domain.Each mode   is required to be mostly compact around a center pulsation   determined along with the decomposition.The VMD to process the input signal is given as follows.
(1) For each mode   , compute the associated analytic signal by means of the Hilbert transform to obtain a unilateral frequency spectrum.
(2) Shift the frequency spectrum of each mode to the baseband by mixing with an exponential tuned to the estimated center frequency, respectively.
(3) Estimate the bandwidth through the  1 Gaussian smoothness of the demodulated signal, that is, the squared  2 -norm of the gradient.As a consequence, the constrained variational problem is given by min where  is the input signal,  is the number of modes,   is the th mode,   is the central frequency of th mode,  is the Dirac distribution,  is the time script, and * denotes convolution.
(4) Both a quadratic penalty term and Lagrangian multipliers  to render the problem unconstrained are used as the reconstruction constraints.The combination of both terms thus benefits greatly from the nice convergence properties of the quadratic penalty at finite weight and the strict enforcement of the constraint by the Lagrangian multiplier .Therefore, the solution to find the optimal   with the inclusion of Lagrange multipliers  and quadratic penalty is given by where  denotes the balancing parameter of the data-fidelity constraint and the Lagrangian multiplier  is a common way of enforcing constraints strictly.
(5) Alternate direction method of multipliers (ADMM) optimization algorithm [37][38][39] is used to solve (2) to produce different decomposed modes and the center frequencies of these modes during each shifting operation.The procedures of these operations are shown in following.
(i) To update the modes   , the subproblem is formulated as (3) which is rewritten as the equivalent minimization problem shown in (4): Subsequently, the solution to find the optimal   in spectral domain is given by (5).Moreover, the two regularized terms are written as half-space integrals over the nonnegative frequencies based on exploiting the Hermitian symmetry of the real signals in the reconstruction fidelity term.As a result, the solution of this quadratic optimization problem is readily found as (6) by forcing the first variation vanish for the positive frequencies which is identified as Wiener filter.The decomposed mode in time domain is eventually obtained by the inverse Fourier transform of the filtered analytic signal: (ii) To update the central frequency   , the other subproblem is formulated as (7).Because the center frequency   does not appear in the reconstruction fidelity term, but only in the bandwidth prior, the relevant problem is thus read as (8): As before, the optimization can also take place in Fourier domain, and the solution to find the optimal   in spectral Shock and Vibration domain is given by (9).Therefore, this quadratic problem is easily solved as (10): which puts the new   at the center of gravity of the corresponding mode's power spectrum.The mean carrier frequency is the frequency of a least squares linear regression to the instantaneous phase observed in the decomposed mode.The complete algorithm of the VMD in detail can be found in [22].
It is applied on a continuous signal () which is defined as The definition of the original discrete-time TEO is then given by The derived version of this operator has been used to separate the frequency and amplitude modulations of signals.
It is logic to conjecture that the extraction of certain information relying on the combined amplitude and frequency demodulations could be carried out directly by this operator [17].When we consider an arbitrary signal () defined as (13), the corresponding discrete-time equivalent of the energy operator is given as (14): where () = () cos(()) and Ω() = () − ( − 1).
According to (14), the TEO extracts both the amplitudemodulation (AM) and frequency-modulation (FM) information of the signal.Although the energy operator has mainly been used to separate the amplitude and frequency modulations of a given signal, the separation of such information is not necessary in the context of machinery fault diagnosis [17].The information of interest for machinery fault detection is the transient nature of the fault impulses resulting from both amplitude and frequency modulations.It is well known that TEO is sensitive to transient impact, where an impact means a signal that is concentrated in a short time interval and at a high frequency band, and valuable means of accentuating the transient fault characteristics relative to the other components of vibration signal such as gear meshing, shaft imbalance vibration, and noise component.However, the TEO gets more sensitive to high noisy peaks than to the true impact under a low SNR or background noise of the obtained signal at high frequencies, and the performance of the TEO as an impact detector degrades rapidly.Moreover, the negative value phenomenon which is nonphysical easily arises in the transformed signal by TEO.In the following subsection, the MTEO is introduced to avoid these drawbacks.

Multiresolution TEO (MTEO).
In paper [35,36], MTEO is proposed to accurately identify the action potentials in the neural signal by tuning TEO to the frequency range of the impacts with the multiresolution parameter  ( ≥ 2).The definition of the discrete-time MTEO is given in the following sense: However, the physical meaning of ( 15) is not given clearly in [35,36].Hereby, we deduce the essential mechanism to deeply uncover (15).We define differentiation operator   which is useful to suppress low frequency interferences and noise, integration operator   which can enhance the signal in the presence of noise or increase signal-to-noise ratio (SNR) due to its smoothing effects, and composite operator Δ as follows: where the   ,   , Δ  represent the  order differentiation operator, integration operator, and composite operator, respectively.It should be noted that the composite operator Δ  is only defined as the two formations in (18).According to (11), a new transform operator Ψ  [⋅] is formed as where and Δ  (Δ  (())) is shown in the following: As a result, ( 15) can be obtained by using (19).It can be concluded from (19) that MTEO is formed by the composite operator which consists of differentiation and integration operators.Therefore, MTEO takes advantages of two original operators for suppressing low frequency interferences and increasing SNR, besides the virtue of TEO.
Another important element that affects the impact detector performance is the smoothing window.When using the -TEO as a tool of enhancing impulses, it only uses three samples to calculate an output value at a time instance.One prominent noise sample can induce a peak at the output and disturb the accurate fault impact.When the SNR is low, such noisy peaks cause serious problems and must be removed by a smoothing window.As suggested by paper [35], the hamming window was suitable for the matched filter of impulses and paper [36] also gave a logical interpretation to the smoothing window using the matched filter theory, and both of them recommended the hamming window with length of 4 + 1 as an optimal window for the -TEO.In addition, after the components obtained from VMD are processed by -TEO, they still contain some noises around the impulses.Therefore, we select the hamming window as the smoothing window to extract the more pure impulses.This optimal smoothing window, together with -TEO, is used to process the faultrelated signal.We can adjust the -TEO to be sensitive to the frequency of the fault-related impacts.When we take the effect of the output window into account, the windowed output of a white Gaussian noise () is where ℎ  () denotes the th coefficient of the smoothing window matched to the -TEO.
According to the central limit theorem,   () is approximated as a Gaussian random variable with mean  2 ∑ ℎ() and variance 3 4 ∑ ℎ 2 ().Therefore, the mean-square value of the windowed noise () is The windowed output can be regarded as a signal of defect-related feature mixed in the background noise   ().However, the value of { 2  ()} could change at different resolution parameters  where some false impacts may be caused by the in-band noise.Therefore, let the window coefficients be (23) to normalize the output noise power to a constant value: For MTEO, all these elements, that is, differentiation, integration and energy operator, and hamming window, are implemented by a simple formula in one step.The main attractive advantages of MTEO include its simplicity, computational efficiency, excellent time resolution, and the leavingout of the band-pass filtering process.As such, it is suited to on-line bearing fault detection in a noisy environment with multiple vibration interference.

Summary of the Proposed Method.
In summary, the procedure of the proposed method can be described briefly as follows: (1) Conduct VMD on the measured raw vibration signal to transform it into several modes.Considering that fault-related resonance bands are identified preferentially by VMD, the number of the resonance bands may be more than one in the vibration signal of early stage defect generally.Therefore, the number of decomposed modes by VMD is set as at least 2, but should not be too large for saving the computing time.
(2) Calculate the kurtosis value of these decomposed modes.
(3) Employ MTEO and smoothing window to enhance the fault-related waveforms of the decomposed modes with a large kurtosis value on four scales with spacing 3. The four scales could cover the analysis frequency range as suggested in paper [35,36].
(4) Observe the enhanced time domain waveform and its frequency spectrum from the four enhanced waveforms to diagnose the faulty bearing.If there are multiresonance bands in the obtained modes, the combination of them will be used to further enhance the fault-related waveform.
Since the VMD could adaptively extract the modes in the resonance bands, the proposed method does not need more prior knowledge.In addition, the MTEO is a nonlinear method.Therefore, it can reveal the nonlinear transient envelope information of machinery systems by combining one or more decomposed modes.Thus the result analyzed by the proposed method is taken for further pure impulse and spectrum analysis to identify the fault-related characteristic and locate the defect position in a machine.The frame of the proposed method is shown in Figure 1.

Simulated Defect-Related Signal with Single Resonance
Frequency Band.To verify the performance of the proposed method, a simulated signal with single resonance frequency band [40,41] is first constructed.In general, the vibration signal of a faulty bearing consists of periodic impulses and white noise.The simulated signal can be computed by the following formula: where   = 1000 Hz is the central frequency of the resonance band and  = 0.02 is the damping ratio. = 1 represents the initial magnitude of a single free vibration;  = 0.02 s is the fault repetition period, and hence the characteristic  frequency of this faulty signal is   = 50 Hz. = 14 stands for the number of impulses; () is the added white noise.Figure 2(a) shows the pure simulated signal with the sampling frequency   = 12000 Hz.It can be seen from the waveform that these resonance impulses derived from the impacts at the localized fault of a machine exist periodically at the interval of  = 0.02 s.However, in the presence of noise, the result may be different.Figure 2(b) shows the simulated signal with the additive white noise resulting in a SNR of 0 dB by using the MATLAB function AWGN(, SNR) in which the power of  defaults to 0 dBW.It can be seen that the noise corruption makes the periodic transient impulses difficult to identify from the waveform.As shown in Figure 2(c), the resonance band in the power spectrum is also completely distorted.As described in Figure 3(a), the Hilbert transform of the noisy signal is obtained.It can be seen that the envelope waveform of fault-related impulses cannot be discovered in the envelope waveform of Figure 3(a), and there still exist a large number of interfering components.It can also be seen from Figure 3(b) that the envelope spectrum cannot show the faulty characteristic frequency.
The proposed method is introduced to analyze the simulated single resonance signal shown in Figure 2(b).First, the VMD is addressed on the signal and five modes are provided.Figure 4(a) shows the power spectra of the five decomposition modes.As shown in Figure 4(b), the five modes converge to their central frequencies where the central frequency of the second mode is 1100 Hz, which approximates the true carrier frequency 1000 Hz of the resonance band.
It can be seen from Figure 4(c) that the second mode has the maximum kurtosis value and the formula of kurtosis is shown in (25).As displayed in Figure 4(d), the waveform of the second mode demonstrates that its fault-related impulses are more apparent to compare with Figure 2(b), but still contaminated by the noise.Then, the second mode is further processed by MTEO with different resolution parameters.It can be seen from the waveform of Figure 4(e) that the faultrelated impulses are more distinct with  = 7.The frequency spectrum of Figure 4(f) shows that the fault-related frequency and its harmonics are very clear.Subsequently, the second mode is also analyzed by TEO.The result is provided in Figure 4(g).We can find that the in-band noise even corrupts the fault-related impulses and the negative values which are nonphysical arise in the transformed signal.Therefore, it is verified that the TEO does not work as well as MTEO.The more detailed explanation for causing the negative values can be seen in the paper [42]: where  is the mean of signal ,  is the standard deviation of , and (⋅) represents the expected value.As a comparison, the simulated single resonance signal is also analyzed by the EMD-based demodulation method [27], which employs the EMD as an adaptive filter followed by the Hilbert envelope analysis on the decomposed IMF that contains the characteristic frequency component.As described in Figures 5(a We can find that the noise component completely corrupts the fault-related impulses in Figure 6(c).As displayed in Figure 6(d), the characteristic frequency in envelope spectrum is not obvious enough by comparing with the proposed method.
In order to further investigate the improved performance of the proposed method in resonance band detection, the fast Kurtogram [2,13] is used to analyze the same simulated signal shown in Figure 2(b).The paving of the fast Kurtogram is shown in Figure 7(a), where an optimal filter with a centre frequency of 5812.5 Hz and a bandwidth of 375 Hz is automatically chosen.The signal filtered by the optimal filter and its envelope spectrum are shown in Figure 7

Simulated Faulty Signal with Double Resonance Frequency
Bands.In this subsection study, the simulated vibration signal with double resonance frequency bands is constructed.The central frequencies of two resonance bands are  1 = 1000 Hz and  2 = 2000 Hz, respectively.The fault characteristic frequency is also   = 50 Hz.The simulated signal can be computed by the following formula: Figure 8(a) shows the pure simulated signal occupying double resonance frequency bands with the sampling frequency   = 12000 Hz.As presented in Figure 8(b), the simulated faulty signal in Figure 8(a) is also mixed with the additive white noise resulting in a SNR of 0 dB by using the MATLAB function AWGN(, SNR) in which the power of  defaults to 0 dBW.It can be seen that the transient feature is almost swallowed by noise signal.The corresponding frequency domain representation of the obtained noise-contaminated signal is illustrated in Figure 8(c).The resonance bands in the power spectrum are also corrupted.As represented in Figure 9(a), the Hilbert transform is used to analyze the noise-contaminated signal.It can be seen that the envelope waveform of fault-related impulses cannot be seen in Figure 9(a).Similarly, as delineated in Figure 9(b), the faulty characteristic frequencies are contaminated by the interference components marked in the rectangles which make some trouble in the defect detection.
The proposed method is employed to process the simulated double resonance noisy signal.To begin with, VMD technique is applied to decompose the noise-contaminated signal and five modes are obtained.Figure 10(a) shows the power spectra of the decomposed components by VMD, in which the components located in the resonance bands are successfully extracted.As shown in Figure 10(b), the central frequencies of the second and third modes are 1067 Hz and 1964 Hz, respectively, which approximate the carrier frequencies 1000 Hz and 2000 Hz assumed in the simulated signal.Figure 10(c) shows the kurtosis curve of the five modes.It can be observed from Figure 10(c) that the third mode has the maximum kurtosis value and the second mode has the second maximum kurtosis value.As displayed in Figures 10(d) and 10(e), the second and third modes demonstrate that the fault-related impulses in the waveform are apparent compared with the result of Figure 8(b).Nevertheless, they are still contaminated by the noise.Then, the two modes are analyzed by TEO.The results are provided in Figures 11(a) and 11(b).It still can be seen that the fault-related impulses are incorporated into in-band noise and the negative values which are nonphysical also arise in the transformed signal.Therefore, the two modes obtained by VMD are further analyzed by MTEO with different resolution parameters, respectively.It can be observed from Figures 11(c) and 11(e) that the fault-related impulses in the waveforms are more pure and distinct with  = 7.As shown in Figures 11(d show that the fault-related frequency and its harmonics are very clear.However, some pulses are still weak as shown in the rectangles of Figures 11(c) and 11(e).The idea that the combination of the second and third modes is designed to produce a new mode is constructed.As delineated in Figure 11(g), the new mode has more apparent impacts than both of the second and third modes.As provided in Figure 11(h), the frequency spectrum of new mode shows the more clearly fault-related frequency and its harmonics.For comparison purpose, the method of EMD-based demodulation is also used to analyze the vibration signal in Figure 8(b).The first four IMFs decomposed by EMD are shown in Figure 12(a).As displayed in Figure 12(b), these IMFs locate in the high frequency ranges and the frequency bands of IMF2 and IMF3 approximate the 1000 Hz and 2000 Hz, respectively.The frequency bands of the rest IMFs are far away from the 1000 Hz shown in Figure 12(c).In addition, it can be seen from Figure 13(a) that there are 13 IMFs shown in the kurtosis curve.Although the IMF9 has the maximum kurtosis value, it locates in the low frequency region and does not contain fault-related information.Therefore, the IMF2 and IMF3 owning the larger kurtosis values in high frequency regions are selected.However, the waveforms of them only show the slight impacts.Both of them are further analyzed by the Hilbert transform.We cannot see the faultrelated impulses from the envelope waveforms in Figures 13(b) and 13(d).Similarly, as described in Figures 13(c) and 13(e), the envelope spectra of the IMF2 and IMF3 cannot display the accurate fault characteristic and its harmonics which are severely contaminated by noise.
The fast Kurtogram is also applied to analyze the same mixed signal with double resonant frequency shown in Figure 8(b).The paving of the fast Kurtogram is plotted in Figure 14(a), where it is indicated that the optimal filter has a centre frequency of 1875 Hz and a bandwidth of 750 Hz.The centre frequency of 1875 Hz deviates from the frequency of 2000 Hz which is one of the original setting natural frequencies.But our proposed method can accurately identify both of the original setting natural frequencies.The signal filtered by the optimal filter and its envelope spectrum show the fault-related signatures in Figure 14(b).It is found that the visual inspection ability of the fast Kurtogram is not as good as that of the proposed method for showing fault characteristic where some interferential components exist.The extracted impulses by the proposed method are more pure than that of the fast Kurtogram.As such, the effectiveness of the VMD technique in detecting multiresonance bands of incipient defective bearing is well demonstrated by using the simulated signal.The validity will be further evaluated using the practical vibration signal in the following section.

Experimental Verification
To verify the effectiveness of proposed method in practical applications for enhancing the fault diagnosis of rotating machines, the experimental data from a defective bearing are analyzed in this section.The experiment is carried out on a deep groove ball bearing.The kinematical parameters and the corresponding fault frequencies of rolling element bearings are listed in Tables 1 and 2, respectively.Experimental data with the sample frequency 97656 Hz, input shaft rate   = 25 H, and under the load of 1200 N are provided by the MFPT (Mechanical Failure Prevention Technology) [43].
A vibration signal measured from the experimental bearing is plotted in Figure 15(a), and its power spectrum is described in Figure 15  observed from the envelope waveform.However, the most fault-related impulses are masked by noise and not evident enough to detect the existence of bearing fault.From the envelope spectrum shown in Figure 15(d), the frequency 78 Hz which approximates the third harmonic of the rotational frequency is quite dominant but the fault frequency of the outer race   = 81.13cannot been seen.This means that the envelope analysis fails to discover the fault characteristics and could cause a misjudgment about the outer race fault bearing in noisy signal.Subsequently, we apply the proposed method to process the vibration signal in Figure 15(a).First, the VMD is addressed on the raw vibration signal and five modes are obtained.Figure 16(a) shows the power spectra of five decomposition modes.The central frequencies of five modes are depicted in Figure 16(b).The kurtosis curve of five modes is also described in Figure 16(c).It can be also seen from Figure 16(c) that the first three modes have the large kurtosis value.As shown in Figure 16(b), the central frequencies of the first three modes are 2666 Hz, 9268 Hz, and 18230 Hz, respectively.Compared with Figures 15(a) and 15(c), the first three modes displayed in Figures 16(d), 16(e), and 16(f) demonstrate that the fault-related impulses in the waveform are more apparent but still contaminated by the noise.Then, the first three modes are processed by MTEO with different resolution parameters.It can be seen from Figures 17(a) and 17(c) that the fault-related impulses of the enhanced waveforms of first two modes are more distinct in  = 10.The frequency spectra of them show that the fault-related frequency and its harmonics are very clear in Figures 17(b) and 17(d).However, the enhanced waveform of the third mode does not illustrate the remarkable faultrelated impulses compared with the first and second modes in Figures 17(a) and 17(c).In addition, the frequency spectrum of Figure 17(f) only shows the fault-related frequency but it cannot display its harmonics.According to the above analysis, we can deduce that this faulty experimental bearing could excite two resonance bands in this experiment and the centers of the resonance bands approximate 2666 Hz and 9268 Hz, respectively.The resonance band with central frequency 18230 Hz may be slightly excited.Although the enhanced waveforms of first two modes have more distinct impulses, Besides, it can also be seen that the fault-related impulses are incorporated into the in-band noise and the negative values which are nonphysical also arise in the transformed signal by TEO.All of these results verify that the method of MTEO can obtain the more pure impulses and suppress the effect of negative value produced by TEO based on the extracted resonance bands.
In addition, the method of EMD-based demodulation is also used to analyze the raw vibration signal in Figure 15(a).There are fifteen IMFs obtained by EMD decomposition.Considering that the frequency ranges of the IMFs after IMF8 which have nothing to do with the bearing fault information is less than 1000 Hz, we show only the time waveforms and spectra of the first eight IMFs in Figure 18.It can also be seen that these IMFs have a heavy mode mixing with each other.Then, the kurtosis curve of all IMFs is given in Figure 19(a).There are two regions which have the larger kurtosis value shown in the rectangles.The high frequency components belong to the first region where the fault-related information is contained.Instead, the second region only has some low frequency components.Moreover, the envelope waveforms and their spectra of IMF3-IMF7 which have more impacts than the rest IMFs are shown in Figures 19(b)-19(k), respectively.We can see from these envelope waveforms that they cannot show the regular impacts.Similarly, the envelope spectra of them cannot also exhibit the accurate frequency of outer race defective bearing and are severely disturbed by the noise.
For further comparison, the fast Kurtogram is applied to analyze the same rolling element fault signal shown in Figure 15(a).The paving of the fast Kurtogram is shown in Figure 20(a), where an optimal filter with a centre frequency of 8138 Hz and a bandwidth of 16276 Hz is automatically defined.The signal filtered by the optimal filter and its envelope spectrum are shown in Figure 20(b).It can be concluded that the impulses extracted by the fast Kurtogram have a lower SNR than that extracted by the proposed method as shown in Figure 17(g) and the envelope spectrum obtained by our proposed method is more clear by comparing Figures 17(b), 17(d), and 17(h) with Figure 20(b).Furthermore, as shown in Figure 21(a), a vibration signal with the more serious impact waveform which is collected from the same test bearing as the signal in Figure 15(a) is also provided by the MFPT.To demonstrate the outstanding performance of the proposed method in identifying the multiresonance bands, the time-frequency representation of the signal in Figure 21(a) which could make a rough guide for the locations of three resonant frequencies is computed by the Short Time Fourier Transform (STFT) [44].As displayed in Figure 21(b), the locations of these resonant frequencies approximate the values identified by the proposed method.Besides, Dr. Eric Bechhoefer also pointed out that the faulty bearing shows a natural frequency about 9000 Hz in [43].Consequently, the proposed method is better than the fast Kurtogram for indicating the existence of the multiresonance bands.
Finally, an outer race defect is discovered in test bearing as shown in Figure 22 which is provided by MFPT [43].It also means that the proposed method is validated and able to more   clearly extract fault characteristics and accurately identify the multiresonance bands than the conventional methods.

Discussions
The effectiveness of the proposed method is verified by simulation cases and engineering applications in Sections 3 and 4, respectively.The proposed method can extract the useful impulse features hidden in the signal.Here, a brief discussion is given on the proposed method.Specifically, from the observation of the analysis results, the time waveforms extracted by the proposed method show the more pure impulses and the frequency spectra of the enhanced waveforms describe the clear fault-related frequency and its harmonics, while the other contrastive methods can merely extract partly periodical impulse features or fail to detect them.The selections of key parameters are also discussed in detail as follows.
For the VMD method, the fault-related resonance bands which are identified preferentially by VMD usually are one or more in vibration signal at early stage defective bearing.Therefore, the number  of the decomposed modes in VMD is set as at least 2.However, the much time will be consumed with increasing of .Considering that the vibration signals at early stage defective bearing only have a few resonance bands, we suggest that the number  of the decomposed modes should be set between 2 and 7.For the MTEO method, the resolution  of MTEO ranges from 1 to 10 with spacing 3 in our paper, which could cover the analysis frequency ranges as suggested in paper [35,36].Therefore, we can directly and easily observe the enhanced time waveform for the only four situations in MTEO.As a result, the resolution  of MTEO is set as 7 in the two simulated cases, and the resolution  is set as 10 in the experimental study when the superior fault-related impulses are displayed.Although the hamming window with length of 4 + 1 is very effective for purifying the noisy impulses as demonstrated in simulation study and experimental verification, some other window functions are also worth studying for achieving the better pure impulses in the future.As a consequence, the proposed method has the following main contributions and advantages: (1) VMD is firstly uncovered to extract the single resonance band or more resonance bands which contain the fault-related impulses.(2) MTEO is also firstly developed to enhance the fault-related impulses.
(3) The physical meaning of MTEO is studied and given in our paper.(4) The idea of combining multiresonance bands is used to further enhance the fault-related impulses.(5) The vibration signal of defective bearing can be transformed to the pure impact waveform in time domain which would make diagnosis of the bearing fault easier for technical personnel.The results of simulation studies and the experimental applications demonstrate the excellent performance of our proposed method in extracting fault characteristics from the strong noisy vibration signals.For future research, the identified natural frequencies by the proposed method is planned for us to use as the prior knowledge of traditional resonance demodulation method and design a more simple and efficient fault diagnosis method.

Conclusion
Aiming at the accurate extracting of the fault-related resonance bands and enhancing the fault-related impulses, our paper proposes a novel method for enhancing fault diagnosis of rotating machines by combining the VMD and MTEO.The proposed method concerns the situation of multiresonance bands that contain the rich fault-related information and nonlinearity of machinery dynamic systems at the same time.The resonance bands which contain the plentiful faultrelated information can be adaptively determined through the VMD method.Moreover, the MTEO in the proposed method which shows the merits of in-band noise suppression and fault-related impulses preserving can expose the remarkable fault-related impulses structure in machinery fault and the computation of MTEO is also quite efficient by comparing with the conventional method.Hence, the proposed method could process the longer data online.In particular, the combination of multiresonance bands can be used further to enhance the fault-related impulses buried in the vibration signal of the early defective bearing.The effectiveness and superiority of the proposed method are verified by simulation studies and the experimental applications.Furthermore, the performance of the proposed method is verified to outperform the traditional Hilbert-based, EMDbased demodulation methods, and fast Kurtogram analysis.

Figure 1 :Figure 2 :Figure 3 :
Figure 1: Flow chart of the proposed scheme for weak faulty signal.

Figure 4 :
Figure 4: Analyzed results of the simulated signal shown in Figure 2(b) by the proposed method: (a) power spectrum obtained by VMD; (b) the number of iterations with central frequencies; (c) kurtosis curve of five decomposed modes; (d) waveform of the second mode; (e) enhanced waveform of the second mode by MTEO; (f) power spectrum of the enhanced waveform by MTEO; (g) enhanced waveform of the second mode by TEO.
(b) which cannot provide any fault-related information to indicate the existence of the fault characteristic frequency.It can be seen that the fast Kurtogram indicates a wrong location for selecting resonance band.

Frequency 3 Figure 7 :
Figure 7: Analyzed results of the simulated signal shown in Figure 2(b) by the fast Kurtogram: (a) Kurtogram of the simulated signal; (b) envelope of the filtered signal and corresponding demodulation spectrum.

Figure 8 :Figure 9 :Figure 10 :Figure 11 :
Figure 8: Simulated signal with double resonance frequency band: (a) time waveform for the pure signal; (b) waveform for the mixed signal; (c) power spectrum for the mixed signal.

FrequencyFigure 14 :Figure 15 :
Figure 14: Analyzed results of the simulated signal shown in Figure 8(b) by the fast Kurtogram: (a) Kurtogram of the simulated signal; (b) envelope of the filtered signal and corresponding demodulation spectrum.

Figure 16 :
Figure 16: Analyzed results of the experimental vibration signal in Figure 15(a) by the proposed method: (a) power spectrum of five modes decomposed by VMD; (b) the number of iterations with central frequencies; (c) kurtosis value of five decomposed modes; (d) waveform of the first mode with the maximum kurtosis value; (e) waveform of the second mode with the third maximum kurtosis value; (f) waveform of the third mode with the second maximum kurtosis value.

Figure 17 :
Figure 17: Enhanced waveform and its power spectrum by the proposed method: (a) enhanced waveform of the first mode by MTEO; (b) power spectrum of the first mode enhanced waveform by MTEO; (c) enhanced waveform of the second mode by MTEO; (d) power spectrum of the second mode enhanced waveform by MTEO; (e) enhanced waveform of the third mode by MTEO; (f) power spectrum of the third mode enhanced waveform by MTEO; (g) enhanced waveform of the combined mode by MTEO; (h) power spectrum of the combined mode enhanced waveform by MTEO; (i) enhanced waveform of the first mode by TEO; (j) enhanced waveform of the second mode by TEO.

Figure 20 : 2 XFigure 21 :
Figure 20: Analyzed results of the simulated signal shown in Figure 15(a) by the fast Kurtogram: (a) Kurtogram of the experimental signal; (b) envelope of the filtered signal and corresponding demodulation spectrum.

Figure 22 :
Figure 22: Picture of outer race fault bearing.

Table 1 :
Parameters of the test rolling bearing.