Tacholess Rotational Speed Estimation via the Integration of Continuous Wavelet Transform and Time-dependent Autoregression

Speed estimation is crucial to monitor the conditions of rotational machinery. Most speed measurements are carried out by installing encoders or tachometers inside the machines. In many cases, such method could be cumbersome or even inaccessible. This paper proposes a vibration-based speed estimation method. The vibration sensors are often cheaper and easier to install than angle encoders. In the proposed method, the continuous wavelet transform (CWT) is used as a preprocessing technique to extract the signal of importance. Then, the time-varying autoregressive (TAR) model is applied to analyze the rotational frequency. Additionally, the paper presents a fast algorithm for implementation. The proposed method is validated by both synthetic and empirical data.


Introduction
The rotational speed estimation of machinery is essential to condition monitoring because of its significant impacts on emission, fuel consumption, noise, dynamics and malfunction diagnosis [1] [2] [3] [4] [5]. In many cases, the direct measurement of rotational speed could be extremely difficult and costly [6]. This is because the encoders or tachometer needs to be positioned on or close to the rotational components and the working environment is harsh [7]. Therefore, the estimation of rotating speed from vibration signal has been received tremendous attentions.
The frequency contents of vibration signals of rotation machines are very complex [1] [8] [9] [10]. Due to the low signal to noise ratio (SNR) and various undesired frequency components, the accurate estimation of rotational speed from vibration is still a great challenge to existing applications. Further, many rotating machines operate at fluctuated speed which produces non-stationary signals [11]. Hence, the speed estimating requires robust and adaptive techniques to deal with the aforementioned challenges.
Numerous achievements have been accomplished in this area of vibration analysis. K.C. Gryllias and I.A. Antoniadis [12] present a complex shifted Morlet wavelets based instantaneous speed estimation algorithm. Its advantage is the free selection of the wavelet's center frequency and bandwidth, making the algorithm adaptive to signal. Nonetheless, a limited number of center frequency selected could be too coarse to identify the speed variation. An automatic speed estimation for gear pair is proposed in [13], combining the singular spectrum decomposition and phase coupling analysis. This method depends solely on the harmonics of the gear-mesh frequency and lacks of noise suppression. H. Lin and K. Ding [14] provides a four stroke engine speed measurement. Its principle is that the half harmonics of the torque of combustion engine can be explicitly recognized via the spectrum autocorrelation. One problem with the method is the highest peak results from the autocorrelation may not always be the half harmonic order.
Parallel, the joint use of filter method and parametric modelling is an alternative way for speed estimation. Its main advantages are the adaptability of filterbank setting and high estimating resolution comes from the modelling [15]. Hence, such concept has great potentials and should be open studied. However, quite few literatures with respect to the concept are found in measuring the rotational speed. One research can be cited is [16]. The paper proposes an instantaneous frequencies estimation based on the combination of complex Morlet wavelet and ESPRIT. By setting appropriate wavelet filterbank, the ESPRIT is able to recover the rotating frequency without interference of undesired components and noise. The restrictions of using this technique are the non-strictly analytic properties of the complex Morlet wavelet and the intensive computation of subspace algorithm [17] [18].
In this paper, we propose a fusion algorithm combing the CWT and TAR to estimate the rotational speed. The CWT works as bandpass filters to suppress the noise and undesired frequency component while retaining the signal of importance. Thus, the TAR built by the wavelet coefficients is adaptable. Additionally, a fast implementation of the algorithm is presented to reduce the computational cost. Compared with [16], the main novelties are the applicability of a wide range of wavelet families and the efficiency in computation. The estimating performance is validated via both synthetic and empirical data.

Methods
The principle of the proposed algorithm is using the wavelet coefficients to build the Modified Yule-Walker (MYW) equation. The wavelet transform is essentially a bandpass filter of uniform shape and varying location and width [19]. By appropriating setting the CWT filters, the model becomes adaptive to the signal. Once the MYW equation is created, we solve it by total least square (TLS) to obtain the model's parameters. Finally, the rotating frequency is estimated via Cadzow's method.
It is worth noting that the TAR is the moving window version of classical autoregressive (AR) model and is detailed in [15]. For simplicity, we assume the described data to be windowed. Thus, the TAR is reduced to AR in this section.
The proposed algorithm is mainly divided into two parts: Autocovariance sequence (ACS) of CWT coefficients (fast implementation); II.
Solve the model via the ACS based MYW equation and estimate the frequency by Cadzow method.

ACS of CWT coefficients
The proposed algorithm uses the ACS to build the MYW equation. In this subsection, the CWT is briefly described first; then, the ACS of wavelet coefficients is introduced.
The CWT of signal f at scale s and time t can be expressed as: where * represents the complex conjugate and  is the mother wavelet. 1 s is the amplitude normalization factor. By modifying the scale parameter, the wavelet is able to construct desired bandpass filters [20].
Summing over a subset of the scales in (1), the filtered signal is then obtained: Such filter has a frequency response given by the sum of the wavelet functions between scales i1 and i2, discriminating frequency components outside its bandwidth.
Next, the ACS of the wavelet coefficients is calculated by: where  .
E denotes the expectation operator and the sequence is assumed to be zero mean. Once the ACS is obtained, the AR model can be solved from the MYW equation: where α is the AR parameter and p represents the AR order. The model is adaptive to the analyzed signal because the ( ) r k  is calculated from the filtered signal.

Fast Implementation
The AR model built directly from (3) is computationally expensive. Compared with the convolution operation in time domain, it is considerable faster to do the calculation in frequency domain. The paper provides a Fourier transform (FT) based algorithm to calculate the ( ) r k  in order to efficiently create the MYW equation. Its principle is that the ( ) r k  is essentially a wavelet transform of the ACS of the raw signal [21]. This can be calculated in Fourier space; the wavelet used is a new wavelet resulting from the autocovariance of the original wavelet. The fast algorithm is illustrated in Fig. 1.
Let  denotes the covariance operator and  denotes the convolution operator respectively. The fast implementation procedures are as follows: Calculate the new wavelet in frequency domain. The new wavelet, termed  , is obtained through the autocovariance of the original wavelet. A simple way of the calculation derives from the fact that many wavelets have symmetric or antisymmetric structure. By combining the correlation property ( ) , the new wavelet is then: Equation (5) and (6) leads to: where and are the Fourier transform of  and respectively. Thus, the new wavelet is obtained.
Taking the m th order derivatives of Gaussian wavelets (DOG) as an example. Its expression in frequency domain is [22]: The new wavelets of the autocovariance of the DOG is then given by: Calculate the autocovariance of raw signal in frequency domain. The autocovariance can also be calculated via FT. For the zero-mean sequence, the autocovariance and correlogram are a Fourier transform pair: where ( ) r k is the autocovariance of the sequence;f is the Fourier transform of the sequence; and  represents the Fourier transform operator.
Calculate the ACS of wavelet coefficients. By multiplying equation (10)

Frequency Estimation from the MYW equation
Many effective algorithms have been developed to solve the AR model [18] [23]. One of the most commonly used techniques is the MYW method; it is often complemented with the TLS. After solving the model, the rotation frequency is then derived from the Cadzow's method. The frequency estimation consists of the following steps: Step Step 2: AR parameters estimation via TLS. The AR parameters p a are then approximated by TLS. This is done by performing the singular value decomposition (SVD) on (12), resulting in a diagonal matrix (  ) of non-zero diagonal elements (singular values): where U: l l  , : ( 1) Next, partition theV corresponding to the p largest singular value: 11 12 where the dimension is 12ˆ: ( 1 ) (16) where ( ) ( ) k k     and can be obtained from the covariance of the signal: Once the PSD is obtained, the rotating frequency is derived from the expectation of the PSD at every time interval.

Results and Discussion
In this section, the proposed algorithm is applied to both synthetic and empirical signals. Further, two other frequency estimation methods are presented for comparison purpose. One is the CWT which estimate the frequency by taking the expectation of the scalogram. In the analysis, the CWT shares the same wavelet parameters with the proposed algorithm. Another method is first conditional spectral moment method, termed FCSM. It is a common application that estimates the instantaneous frequency as the first conditional spectral moment of the time-frequency distribution of the signal [25].

Synthetic Signal
The synthetic signal consists a non-linear chirp with "instantaneous" frequency that varies sinusoidally between 120Hz to 280Hz. Also, a white Gaussian noise is embedded into the signal with a significant SNR 10dB. The sampling frequency is set to be 5000Hz. The signal is depicted in Fig. 2.
The wavelet is the generalized Morse wavelets with The wavelet filterbank is shown in Fig. 3 that has -3dB bandwidth between 100Hz to 300Hz. The rectangular window has 40 samples and steps forward every 20 samples. Fig. 2 The synthetic signal. a. signal waveform. b. signal spectrum.      The PSD and estimated results are presented in Fig. 4 and Fig. 5 respectively. Also, the estimation of CWT and FCSM are depicted in Fig. 6 and Fig. 7. The reference frequency is marked in red dashed line. Obviously, the proposed method demonstrates high estimating accuracy compared with the CWT and FCSM. The use of six wavelets in CWT lacks of enough 'scale resolution', leading to the drastically large variance. The FCSM presents large fluctuations in comparison to the proposed method because of its non-adaptability, being corrupted by the noise. The mean and median absolute errors are shown in Fig. 8.
The wavelet filter makes the proposed algorithm adaptive to the signal, suppressing undesired frequency components out side its bandwidth. Further, the AR model offers a high frequency resolution that enables the high accuracy in estimation.

Aircraft Engine Data
In this case, the data is the aircraft engine vibration provided by [27]. The signal represents a slow acceleration progress from idle to full power as displayed in Fig. 9. An accelerometer and a tachometer are installed to record the vibration signal and the angular speed respectively. The vibration signal is initially sampled at 50K Hz and then down sampled to 31.2K Hz.
The generalized Morse wavelet ( 3   and 20   ) is used to set the filters. Owing to the engine characteristics, the wavelets are set to cover the frequency band 175Hz to 260Hz. Hence, the undesired components and noise outside the band are suppressed. The wavelet filters are presented in Fig. 10. The window length is 4000 samples with 50% overlap. Fig. 9 The vibration signal. a. signal waveform. b. signal spectrum.     The estimated result is depicted in Fig. 11. During the idle stage, the estimated frequency remains at 183.5Hz with a mean error of 2.8Hz. This is actually an expected result since the tachometer is not mounted on the main shaft where the accelerometer is closed to. It is also seen that the sudden increase in rotational speed around 45s strains the estimated accuracy due to the smearing of harmonics. Due to the increased SNR when the engine reaches to its full power, the estimated frequency becomes a bit choppy but neglectable. Overall, the proposed method effectively recovers the rotating speed from vibration signal with high accuracy. Fig. 12 illustrates the CWT estimation in which the engine speed can hardly be recognized. This is mainly due to the four wavelet scales are too coarse to identify the frequency content. The FCSM is presented in Fig. 13. In this case, a bandpass filter is applied first with passband range [175Hz 260Hz]. This is because the raw estimation of FCSM is significantly corrupted by the strong background noise, making the comparison meaningless. Its estimating accuracy is inferior to the proposed method. Two main distortions are pointed out in Fig. 13. The mean and median absolute error of the three methods are presented in Fig. 14.

Conclusion
The paper presents a vibration signal based rotational speed estimation method that integrates CWT and AR model. By combining the advantages of adaptive filter and parametric modeling, the speed frequency can be recovered with good accuracy. It should be noting that the direct measurement approaches (e.g. tachometer) is still superior in accuracy.
Nonetheless, the installations of such devices required to be on the shaft of interest which can be unfeasible. On the other hand, the vibration sensors are often cheaper and easier to instrument than tachometers. Hence, the proposed method is an alternative way to effectively measure the rotational speed.