Control of phases by ESPRIT and WLSE methods for the early detection of gear cracks

– The early detection of gear faults remains a major problem, especially when the gears are subjected to non stationary phenomena due to defects. In industrial applications, the crack of tooth is a default very diﬃcult to detect whether using the time descriptors or the frequency analysis. In this work and based on a numerical model, we prove that the crack default aﬀects directly the phase of the frequency component of the defective wheel (frequency modulation). To properly estimate the phases, we suggest two high-resolution techniques (Estimation of Signal Parameters via Rotational Invariance Techniques ESPRIT with a sliding window and Weighted Least Squares Estimator WLSE). The results of both methods are compared to the phase obtained by Hilbert transform. The three techniques are then applied on a multiplicative signal with a frequency modulation to show the inﬂuence of the amplitude modulation on the quality of phase estimation. We note that the ESPRIT method is much better in the estimation of frequencies while WLSE shows much eﬃciency in the estimation of phases if we keep the frequencies almost stables.


Introduction
In industrial applications, gear is considered among one of the most important elements as critical mechanical components in rotating machinery.Its improper installation, the overloading, the fluctuations of the load, the velocity and bad lubricants may lead to the generation of noise and damage to teeth.Several types of gear monitoring techniques have been developed to improve the precision of fault identification.Vibration analysis remains the most recognized tool for diagnosis of rotating machines [1].Among these techniques, time-series averaging and amplitude demodulation [2], Cepstral analysis [3], time-frequency analysis [4], wavelet, empirical mode decomposition (EMD) and local mode decomposition (LMD) filtering [5,6] have been suggested for structural health monitoring (SHM).These techniques give good indications for common faults (unbalance, eccentricity, misalignment, peeling. . . ) [7].However, they are not always suitable with regard to non-stationary phenomena [8], when frequency and amplitude modulations a Corresponding author: marc.thomas@etsmtl.caare combined.Crack fault remains one of the most serious defect and the more difficult to detect.In order to understand the dynamic behavior of a cracked structure and its vibratory signature, a numerical model for breathing crack of structures has been developed [9].
This paper is aimed to the improvement of signal processing methods to perform an early diagnosis of faults in mechanical systems using vibration signals.For processing the signals, the techniques of high-resolution ESPRIT [10] and WLSE [11] which allow for estimating the frequency components and their energies from vibratory signals are used in this paper.A sliding window is introduced into the ESPRIT method for monitoring the instantaneous phase variation due to cracks.The results are then compared with those obtained from WLSE and the instantaneous phase as calculated by the Hilbert transform.Since amplitude and frequency modulations appear in a breathing crack, a multiplicative signal with frequency-modulation is first used to analyze the influence of the amplitude modulation phenomenon on the quality of phase estimation.In Section 2, the numerical model [9] of a cracked structure is presented with the resultant vibratory signals.The formulation of the data model is then shown in Section 3. The high-resolution ESPRIT technique with the sliding window is discussed in Section 4. The WLSE technique for the estimation of frequencies and phases is discussed in Section 5.The estimation of phases by the Hilbert method is presented in Section 6.A comparison between the ESPRIT, WLSE and the Hilbert technique is then performed for estimating the phases of numerical models in Section 7. The results are finally discussed and conclusions are given in Section 8.

Tooth crack modeling
Literature shows some studies aimed to the diagnosis of cracks by exploiting the phase variations.The works of McFadden [12] propose to control the estimated phase of the meshing frequency by using the Fourier Transform.
To show the effects of the crack on the stiffness and therefore phases, we will rely in this study on the idea that the tooth crack presents a varying opening during the meshing [13] under the varying load.To understand and develop a model for the stiffness variation due to the crack, the tooth is considered as a beam with a breathing crack [9].The tooth sections are 2.88 × 20 mm with a length of 6 mm.The application of a harmonic force on the tooth generates oscillations with tension stresses and the opening of the crack varies.Consequently, the stiffness varies as a function of time, and hence, the modal parameters.This stiffness variation under a varying load has been validated by measurements and the results can be found in [14].From Equation (1), the variation in stiffness may then be expressed as: where -T is the period of gear meshing (7.5E-3 sec); k 0 represents the mean stiffness for acracked tooth without a tension stress (1.1E8 N.m −1 ); δk (1.25E7 N.m −1 ) is the amplitude of the stiffness variation due to the crack (11.4%).
The values of the cracked tooth stiffness were extracted from finite element simulations [15].By assuming a one degree of freedom system, the equation of motion gives a non stationary behavior after introducing the non-linear variation of stiffness, as follows: The variation of stiffness leads to a variation of the natural frequency and critical damping, and thus of the damping rate.Figure 1 shows the variation of the stiffness, the natural frequency and the damping rate in function of time.
These variations due to the presence of a crack have thus an influence on the modal properties of the tooth,  and on the vibratory responses.The non linearity generates harmonics of the excitation frequency.The proposed method consists in exciting harmonically the structure at a frequency equal to half of its natural frequency.The amplitude of the second harmonic is then amplified by the coincidence with the resonance (Fig. 2).This method results in amplitude and frequency modulations [9].
The crack generates a variation in the instantaneous phase and thus a variation of the natural frequency of the tooth.Figure 3 shows the phase shift of the acceleration signal between the cracked tooth without and with the tension stress.

Data model formulation
This section is aimed to define the problem for estimating the frequencies of d complex sinusoids S k (t) of the form: S k (t) = α k e j(w k t+ϕ k ) , where α k is the real amplitude of the sinwave, ω k is the unknown frequency, and ϕ k is the initial phase.We assume that N samples are available from a noisy measurement z(t) defined as: The vector b(t) presents a complex noise, zero mean and Gaussian random vector with where I m is a (m × m) identity matrix.As it intends to apply a subspace approach [16] in the vibration domain whereas only one sensor is used, we define a data vector y(t) viewed as a windowed partition of the whole data set acquired, where y (t) = (z (t) z(t + 1) . . .z(t + m − 1)) T , with m the window length.
The signal may hence be represented accordingly with the so-called matrix form: where

. , b(t + L)]
are independent whereas the columns are mutually correlated (with respect to the window length).We also point out that the matrix A could be re-written in a centro-symmetric manner without loss of generality.The covariance matrix of the observed data is usually estimated by means of time averaging leading to R (t) = y (k) y H (k) when the complex random vectors are reputed proper.We suggest in this section a classical rotational invariance technique applied to the model (4) supposed to have a better performance for the estimation of frequencies and phases.

ESPRIT approach
The ESPRIT approach is based on the exploitation of the subspace signal of the autocorrelation matrix R(t).This could be obtained by an Eigen or Singular value decomposition of R(t).When using this technique, we have: In the same way, Σb contains the (2m − d) remaining eigenvalues and the columns of Êb are the associated eigenvectors.A rough estimation of the noise variance is given by σ2 = 1 2m−d Trace Σb .In this position, we can now rewrite the system (5) as follows: Ês↓ Φ = Ês↑ (7) The Φ matrix can be estimated by the Least Squares approach and the d frequencies are next achieved by taking the argument of the eigenvalues of Φ.On other hand, the matrix Ŝ can be estimated as follows: To well estimate the instantaneous phases by this technique, we propose to use a sliding window with a step of Δt along the signal (Fig. 4).For each window, we firstly estimate the desired frequency and then its phase by: where δ k = α k e jϕ k represent the elements of the first column of the matrix Ŝ given on Equation ( 8).This will allow for controlling the slightest variation in the phase.

WLSE approach
In this section, we present the estimation of the instantaneous phases by the Weighted Least Square Estimator (WLSE).A frequency estimation algorithm is also proposed to track the frequencies when the frequency is varying.We assumed that the processed signal is given by: (10) where Ē is the amplitude, ω is the angular frequency, ∅ is the initial phase, E d = Ē cos(∅), and E q = Ē sin(∅).Note that E d , E q and ∅ are constant under the steady state.
Expressing (10) in the matrix form, we obtain: where x The cost function is chosen such that: where λ ∈ (0, 1) is the forgetting factor.The solution x (t i ) that minimizes the cost function J [x (t i )] is obtained by following the least-squares algorithm such that [17]: 16) where x (t 0 ) = 0, P (t 0 ) = π 0 I ∈ R 2×2 and π 0 > 0 is the initial covariance constant.From Equations ( 13) and ( 16), we obtain the instantaneous phase such that where atan2 is the arctangent function being capable of distinguishing 45 • from 225 • .The proposed phase angle estimation algorithm can be extended to estimate the frequency ω, when it is varying.
If ω = ω then ∅ is no more constant, but varies such that ∅ ( From Equation ( 21), one can recognize that if Δ ∅ = 0, then there is a frequency estimation error.The basic idea for updating ω is to employ a PI regulator such that Δ ∅ is regulated to a zero value.Then, we claim that the output of the PI regulator can be used for ω.Note that ω is utilized for Equation ( 12) along with the WLSE algorithm.To avoid a possible instability of the WLSE algorithm, the PI gain of this frequency estimation algorithm must not be large.Note that ∅0 is an initial phase angle.The proof of the convergence of the frequency estimator was presented in [18].Additionally, to avoid a possible erroneous behavior caused by a large amount of phase angle variation when the covariance is reset, ω is kept a constant during N ω τ time after the covariance resetting.According to simulation study, N ω ∈ [2, 20] is recommended.In the numerical applications, we will use the data model developed in Section 3, where the matrix H(t) will be replaced by the Vandermande matrix A and the E s (t) matrix by y(t) and the estimated matrix x(t) will wear the desired phases of our signal that they can be separated using the arctang function.

The Hilbert transform
Hilbert transforms are essential in understanding envelop methods.Let x(t) denote a real signal and x(t) its Hilbert transform.The analytic signal of x (t) as proposed by [19] is given by: The amplitude modulation (AM), phase modulation (PM) and frequency modulation (FM) are obtained by: where denotes the time derivative and tan −1 (.) is the inverse of the tangent function, which gives the phase values in the range [−π, +π].
In our study, we are only interested in the formula (24) for the estimation of the instantaneous phases.

Numerical analysis
In this section, we present some theoretical numerical simulations to illustrate the performance of the techniques WLSE and ESPRIT with a sliding window for the estimation of phases and detecting cracks.

Phase analysis of a multiplicative signal with amplitude and frequency modulation
In order to investigate the efficiency of methods (ES-PRIT, WLSE and Hilbert) when amplitude and frequency modulations are both present, a multiplicative signal which has the same spectral shape than gear signals, has been generated.Consider a multiplicative signal as follows: where τ e , τ r1 and τ r2 represent the meshing period and the rotational periods of the two wheels, respectively.G e (t), G r1 (t) and G r2 (t) are harmonic signals [3] representing the meshing signal and the modulation caused by the two wheels, respectively.∅ 1 (t) is a square wave which has the same pulsation than the defective pinion G r1 (t) aimed for representing the phase variation due to the coincidence with the defected tooth at each revolution.This phase ∅ 1 (t) has been assumed as integrated into G e (t), G r1 (t) and G r2 (t) because these three components are affected by the defect tooth of G r1 (Eq.( 26)).
The variation of the phase depends on β = T 0 /T r1 , where T 0 is the opening time of the crack and T r1 presents the period of the pinion.Figure 5a shows the instantaneous phase of the multiplicative signal, estimated by Hilbert where the amplitudes of G r1 and G r2 are filtered around the gear frequency (f e = 450 Hz) in order to remove the frequencies generated by the amplitude modulations.The peaks (with low amplitudes) represent the sudden variation of the phase when there is coincidence with the defected tooth at each revolution (1/F 1 ).However the duration of the squared wave ∅ 1 (t) is difficult to distinguish.Figure 5b shows the frequency spectrum of (Fig. 5a) where we can distinguish the frequency of the defective gear (F 1 = 30 Hz) and its harmonics.These results show that the Hilbert technique gives a good estimation of the variation of the phase, which is corresponding to the period of the pinion, when the signal is not modulated in amplitude.
Figure 5c shows the instantaneous phase (with high amplitudes) when there is an amplitude modulation.When the signal is modulated in amplitude, it is more difficult to identify the square variation at each period (1/F 1 ) corresponding to the revolution, because it is difficult to distinguish the frequencies coming either from amplitude or frequency modulations.Figure 5d presents the FFT of the phase of the signal modulated in amplitude.In this spectrum, the coincidence frequency (6 Hz = f e /PPCM = 450/(3 × 5 × 5)) and its harmonics appear very clearly because ∅ 1 (t) affects both gears, but the frequency F 1 and its harmonics have low amplitudes.The frequency of coincidence (6 Hz) can also be noticed in Figure 5b when there is no amplitude modulation, but with very low amplitude.Figure 6a presents the results obtained by ESPRIT with a sliding window when amplitude modulation is present.The squared wave can be identified at each revolution, but not its duration.Figure 6b shows the variation of the phase obtained by WLSE where it is obvious that the duration of the squared phase variation injected into the Equation (26) (red curve) is the same as that the phase variation estimated by WLSE (blue curve).In this case, the main frequencies have been first computed by FFT.
Despite the existence of the amplitude modulation in this case, ESPRIT and WLSE clearly illustrate the variation of the instantaneous phase because in these two algorithms, the estimate of the phase of each frequency component is independent of all other frequencies components.The disadvantage of WLSE relies on the fact that when the estimation of frequencies is required, the phases must be stable that is not the case, here.Otherwise it will create small fluctuations in the estimated frequencies.Figure 7 shows the instantaneous frequency variation when computed from WLSE instead of FFT which leads to the disruption of phases estimated even with the use of a PI regulator as shown in Figure 8. Consequently this technique (WLSE) would not be useful for estimating the frequencies if the phases change at a high rate.In this case, it is preferable to apply the FFT for estimating the frequencies before to use WLSE in order to estimate the phase as made in Figure 6b.

Phase analysis of a cracked tooth model
By applying the principle to the crack model, the tooth has been harmonically excited to half of its natural frequency (1325 Hz).The second harmonic which corresponds to the natural frequency is thus subjected to  amplitude modulation and the variation of natural frequencies between 2650 to 2494 Hz produces frequency modulations, as described in Figure 2.
Figure 9 presents the behavior of phases obtained through Hilbert transform (Eq.( 24)) for the case of a cracked tooth without tension stress (Fig. 9a) and its FFT to see the amplitudes (Fig. 9b).It can be noticed the frequencies of (1325, 2650 Hz) corresponding to the excitation and the resonance respectively.When a tension stress is applied, the FFT of the defect tooth (Fig. 9d) shows several harmonics to the excitation frequency, with higher amplitudes and the phase is harmonically   modulated (Fig. 9c).The comb function in frequency is more characteristic of periodic shocks than modulation.
Figures 10 and 11 show the phases as computed by ESPRIT with a sliding window and its FFT respectively for the cracked tooth without stress.Figure 10 shows the phase with constant amplitude and Figure 11 only exhibits the mean natural frequency (2650 Hz).For the case of cracked tooth with stress, Figure 12 shows the phase with a large variation in amplitude and Figure 13 shows the amplitude modulation around the resonance frequency due to the opened crack.Figure 11 exhibits  only the natural frequency when there is crack without tension stress while the natural frequency is modulated by the excitation frequency (Fig. 13) when a tension stress is present.In this last case, the amplitude at the natural frequency decreases with the crack depth while it increases at its modulation frequency as it is generally the case when investigating a phase modulation.
Figures 14 and 15 show that WLSE provides similar results as obtained by ESPRIT with sliding window for   the cracked tooth without stress, except the 4th harmonic of the excitation but with small amplitude.
Regarding the defective tooth, Figure 16 shows the variation in amplitude of the phase estimated by WLSE and Figure 17 reveals a spectrum similar to those of Hilbert Transform (Fig. 9d) with many harmonics of the excitation frequency (1325 Hz) with large amplitudes.
The modulation frequencies are not present.This may be due to the choice of the parameter λ, since when increasing λ close to 1, there will be several harmonics and vice versa.

Conclusion
The objective of this study is to investigate the performance of Hilbert, ESPRIT and WLSE techniques for detecting cracks in gear teeth.It is proved that cracks affect the stiffness of the tooth, which produces non-stationarity in the mechanical system.Since the non-linearity generates harmonics of the excitation frequency, a natural amplification is obtained by exciting the tooth by a harmonic force at half of its natural frequency.This gives amplitude and frequency modulations at the second harmonic.By generating an analytical signal with and without amplitude modulation, the application of Hilbert transform for the estimation of the instantaneous phase has proven reliability only if the signal was not modulated in amplitude because it is difficult to distinguish the frequencies coming either from amplitude or frequency modulations.WLSE and ESPRIT methods revealed in this case to be efficient even if the signal is modulated in amplitude and frequency because in these two algorithms, the estimate of the phase of a frequency component is independent of all other frequencies.However, it is well known that the transmission error in gears produces amplitude modulation.Consequently, it is not possible to practically separate these two modulations.The advantage of the WLSE method is that it well estimates the duration of the phase variations.However, large phase variations degrade the quality of the frequency estimation which disrupts the estimation of phases.To get a good estimate of the phases by WLSE, we must previously use the FFT method to define the frequencies of our system and then kept them constant during the estimation of phases.However, the choice of parameters (IP and λ) influences and degrades the performance of the estimation of the signal components.ESPRIT technique with the sliding window gives a good result for the frequency estimation.ESPRIT has been found able to distinguish the variation of the instantaneous phases, due to the variation in opening of the crack, despite the presence of the amplitude modulation.The amplitude modulation can be easily identified from the FFT of the phase estimated by ESPRIT in the case of cracked tooth.These three methods should be considered in Structural Health Monitoring (SHM) of gears.However, ESPRIT with sliding window remains the best performing among these three techniques for identifying the signal frequencies and distinguishing the modulation frequencies due to a defect.

Fig. 1 .
Fig. 1.Variations of the dynamic properties of a cracked tooth.

Fig. 3 .
Fig. 3. Responses of a cracked tooth with and without a tension stress.
) with D = diag(e jω1 , . . ., e jω d ) and for any matrix T , the subscript .↓ means the elimination of the last row of T and the subscript .↑ means the elimination of the first row of T .To improve the calculation time of this algorithm, Haard[10] suggests a basic selection matrix based on Kronecker product for eliminating the appropriate rows.Let R = Ês Σs ÊH s + Êb Σb ÊH b (6) be an Eigen-decomposition of R where Σs is a diagonal matrix which carries the d largest eigenvalues, and columns (2m × d) of Ês correspond to their eigenvectors.

Fig. 4 .
Fig. 4. Sliding window to estimate the instantaneous phase of the signal.

Fig. 5 .
Fig. 5. Phase estimation by Hilbert for signals: (a) without amplitude modulation.(b) FFT of the phase of (a).(c) With amplitude modulation.(d) FFT of the phase of (c).

Fig. 6 .
Fig. 6.(a) Phase estimation by ESPRIT with sliding window.(b) Phase estimation by WLSE, of a signal modulated in amplitude and frequency.

Fig. 9 .
Fig. 9. Instantaneous phases by Hilbert Transform: (a) and (b) cracked tooth without tension stress and its FFT.(c) and (d) cracked tooth with tension stress and its FFT.

Fig. 13 .
Fig. 13.FFT of the phase estimated by ESPRIT (cracked tooth with stress).

Fig. 17 .
Fig. 17.FFT of the phase estimated by WLSE (cracked tooth with tension stress).