Planetary Gearbox Fault diagnosis via Joint Amplitude and Frequency Demodulation Analysis Based on Variational Mode Decomposition

Zhipeng Feng 1 ID , Dong Zhang 1 ID and Ming J. Zuo 2,3,* 1 School of Mechanical Engineering, University of Science and Technology Beijing, Beijing 100083, China; fengzp@ustb.edu.cn (Z.F.); s20150504@xs.ustb.edu.cn (D.Z.) 2 School of Mechatronics Engineering, University of Electronic Science and Technology of China, Chengdu 611731, China 3 Department of Mechanical Engineering, University of Alberta, Edmonton, AB T6G 1H9, Canada * Correspondence: ming.zuo@ualberta.ca; Tel.: +1-780-492-4466


Introduction
Planetary gearboxes are widely employed in many engineering systems, such as helicopters and wind turbines, for their unique merits of high load bearing capacity and large transmission ratio in a compact structure [1][2][3].Since planetary gearboxes are usually key units in drivetrains, it is vital to maintain them working reliably and steadily.Once a fault occurs to planetary gearboxes, it may result in transmission deficiency, or even lead to shutdown of entire machinery.Therefore, planetary gearbox fault diagnosis plays an important role.
To date, researchers have made important contributions to planetary gearbox fault diagnosis.This includes fault vibration mechanism, fault feature extraction, fault pattern identification, and fault prognosis, reported in [4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20], to name a few examples.Through dynamics modeling, Chaari et al. [4,5] analyzed the effect of gear faults (tooth pitting and crack) on the meshing stiffness, and further examined the effects of gear eccentricity and tooth profile error on vibration responses.Mark and Hines [6,7] studied the effects of non-uniform planet loading due to gear imperfections and the modulation effects of planet carrier torque on vibration spectral sidebands.Inalpolat and Kahraman [8,9] investigated the gearbox configuration (such as the number of planets, the planet position phasing, and the number of gear teeth) and the gear manufacturing errors on vibration spectral sidebands.In order to extract fault signatures from vibration signals, McFadden [10,11] studied the sideband asymmetry of planetary gearbox vibrations, and extended the time domain averaging method to planetary gearbox vibration signal analysis.Samuel and Pines [12,13] proposed a method to separate the planet and sun gear vibration signals, and a constrained adaptive wavelet lifting method to analyze individual tooth mesh waveforms.Barszcza and Randall [14] applied the spectral kurtosis method to detect the ring gear tooth crack in a planetary gearbox.Lei et al. [15] improved the adaptive stochastic resonance method and applied it to extract the weak fault symptoms of a planetary gearbox.Elasha et al. [16] analyzed vibration and acoustic emission data using adaptive filter, spectral kurtosis and envelope spectrum, and found that acoustic emission more susceptible to planetary bearing defects through a comparison study.In order to classify planetary gearbox fault patterns, Qu et al. [17] presented a support vector machine based feature selection method, and applied it to planet gear pitting degree classification.Zhao et al. [18] proposed to preserve the ordinal information by ordinal ranking, and improved the accuracy of ordinal ranking by a correlation coefficients based feature selection method.Lei et al. [19] identified the running condition of a planetary gearbox using the multiclass relevance vector machine, with accumulative amplitudes of the planet carrier orders and energy ratio based on difference spectra as indices.As for planetary gearbox fault prognosis, Djeziri et al. [20] proposed a fault prognosis and remaining useful life prediction method for wind turbine systems based on a physical model, data clustering and geolocation principal.These research results have enriched the literature on planetary gearbox fault diagnosis.
Planetary gearbox vibration signals have strong modulation features, due to the amplitude modulation (AM) and frequency modulation (FM) effects of gear fault, as well as the AM effect of time-varying vibration transfer paths, on the gear meshing vibration.In the time domain, they are a product of the multiple AM terms and the FM signal.According to the convolution property of Fourier transform, in the frequency domain, their Fourier spectra amount to the convolution between the respective Fourier spectrum of each AM term and the FM signal.Meanwhile, the Fourier spectrum of an FM signal involves a Bessel series expansion of infinite orders.Therefore, planetary gearbox vibration signal spectra exhibit a very intricate sideband structure.
Planetary gearbox fault diagnosis often relies on effective detection of the gear fault characteristic frequency and monitoring its magnitude change.Through conventional Fourier spectrum analysis, one has to identify the sideband spacing first, and then conduct fault diagnosis by matching it with the gear fault characteristic frequency.However, it is not easy to pinpoint the gear fault characteristic frequency via sideband analysis.Because of multiple modulations and their coupling effects, planetary gearbox vibration signals feature a highly complex spectral structure: sideband clusters around the carrier frequency (meshing frequency or its harmonics) and even sidebands within a sideband cluster.The multiple modulation sources may result in different sideband spacing values.The spacing between two adjacent sideband peaks may not exactly equal the gear fault characteristic frequency, instead may equal the sum or difference of the gear fault characteristic frequency and the frequency of time-varying vibration transfer paths.This leads to difficulty in planetary gearbox fault diagnosis through conventional Fourier spectrum analysis.
The modulating frequency of both the AM and the FM parts links to the gear fault characteristic frequency.If we can separate the AM and the FM parts, the Fourier spectrum of AM part or FM instantaneous frequency alone can directly reveal the gear fault characteristic frequency, thus avoiding the complex sidebands inherent in Fourier spectrum of the whole AM-FM signal.This inspires the idea of fault characteristic frequency identification through amplitude and frequency demodulation analysis.
However, demodulation analysis needs proper selection of a sensitive AM-FM component containing gear fault information.Most importantly, frequency demodulation analysis requires effective mono-component decomposition for accurate instantaneous frequency estimation.Recently, Dragomiretskiy and Zosso [21] proposed the variational mode decomposition (VMD).This method can non-recursively decompose a complex multi-component signal into constituent AM-FM components, and it is robust to noise.Wang and his colleagues [22][23][24][25] discovered the wavelet packet like frequency band decomposition property of VMD based on fractal Gaussian noise simulations, extended VMD from real-valued to complex-valued signals, and applied them to extract the fundamental, sub-harmonics, super-harmonics, and impacts of rotor-stator rubbing fault, as well as the impulsive components of rolling bearing fault.An and Zeng [26] utilized VMD to analyze the nonstationary pressure fluctuation signal of a hydraulic turbine draft tube.They showed that VMD is better than EMD in suppressing mode mixing and improving time-frequency readability.An and Yang [27] developed a denoising method based on VMD and approximate entropy, they demonstrated its better performance than wavelet transform, and applied it to hydropower unit vibration analysis.Tang et al. [28] proposed to decompose signals into multiple components, and then use independent component analysis to solve the underdetermined blind source separation problem.They detected the compound fault of rolling bearings with the proposed method.Lv et al. [29] extracted features from the VMD based time-frequency representation, and identified rolling bearing faults.Yi et al. [30] improved the robustness of VMD to sampling and noise via particle swarm optimization, and utilized the proposed method to detect rolling bearing faults.Liu et al. [31] determined the number of modes in VMD following a correlation criterion, and presented a time-frequency analysis method based on VMD and energy separation.They extracted the rotor-stator rub and oil-whirl fault feature using the proposed method.Mahgoun et al. [32] used VMD to detect defect impulses under variable speeds and loads, and verified the feasibility via analysis of gear transmission dynamics simulated signals.Jiang et al. [33] incorporated bi-dimensional VMD with discriminant diffusion maps to analyze the vibration signals of coal cutters, and detected the rolling bearing cracks.Yan et al. [34] applied genetic algorithm to optimize the balancing parameter and the number of modes, and combined this improved VMD with 1.5D envelope spectrum to extract compound fault signature of rolling bearings and gearboxes.Zhang et al. [35] decomposed cutting force signals into constituent components via VMD and wavelet packet decomposition, and detected milling chatter based on the energy entropy of constituent components.Zhu et al. [36] proposed an adaptive version of VMD, by optimizing the number of modes and data fidelity constraint with kurtosis as an optimization index via artificial fish swarm algorithm, and they applied the method to rolling bearing fault detection.These studies have demonstrated that VMD can effectively decompose complicated multi-component signals into AM-FM components, and suppress noise interferences, thus enhancing fault features manifested by modulations.Therefore, VMD provides an effective approach to demodulation analysis of planetary gearbox vibration signals.
In this paper, two major contributions are made.Firstly, considering the complex AM-FM features of planetary gearbox vibration signals, and exploiting the merits of VMD in AM-FM component decomposition of complex modulated signals, we propose a joint amplitude and frequency demodulation analysis idea for planetary gearbox fault diagnosis, thereby to overcome the difficulty with conventional Fourier spectrum analysis due to sideband complexity, and reveal the gear fault feature in a more effective way.Secondly, in order for effective implementation of the proposed demodulation analysis, we propose solutions to address the issues existing with application of VMD in planetary gearbox vibration signal demodulation analysis.These solutions can also be extended to applications of VMD in other signal analysis tasks.(1) The number of intrinsic mode functions (IMFs) to be separated is a key parameter for VMD.We propose a criterion to determine the number according to the spectral characteristics of planetary gearbox vibration signals.(2) Among those IMFs obtained from VMD, which IMF contains the major information about and is sensitive to gear fault is also an important issue.We propose a criterion to select the sensitive IMF for further demodulation analysis, based on the modulation characteristics of gear vibration signals and the wavelet packet like decomposition property of VMD.(3) Accurate instantaneous frequency estimation is necessary for effective frequency demodulation analysis, but conventional Hilbert transform (HT) based approach is subject to the constraint imposed by Bedrosian and Nuttall theorem [37].We preprocess the IMFs with empirical AM-FM decomposition by exploiting its capability to separate the carrier signal from amplitude envelope of an AM-FM signal, to satisfy the non-overlap requirement of amplitude envelope and carrier signal spectra by Bedrosian and Nuttall theorem.
Hereafter, this paper is organized as follows.Sections 2 and 3 introduce the principle of VMD and empirical AM-FM decomposition, respectively.Section 4 proposes a joint amplitude and frequency demodulation analysis method.Section 5 illustrates the method using a simulated planetary gearbox signal.Section 6 further validates the method by analyzing the lab experimental datasets of a planetary gearbox.Section 7 draws conclusions.

Variational Mode Decomposition
This section summarizes the principle of VMD based on [21].VMD is an entirely non-recursive decomposition method, and it extracts the constituent AM-FM components of a complex multi-component signal adaptively and concurrently.It defines IMFs as explicit AM-FM models, and relates the parameters of AM-FM models to the bandwidth of IMFs.According to the narrow-band property of IMFs, the AM-FM parameters can be found by minimizing the bandwidth, thus obtaining IMFs.This method has good merits over other available mode decomposition methods, such as theoretical rationale and robustness to noise and sampling.
In essence, IMFs are AM-FM signals where A k (t) is the amplitude envelope, φ k (t) is the instantaneous phase, and both the amplitude envelope A k (t) and the instantaneous frequency ω k (t) = dφ k (t) dt vary much slower than the instantaneous phase φ k (t).Such IMFs of AM-FM nature have a limited bandwidth.
VMD decomposes a signal x(t) into an ensemble of IMFs u k (t) that are band-limited about their respective center frequency ω k , while reconstructing the signal optimally.It iteratively updates each IMF u k (t) in the frequency domain, and then estimates the center frequency ω k as the center of gravity of the IMF power spectrum.
Motivated by the narrow-band properties of the AM-FM IMF definition, each IMF is assumed to be mostly compact around a center frequency ω k , i.e. it has specific sparsity properties.The sparsity prior of each IMF is described by its bandwidth.For each IMF u k (t), in order to assess its bandwidth, its analytic signal is firstly computed by means of Hilbert transform to obtain a spectrum of unilateral non-negative frequency.Then, its spectrum is shifted to baseband by multiplying with an exponential harmonic tuned to the respective center frequency.Finally, the bandwidth can be estimated through the squared l 2 norm of the gradient.The resulting constrained variational optimization problem is min , s.t.
where δ(•) is the Dirac delta function, the symbol * denotes the convolution operator, and K is the number of IMFs to be extracted.
To render the optimization problem, Equation (2), into an unconstrained form, a quadratic penalty term and a Lagrangian multiplier are introduced, for quick convergence and strict enforcement of the constraint.Then, the objective function to be minimized becomes an augmented Lagrangian where λ(t) is the Lagrange multiplier, α is the balancing parameter of the data-fidelity constraint, and •, • stands for inner product.
The solution to the minimization problem, Equation ( 2), can now be found as the saddle point of the augmented Lagrangian, Equation (3), in a sequence of iterative sub-optimizations [21].
Each IMF u k (t) can be updated as a solution to a minimization problem equivalent to Equation (3) In the frequency domain, the solution to Equation (4) can be found as [19] ûk (ω) = Then, the IMF in the time domain u k (t) can be obtained by inverse Fourier transforming Equation (5) and taking the real part.The center frequency ω k associated with each IMF u k (t) can also be updated as a solution to a minimization problem equivalent to Equation ( 3) It can be found as the center of gravity of the associated IMF's power spectrum [21] The complete algorithm of VMD is summarized as follows: t), n be 0, and predefine convergence threshold ε and the number of IMFs K to be separated.
2. Update each IMF u k (t) and its associated center frequency ω k , for k = 1 : K and all ω ≥ 0 ûn+1 where τ is the Lagrangian multiplier update parameter.4. Check the convergence condition , and terminate the decomposition.Otherwise, let n = n + 1, return to Step 2.
VMD decomposes a complex signal into a specific number of IMFs.These IMFs are AM-FM in nature, thus we can estimate their amplitude envelope and instantaneous frequency.More importantly, the non-recursive and concurrent decomposition nature of VMD effectively avoids the shortcoming of recursive decomposition algorithms, such as sensitivity to noise and sampling, over-shooting or under-shooting of upper and lower envelope via interpolation to extrema.

Empirical AM-FM Decomposition
To estimate instantaneous frequency accurately via commonly used Hilbert transform, it requires the signal to be mono-component and satisfy Bedrosian and Nuttall theorem, i.e., the spectra of carrier signal and amplitude envelope do not overlap in frequency domain [37].However, the IMFs obtained from VMD do not satisfy this condition automatically.In this case, amplitude envelope will contaminate the carrier signal, and instantaneous frequency estimation via Hilbert transform is subject to influences of amplitude envelope, leading to occasional negative frequency values.
To address this issue, we preprocess the IMFs with empirical AM-FM decomposition.For any AM-FM component, its instantaneous frequency is essentially contained in the carrier signal.Empirical AM-FM decomposition can separate an IMF of AM-FM nature into corresponding amplitude envelope and carrier signal, thus removing the effect of AM term on instantaneous frequency estimation, and fulfilling the requirement by Bedrosian and Nuttall theorem.
Empirical AM-FM decomposition is based on iterative applications of cubic spline fitting to the local maxima of the signal absolute value, and separates the AM and FM parts of any IMF signal uniquely but empirically through a normalization scheme [38]: 1. Initialization: Given an IMF u(t), set iteration index i = 1, residual signal r 0 (t) = u(t).
2. For the residual signal r i−1 (t), identify all the local maxima of its absolute value.
3. Fit to the local maxima using a cubic spline, obtaining the empirical envelope e i (t).
4. Normalize the residual signal r i−1 (t) using the empirical envelope e i (t), obtaining the normalized residual r i (t) = (t) .Using the absolute value data fitting, the normalized signal is guaranteed symmetric about the zero mean.
5. Check the stop criterion, if r i (t) ≤ 1 for all t, terminate the normalization, and designate the empirical FM part c(t) = r i (t) = cos[φ(t)], which is a pure FM function with unity amplitude, and the AM part (amplitude envelope)

go to Step 1 and repeat
Steps 2-5.
After empirical AM-FM decomposition, the IMF can be written as It separates any IMF empirically and uniquely into the corresponding amplitude envelope (AM) and the carrier signal (FM) parts.The obtained normalized carrier signal c(t) has unity amplitude, thus satisfying the requirement by Bedrosian and Nuttall theorem automatically.

Procedure of Demodulation Analysis
Planetary gearbox vibration signals have strong AM-FM features.Gear faults are manifested by the gear fault characteristic frequency that modulates the gear meshing vibration in terms of both amplitude and frequency.Amplitude and frequency demodulation analysis is effective to identify the modulating frequency of both the AM and the FM parts, thus providing an approach to extraction of gear fault features.
Amplitude demodulation needs properly selecting a constituent component sensitive to the gear fault, and frequency demodulation requires fine mono-component decomposition for accurate estimation of instantaneous frequency.However, planetary gearbox vibration signals are usually composed of complex multiple AM-FM components.It is difficult to reveal the modulation characteristics from raw signals, due to their multi-component nature.
VMD is effective in decomposing complex signals into constituent AM-FM components, and empirical AM-FM decomposition is capable to separate the carrier signal from any AM-FM component.A fusion of these two methods provides an effective solution to estimation of instantaneous frequency and amplitude envelope, and thereby offering insight into the modulation nature of planetary gearbox vibration signals.In our proposed method, we exploit the advantages of VMD and empirical AM-FM decomposition to jointly amplitude and frequency demodulate AM-FM signals for planetary gearbox fault diagnosis.The procedure is detailed below.
1. Decompose the signal into some IMFs via VMD.The number of IMFs to be separated K is a key parameter for VMD in real applications.We propose a criterion to determine the number according to the AM-FM nature and the spectral characteristics of planetary gearbox vibration signals.Planetary gearbox vibration signals can be modeled as a sum of some AM-FM signals considering their modulation nature (see Appendix A).Each of these AM-FM signals has a specific center frequency equal to one of the gear meshing frequency harmonics.Conversely, each gear meshing frequency harmonic corresponds to an AM-FM component (an IMF equivalently).Therefore, the gear meshing frequency harmonic orders covered in the effective frequency band of a signal, i.e. 0 through half the sampling frequency, roughly amounts to the number of IMF center frequencies.Accordingly, we can estimate the number of IMFs K = F s 2 f m , where • is the floor operator, F s is the sampling frequency, and f m the gear meshing frequency.
2. Separate the amplitude envelope a(t) and the carrier signal c(t) of each IMF via empirical AM-FM decomposition.For each IMF, calculate its instantaneous frequency based on the local phase derivative of carrier signal c(t) using Hilbert transform where H(•) stands for Hilbert transform.This instantaneous frequency estimation approach is called normalized Hilbert transform (NHT).
3. Choose the IMF with an instantaneous frequency fluctuating around the gear meshing frequency or its harmonics as a sensitive mono-component for further demodulation analysis, because gear fault information are mainly carried by the meshing vibration or its harmonics.If more than one IMF satisfies this condition, choose the IMF with a higher center frequency, since gear fault induced vibration signals feature high frequency impulses.
4. Fourier transform the amplitude envelope and the instantaneous frequency of the selected sensitive IMF, thus obtaining the amplitude and frequency demodulated spectra, i.e. envelope spectrum and Fourier spectrum of instantaneous frequency.5. Diagnose the gear fault by matching peaks in the demodulated spectra with the gear fault characteristic frequencies.
The envelope spectrum and the Fourier spectrum of instantaneous frequency are free from complex sidebands.Therefore, the proposed amplitude and frequency demodulation analysis is expected to be able to identify the modulating frequency related to gear fault characteristic frequency more effectively.

Numerical Simulation Evaluation
In this section, we illustrate the proposed demodulation analysis with a numerical simulated vibration signal of planetary gearboxes.Without loss of generality, we consider a signal synthetic of three AM-FM components centered at the first three gear meshing frequency harmonics respectively.According to the planetary gearbox fault vibration signal model in [39], the simulated signal is generated by where f t , f d and f m are the vibration transfer path varying frequency, the gear fault characteristic frequency and the gear meshing frequency respectively, and n(t) is a white Gaussian noise at a signal-to-noise-ratio of 0 dB mimicking the background interferences.Table 1 lists the parameter values.Figure 1 shows the signal waveform, its Fourier spectrum and envelope spectrum.The effective frequency band of this simulated signal is 0-650 Hz (F s /2 = 650 Hz), and covers up to the third harmonics of gear meshing frequency f m = 182 Hz.According to the proposed criterion for selecting the number of IMFs to be separated, we set the number K = 3 in the VMD.
Appl.Sci.2017, 7, 775 9 of 18 Figure 2 shows the three IMFs that are obtained by VMD and their instantaneous frequencies.For each IMF, its instantaneous frequency is calculated via both NHT based on empirical AM-FM decomposition and conventional Hilbert transform (HT) approaches.As shown in Figure 2b,d,f, negative values exist with HT results.This phenomenon does not comply with the physical meaning of frequency.On the contrary, NHT results are free from the presence of misleading negative values.This demonstrates that empirical AM-FM decomposition can improve instantaneous frequency estimation, thus being helpful for further frequency demodulation analysis.In the Fourier spectrum, (Figure 1b), the gear meshing frequency, its harmonics up to third order, and their associated sidebands dominate the signal.In the envelope spectrum, which is obtained by Fourier transforming the amplitude envelope of raw signal (Figure 1c), the gear fault frequency, transfer path frequency, and their sum and difference combinations show up in a pronounced magnitude.
Figure 2 shows the three IMFs that are obtained by VMD and their instantaneous frequencies.For each IMF, its instantaneous frequency is calculated via both NHT based on empirical AM-FM decomposition and conventional Hilbert transform (HT) approaches.As shown in Figure 2b,d,f, negative values exist with HT results.This phenomenon does not comply with the physical meaning of frequency.On the contrary, NHT results are free from the presence of misleading negative values.This demonstrates that empirical AM-FM decomposition can improve instantaneous frequency estimation, thus being helpful for further frequency demodulation analysis.For each IMF, its center frequency increases with the order of IMFs, and its instantaneous frequency fluctuates around the gear meshing frequency or its harmonic.These characteristics are consistent with both the numerical simulation settings and the findings from the Fourier spectrum.They also match well with the wavelet packet like property of VMD.Based on these findings, the three IMFs contain the major information about gear fault.
According to the proposed sensitive IMF selection criterion, we choose IMF1 as an informative one, and calculate its amplitude envelope and instantaneous frequency following Step 5 of empirical AM-FM decomposition in Section 3 and Step 2 of proposed analysis procedure in Section 4, respectively, for further amplitude and frequency demodulation analysis.Figure 3 shows its envelope spectrum and the Fourier spectrum of instantaneous frequency.The envelope spectrum (Figure 3a) shows a simpler structure than the Fourier spectrum (Figure 1b), since it is free of involute sidebands.Peaks appear at the time-varying transfer path frequency f t , the gear fault characteristic frequency f d , as well as their sum and difference combination f d ± f t .This feature is similar to the envelope spectrum of raw signal (Figure 1c), and shows that the proposed amplitude demodulation analysis of selected IMF is equally effective as conventional envelope spectrum of raw signal.The Fourier spectrum of instantaneous frequency (Figure 3b) exhibits a simpler structure than the envelope spectrum of both raw signal and selected IMF (Figures 1c and 3a).Only one pronounced peak emerges, exactly corresponding to the gear fault characteristic frequency f d .The proposed frequency demodulation analysis outperforms both Fourier and envelope spectra.It not only overcomes the intricate sideband difficulty inherent with Fourier spectrum, but also removes the irrelevant transfer path frequency existing in envelope spectrum.Most of all, these findings are consistent with the gear fault vibration symptoms expected in amplitude and frequency demodulated spectra [40], demonstrating the effectiveness of proposed amplitude and frequency demodulation analysis method.
respectively, for further amplitude and frequency demodulation analysis.Figure 3 shows its envelope spectrum and the Fourier spectrum of instantaneous frequency.The envelope spectrum (Figure 3a) shows a simpler structure than the Fourier spectrum (Figure 1b), since it is free of involute sidebands.Peaks appear at the time-varying transfer path frequency t f , the gear fault characteristic frequency d f , as well as their sum and difference combination  d t f f .This feature is similar to the envelope spectrum of raw signal (Figure 1c), and shows that the proposed amplitude demodulation analysis of selected IMF is equally effective as conventional envelope spectrum of raw signal.The Fourier spectrum of instantaneous frequency (Figure 3b) exhibits a simpler structure than the envelope spectrum of both raw signal and selected IMF (Figures 1c and 3a).Only one pronounced peak emerges, exactly corresponding to the gear fault characteristic frequency d f .The proposed frequency demodulation analysis outperforms both Fourier and envelope spectra.It not only overcomes the intricate sideband difficulty inherent with Fourier spectrum, but also removes the irrelevant transfer path frequency existing in envelope spectrum.Most of all, these findings are consistent with the gear fault vibration symptoms expected in amplitude and frequency demodulated spectra [40], demonstrating the effectiveness of proposed amplitude and frequency demodulation analysis method.

Experimental Settings
Figure 4 shows the planetary gearbox experimental system (University of Science and Technology Beijing, China), which consists of a drive motor, an encoder, a one-stage planetary gearbox and a magnetic powder brake.Table 2 lists the configuration parameters of the planetary gearbox.During experiments, the rotating frequency of the drive motor connecting to the sun gear shaft is set to a constant speed of 24.94 Hz.Given the planetary gearbox configuration parameters and the sun gear shaft rotating frequency, the characteristic frequencies are calculated, as listed in Table 3.A load of 18.75 Nm is applied to the planet carrier shaft by the magnetic powder brake.Accelerometers are fixed to the planetary gearbox casing to measure vibrations, and the signal is collected at a frequency of 20,480 Hz.In this section, the signal from the accelerometer at the top of the gearbox casing will be analyzed.To simulate a gear fault, localized chipping is introduced to one tooth of the sun, one planet and ring gear, respectively.Figure 5 shows the gear fault photos.

Experimental Settings
Figure 4 shows the planetary gearbox experimental system (University of Science and Technology Beijing, China), which consists of a drive motor, an encoder, a one-stage planetary gearbox and a magnetic powder brake.Table 2 lists the configuration parameters of the planetary gearbox.During experiments, the rotating frequency of the drive motor connecting to the sun gear shaft is set to a constant speed of 24.94 Hz.Given the planetary gearbox configuration parameters and the sun gear shaft rotating frequency, the characteristic frequencies are calculated, as listed in Table 3.A load of 18.75 Nm is applied to the planet carrier shaft by the magnetic powder brake.Accelerometers are fixed to the planetary gearbox casing to measure vibrations, and the signal is collected at a frequency of 20,480 Hz.In this section, the signal from the accelerometer at the top of the gearbox casing will be analyzed.To simulate a gear fault, localized chipping is introduced to one tooth of the sun, one planet and ring gear, respectively.Figure 5 shows the gear fault photos.

Signal Analysis
In order to alleviate the computational burden in VMD, we resample all the experimental signals at a new sampling frequency of 3200 Hz.The new effective frequency band (0-1600 Hz) covers the first five harmonics of the gear meshing frequency.The new resampled signal still contains the gear fault information, because gear fault information is mainly carried by the gear meshing vibration and its harmonics.According to the proposed criterion, we set the number of IMFs to be separated K = 5, in the following analyses.

Sun Gear Fault
Figure 6 shows the sun gear fault signal, its Fourier spectrum and envelope spectrum.Figure 7 displays the five IMFs output from VMD and their respective instantaneous frequency waveform.Among them, IMF3 has a center frequency of 1136 Hz corresponding to four times the gear meshing frequency, and its instantaneous frequency fluctuates around 1136 Hz, as shown in Figure 7f.According to the sensitive IMF selection criterion, we choose IMF3 for further amplitude and frequency demodulation analysis.
In order to alleviate the computational burden in VMD, we resample all the experimental signals at a new sampling frequency of 3200 Hz.The new effective frequency band (0-1600 Hz) covers the first five harmonics of the gear meshing frequency.The new resampled signal still contains the gear fault information, because gear fault information is mainly carried by the gear meshing vibration and its harmonics.According to the proposed criterion, we set the number of IMFs to be separated K = 5, in the following analyses.

Sun Gear Fault
Figure 6 shows the sun gear fault signal, its Fourier spectrum and envelope spectrum.Figure 7 displays the five IMFs output from VMD and their respective instantaneous frequency waveform.Among them, IMF3 has a center frequency of 1136 Hz corresponding to four times the gear meshing frequency, and its instantaneous frequency fluctuates around 1136 Hz, as shown in Figure 7f.According to the sensitive IMF selection criterion, we choose IMF3 for further amplitude and frequency demodulation analysis.In order to alleviate the computational burden in VMD, we resample all the experimental signals at a new sampling frequency of 3200 Hz.The new effective frequency band (0-1600 Hz) covers the first five harmonics of the gear meshing frequency.The new resampled signal still contains the gear fault information, because gear fault information is mainly carried by the gear meshing vibration and its harmonics.According to the proposed criterion, we set the number of IMFs to be separated K = 5, in the following analyses.

Sun Gear Fault
Figure 6 shows the sun gear fault signal, its Fourier spectrum and envelope spectrum.Figure 7 displays the five IMFs output from VMD and their respective instantaneous frequency waveform.Among them, IMF3 has a center frequency of 1136 Hz corresponding to four times the gear meshing frequency, and its instantaneous frequency fluctuates around 1136 Hz, as shown in Figure 7f.According to the sensitive IMF selection criterion, we choose IMF3 for further amplitude and frequency demodulation analysis.(e) (f) Figure 8 shows the envelope spectrum and the Fourier spectrum of instantaneous frequency of IMF3.The envelope spectrum (Figure 8a) presents a similar structure to that of raw signal (Figure 6c).The sun gear fault characteristic frequency s f appears with a prominent magnitude, indicating the sun gear fault.In addition, the fractional harmonics of sun gear fault characteristic frequency, s 1/3f and s 2/3f , also emerge with pronounced magnitudes, further confirming the occurrence of sun gear fault.For a real planetary gearbox, it is impossible for all the planet gears to be perfectly identical.
The inevitable difference between planet gears will generate three different impulse trains when the faulty sun gear meshes with the three planet gears, thus leading to presence of the fractional harmonics.Meanwhile, the planet carrier rotating frequency harmonics c (1~3)f , as well as their sum and difference combinations with the fractional harmonics of sun gear fault characteristic frequency ± s c /3 m f nf also have high magnitudes.This is because the sun gear fault will result in an uneven load distribution among planet gears, and this uneven load distribution will magnify the modulation effect of planet carrier rotation on gear meshing vibrations.The Fourier spectrum of instantaneous frequency (Figure 8b) shows a simpler structure than both the Fourier and envelope spectra (Figure 6b,c).Only four prominent peaks appear at the sun gear fault characteristic frequency s f , as well as its sum combination with the planet carrier rotating frequency +  Figure 8 shows the envelope spectrum and the Fourier spectrum of instantaneous frequency of IMF3.The envelope spectrum (Figure 8a) presents a similar structure to that of raw signal (Figure 6c).The sun gear fault characteristic frequency f s appears with a prominent magnitude, indicating the sun gear fault.In addition, the fractional harmonics of sun gear fault characteristic frequency, 1/3 f s and 2/3 f s , also emerge with pronounced magnitudes, further confirming the occurrence of sun gear fault.For a real planetary gearbox, it is impossible for all the planet gears to be perfectly identical.The inevitable difference between planet gears will generate three different impulse trains when the faulty sun gear meshes with the three planet gears, thus leading to presence of the fractional harmonics.Meanwhile, the planet carrier rotating frequency harmonics (1~3) f c , as well as their sum and difference combinations with the fractional harmonics of sun gear fault characteristic frequency m/3 f s ± n f c also have high magnitudes.This is because the sun gear fault will result in an uneven load distribution among planet gears, and this uneven load distribution will magnify the modulation effect of planet carrier rotation on gear meshing vibrations.
The Fourier spectrum of instantaneous frequency (Figure 8b) shows a simpler structure than both the Fourier and envelope spectra (Figure 6b,c).Only four prominent peaks appear at the sun gear fault characteristic frequency f s , as well as its sum combination with the planet carrier rotating frequency f s + f c , twice the sun gear rotating frequency 2 f s (r) and the third harmonic of planet carrier rotating frequency 3 f c .This implies existence of the sun gear fault, consistent with the actual experimental setting.
A comparison study again demonstrates the outperformance of proposed joint amplitude and frequency demodulation analysis method over conventional Fourier and envelope spectra.The proposed amplitude demodulation analysis of sensitive IMF avoids the difficulty due to intricate sidebands in Fourier spectrum, and is equally effective as the envelope spectrum of raw signal (see Figure 8a vs. Figure 6b,c).The proposed frequency demodulation analysis of sensitive IMF is free from intricate sidebands, further suppresses the irrelevant transfer path effect, and highlights gear fault features (see Figure 8b vs. Figure 6b,c).
Appl.Sci.2017, 7, 775 14 of 18 A comparison study again demonstrates the outperformance of proposed joint amplitude and frequency demodulation analysis method over conventional Fourier and envelope spectra.The proposed amplitude demodulation analysis of sensitive IMF avoids the difficulty due to intricate sidebands in Fourier spectrum, and is equally effective as the envelope spectrum of raw signal (see Figure 8a vs. Figures 6b,c).The proposed frequency demodulation analysis of sensitive IMF is free from intricate sidebands, further suppresses the irrelevant transfer path effect, and highlights gear fault features (see Figure 8b vs. Figures 6b,c).

Planet Gear Fault
Figure 9 shows the planet gear fault signal and its Fourier spectrum.For this case, we choose IMF5 for demodulation analysis based on the sensitive IMF selection criterion, because its center frequency 1422 Hz is the highest and its instantaneous frequency ripples around the fifth harmonic of gear meshing frequency.Figure 10 shows the amplitude and frequency demodulation analysis result.In the envelope spectrum (Figure 10a), the planet carrier rotating frequency and its second harmonic mf nf , are also prominent.These findings accord with theoretical expectations on the planet gear fault symptoms, and indicate the planet gear fault.The planet gear fault will result in an uneven load distribution among planet gears and thereby will magnify the AM effect of planet carrier rotation on gear meshing vibrations.Moreover, the AM effect of planet gear fault is modulated by the time-varying transfer path due to planet carrier rotation.These factors lead to presence of the planet gear fault characteristic frequency harmonics, the planet carrier rotating frequency harmonics, and their sum and difference combinations, in the envelope spectrum of planet gear fault case.
In the Fourier spectrum of instantaneous frequency (Figure 10b), peaks also appear at the planet gear fault characteristic frequency harmonics

Planet Gear Fault
Figure 9 shows the planet gear fault signal and its Fourier spectrum.For this case, we choose IMF5 for demodulation analysis based on the sensitive IMF selection criterion, because its center frequency 1422 Hz is the highest and its instantaneous frequency ripples around the fifth harmonic of gear meshing frequency.
Appl.Sci.2017, 7, 775 14 of 18 A comparison study again demonstrates the outperformance of proposed joint amplitude and frequency demodulation analysis method over conventional Fourier and envelope spectra.The proposed amplitude demodulation analysis of sensitive IMF avoids the difficulty due to intricate sidebands in Fourier spectrum, and is equally effective as the envelope spectrum of raw signal (see Figure 8a vs. Figures 6b,c).The proposed frequency demodulation analysis of sensitive IMF is free from intricate sidebands, further suppresses the irrelevant transfer path effect, and highlights gear fault features (see Figure 8b vs. Figures 6b,c).

Planet Gear Fault
Figure 9 shows the planet gear fault signal and its Fourier spectrum.For this case, we choose IMF5 for demodulation analysis based on the sensitive IMF selection criterion, because its center frequency 1422 Hz is the highest and its instantaneous frequency ripples around the fifth harmonic of gear meshing frequency.Figure 10 shows the amplitude and frequency demodulation analysis result.In the envelope spectrum (Figure 10a), the planet carrier rotating frequency and its second harmonic mf nf , are also prominent.These findings accord with theoretical expectations on the planet gear fault symptoms, and indicate the planet gear fault.The planet gear fault will result in an uneven load distribution among planet gears and thereby will magnify the AM effect of planet carrier rotation on gear meshing vibrations.Moreover, the AM effect of planet gear fault is modulated by the time-varying transfer path due to planet carrier rotation.These factors lead to presence of the planet gear fault characteristic frequency harmonics, the planet carrier rotating frequency harmonics, and their sum and difference combinations, in the envelope spectrum of planet gear fault case.
In the Fourier spectrum of instantaneous frequency (Figure 10b), peaks also appear at the planet gear fault characteristic frequency harmonics Figure 10 shows the amplitude and frequency demodulation analysis result.In the envelope spectrum (Figure 10a), the planet carrier rotating frequency and its second harmonic (1~2) f c are dominant, the planet gear fault characteristic frequency harmonics (1~3) f p , as well as their difference combinations with the planet carrier rotating frequency m f p − n f c , are also prominent.These findings accord with theoretical expectations on the planet gear fault symptoms, and indicate the planet gear fault.The planet gear fault will result in an uneven load distribution among planet gears and thereby will magnify the AM effect of planet carrier rotation on gear meshing vibrations.Moreover, the AM effect of planet gear fault is modulated by the time-varying transfer path due to planet carrier rotation.These factors lead to presence of the planet gear fault characteristic frequency harmonics, the planet carrier rotating frequency harmonics, and their sum and difference combinations, in the envelope spectrum of planet gear fault case.
In the Fourier spectrum of instantaneous frequency (Figure 10b), peaks also appear at the planet gear fault characteristic frequency harmonics m f p , the planet carrier rotating frequency harmonics n f c , and their difference combinations m f p ± n f c .Although it has a more complex structure than theoretical expectations, considering the features revealed in the envelope spectrum, a comprehensive analysis of amplitude and frequency demodulated spectra can still confirm the planet gear fault.Figure 12 displays the amplitude and frequency demodulation analysis result of IMF2.In the envelope spectrum (Figure 12a), the ring gear fault characteristic frequency and its second harmonic r (1~2)f are present with pronounced magnitudes.Meanwhile, prominent peaks show up at fractional harmonics of the ring gear fault characteristic frequency r (1~8)/3f .Similar to the sun gear fault case, it is the difference among the three planet gears that leads to three different impulse trains when the faulty ring gear tooth meshes with planet gears and therefore presence of the ring gear fault characteristic frequency fractional harmonics.In the Fourier spectrum of instantaneous frequency (Figure 12b), peaks also appear at the ring gear fault characteristic frequency harmonics r (1~3)f , and its fractional harmonics r (1~8)/3f .These findings confirm the ring gear fault.Figure 12 displays the amplitude and frequency demodulation analysis result of IMF2.In the envelope spectrum (Figure 12a), the ring gear fault characteristic frequency and its second harmonic r (1~2)f are present with pronounced magnitudes.Meanwhile, prominent peaks show up at fractional harmonics of the ring gear fault characteristic frequency r (1~8)/3f .Similar to the sun gear fault case, it is the difference among the three planet gears that leads to three different impulse trains when the faulty ring gear tooth meshes with planet gears and therefore presence of the ring gear fault characteristic frequency fractional harmonics.In the Fourier spectrum of instantaneous frequency (Figure 12b), peaks also appear at the ring gear fault characteristic frequency harmonics  Figure 12 displays the amplitude and frequency demodulation analysis result of IMF2.In the envelope spectrum (Figure 12a), the ring gear fault characteristic frequency and its second harmonic (1~2) f r are present with pronounced magnitudes.Meanwhile, prominent peaks show up at fractional harmonics of the ring gear fault characteristic frequency (1~8)/3 f r .Similar to the sun gear fault case, it is the difference among the three planet gears that leads to three different impulse trains when the faulty ring gear tooth meshes with planet gears and therefore presence of the ring gear fault characteristic frequency fractional harmonics.In the Fourier spectrum of instantaneous frequency (Figure 12b), peaks also appear at the ring gear fault characteristic frequency harmonics (1~3) f r , and its fractional harmonics (1~8)/3 f r .These findings confirm the ring gear fault.Figure 12 displays the amplitude and frequency demodulation analysis result of IMF2.In the envelope spectrum (Figure 12a), the ring gear fault characteristic frequency and its second harmonic r (1~2)f are present with pronounced magnitudes.Meanwhile, prominent peaks show up at fractional harmonics of the ring gear fault characteristic frequency r (1~8)/3f .Similar to the sun gear fault case, it is the difference among the three planet gears that leads to three different impulse trains when the faulty ring gear tooth meshes with planet gears and therefore presence of the ring gear fault characteristic frequency fractional harmonics.In the Fourier spectrum of instantaneous frequency (Figure 12b), peaks also appear at the ring gear fault characteristic frequency harmonics

Summary and Conclusions
Planetary gearbox vibration signals have complex and multiple AM-FM features, and the modulating frequency links to the gear fault.Thus, amplitude and frequency demodulation are expected to be effective in gear fault feature extraction.VMD is effective in decomposing complex modulated signals into their constituent AM-FM components, thus providing a potential approach to demodulation analysis.In this paper, considering the modulation feature of planetary gearbox vibration signals and the merits of VMD, we propose an amplitude and frequency demodulation analysis method for planetary gearbox fault diagnosis.In order to make VMD efficient in real complex signal analysis, the number of IMFs to be separated can be set as how many meshing frequency harmonic orders is covered in the effective frequency band from 0 to half sampling frequency, according to the spectral distribution characteristics around carrier frequencies of modulated signals.Furthermore, to enable accurate instantaneous frequency estimation via Hilbert transform, IMFs can be preprocessed via empirical AM-FM decomposition, thus avoiding contamination of AM on the carrier signal, particularly negative frequency values.To extract gear fault feature effectively, sensitive IMF can be selected as the one with a center frequency equal to meshing frequency or its harmonics, based on the wavelet packet like property of VMD and the modulation characteristics of gear vibrations.In the amplitude and frequency demodulated spectra, gear fault feature can be identified by matching present peaks with theoretical gear characteristic frequencies.We validate the proposed method through both numerical simulated and lab experimental signal analyses.Localized fault on all the three types of gears (sun, planet and ring) are successfully diagnosed.are the AM and FM functions, respectively; A kn > 0 and B kl > 0 are the magnitudes of AM and FM, respectively; c is a dimensionless constant depending on signal amplitude; θ, φ and ϕ are the initial phases of a signal component, AM and FM, respectively; f m is the gear meshing frequency; f d is the gear fault characteristic frequency (as listed in Table A1); and N and L are the highest order of AM and FM, respectively.In Equation (A2), when n = 0, A k0 ≡ 1 and φ k0 ≡ 0.
The function h(t) in Equation (A1) describes the effect of vibration transfer paths on the sensor signal when the gear fault induced vibration propagates from the meshing location of the faulty gear tooth to the sensor.Table A1 lists the function h(t) for both localized and distributed fault cases on each type of gear [39,40].
s cf f , twice the sun gear rotating frequency( ) r s 2 fand the third harmonic of planet carrier rotating frequency c 3 f .This implies existence of the sun gear fault, consistent with the actual experimental setting.

c(
1~2)f are dominant, the planet gear fault characteristic frequency harmonics p (1~3)f , as well as their difference combinations with the planet carrier rotating frequency  p c

p
mf , the planet carrier rotating frequency harmonics c nf , and their difference combinations  p cmf nf .Although it has a more complex structure than theoretical expectations, considering the features revealed in the envelope spectrum, a
planet gear fault characteristic frequency harmonics p (1~3)f , as well as their difference combinations with the planet carrier rotating frequency  p c

Figure 10 .Figure 11 .
Figure 10.Demodulation analysis of IMF5: (a) Envelope spectrum; and (b) Fourier spectrum of instantaneous frequency.6.2.3.Ring Gear FaultFigure11exhibits the ring gear fault signal and its Fourier spectrum.According to the sensitive IMF selection criterion, we choose IMF2 for demodulation analysis, because it has a center frequency of 284 Hz corresponding to the gear meshing frequency.

Figure 10 .Figure 11 .
Figure 10.Demodulation analysis of IMF5: (a) Envelope spectrum; and (b) Fourier spectrum of instantaneous frequency.6.2.3.Ring Gear FaultFigure11exhibits the ring gear fault signal and its Fourier spectrum.According to the sensitive IMF selection criterion, we choose IMF2 for demodulation analysis, because it has a center frequency of 284 Hz corresponding to the gear meshing frequency.

Table 1 .
Parameters in simulated signal.

Table 2 .
Configuration parameters of planetary gearbox.