## 1 Introduction  
Gearbox fault detection is significant for the security and reliability of mechanical equipment [1]. Faults caused by gear failures account for 10% of the malfunctions in rotating machines and 80% of those in mechanical transmission systems [2]. Therefore, it is necessary to effectively detect faults in a gearbox in order to prevent accidents and guarantee its efficient operation [3]. Model-based methods in dynamic systems have been continuously proposed for gear fault diagnosis [4; 5; 6]. Nowadays, vibration-based analysis approaches are widely used to identify tooth faults due to their effectiveness and ease of measurement [2]. However, the collected vibration signals of a multi-stage gearbox have strong modulation signatures [7; 8] due to the existence of meshing frequencies, multi-harmonics, and coupled frequencies generated by modulation. Thus, in order to extract important signatures from raw vibration signals, various signal processing techniques have been utilized in gear defect detection, such as wavelet packet transform (WPT) [9] and empirical mode decomposition (EMD) [10]. Nevertheless, WPT alone is not suitable for fault diagnosis in a multi-stage gearbox because many close or identical frequency components occur in the signals [9]. In addition, EMD is an adaptive signal decomposition approach and has been extensively used to decompose a multi-component signal into several sub-signals, but it still has a major drawback of mode mixing [11]. Recently, variational mode decomposition (VMD) has been pioneered by [PERSON] and [PERSON] [12; 13], which is a non-recursive signal processing technique with a firm theoretical foundation compared with EMD. Meanwhile, the performances of VMD on non-stationary signal analysis have been also evaluated and compared with other adaptive signal processing methods. For instance, the equivalent filter bank property of the VMD has been explored, and the results demonstrated that VMD can be considered as a WPT or a generalized short time Fourier transform (STFT) with a varying window width [14]. Moreover, it has been demonstrated that VMD can better extract the bearing fault signatures of acoustic signals in comparison with EMD [15]. A hybrid defect detecting technique has been also proposed for denoising and non-stationary signature extraction, which joins VMD and the majoriation-minimization-based total variation denoising (TV-MM) technique [16]. Because of these merits of VMD, it has been well-utilized for detecting rubbing faults in a rotor system [13]. However, VMD still has some limitations regarding signal demodulation analysis. For example, VMD cannot directly identify the numerical value of demodulated signatures. Moreover, results using the above-mentioned signal decomposition methods alone are usually unsatisfactory for detecting tooth defects, especially in a multistage gearbox. Demodulation analysis is also a popular technique for detecting defects in rotating machines [17]. Some new works on demodulation analysis for fault diagnosis in bearings and gears have been developed; for example, a modulation signal bispectrum (MSB) and the modulation signal bispectrum-based sideband estimator (MSB-SE) method [18; 19]. One advantage of some demodulation techniques is that the numerical modulated information can be directly detected in the time domain [20] or via power spectral density in the spectral domain [21]. Furthermore, complex demodulated signals can be investigated using joint time-frequency or time-scale analysis techniques [22], but the meaning of the time-frequency plot is still hard to identify for non-experts. Consequently, the cyclostationary method seems to be a promising technique to identify periodic modulation components [23], which is very effective when the analyzed signal contains various components for complex machinery [24]. Modulation intensity distribution (MID) is regarded as a generalization of the spectral correlation density and can be applied to demodulation analysis of a cyclostationary signal [17].