Extraction method of composite fault features of gear transmission system based on demodulation of NMD and Teager energy operators

In order to effectively identify and extract the composite fault characteristics of the gear transmission system, a composite fault diagnosis method combining nonlinear mode decomposition (NMD) and Teager energy operator demodulation is proposed. Because the envelope demodulation of Hilbert transform has the disadvantages of large amount of calculation and end effect, it uses Teager energy operator to solve the problem of large amount of calculation, and NMD solves the problem that the fault signal features is not easy to extract under the mode aliasing. First, the NMD method is used to decompose the fault simulation signal, and the nonlinear modal component with practical physical significance is obtained. Secondly, the Teager energy operator demodulation is carried out for the nonlinear modal components, and the demodulation results are analyzed to verify the feasibility of the method. Then, the method is applied to the composite fault diagnosis of gear pitting wear in gear transmission system, and the characteristic frequency obtained from the test data is compared with the calculated characteristic frequency. The comparative analysis shows that the method can separate the fault characteristic frequencies of large and small gears. The comparative analysis with EMD and EEMD methods in simulation signal analysis and experimental research shows that this method is superior.


Introduction
Gearbox is an integral part of most mechanical products. In order to prevent accidental failure and possible life loss, the fault diagnosis of gearbox has received great attention. 1,2 The vibration signal of gearbox is usually complex, non-stationary and nonlinear. 3 When a composite fault occurs, multiple faults will interfere with each other and be coupled with each other. 4 In recent years, time-frequency analysis methods have been widely used in composite fault diagnosis. Wang et al. 5 proposed a composite fault diagnosis method for Gearbox Based on improved integrated local mean decomposition. Through the analysis of simulation and experimental signals, it can be seen that the method effectively overcomes the mode mixing phenomenon and reduces the reconstruction error. At the same time, the method avoids the appearance of pseudo components and reduces the calculation amount. It is verified that the improved integrated local mean decomposition method is an effective method to extract composite fault features. Isham et al. 6 proposed the VMDEA and ELM combined fan gearbox fault diagnosis method. This study provides a good solution for determining optimized VMD parameters to decompose vibration signals, and also provides a more effective fault diagnosis method for wind turbine gearboxes. Li et al. 7 proposed a gearbox composite fault diagnosis method based on SWD-AVDIF. By analyzing the simulation signals and experimental signals of the gearbox 1 School of Automobile and Transportation, Shenyang Ligong University, Shenyang, China 2 composite fault, the results showed that the method could effectively separate the composite fault signals and strengthen the fault characteristics. Yang et al. 8 put forward a gear fault diagnosis method combining EEMD and Choi Williams. Through the analysis of simulation and experimental signals, the adaptability and effectiveness of this method in gearbox fault diagnosis were verified. The traditional time-domain or frequency-domain analysis method is only suitable for linear stationary signals. The vibration signals of composite faults are mostly non-linear and non-stationary signals. Using only time-domain or frequency-domain analysis methods often fails to effectively extract fault features. 9 Therefore, a time-frequency analysis method should be used in feature extraction of composite faults. Windowed Fourier transform (WFT) and wavelet transform (WT) are more commonly used timefrequency analysis methods in signal processing. Both are suitable for processing nonlinear, non-stationary signals. 10,11 The difference between them is that they analyze the form of the components that appear in the signal. WFT has linear frequency resolution, that is, it is distinguished according to the frequency difference of each component. WT has logarithmic frequency resolution and is distinguished according to the component frequency ratio. According to Heisenberg's uncertainty principle, for any kind of time-frequency transformation, the frequency resolution and time resolution cannot reach the optimum at the same time. The time resolution of the WT will increase as the frequency increases. The time resolution of the WFT is the same at all frequencies. So we merge the two, Using WT to get high temporal resolution, Use WFT for high frequency resolution. Huang first proposed a new timefrequency analysis method-empirical mode decomposition (EMD). It can decompose complex signals into finite intrinsic mode functions (IMF), which are suitable for nonlinear and non-stationary signals. 12 Since this method came out, it is widely used in fault diagnosis of mechanical equipment. 13,14 Wu and Huang 15 also proposed an EEMD method, adding white noise to the signal and then performing EMD decomposition. However, the signals decomposed by EMD and EEMD are prone to produce pseudo components. The non-sinusoidal complex components generated by the same excitation source will be decomposed into multiple simple components, which will increase the number of components. In addition, they are sensitive to noise and their robustness is not ideal. In view of their disadvantages, Iatsenko et al. 16 proposed a new timefrequency analysis method-nonlinear mode decomposition (NMD), which can decompose the given signal into oscillation with practical physical significance and decompose it into nonlinear mode component (NM). In the process of decomposition, the advantages of WT and WFT are integrated, the accuracy of extraction and reconstruction of main signals is improved, and the noise robustness is very strong. At present, the envelope demodulation method used in many domestic and foreign literatures of fault diagnosis is Hilbert transform envelope demodulation. 17,18 But there are also endpoint effects in Hilbert transform, which will affect the whole signal analysis results. 19 Hilbert has a lot of calculations, the calculation of the Teager energy operator is relatively small. The calculation amount of Teager energy operator demodulation is mainly 4N times real number multiplication, equivalent to N times complex multiplication; The calculation amount of Teager energy operator demodulation is about 1/logN of Hilbert demodulation method. 20 When N = 1024, it is 1/10. Teager energy operator can enhance the transient characteristics of the signal and is suitable for detecting the impact component in the signal. This method has high time resolution, and it has good adaptive ability to the instantaneous change of signal. It has low computational complexity and high efficiency. 21 In this paper, a composite fault feature extraction method combining nonlinear mode decomposition (NMD) and Teager energy operator is proposed. The composite fault signal is decomposed by NMD, and then the obtained nonlinear modal component (NM) is demodulated by Teager energy operator, and the fault type is determined by comparing the obtained fault characteristic frequency with the calculated fault characteristic frequency.
The second part of this paper introduces the basic theory of the proposed method. The third part analyzes the simulation signal. The fourth part compares and analyzes the effectiveness of the three methods for composite fault feature extraction. The fifth part expounds the conclusion of this paper.

Basic theory of NMD
Under a certain excitation, nonlinear system will produce response, which can be regarded as NM component. It is assumed that the signal consists of a series of NM components and noise.
Among them, h t ð Þ is noise, and c i (t) is NM component. c i (t) can be expressed as a series of FM and AM harmonics: Where v u t ð Þ ½ period is 2 p. r h (t) is the h-th harmonic of c i (t). r 1 (t) is the main harmonic. The rest are subharmonics (h = 2, 3, 4 .).
The purpose of NMD decomposition is to decompose the NM signal. The characteristic parameters of each harmonic r h (t) are obtained. Amplitude is A h (t), phase is u h (t), frequency is v h t ð Þ = u 0 h (t). To construct WT and WFT: (1) Wavelet transform the original signal. Given the signal s(t), WT can be expressed as: In the formula, v is the frequency, m and e are the other forms of time t in the time and frequency domains respectively; v c is the peak frequency of wavelet; c t is a normal wavelet basis function;ĉ and c Ã are the Fourier transform and complex conjugate function of c(t) respectively.
(2) Extraction time-frequency ridge v p A curve consisting of a series of ridge points is called a ridge curve. Represents the peak frequency of each moment. In the paper, the ridge extraction method proposed by Iatsenko et al. 22 is adopted. The ridge r 1 (t) of the main harmonic v p (t) is a sequence of peaks that maximizes the function path. For wavelet transform (WT), the path function can be expressed as: Where, G s (v, t) is the formula represented by equation (4). v p (t n ) is the expression of ridge v p when t = t n , t n is the n-th peak; \ f(t) . and std f(t) ½ represent the average and standard deviation of f(t) respectively.

(3) Main harmonic reconstruction
After the ridgeline curve is extracted, the main harmonic r 1 (t) needs to be reconstructed. There are two methods to reconstruct the main harmonic: the direct method and the ridge method. Because the ridge method has good robustness. Therefore, we first use ridge method to reconstruct the main harmonic. For wavelet transform, the ridge method can be expressed as: (4) WT or WFT selection criteria The paper combines the advantages of WT and WFT. It enables NMD to adaptively select the most suitable type of time-frequency transformation. The empirical formula can be expressed as: When R \ 1, WFT is used, and WT is used otherwise.

(5) Determination of reconstruction method
Since there are two reconstruction methods, in order to make the reconstructed harmonics more accurate and reasonable. According to different time-frequency analysis methods and two reconstruction methods, the final reconstruction method is determined.
(1) When WFT is the optimal time-frequency transform, the given signal G(s) is transformed by WFT: Where, g(t) is the Gaussian window function. The timefrequency ridge can be extracted according to equation (5). For WFT, according to the difference between the two, it is only necessary to change the logarithmic frequency scale ln v p in equation (5) to v p . Then the para- , respectively, using the direct method and ridge line method refactoring (d represent direct method, r represent direct method): Where, d is the discrete correction coefficient. First, calculate the weighted averageÃ (d, r) ,ũ (d, r) , v (d, r) . The reconstruction method is determined according to the inconsistency of the reconstruction method. For harmonic parameters, the inconsistency can be expressed as: Where (2) According to formula (11), if e d A \ e r A , the direct method is used, otherwise the ridge method is used. According to the above, WT is set as the best, the direct method is used to reconstruct the harmonic parameters.
(6) Extraction of sub-harmonics After obtaining the main harmonic r 1 (t), its instantaneous frequency is v 1 ð Þ (t). For the sub-harmonic r h (t), its ridge line v (h) p is the closest to the peak sequence of hv 1 (t) in the direction of amplitude increase. After finding the sub-harmonic ridge line v (h) p , the direct method and the ridge method are used to reconstruct the subharmonics r h (t). Then, the optimal sub-harmonic parameter is determined by equation (11).

(7) Identification of true harmonics
The amplitude, phase, and frequency of true harmonics should be related to the main harmonic. After extracting the primary harmonic r 1 (t) and the secondary harmonic r h (t), we need to determine if they are true harmonics. By using the method of time shift instead of data, the independent zero hypothesis test of the main harmonic r 1 (t) and the sub-harmonic r h (t) is carried out, so as to identify the authenticity of the harmonic. 23 The alternative parameter of the main harmonic r 1 (t) can be expressed as: Where, N is the sampling length, f s is the sampling frequency, d is the number of time shift points of the substitute data and M is the maximum time shift of the substitute data.
Alternative parameter process of subharmonic r h (t): (1) Move its signal forward in the time domain by (2) Use step (6) to select the optimal substitute har- Introduce the statistical feature q_(A,v)^((h)) to measure the degree of consistency between the primary and secondary harmonics. Its dependence can be expressed as: Consistency coefficient can be expressed as: The parameter v (h) A, u, v is the weight of each consistency coefficient q (h) A, u, v . r (h) 0 is the consistency coefficient when DT 0 = 0.
(4) In order to reduce the error caused by harmonic interference, the previous true harmonic must be subtracted and then harmonic identification (it is also applicable to the main harmonic. When identifying the sub-harmonic r h (t), the main harmonic r 1 (t) should be subtracted).
The extracted true harmonics are added to obtain the NM component, which is then subtracted from the original signal. Until there is only noise in the remaining signal. The surrogate test method with Fourier transform is used to detect whether the residual signal contains noise. The statistical characteristic quantity D of the alternative data is expressed as: Where, Q½Â(e) and Q½v(e) are the discrete spectral entropy of amplitude and frequency respectively. The inspection steps are as follows: (1) Calculate the statistical characteristic quantity D 0 (a A , a v ). (2) Create N s Fourier transform alternatives, and calculate the corresponding D s = 1, 2, ..., N s (a A , a v ). (3) If the confidence level of D s . D 0 is greater than or equal to the preset confidence level, the noise test zero hypothesis is rejected. D(1,1), D(0,1), and D (1,0) are used as discriminant statistics to test three different ones (a A , a v ). If at least one of them negates the zero hypothesis, it can be considered that the signal is not noise, and the test should be continued. The NMD flowchart is as follows (Figure 1):

Basic theory of Teager energy operator
The Teager energy operator is a non-linear operator, and the energy operator method is developed on the signal form of equation (17): An energy operator is defined as follows: _ x and € x are the first reciprocal and second derivative of x(t), respectively. There are: In general, the variation of the modulated signal is much slower than that of the carrier signal, so the amplitude-frequency value can be treated approximately as a constant. Thus there are: In the same way: The instantaneous amplitude and phase of the modulated signal are: Therefore, this method can be used to demodulate and analyze the signal. In this paper, NMD method and Teager energy operator demodulation method are combined to perform NMD decomposition for composite fault signals, and then Teager energy operator demodulation for decomposed signals. At the same time, it is compared with EMD and EEMD. The specific fault diagnosis process is shown in Figure 2.

Simulation analysis
This paper compares NMD with EMD and EEMD to verify the advantages of NMD. In the paper, the sum of two Am and FM signals and one random signal is selected for simulation analysis. The simulated signals and components are as follows: Among them, randn is random noise of N [0,1], the sampling time is 1 s, and the sampling frequency is 1024 Hz. The simulation signal and each component is shown in the Figure 3: Next, the NMD, EMD, and EEMD decompositions are performed on the simulation signals. Their decomposition results are shown in Figures 4 to 6 respectively.
It can be seen from the decomposition results that NMD only decomposes the simulated signal into two NM components. There are no redundant other components. The decomposed NM1 and NM2 correspond to the simulation signals X1 and X2 respectively, and are basically consistent. However, after EMD decomposition, we can actually get more components, six of which are IMF1-IMF6, among which IMF4-IMF6 is basically meaningless. After EEMD decomposition, more components can be obtained, the first six are IMF1-IMF6, among which IMF4-IMF6 is basically a pseudo component, which has no practical significance. Therefore, from the perspective of the decomposition effect of simulation signals, the decomposition effect of NMD is indeed better than that of EMD and EEMD. Next, the simulated components obtained from NMD, EMD and EEMD decomposition are demodulated by Teager energy operator respectively, obtain the NM simulation component demodulation spectrum (as shown in Figure 7), the IMF simulation component demodulation spectrum after EMD decomposition (as shown in Figure 8) and the IMF simulation component demodulation spectrum after EEMD decomposition (as shown in Figure 9). From the demodulation results of the NM simulation component and the demodulation results of the IMF simulation component decomposed by EMD and EEMD, it can be seen that the NM simulation component can be extracted to 1 and 2 times the characteristic frequency of the simulated signal fault, respectively. The IMF1 simulation component of EMD decomposition can be extracted to 1 frequency and 2 frequency, but the IMF2 simulation component of EMD decomposition can only be extracted to 1 frequency. The IMF1 simulation component of EEMD decomposition can also extract 1 and 2 times of frequency, but the IMF2 simulation component of EEMD decomposition can only extract 1 times of frequency. Although the demodulation results of EMD and EEMD simulation signal components can also get fault characteristic frequency, the number of peak points is less, and the phenomenon of mode aliasing occurs in the demodulation spectrum of EMD and EEMD simulation components. Therefore, the demodulation results show that the decomposition effect of NMD is indeed better than EMD and EEMD, and it can effectively suppress the phenomenon of modal aliasing, the decomposition and demodulation results of the simulated signal illustrate the advantages of the NMD method. Therefore, this paper uses the decomposition method of NMD.
Experimental analysis of gearbox composite fault diagnosis NMD, EMD, and EEMD methods are compared and applied to the complex fault diagnosis of gear-gear transmission system respectively, so as to verify the  effectiveness of this method in practical application. The equipment model used in the experiment is HFXC-I. The experimental equipment is shown in Figure 10. By manual setting, it is set up as a composite failure of pitting of big gear and wear of small gear. The acceleration sensor is used in the experiment, and the experiment is carried out under the load of 0.2 A, the sampling time is 2 S, the sampling frequency is 5120 Hz, and the measured speed is n = 825 r/min. The number of teeth of the small gear and the large gear is Z 1 = 55 and Z 2 = 75, The pinion shaft is the driving shaft, By calculation, the fault characteristic frequency of the large gear is f 0 =10.08 Hz, The fault characteristic frequency of the pinion is f r =13.75 Hz. The calculation formula of fault characteristic frequency is shown in equations (26) and (27).
The original composite fault signals were decomposed by EMD, EEMD and NMD methods, and the first two signals obtained by the decomposition of EMD and EEMD were taken, namely, IMF1 and IMF2, The results are shown in Figures 11 and 12, NMD decomposition results in two NM components as shown in Figure 13. Then the signals IMF1 and IMF2 obtained from the decomposition of EMD and EEMD were demodulated by Teager energy operator respectively and the demodulation results of the EMD decomposed signal is shown in Figure 14, and the demodulation results of the EEMD decomposed signal is shown in Figure 15. Then the two NM components obtained by NMD decomposition are demodulated by Teager energy operator, and the demodulation results is shown in Figure 16.
As shown in Figure 14(a), the demodulation results of the IMF1 component after EMD decomposition only extracted the 1-fold, 2-fold, and 3-fold frequency of the characteristic frequency of the large gear. As shown in Figure 15(a), the demodulation spectrum of IMF1 component decomposed by EEMD extracts the 1-fold frequency, 2-fold frequency, and 3-fold frequency of large gear fault feature frequency, and the 2fold frequency of small gear fault feature frequency. As shown in Figure 16  can extract the 1-fold, 2-fold, and 3-fold frequency of the fault feature frequency of large gear. Although EEMD can also extract the 1-fold, 2-fold, and 3-fold frequency of the fault feature frequency of large gear, it also extracts the 2-fold frequency of the fault feature frequency of small gear, which shows that EEMD cannot extract the fault feature frequency of large gear alone. As shown in Figure 14(b), the demodulation results of imf2 component after EMD decomposition extract the 1-fold and 2-fold frequency of pinion fault feature frequency, and extract the 1-fold, 2-fold, and 3fold frequency of large gear fault feature frequency. As shown in Figure 15(b), the demodulation spectrum of IMF2 component after EEMD decomposition only extracted the 1-fold, 2-fold, and 3-fold the fault characteristic frequency of pinion. As shown in Figure 16(b), the NM2 demodulation spectrum only extracted the 1fold, 2-fold, and 3-fold frequency of the characteristic frequency of the pinion. The demodulation results show that the EEMD and NMD methods can extract the fault characteristic frequency of the pinion. The effect is better than the EMD method. Therefore, through the comparative analysis of the demodulation results, we can know that the EMD method can extract the fault feature frequency of the large gear separately, but it can't extract the fault feature frequency of the small gear separately. The EEMD method can extract the fault characteristic frequency of  the small gear, but it cannot extract the fault characteristic frequency of the large gear. It shows that the EMD and EEMD methods cannot effectively extract the fault characteristic frequency of the large and small gears separately. NMD method can not only extract the fault feature frequency of small gear, but also the fault feature frequency of big gear. The superiority of NMD method is illustrated. The NMD method can separate the fault characteristic frequency of the gears separately. It shows that the NMD method is better than the EMD and EEMD methods. Therefore, NMD method is suitable for composite fault diagnosis of gear transmission system, and the effect is good.

Conclusion
The structure of gear transmission system is complex, and it is easy to be affected by noise and other external factors, so it is not easy to extract fault feature frequency and identify fault type. Therefore, this paper proposes a method of Composite Fault Feature Extraction Based on the combination of nonlinear mode decomposition (NMD) and Teager energy operator.
1. Through comparative analysis with EMD and EEMD methods, it can be seen that the composite fault feature extraction method combining the nonlinear mode decomposition (NMD) and  Teager energy operator proposed in this paper effectively suppresses the modal aliasing phenomenon. 2. Through the analysis of simulation signals and experimental data, we can see that the method proposed in this paper can effectively extract the composite fault characteristics and accurately identify the fault type. The validity of the proposed method is verified.
The proposed method has strong noise robustness and can effectively eliminate modal aliasing. This paper can effectively identify and extract the composite fault features of gear pitting wear. However, the correctness and accuracy of composite fault diagnosis of other mechanical equipment need to be further explored.