# An Instantaneous Engine Speed Estimation Method Using Multiple Matching Synchrosqueezing Transform

Authors: Youyong Liu; He Wen; Zongying Ding; Lei Xu; Hongjiang Chen; Janusz Smulko
Journal: Journal of Sensors
 (2021)

---

## Abstract

Instantaneous rotational speed measurement of the engine is crucial in routine inspection and maintenance of an automobile engine. Since the contact measurement of rotational speed is not always available, the vibration measurement has been used for noncontact rotational speed estimation methods. Unfortunately, the accuracy of the noncontact estimation methods by analyzing engine vibration frequency is not satisfactory due to the influence of noise and the strong nonstationary characteristic of the vibration signal. To overcome these problems, based on the multiple matching synchrosqueezing transform (MSST) (MMSST, improved MSST with multiple squeeze operations), a novel noncontact method is proposed to accurately estimate the instantaneous rotational speed of automobile engine in this paper. Firstly, a MMSST is proposed to process the vibration signal to obtain a concentrated time-frequency (TF) representation. Secondly, the instantaneous frequency (IF) detection algorithm is employed to extract the fundamental frequency from the TF result. Finally, the rotational speed of the engine is calculated according to the relationship between the fundamental frequency and rotational speed. Results from numerical simulations and test on real engine have proven that the proposed method can obtain much higher frequency resolution and more precise IF estimation of the engine vibration signal and more accurate rotational speed estimation result compared with the MSST method. Furthermore, the proposed method is verified to have a stronger noise robustness and can provide satisfactory estimation results for engine vibration signal containing nonlinear frequency-modulated components.


---

## Body

## 1. Introduction

The engine is an essential part of an automobile, and its quality can comprehensively reflect the performance of the automobile’s dynamic. The rotational speed of the engine is a key parameter in the automobile safety inspection and emission test. For example, when the automobile is idling, the unstable engine speed indicates that the engine needs to be repaired to reduce exhaust emissions. Therefore, it is important to measure the engine speed stably and accurately. At present, the measurement methods of the engine rotational speed can be divided into two types (i.e., the contact-type and noncontact-type). The contact-type measurement methods need to install a toothed or notched encoder disc which rotates synchronously, at the end of the crankshaft of the automobile [1–4]. As the crankshaft rotates, the tooth or notch in the encoder disc triggers the sensor to generate a series of pulses, and the rotational speed is obtained by calculating the number of pulses between two consecutive pulses, or within a fixed time [5, 6]. However, these methods are time consuming and have poor portability due to the encoder being not always accessible, and the extremely high sampling frequency is required. Thus, the contact-type measurement methods are not suitable for the routine inspection of the automobile and other occasions that need to obtain the engine rotational speed of vehicles quickly and accurately.In general, a normally functioning engine generates a regular vibration which is closely related to the rotational speed [7]. The noncontact-type measurement methods convert the problem of estimating the rotational speed into that of estimating the frequency of a particular spectral component, which means noncontact methods can be realized by analyzing engine vibration frequency; it overcomes the aforesaid disadvantages. Lin and Ding [8] and Shan et al. [9] use discrete spectrum correction theory [10] to improve the precision of fast Fourier transform (FFT) analysis and get a satisfactory estimation result of rotational speed. Still, the accuracy and calculation burden depend on the characteristics of the chosen window function and the number of points involved in correction.Recently, time-frequency (TF) analysis (TFA) methods have attracted considerable attention in instantaneous frequency (IF) estimation for instantaneous rotational speed estimation [11–13], as they can provide a powerful tool for the analysis of time series signals and produce more reliable and more robust results [14]. Because the IF of nonstationary signal is the ridge point in the TF representation, the speed-related frequency can be estimated by detecting the ridge of the TF representation of vibration signals [15–23]. Urbanek et al. [24] utilize phase demodulation and short-time Fourier transform (STFT) to obtain the IF that is related to the speed. STFT and continuous wavelet transform [25, 26] (CWT) are two of the most commonly employed TFA techniques that map the TF distribution of a nonstationary signal by its inner product with a TF atom. However, the TF resolution is limited by the Heisenberg uncertainty principle, i.e., both time and frequency resolutions cannot become arbitrary small at the same time. The Wigner-Ville distribution [27] (WVD) is an effective approach to solve the contradiction between time localization and frequency resolution. However, WVD can generate nonnegligible cross-terms in analyzing multicomponent signals. To address these issues, a new trend in the development of TFA techniques is to characterize the nonlinear features of nonstationary signal by using postprocessing procedures, such as reassignment method [28] (RM) and synchrosqueezing transform [29, 30] (SST). The RM technique reassigns the diffused TF energy into the IF trajectory in both time and frequency directions to obtain a sharper TF result; the drawback of RM is that the TF result cannot be used for signal reconstruction.In contrast, since SST reassigns the TF coefficients in the frequency direction, it retains the ability to reconstruct signals from the TF result. It has been proved that the SST result is equivalent to the ideal TF representation when addressing a purely harmonic signal [31]. However, when dealing with time-varying signals, e.g., chirp signals [32] or nonlinear frequency-modulated signals, SST cannot generate a concentrated TF result. When the automobile is in the accelerating stage, the IF of the vibration signal of engine changes rapidly and has strong nonlinearity [33]. Therefore, SST is not suitable for analyzing the vibration signal of engine. Wang et al. [34] proposed a matching SST (MSST) approach to analyze the time-varying signal. It constructs a matching IF estimator according to the law of frequency-modulated signal and reassigns the TF coefficients in the frequency direction. For chirp signals or nonstationary signals that can be locally approximated as chirp signals, MSST can obtain a concentrated TF result. However, if the limitation described above is not met, i.e., nonlinear frequency-modulated signals, the TF result obtained through MSST will still suffer from the diffused TF energy problem. Multisynchrosqueezing chirplet transform [15] (MSSCT) is proposed to overcome this problem. It shows its effectiveness in generating a highly energy-concentrated TF result. However, a parameter that well matches the real chirp rate of analyzed signal needs to be determined in advance, which will influence the performance of MSSCT approach. The parameter estimation [35] with high precision is very difficult, especially for addressing nonlinear frequency-modulated engine vibration signal.In this paper, a multiple matching synchrosqueezing transform (MMSST) is proposed to obtain a more concentrated TF result and improve the IF estimation accuracy of nonlinear strong frequency-modulated components. By introducing multiple squeeze operations in MSST, MMSST iterates on the IF estimator of MSST to obtain a more accurate IF estimator which is close to the signal true IF and generates the concentrated TF result by reassigning the TF coefficients in the frequency direction. In order to achieve noncontact rotational speed measurement of engine, an accurate noncontact engine rotational speed estimation method that consists of a proposed MMSST and the IF detection algorithm [31] is proposed in this paper. The IF detection algorithm is employed to extract the fundamental frequency which is related to the rotational speed. Compared with the contact measurement methods, the proposed method does not require special device such as encoder, and the only hardware needed is a one-axis acceleration sensor used for acquiring engine block vibration signals which is easy to install, and a commonly signal pickup assembly. In addition, due to the proposed MMSST can provide a sharper TF result than MSST and do not need to predetermine a complex parameter like MSSCT in processing nonlinear strong frequency-modulated signals, the proposed method can accurately and quickly estimate the rotational speed when the engine is working in a steady state or under conditions of large fluctuations in speed. The two advantages make the proposed method suitable for routine inspection of the automobile and other occasions that need to obtain the engine rotational speed of vehicles quickly and accurately.


## 2. Method

### 2.1. Estimation Principle

When an engine is running, its power is transmitted and transformed through the crank and connecting rod mechanism. The excitation loads acting on the crankshaft of an internal combustion engine are mainly derived from the driving torque due to gas forces generated within the cylinders. In a single cylinder, the driving torque produced by the gas pressure can be expressed as [8, 36]
(1)Mg=M0+∑v∞Mvsinvωt+φv,where M0 represents the average torque generated by the gas pressure, v denotes harmonic order, Mv and φv represent the initial phase and amplitude of the vth harmonic order of the torque waveform, respectively, and ω is crankshaft angular velocity. Since each cylinder fires once every two engine revolutions in a four-stroke engine, the rate of driving torque change is half of the rotational angular velocity of the crankshaft, there is v=1/2,1,3/2,⋯, and 1/2 is the fundamental harmonic order.For a four-stroke engine withτ cylinders, according the knowledge of engine dynamics, when harmonic order v=τ/2,τ,3τ/2,⋯, the driving torques of each cylinder are in the same direction, and these orders are called main harmonic orders of the engine. In the spectrum of the engine vibration signal, the amplitude of these main harmonic orders and its frequency doubling components are prominent, and the first main harmonic order v=τ/2 has the largest amplitude. Therefore, the engine rotational speed can be worked out according to the first main harmonic order when the number of cylinders is known. Supposed that the frequency of the first main harmonic order component is f0, the engine speed is
(2)n=60ω2π=602f0τ=120f0τ.For an automobile engine with stroke numberi and cylinder number τ, the mathematical relationship between rotational speed n and the fundamental frequency f0 of total torque is [37]
(3)n=60i2f0τ.As seen from Equation (3), the key to engine speed estimation methods is to accurately estimate the fundamental frequency f0 of the vibration signal. The proposed method in this paper utilizes a novel TF analysis method to get a TF result with higher frequency resolution and uses the IF detection algorithm to look for the fundamental frequency f0.


### 2.2. Multiple Matching Synchrosqueezing Transform

#### 2.2.1. Matching Synchrosqueezing Transform (MSST)

The STFT ofx∈L2R concerning a real-valued and even function window g∈L2R is represented as follows:
(4)Gt,ω=∫−∞+∞xugu−te−jωu−tdu,in which gu−t denotes the moving window along with the time axis, xu is the analyzed signal, and Gt,ω2is the spectrogram of the analyzed signal.Restricted by the Heisenberg uncertainty principle, the TF result is blurry around the IF position, which means the TF resolution along with time and frequency directions cannot become arbitrarily small at the same time. The SST method is proposed to serve as a postprocessing method of STFT; it squeezes these diffused TF coefficients into the correct IF position and its formula is given as follows:(5)Tt,ξ=∫−∞∞Gt,ωδξ−ω~t,ωdω,where ∫‐∞∞δξ−ω~t,ωdω denotes the synchrosqueezing operator and ω~t,ω is the two-dimensional IF estimator defined as follows:
(6)ω~t,ω=∂tGt,ωjGt,ωGt,ω≠0,where ∂t represents the derivative with respect to t.The unbiased IF estimator of SST in Equation (6) is deduced on the hypothesis that the analyzed signal can only be a pure harmonic form. In other words, to generate a concentrated TFR, the SST requires to be established on the hypothesis that the analyzed signal should be weakly time varying. When addressing the strong frequency-modulated signal, the error between the estimated IF and true IF will become large, which eventually results in a blurry SST representation. To overcome this problem, MSST constructs a matching IF (MIF) estimator that can precisely characterize the time-varying IF with second order, e.g., chirp signal.In the MSST method, the group delay (GD) estimatort~t,ω and chirp rate estimator c~t,ω are introduced; they are defined as follows:
(7)t~t,ω=t+j∂ωGt,ωGt,ω,(8)c~t,ω=∂tω~t,ω∂tt~t,ω,where ∂ω represents the derivative with respect to ω. Combined with the IF estimator in Equation (6), a more precise IF estimator called MIF estimator is defined in MSST and we can formulate it as
(9)ω~mt,ω=ω~t,ω+c~t,ωt−t~t,ω.Because the MIF estimator incorporates the information of the original IF estimator in (5), the GD estimator in (7) and the chirp-rate estimator in (8), compared with the IF estimator in SST, can entirely reduce the IF estimation error of analyzing strong frequency-modulated signal. MSST uses the MIF estimator to replace the IF estimator in SST, which is defined as follows:
(10)Tmt,ξ=∫−∞∞Gt,ωδξ−ω~mt,ωdω,which means that the MSST reassigns the TF coefficients form the t,ω plane to another t,ω~mt,ω plane. Since the coefficients are only assigned in the frequency domain, the reconstruction is still achievable by summing up the synchrosqueezing coefficients in the frequency domain as follows:
(11)xt=12πg0∫ξ−φ′t≤dsTmt,ξdξ.In MSST, the MIF estimator is constructed based on the IF structure of chirp signal, which solves the defect of the IF estimator in SST processing time-varying signal, and improves the resolution of the TF result. For chirp signal or signals that can be locally approximated as chirp signal, MSST can get a more concentrated TF result than SST. However, for nonlinear frequency-modulated signals that do not meet the above requirements, i.e., signals that have strongly time-varying IF, with high order, MSST will still obtain a blurry TF result like SST. The engine vibration signal in real world usually owns a strongly time-varying IF and cannot be locally approximated as chirp signal while the automobile is speeding up or slowing down. Therefore, a novel TFA method with good performance in processing strongly time-varying signal is needed to analyze engine vibration signal for rotational speed estimation.


#### 2.2.2. Multiple Matching Synchrosqueezing Transform (MMSST)

MMSST can be viewed as a right combination of MSST and multiple squeezing technique. Although MSST cannot get an ideal TF result for strong frequency-modulated signals that have time-varying IF with high order, it can obtain an unbiased MIF estimator for chirp signals. Based on the result in (10), a multiple squeeze operation is applied to MSST iteratively to generate a new TF representation with high-frequency resolution, that is,
(12)Tm2t,η=∫−∞∞Tm1t,ωδη−ω~mt,ωdω,Tm3t,η=∫−∞∞Tm2t,ωδη−ω~mt,ωdω,⋯⋯TmNt,η=∫−∞∞TmN−1t,ωδη−ω~mt,ωdω,where Tm1t,ω is the MSST result Tmt,ξ in (9) and NN≥2 denotes the number of iterations, considering the case where N is equal to 2. Then, we can get
(13)Tm2t,η=∫−∞∞Tm1t,ωδη−ω~mt,ωdω=∫−∞∞Tmt,ξδη−ω~mt,ξdξ=∫−∞∞∫−∞∞Gt,ωδξ−ω~mt,ωdωδη−ω~mt,ξdξ=∫−∞∞Gt,ωδη−ω~mt,ω~mt,ωdω.It can be seen that the MMSST (N=2) constructs a novel IF estimator ω~mt,ω~mt,ω. If we further substitute TmN−1t,η into TmNt,η, the IF estimator of MMSST with more iterations can be calculated, for instance, ω~mt,ω~mt,ω is the IF estimator of MMSST (N=3). When iterations is N, the IF estimator is ω~mNt,ω=ω~mt,ω~mNt,ω, where ω~m1t,ω is the IF estimator ω~mt,ω in (9), and the expression of MMSST can be rewritten as
(14)TmNt,η=∫−∞∞Tmt,ωδη−ω~mNt,ωdω.According to the above calculation, the proposed MMSST is easy to implement by executing multiple squeeze operations. The entire procedure of MMSST is summarized in Algorithm1.Algorithm 1: MMSST algorithm.
1: Choose the window function g and iteration numberN;2: Calculate MSSTTmt,ω form (9), and the MIF estimate ω~mt,ω from (10);3: Letω~m1t,ω←ω~mt,ω4: ifN>15: forn=2:N6:ξ←ω~mn‐1t,ω7:ω~mnt,ω←ω~mn‐1t,ξ8: end for9: end if10: LetTmNt,η←0 for all t and η11:η←ω~mNt,ω12:TmNt,η←TmNt,η+Tmt,ω13: OutputTmNt,ηThe multiple squeeze operations can entirely reduce MIF estimation error in analyzing nonlinear strong frequency-modulated signal and improve the resolution of TF representation. Since MMSST reassigns the coefficients only in the frequency direction and no information missing, the signal reconstruction is still achievable by summing the synchrosqueezing coefficients in the frequency domain as follows:(15)xt=12πg0∫ξ−φ′t≤dsTmNt,ηdη


### 2.3. The Instantaneous Frequency (IF) Detection Algorithm

An effective detection algorithm is described in [31], which can provide a reliable IF estimation result. Knowing the number K of modes, this algorithm calculates the local minimum value of the function
(16)Eϕ=∑k=1K−∫−∞+∞TFRt,ϕkt2dt+∫−∞+∞λ·ϕ′t2+βϕkt2dt,where ∑k=1Kt,ϕkt is the estimation of IF trajectories in the TF plane, and λ and β are two parameters to adjust the level of regularization. Equation (16) suggests that both the energy of TF coefficients and the local smoothness of the detected trajectory should be considered. The practical implementation of the IF detection method adopted in this paper is that the entire TF plane is suggested to be first separated into several segments. The starting points are determined from the maximum value in of these individual TF segments, instead of the global TF plane. Then, the forward and backward procedures are executed to search the IF trajectories. Eventually, the best IF trajectory is selected from these trajectories according to the TF energy. Given the TF representation TmNt,η, it is separated into F segments and there are K modes needed to be detected. In the engine speed estimation, only the fundamental frequency of vibration signal needs to be obtained. Since the fundamental frequency has the largest amplitude on the TF plane, for simplicity, it is only necessary to set the parameter K to 1.


### 2.4. The Proposed Method

Combining the proposed MMSST and the IF detection algorithm, the proposed MMSST-based method for rotational speed estimation of an automobile engine in this paper consists of the following steps:(1)
Data acquisition and preprocessing: collect the vibration signals of the engine cylinder through the accelerator sensor, and a low-pass filter is applied to restrain the high-frequency interference(2)
TF representation calculation: use the proposed MMSST to obtain the TF representation of the vibration signal, and the implementation of the MMSST is described in detail in Algorithm1(3)
Rotational speed estimation: utilize the IF detection algorithm to extract the fundamental frequency of the vibration signal from the TF representation obtained by step (2), and then, the rotational speed of engine is calculated through Equation (3)


## 2.1. Estimation Principle

When an engine is running, its power is transmitted and transformed through the crank and connecting rod mechanism. The excitation loads acting on the crankshaft of an internal combustion engine are mainly derived from the driving torque due to gas forces generated within the cylinders. In a single cylinder, the driving torque produced by the gas pressure can be expressed as [8, 36]
(1)Mg=M0+∑v∞Mvsinvωt+φv,where M0 represents the average torque generated by the gas pressure, v denotes harmonic order, Mv and φv represent the initial phase and amplitude of the vth harmonic order of the torque waveform, respectively, and ω is crankshaft angular velocity. Since each cylinder fires once every two engine revolutions in a four-stroke engine, the rate of driving torque change is half of the rotational angular velocity of the crankshaft, there is v=1/2,1,3/2,⋯, and 1/2 is the fundamental harmonic order.For a four-stroke engine withτ cylinders, according the knowledge of engine dynamics, when harmonic order v=τ/2,τ,3τ/2,⋯, the driving torques of each cylinder are in the same direction, and these orders are called main harmonic orders of the engine. In the spectrum of the engine vibration signal, the amplitude of these main harmonic orders and its frequency doubling components are prominent, and the first main harmonic order v=τ/2 has the largest amplitude. Therefore, the engine rotational speed can be worked out according to the first main harmonic order when the number of cylinders is known. Supposed that the frequency of the first main harmonic order component is f0, the engine speed is
(2)n=60ω2π=602f0τ=120f0τ.For an automobile engine with stroke numberi and cylinder number τ, the mathematical relationship between rotational speed n and the fundamental frequency f0 of total torque is [37]
(3)n=60i2f0τ.As seen from Equation (3), the key to engine speed estimation methods is to accurately estimate the fundamental frequency f0 of the vibration signal. The proposed method in this paper utilizes a novel TF analysis method to get a TF result with higher frequency resolution and uses the IF detection algorithm to look for the fundamental frequency f0.


## 2.2. Multiple Matching Synchrosqueezing Transform

### 2.2.1. Matching Synchrosqueezing Transform (MSST)

The STFT ofx∈L2R concerning a real-valued and even function window g∈L2R is represented as follows:
(4)Gt,ω=∫−∞+∞xugu−te−jωu−tdu,in which gu−t denotes the moving window along with the time axis, xu is the analyzed signal, and Gt,ω2is the spectrogram of the analyzed signal.Restricted by the Heisenberg uncertainty principle, the TF result is blurry around the IF position, which means the TF resolution along with time and frequency directions cannot become arbitrarily small at the same time. The SST method is proposed to serve as a postprocessing method of STFT; it squeezes these diffused TF coefficients into the correct IF position and its formula is given as follows:(5)Tt,ξ=∫−∞∞Gt,ωδξ−ω~t,ωdω,where ∫‐∞∞δξ−ω~t,ωdω denotes the synchrosqueezing operator and ω~t,ω is the two-dimensional IF estimator defined as follows:
(6)ω~t,ω=∂tGt,ωjGt,ωGt,ω≠0,where ∂t represents the derivative with respect to t.The unbiased IF estimator of SST in Equation (6) is deduced on the hypothesis that the analyzed signal can only be a pure harmonic form. In other words, to generate a concentrated TFR, the SST requires to be established on the hypothesis that the analyzed signal should be weakly time varying. When addressing the strong frequency-modulated signal, the error between the estimated IF and true IF will become large, which eventually results in a blurry SST representation. To overcome this problem, MSST constructs a matching IF (MIF) estimator that can precisely characterize the time-varying IF with second order, e.g., chirp signal.In the MSST method, the group delay (GD) estimatort~t,ω and chirp rate estimator c~t,ω are introduced; they are defined as follows:
(7)t~t,ω=t+j∂ωGt,ωGt,ω,(8)c~t,ω=∂tω~t,ω∂tt~t,ω,where ∂ω represents the derivative with respect to ω. Combined with the IF estimator in Equation (6), a more precise IF estimator called MIF estimator is defined in MSST and we can formulate it as
(9)ω~mt,ω=ω~t,ω+c~t,ωt−t~t,ω.Because the MIF estimator incorporates the information of the original IF estimator in (5), the GD estimator in (7) and the chirp-rate estimator in (8), compared with the IF estimator in SST, can entirely reduce the IF estimation error of analyzing strong frequency-modulated signal. MSST uses the MIF estimator to replace the IF estimator in SST, which is defined as follows:
(10)Tmt,ξ=∫−∞∞Gt,ωδξ−ω~mt,ωdω,which means that the MSST reassigns the TF coefficients form the t,ω plane to another t,ω~mt,ω plane. Since the coefficients are only assigned in the frequency domain, the reconstruction is still achievable by summing up the synchrosqueezing coefficients in the frequency domain as follows:
(11)xt=12πg0∫ξ−φ′t≤dsTmt,ξdξ.In MSST, the MIF estimator is constructed based on the IF structure of chirp signal, which solves the defect of the IF estimator in SST processing time-varying signal, and improves the resolution of the TF result. For chirp signal or signals that can be locally approximated as chirp signal, MSST can get a more concentrated TF result than SST. However, for nonlinear frequency-modulated signals that do not meet the above requirements, i.e., signals that have strongly time-varying IF, with high order, MSST will still obtain a blurry TF result like SST. The engine vibration signal in real world usually owns a strongly time-varying IF and cannot be locally approximated as chirp signal while the automobile is speeding up or slowing down. Therefore, a novel TFA method with good performance in processing strongly time-varying signal is needed to analyze engine vibration signal for rotational speed estimation.


### 2.2.2. Multiple Matching Synchrosqueezing Transform (MMSST)

MMSST can be viewed as a right combination of MSST and multiple squeezing technique. Although MSST cannot get an ideal TF result for strong frequency-modulated signals that have time-varying IF with high order, it can obtain an unbiased MIF estimator for chirp signals. Based on the result in (10), a multiple squeeze operation is applied to MSST iteratively to generate a new TF representation with high-frequency resolution, that is,
(12)Tm2t,η=∫−∞∞Tm1t,ωδη−ω~mt,ωdω,Tm3t,η=∫−∞∞Tm2t,ωδη−ω~mt,ωdω,⋯⋯TmNt,η=∫−∞∞TmN−1t,ωδη−ω~mt,ωdω,where Tm1t,ω is the MSST result Tmt,ξ in (9) and NN≥2 denotes the number of iterations, considering the case where N is equal to 2. Then, we can get
(13)Tm2t,η=∫−∞∞Tm1t,ωδη−ω~mt,ωdω=∫−∞∞Tmt,ξδη−ω~mt,ξdξ=∫−∞∞∫−∞∞Gt,ωδξ−ω~mt,ωdωδη−ω~mt,ξdξ=∫−∞∞Gt,ωδη−ω~mt,ω~mt,ωdω.It can be seen that the MMSST (N=2) constructs a novel IF estimator ω~mt,ω~mt,ω. If we further substitute TmN−1t,η into TmNt,η, the IF estimator of MMSST with more iterations can be calculated, for instance, ω~mt,ω~mt,ω is the IF estimator of MMSST (N=3). When iterations is N, the IF estimator is ω~mNt,ω=ω~mt,ω~mNt,ω, where ω~m1t,ω is the IF estimator ω~mt,ω in (9), and the expression of MMSST can be rewritten as
(14)TmNt,η=∫−∞∞Tmt,ωδη−ω~mNt,ωdω.According to the above calculation, the proposed MMSST is easy to implement by executing multiple squeeze operations. The entire procedure of MMSST is summarized in Algorithm1.Algorithm 1: MMSST algorithm.
1: Choose the window function g and iteration numberN;2: Calculate MSSTTmt,ω form (9), and the MIF estimate ω~mt,ω from (10);3: Letω~m1t,ω←ω~mt,ω4: ifN>15: forn=2:N6:ξ←ω~mn‐1t,ω7:ω~mnt,ω←ω~mn‐1t,ξ8: end for9: end if10: LetTmNt,η←0 for all t and η11:η←ω~mNt,ω12:TmNt,η←TmNt,η+Tmt,ω13: OutputTmNt,ηThe multiple squeeze operations can entirely reduce MIF estimation error in analyzing nonlinear strong frequency-modulated signal and improve the resolution of TF representation. Since MMSST reassigns the coefficients only in the frequency direction and no information missing, the signal reconstruction is still achievable by summing the synchrosqueezing coefficients in the frequency domain as follows:(15)xt=12πg0∫ξ−φ′t≤dsTmNt,ηdη


## 2.2.1. Matching Synchrosqueezing Transform (MSST)

The STFT ofx∈L2R concerning a real-valued and even function window g∈L2R is represented as follows:
(4)Gt,ω=∫−∞+∞xugu−te−jωu−tdu,in which gu−t denotes the moving window along with the time axis, xu is the analyzed signal, and Gt,ω2is the spectrogram of the analyzed signal.Restricted by the Heisenberg uncertainty principle, the TF result is blurry around the IF position, which means the TF resolution along with time and frequency directions cannot become arbitrarily small at the same time. The SST method is proposed to serve as a postprocessing method of STFT; it squeezes these diffused TF coefficients into the correct IF position and its formula is given as follows:(5)Tt,ξ=∫−∞∞Gt,ωδξ−ω~t,ωdω,where ∫‐∞∞δξ−ω~t,ωdω denotes the synchrosqueezing operator and ω~t,ω is the two-dimensional IF estimator defined as follows:
(6)ω~t,ω=∂tGt,ωjGt,ωGt,ω≠0,where ∂t represents the derivative with respect to t.The unbiased IF estimator of SST in Equation (6) is deduced on the hypothesis that the analyzed signal can only be a pure harmonic form. In other words, to generate a concentrated TFR, the SST requires to be established on the hypothesis that the analyzed signal should be weakly time varying. When addressing the strong frequency-modulated signal, the error between the estimated IF and true IF will become large, which eventually results in a blurry SST representation. To overcome this problem, MSST constructs a matching IF (MIF) estimator that can precisely characterize the time-varying IF with second order, e.g., chirp signal.In the MSST method, the group delay (GD) estimatort~t,ω and chirp rate estimator c~t,ω are introduced; they are defined as follows:
(7)t~t,ω=t+j∂ωGt,ωGt,ω,(8)c~t,ω=∂tω~t,ω∂tt~t,ω,where ∂ω represents the derivative with respect to ω. Combined with the IF estimator in Equation (6), a more precise IF estimator called MIF estimator is defined in MSST and we can formulate it as
(9)ω~mt,ω=ω~t,ω+c~t,ωt−t~t,ω.Because the MIF estimator incorporates the information of the original IF estimator in (5), the GD estimator in (7) and the chirp-rate estimator in (8), compared with the IF estimator in SST, can entirely reduce the IF estimation error of analyzing strong frequency-modulated signal. MSST uses the MIF estimator to replace the IF estimator in SST, which is defined as follows:
(10)Tmt,ξ=∫−∞∞Gt,ωδξ−ω~mt,ωdω,which means that the MSST reassigns the TF coefficients form the t,ω plane to another t,ω~mt,ω plane. Since the coefficients are only assigned in the frequency domain, the reconstruction is still achievable by summing up the synchrosqueezing coefficients in the frequency domain as follows:
(11)xt=12πg0∫ξ−φ′t≤dsTmt,ξdξ.In MSST, the MIF estimator is constructed based on the IF structure of chirp signal, which solves the defect of the IF estimator in SST processing time-varying signal, and improves the resolution of the TF result. For chirp signal or signals that can be locally approximated as chirp signal, MSST can get a more concentrated TF result than SST. However, for nonlinear frequency-modulated signals that do not meet the above requirements, i.e., signals that have strongly time-varying IF, with high order, MSST will still obtain a blurry TF result like SST. The engine vibration signal in real world usually owns a strongly time-varying IF and cannot be locally approximated as chirp signal while the automobile is speeding up or slowing down. Therefore, a novel TFA method with good performance in processing strongly time-varying signal is needed to analyze engine vibration signal for rotational speed estimation.


## 2.2.2. Multiple Matching Synchrosqueezing Transform (MMSST)

MMSST can be viewed as a right combination of MSST and multiple squeezing technique. Although MSST cannot get an ideal TF result for strong frequency-modulated signals that have time-varying IF with high order, it can obtain an unbiased MIF estimator for chirp signals. Based on the result in (10), a multiple squeeze operation is applied to MSST iteratively to generate a new TF representation with high-frequency resolution, that is,
(12)Tm2t,η=∫−∞∞Tm1t,ωδη−ω~mt,ωdω,Tm3t,η=∫−∞∞Tm2t,ωδη−ω~mt,ωdω,⋯⋯TmNt,η=∫−∞∞TmN−1t,ωδη−ω~mt,ωdω,where Tm1t,ω is the MSST result Tmt,ξ in (9) and NN≥2 denotes the number of iterations, considering the case where N is equal to 2. Then, we can get
(13)Tm2t,η=∫−∞∞Tm1t,ωδη−ω~mt,ωdω=∫−∞∞Tmt,ξδη−ω~mt,ξdξ=∫−∞∞∫−∞∞Gt,ωδξ−ω~mt,ωdωδη−ω~mt,ξdξ=∫−∞∞Gt,ωδη−ω~mt,ω~mt,ωdω.It can be seen that the MMSST (N=2) constructs a novel IF estimator ω~mt,ω~mt,ω. If we further substitute TmN−1t,η into TmNt,η, the IF estimator of MMSST with more iterations can be calculated, for instance, ω~mt,ω~mt,ω is the IF estimator of MMSST (N=3). When iterations is N, the IF estimator is ω~mNt,ω=ω~mt,ω~mNt,ω, where ω~m1t,ω is the IF estimator ω~mt,ω in (9), and the expression of MMSST can be rewritten as
(14)TmNt,η=∫−∞∞Tmt,ωδη−ω~mNt,ωdω.According to the above calculation, the proposed MMSST is easy to implement by executing multiple squeeze operations. The entire procedure of MMSST is summarized in Algorithm1.Algorithm 1: MMSST algorithm.
1: Choose the window function g and iteration numberN;2: Calculate MSSTTmt,ω form (9), and the MIF estimate ω~mt,ω from (10);3: Letω~m1t,ω←ω~mt,ω4: ifN>15: forn=2:N6:ξ←ω~mn‐1t,ω7:ω~mnt,ω←ω~mn‐1t,ξ8: end for9: end if10: LetTmNt,η←0 for all t and η11:η←ω~mNt,ω12:TmNt,η←TmNt,η+Tmt,ω13: OutputTmNt,ηThe multiple squeeze operations can entirely reduce MIF estimation error in analyzing nonlinear strong frequency-modulated signal and improve the resolution of TF representation. Since MMSST reassigns the coefficients only in the frequency direction and no information missing, the signal reconstruction is still achievable by summing the synchrosqueezing coefficients in the frequency domain as follows:(15)xt=12πg0∫ξ−φ′t≤dsTmNt,ηdη


## 2.3. The Instantaneous Frequency (IF) Detection Algorithm

An effective detection algorithm is described in [31], which can provide a reliable IF estimation result. Knowing the number K of modes, this algorithm calculates the local minimum value of the function
(16)Eϕ=∑k=1K−∫−∞+∞TFRt,ϕkt2dt+∫−∞+∞λ·ϕ′t2+βϕkt2dt,where ∑k=1Kt,ϕkt is the estimation of IF trajectories in the TF plane, and λ and β are two parameters to adjust the level of regularization. Equation (16) suggests that both the energy of TF coefficients and the local smoothness of the detected trajectory should be considered. The practical implementation of the IF detection method adopted in this paper is that the entire TF plane is suggested to be first separated into several segments. The starting points are determined from the maximum value in of these individual TF segments, instead of the global TF plane. Then, the forward and backward procedures are executed to search the IF trajectories. Eventually, the best IF trajectory is selected from these trajectories according to the TF energy. Given the TF representation TmNt,η, it is separated into F segments and there are K modes needed to be detected. In the engine speed estimation, only the fundamental frequency of vibration signal needs to be obtained. Since the fundamental frequency has the largest amplitude on the TF plane, for simplicity, it is only necessary to set the parameter K to 1.


## 2.4. The Proposed Method

Combining the proposed MMSST and the IF detection algorithm, the proposed MMSST-based method for rotational speed estimation of an automobile engine in this paper consists of the following steps:(1)
Data acquisition and preprocessing: collect the vibration signals of the engine cylinder through the accelerator sensor, and a low-pass filter is applied to restrain the high-frequency interference(2)
TF representation calculation: use the proposed MMSST to obtain the TF representation of the vibration signal, and the implementation of the MMSST is described in detail in Algorithm1(3)
Rotational speed estimation: utilize the IF detection algorithm to extract the fundamental frequency of the vibration signal from the TF representation obtained by step (2), and then, the rotational speed of engine is calculated through Equation (3)


## 3. Numerical Validation

In this section, we focus on the comparisons between the proposed TFA method and other TFA methods in addressing complex signals, for instance, noisy signal and strongly time-varying signal. The comparisons mainly focus on the TF energy concentration and the performance of the IF estimation. According to the characteristics of speed variating engine, the simulation signal is modeled as(17)xt=e−0.2tsin2π20t+2sin4.4t+nt,where nt denotes the white Gaussian noise. For numerical signal, the signal-to-noise ratio (SNR) is set to 30 dB, the sampling frequency is 100 Hz, the time duration is 4 s, and the iterations in multiple squeeze of MMSST is 4.


### 3.1. Comparison Test of TF Energy Concentration

The TF representations generated by SST, MSST, and MMSST are shown in Figures1(a)–1(f), and the local zoom of result follows in the right side. For the result of SST shown in Figure 1(a), it is evident that the energy of TF representation smears around the IF trajectory and the frequency resolution is low. It is because the signal defined in (17) is a sinusoidal frequency-modulated signal that cannot be approximated locally as a pure harmonic signal, decreasing the IF estimator accuracy of the SST. As a result, the TF coefficients cannot reassign into the true IF position as it can be seen in Figure 1(c). Compared with SST, MSST can improve the resolution and obtain a more concentrated TF result when addressing the time-varying signal. However, MSST cannot generate the significantly concentrated result for the nonlinear varying part that cannot approximate locally as chirp signal as mentioned above. It can be observed that the MSST generates a TF representation as blurry as SST in Figure 1(d). This shortcoming limits the application scope of MSST method. Therefore, SST and MSST are not suitable for analyzing the strong nonlinear frequency-modulated signals. As shown in Figures 1(e) and 1(f), the MMSST approach generates a TF result with high resolution, and the energy of TF result concentrated on the true IF trajectory. We can conclude that the MMSST can characterize correctly the strongly time-varying feature, leading to the result closer to an ideal TF representation.Figure 1
TF result of signal (17) by the method: (a) SST, (b) local zoom on SST, (c) MSST, (d) local zoom on MSST, (e) MMSST, and (f) local zoom on MMSST.
(a)(b)(c)(d)(e)(f)To evaluate the energy concentration of different TFA methods quantitatively, in this paper, the Renyi entropy [38] is adopted to evaluate the spectrum concentration. The Renyi entropy of order α for TF result is defined as
(18)Rα=11‐αlog2∫∫TFRt,ωαdtdω∫∫TFRt,ωdtdω,where the order is usually set to 3, TFRt,ω denotes the TF result by different TFA methods. The Renyi entropies of these TF results are calculated and listed in Table 1, where a lower Renyi entropy value represents a more concentrated TF representation. It is shown that MMSST provides the most concentrated result among all TFA methods. To verify the performance of the MMSST, the Renyi entropies of TF representations generated by different TFA methods under different SNR are tested. The results are shown in Figure 2. In the condition of either of high or low SNR, MMSST derives the most concentrated TF representation.Table 1
The Renyi entropy of TF result by three TFA methods (SNR=30dB).
TFASSTMSSTMMSSTRenyi entropy10.95279.8778.7512Figure 2
The Renyi entropy of TF result under different SNR.


### 3.2. Comparison Test of IF Estimation

The IF of the signal is an important feature. It includes essential information about the analyzed object, for example, the rotational speed of an engine. Therefore, we have tested the IF estimation performance of the MMSST method under different noise levels. The signal (17) polluted by Gaussian white noise is employed to analyze, and the SST and MSST methods are selected as the comparison reference. Based on TF results, we can obtain the IF of the analyzed signal by utilizing the IF detection algorithm, and the results can be evaluated by the mean relative error [15] (MRE), which is calculated by
(19)MRE=1NlIF~−IFIF1,where Nl is the length of the IF sequence, •1 denotes the l1-norm, and IF~ and IF represent the estimated IF and the original true IF, respectively. The following experiment is carried out 20 times, and the average performance is recorded as the final result.Figure3 shows that the calculated errors of an estimated IF by the three TFA methods. From the results, it is clear that the MMSST approach achieves more accurate IF estimation than SST. Compared with MSST, MMSST has a similar performance in the condition of high SNR, and the IF estimation error is kept at a stable value. But it can be found that MMSST can get a smaller MRE; in the condition of low SNR, the MMSST provides an estimated IF value closer to the original IF. It means that the MMSST is more robust against noise and can characterize more accurately the strongly time-varying signal.Figure 3
The MRE of the detected IF under different SNR.


### 3.3. Iteration Number Analysis

Since MMSST introduces multiple squeeze operations on MSST to improve the frequency resolution of TF representation, the running time is bound to increase with the increase of iteration times. The efficiency of TFA method is essential in a real-time application. Therefore, we must take both a frequency resolution and a time of running into account to find a suitable iteration number for MMSST.The simulated signal (17) merged in the Gaussian white noise with the SNR=15dB is employed to analyze the impact of iteration number on MMSST method. Figure 4 shows the change of the Renyi entropy of the TF representations generated by MMSST and the running time as the number of iterations increases. It illuminates that the rise in the iterations leads to a more concentrated time-frequency expression, and it also increases the running time. But after the fourth iteration, the Renyi entropy indicates to be stable, which means that when the number of iterations is greater than 4, instead of improving the frequency resolution of TF result, MMSST only increases the running time. Therefore, in this paper, if not mentioned, the number of iterations is set to 4.Figure 4
The Renyi entropy of TFR and the running time concerning the number of iterations (SNR=15dB).


## 3.1. Comparison Test of TF Energy Concentration

The TF representations generated by SST, MSST, and MMSST are shown in Figures1(a)–1(f), and the local zoom of result follows in the right side. For the result of SST shown in Figure 1(a), it is evident that the energy of TF representation smears around the IF trajectory and the frequency resolution is low. It is because the signal defined in (17) is a sinusoidal frequency-modulated signal that cannot be approximated locally as a pure harmonic signal, decreasing the IF estimator accuracy of the SST. As a result, the TF coefficients cannot reassign into the true IF position as it can be seen in Figure 1(c). Compared with SST, MSST can improve the resolution and obtain a more concentrated TF result when addressing the time-varying signal. However, MSST cannot generate the significantly concentrated result for the nonlinear varying part that cannot approximate locally as chirp signal as mentioned above. It can be observed that the MSST generates a TF representation as blurry as SST in Figure 1(d). This shortcoming limits the application scope of MSST method. Therefore, SST and MSST are not suitable for analyzing the strong nonlinear frequency-modulated signals. As shown in Figures 1(e) and 1(f), the MMSST approach generates a TF result with high resolution, and the energy of TF result concentrated on the true IF trajectory. We can conclude that the MMSST can characterize correctly the strongly time-varying feature, leading to the result closer to an ideal TF representation.Figure 1
TF result of signal (17) by the method: (a) SST, (b) local zoom on SST, (c) MSST, (d) local zoom on MSST, (e) MMSST, and (f) local zoom on MMSST.
(a)(b)(c)(d)(e)(f)To evaluate the energy concentration of different TFA methods quantitatively, in this paper, the Renyi entropy [38] is adopted to evaluate the spectrum concentration. The Renyi entropy of order α for TF result is defined as
(18)Rα=11‐αlog2∫∫TFRt,ωαdtdω∫∫TFRt,ωdtdω,where the order is usually set to 3, TFRt,ω denotes the TF result by different TFA methods. The Renyi entropies of these TF results are calculated and listed in Table 1, where a lower Renyi entropy value represents a more concentrated TF representation. It is shown that MMSST provides the most concentrated result among all TFA methods. To verify the performance of the MMSST, the Renyi entropies of TF representations generated by different TFA methods under different SNR are tested. The results are shown in Figure 2. In the condition of either of high or low SNR, MMSST derives the most concentrated TF representation.Table 1
The Renyi entropy of TF result by three TFA methods (SNR=30dB).
TFASSTMSSTMMSSTRenyi entropy10.95279.8778.7512Figure 2
The Renyi entropy of TF result under different SNR.


## 3.2. Comparison Test of IF Estimation

The IF of the signal is an important feature. It includes essential information about the analyzed object, for example, the rotational speed of an engine. Therefore, we have tested the IF estimation performance of the MMSST method under different noise levels. The signal (17) polluted by Gaussian white noise is employed to analyze, and the SST and MSST methods are selected as the comparison reference. Based on TF results, we can obtain the IF of the analyzed signal by utilizing the IF detection algorithm, and the results can be evaluated by the mean relative error [15] (MRE), which is calculated by
(19)MRE=1NlIF~−IFIF1,where Nl is the length of the IF sequence, •1 denotes the l1-norm, and IF~ and IF represent the estimated IF and the original true IF, respectively. The following experiment is carried out 20 times, and the average performance is recorded as the final result.Figure3 shows that the calculated errors of an estimated IF by the three TFA methods. From the results, it is clear that the MMSST approach achieves more accurate IF estimation than SST. Compared with MSST, MMSST has a similar performance in the condition of high SNR, and the IF estimation error is kept at a stable value. But it can be found that MMSST can get a smaller MRE; in the condition of low SNR, the MMSST provides an estimated IF value closer to the original IF. It means that the MMSST is more robust against noise and can characterize more accurately the strongly time-varying signal.Figure 3
The MRE of the detected IF under different SNR.


## 3.3. Iteration Number Analysis

Since MMSST introduces multiple squeeze operations on MSST to improve the frequency resolution of TF representation, the running time is bound to increase with the increase of iteration times. The efficiency of TFA method is essential in a real-time application. Therefore, we must take both a frequency resolution and a time of running into account to find a suitable iteration number for MMSST.The simulated signal (17) merged in the Gaussian white noise with the SNR=15dB is employed to analyze the impact of iteration number on MMSST method. Figure 4 shows the change of the Renyi entropy of the TF representations generated by MMSST and the running time as the number of iterations increases. It illuminates that the rise in the iterations leads to a more concentrated time-frequency expression, and it also increases the running time. But after the fourth iteration, the Renyi entropy indicates to be stable, which means that when the number of iterations is greater than 4, instead of improving the frequency resolution of TF result, MMSST only increases the running time. Therefore, in this paper, if not mentioned, the number of iterations is set to 4.Figure 4
The Renyi entropy of TFR and the running time concerning the number of iterations (SNR=15dB).


## 4. Experimental Validation

In this section, to evaluate the accuracy of MMSST-based method in engine rotational speed estimation, a test experiment is carried out in an automobile inspection station, and the tested automobile is randomly selected The tested engine is a four-cylinder four-stroke diesel engine. The acceleration sensor is assembled in the form of a magnetic suction probe, and it is attached on the engine hood of the tested automobile, and the sample rate is set to 512 Hz.The TF representations of the engine vibration signal generated by SST, MSST, and MMSST are shown in Figures5(a), 5(c), and 5(e). The relative error of the rotational speed estimation method are provided in Figure 6. MMSST provides the most concentrated TF result. The Renyi entropies of TF results, generated by different methods, are calculated and listed in Table 2 to evaluate the performance of these methods. As the engine vibration signal contains noise and harmonic interference, and the frequency range is wide, the Renyi entropy is relatively large, but the MMSST gives the lowest Renyi entropy among these TFA methods.Figure 5
The TF results of engine vibration signal and the test results of rotational speed estimation for the applied methods: (a) SST and (b) SST-based estimation result, (c) MSST and (d) MSST-based estimation result, and (e) MMSST and (f) MMSST-based estimation result.
(a)(b)(c)(d)(e)(f)Figure 6
The relative error of the rotational speed estimation method: (a) SST-based estimation result, (b) MSST-based estimation result, and (c) MMSST-based estimation result.
(a)(b)(c)Table 2
The Renyi entropy and MRE of three TFA-based methods.
MethodsSST-basedMSST-basedMMSST-basedRenyi entropy17.910418.025917.6883MRE0.00710.00770.0051Based on the TF result, we can obtain the fundamental frequency (related to the engine rotational speed) by using the IF detection algorithm and calculate the rotational speed according to the relationship between speed and frequency. The results of rotational speed estimation by the proposed method, SST-based and MSST-based, are displayed in Figures5(b), 5(d), and 5(f). In Figure 5(b), the SST-based method produces a large estimation error in about 50-65 seconds when the engine in accelerating stage and the local features are shown in the blue area. As the rotational speed changes rapidly, so does the rotational speed-related fundamental frequency. According to the previous description, SST fails to obtain a concentrated TF result for time-varying signals, resulting in substantial estimation errors. In Figure 5(d), due to the fact that the vibration signal cannot be approximated to a chirp signal locally, the MSST-based method can reduce the estimation error that occurs in SST-based method when the engine is in the accelerating stage, but it cannot eliminate the error. Meanwhile, it also produces an error in the idling stage, and the local features are displayed in a blue and a grey area, respectively. Since the TF result generated by MSST is not concentrated, the MSST-based method leads to extract a wrong fundamental frequency, related to rotational speed. As shown in Figure 5(f), the proposed method can provide a reliable estimation result, and the result is well matched to the right rotational speed. For the engine vibration signal, MMSST can give a full TF result, which provides a precise way to determine the fundamental frequency that is related to the speed.In the end, the relative error of the rotational speed estimation results obtained by SST-based, MSST-based, and the proposed method are presented in Figure2, and the corresponding MRE is calculated and listed in Table 2. It is obvious that the results of SST-based and MSST-based methods show a large deviation, while the proposed method provides a correct and robust estimation result, which illuminates that the proposed method is feasible the accurate estimation of engine rotational speed.


## 5. Conclusions

In this paper, an accurate method for estimating the automobile engine rotational speed based on vibration signals has been proposed. The proposed method consists of a proposed MMSST and the IF detection algorithm. The proposed MMSST is an improvement of the MSST, which overcomes the shortcoming that MSST cannot provide an accurate IF estimator for nonlinear strong frequency-modulated signals by introducing multiple squeeze operations in MSST. Therefore, MMSST could obtain a TF representation with higher frequency resolution for automobile engine signals. The IF detection algorithm is used to extract the IF from the TF representation. As compared to the commonly used contact methods for measuring the engine rotational speed, the proposed method does not require any installation of encoder disks and only needs to use the proposed algorithm to analyze the engine vibration signal.Results from numerical simulations and experiment on real engine have proven that the proposed method could provide accurate IF estimation results and more concentrated TF result. The numerical simulation also shows that the proposed MMSST can obtain a more robust TF result even in a high noise environment compared with MSST and SST. In addition, the results of experiment on real engine show that the proposed method has the lowest MRE value at 0.0051, compared with 0.0077 in MSST-based method and 0.0071 in SST-based method. Meanwhile, the proposed method can accurately estimate the speed when the engine is working in a steady state or under conditions of large fluctuations in speed, which verifies the proposed MMSST can characterize correctly the strongly time-varying feature of considered signal. The proposed method does not require complex sensor equipment and can conveniently get the accurate rotational speed estimation results; it is suitable for routine inspection and maintenance of an automobile engine.



---

*Source: 6650432-2021-04-09.xml*