Parameterized Instantaneous Frequency Estimation Method for Vibration Signal with Nonlinear Frequency Modulation

Abstract

The vibration signal from the rotatory machinery condition monitoring under time-varying speed is usually amplitude-modulated (AM) and frequency-modulated (FM). It is important to efficiently and accurately estimate the instantaneous frequency (IF) of the vibration signal. In this paper, a novel parameterized IF estimation method is proposed. The method employs a high-order polynomial function to approximate the nonlinear IF and subsequently constructs overdetermined systems of linear equations by calculating the Fourier transform of the derivative of the signal. The IF can be estimated by using least squares estimation to solve the equations. The proposed method has high computational efficiency because it can obtain the estimation of IF over a period of time simultaneously; it differs from traditional time-frequency analysis methods that need to calculate the IF at each point in the time axis. It is demonstrated that the proposed method is not only particularly powerful for the nonlinear FM mono-component signal but also applicable to the multi-component signal constructed by multiple harmonics. The numerical simulation validates the effectiveness of the proposed method, and the experiment’s results show that the method is suitable for the IF estimation of the vibration signal from the varying-speed rotor system.

1. Introduction

The condition monitoring of rotating machinery is important to ensure the smooth operation of machines. It helps develop effective condition-based maintenance strategies to shorten unnecessary downtime and reduce labor costs of inspection, in recognizing defects before catastrophic faults to enhance the availability of systems [1]. The vibration measurement method as reliable, effective, and non-intrusive technology has been widely used in rotating machinery conditions monitoring during startup, shutdown, and normal operations [2]. However, the rotating machinery usually operates in nonstationary conditions [3], such as variable rotating speed and unsteady load. It results in the nonstationary vibration signal, which can be modeled as a superposition of several amplitude-modulated (AM) and frequency-modulated (FM) components [4]. The instantaneous frequency (IF) and instantaneous amplitude (IA) as time-frequency features often reveal important physical characteristics of the nonstationary vibration signal.

The time-frequency analysis (TFA) technology is considered one of the most powerful methods to reveal the time-varying characteristics of nonstationary signals [5] as it can map a one-dimensional nonstationary signal into a two-dimensional time-frequency representation (TFR) [6]. Some traditional TFA methods, such as the short-time Fourier transform (STFT) [7] and the continuous wavelet transform (CWT) [8], have been widely used to extract the time-varying features for nonstationary signals [9]. The STFT uses a fixed-length sliding window, which is only available for quasi-stationary signals; the CWT is a powerful method to obtain time and frequency resolution simultaneously in the TFR [10,11], which uses a longer window to analyze the lower frequency that provides a higher frequency resolution and a lower time resolution, and a shorter window to analyze the higher frequency that provides a lower frequency resolution and a higher time resolution [12]. However, as typical linear TFA methods, both STFT and CWT often generate a blurry TFR and fail to extract nonlinear time-varying features for the nonstationary signal, especially when the law of the IF is nonlinear [13].

In recent years, a lot of advanced TFA methods have been put forward to improve the performance of traditional TFA methods. Daubechies et al. [14] proposed the synchrosqueezing transform (SST) based on CWT as a TFA postprocessing method, which can generate high energy-concentrated TFR for a signal with weakly AM-FM components. To obtain more accurate estimates of IFs for strongly AM-FM signals, Meignen et al. [15,16,17] proposed the second-order STFT-based SST (FSST2) and high-order STFT-based SST (FSSTN). However, FSST2 is only demonstrated to generate the high energy-concentrated TFR for linear chirps with Gaussian modulated amplitudes [17], and FSSTN is sensitive to noise and has high computational complexity [18,19]. To improve the capacity of the chirplet transform [20] for analyzing nonstationary signals with nonlinear laws of IFs, the general parameterized time-frequency (GPTF) transform [12,21] was proposed by introducing the generalized kernel function based on the frequency rotation operator and frequency shift operator. By replacing the different kernel functions, the GPTF transform has different variants to match different nonlinear laws of IFs for signals, such as the polynomial chirplet transform (PCT) [22], the spline chirplet transform (SCT) [23], and the generalized warblet transform (GWT) [24]. However, these variants are mainly focused on mono-component nonstationary signals and are less capable of analyzing multi-component signals.

In the present work, we proposed a novel TFA method that adopts the high-order polynomial to approximate the IF of signals with nonlinear frequency modulation, which could rapidly and accurately estimate the IF, especially for signals with nonlinear laws of IFs, named the parameterized IF estimation method (PIFEM). By calculating the Fourier transform of the derivative of the signal, the overdetermined systems of linear equations can be constructed and then solved by the least squares estimation method to obtain the estimated IF. The proposed method has high computational efficiency because it can obtain the estimate of IF over a period of time simultaneously, which differs from traditional time-frequency analysis methods that need to calculate the IF at each time instant. The PIFEM is particularly powerful for nonlinear FM mono-component signals and is also suitable for multi-component signals constructed by multiple harmonics.

The rest of this paper is organized as follows: Section 2 elaborates the principle of the PIFEM algorithm and develops an iterative procedure for the mono-component nonstationary signal, Section 3 verifies the effectiveness of the PIFEM algorithm by the numerical simulation, and the experiment results in Section 4 show that the PIFEM can be used to estimate the IF of the vibration signal from the rotor system under time-varying speed. Finally, the conclusion is drawn in Section 5.

2. Parameterized Instantaneous Frequency Estimation Method

2.1. Mono-Component Signal

Let us firstly consider a mono-component nonstationary signal defined as follows:

$$s\left(t\right)=A\left(t\right)e^{j\phi \left(t\right)}=e^{\ln A\left(t\right)+j\phi \left(t\right)}$$

(1)

where $A\left(t\right)>0$ is the IA and $\phi \left(t\right)$ is the instantaneous phase (IP). Then, the derivative of $s\left(t\right)$ can be defined as follows:

$$s^{\text{'}}\left(t\right)=\left[{\left[\ln \left(A\right)\right]}^{\text{'}}\left(t\right){+j\phi }^{\text{'}}\left(t\right)\right]s\left(t\right)$$

(2)

where ${\phi }^{\text{'}}\left(t\right)>0$ is the IF. If ${\left[\ln \left(A\right)\right]}^{\text{'}}\left(t\right){+j\phi }^{\text{'}}\left(t\right)$ can be approximated by the Nth-order polynomial defined as follows:

$${\left[\ln \left(A\right)\right]}^{\text{'}}\left(t\right){+j\phi }^{\text{'}}\left(t\right)={\sum }_{n=0}^{N-1}\left({\alpha }_n{+j\beta }_n\right)t^n.$$

(3)

Subsequently, the Fourier transform of $s^{\text{'}}\left(t\right)$ can be calculated as follows:

$$V_{s^{\text{'}}}\left(\omega \right){=j\omega V}_s\left(\omega \right)={\sum }_{n=0}^{N-1}\left({\alpha }_n{+j\beta }_n\right)V_s^{t^n}\left(\omega \right),$$

(4)

where $V_s^{t^n}\left(\omega \right)$ denotes the Fourier transforms of $t^{n}s\left(t\right)$.

Because Equation (4) is always valid for every frequency along the whole frequency axis, Equation (4) can be then considered as a system of equations in a form of a matrix as follows:

$$\mathit{y}=\mathit{A}\mathit{x}$$

(5)

with

$$\{\begin{array}{l}{\mathit{y}=j\omega V}_s\left(\omega \right)\in C^{M\times 1}\hfill \\ \begin{array}{l}\mathit{A}=\left[V_s^{t^0}\left(\omega \right){,V}_s^{t^1}\left(\omega \right){,\dots ,V}_s^{t^{N-1}}\left(\omega \right)\right]\in C^{M\times N}\hfill \\ \mathit{x}={\left[{\alpha }_0{+j\beta }_0{,\alpha }_1{+j\beta }_1{,\dots ,\alpha }_{N-1}{+j\beta }_{N-1}\right]}^T\in C^{N\times 1}\hfill \end{array}\hfill \end{array}$$

(6)

where $M$ is the bandwidth of the frequency axis. Therefore, $\mathit{y}$ is an M-vector called observation vector, $\mathit{A}$ is an $M\times N$ matrix called coefficient matrix, and $x$ is the N-vector of unknown parameters. When $M> N$, Equation (5) becomes an overdetermined system of linear equations, and the least squares estimation of $x$ can be obtained by

$$\hat{\mathit{x}}={\left({\mathit{A}}^{\mathrm{H}}\mathit{A}\right)}^{-1}{\mathit{A}}^{\mathrm{H}}\mathit{y}$$

(7)

where ${\mathit{A}}^{\mathrm{H}}$ is the conjugate transpose matrix of $\mathit{A}$. According to Equation (3), the imaginary part of $x$ is the parameter of the estimated IF of $s\left(t\right)$. Then, the IF can be estimated by

$$\widehat{\mathit{IF}}=\mathrm{Im}\left(\mathit{T}\hat{\mathit{x}}\right)$$

(8)

with

$$\mathit{T}=\left[t^0{,t}^1{,\dots ,t}^{N-1}\right]\in {\mathrm{C}}^{L\times N}$$

(9)

where $\mathrm{Im}(·)$ denotes taking the imaginary part and $L$ is the length of $s\left(t\right)$. Thus, the IF can be estimated over the whole time simultaneously if the IF can be well approximated by a polynomial.

Since the desired signal is always contaminated by noise in real applications, we develop an iterative procedure to estimate the IF for mono-component signals with noise. The basic principle of the procedure is performing the frequency demodulation and the band-pass filtering simultaneously in every iteration. The frequency demodulation operation can gradually change a wideband mono-component signal into a narrow band signal, and then the noise can be filtered out by using a band-pass filter with a progressively narrower pass bandwidth. The bandwidth of the filter can be simply updated by the maximum and minimum of the estimated IF in the last iteration. When more noise is filtered out, the estimated IF is more accurate. Finally, the wideband mono-component signal will be demodulated to an approximate single-tone signal. The iteration can be stopped when the estimated IF is no longer changing. The mean absolute percentage error (MAPE) can be used as the break condition, which is defined as follows:

$$\zeta =\frac{1}{L}{\sum }_{l=1}^L\left|\frac{{\widehat{\mathit{IF}}}_i\left(l\right)-{\widehat{\mathit{IF}}}_{i-1}\left(l\right)}{{\widehat{\mathit{IF}}}_i\left(l\right)}\right|<\epsilon .$$

(10)

where L is the sampling number of the signal, and ${\widehat{IF}}_i$ and ${\widehat{IF}}_{i-1}$ are the IF estimated by the current and last iteration, respectively.

The iterative procedure is defined as Algorithm 1. The steps of Algorithm 1 are given as follows:

Algorithm 1

Input: the signal $s\left(t\right)$, the sampling frequency $f_s$, the maximum iteration $I_{\max }$, the error threshold $\epsilon$, and the order of the polynomial $N$.  Initialization: let $i=1$, $s_1=s\left(t\right)$; the initial band-pass filter ${\mathit{BF}}_1$ is the whole frequency axis, and the initial ${\widehat{\mathit{IF}}}_0$ is all zeros.  Loop: while ${i< I}_{\max }$ and $\zeta <\epsilon$  1. Constructing the observation vector and coefficient matrix according to Equation (6) for $s_i$.  2. Filtering out the noise by ${\mathit{BF}}_i$.  3. Estimating the parameters of the ${\widehat{\mathit{IF}}}_i$ according to Equation (7).  4. Calculating the ${\widehat{\mathit{IF}}}_i$ according to Equations (8) and (9).  5. Constructing the phase demodulation function ${\phi }_i$ by ${\widehat{\mathit{IF}}}_i$.  6. Updating the $s_{i+1}$ by demodulating $s_i$ by ${\phi }_i$.  7. Updating the ${\mathit{BF}}_{i+1}$.  8. Calculating the break condition $\zeta$ according to Equation (10), and then $i=i+1$.  End of Loop

End

It is worth noting that, for every iteration, it just needs to execute N Fast Fourier Transform (FFT) operations, and the computational complexity of the least squares estimation in Equation (7) is O(${\mathit{MN}}^2$) for ${\mathit{A}}^{\mathrm{H}}\mathit{A}$, O(MN) for ${\mathit{A}}^{\mathrm{H}}\mathit{y}$, and O($N^3$) for the matrix inversion.

2.2. Multi-Component Signal

The Algorithm 1 above can accurately estimate the nonlinear IF for mono-component signals, but it cannot be directly used to deal with general multi-component nonstationary signals because it is usually difficult to calculate the Fourier transform for each component individually. However, let us consider a special kind of multi-component signal that is constructed by multiple nonstationary harmonic components and defined as follows:

$$s\left(t\right)={\sum }_{k=1}^K{s}_k\left(t\right)={{\displaystyle \sum }}_{k=1}^K{a}_k{A}_0\left(t\right)e^{{\mathit{jp}}_k{\phi }_0\left(t\right)},$$

(11)

with

$$I_k\left(t\right)={\ln A}_0\left(t\right){+\mathit{jp}}_k{\phi }_0\left(t\right)$$

(12)

and

$$I_k^{\text{'}}\left(t\right)={\left[\ln \left(A_0\right)\right]}^{\text{'}}\left(t\right){+\mathit{jp}}_k{\phi }_0^{\text{'}}\left(t\right)={{\displaystyle \sum }}_{n=0}^{N-1}\left({\alpha }_n{+\mathit{jp}}_k{\beta }_n\right)t^n$$

(13)

where $a_k>0$, $p_k>0$, $a_1{=p}_1=1$ and $I_k^{\text{'}}\left(t\right)$ can be represented by an Nth-order polynomial. So, $A_0\left(t\right)$, ${\phi }_0\left(t\right)$, and ${\phi }_0^{\text{'}}\left(t\right)$ are the IA, IP, and IF of the fundamental harmonic component respectively. Then, Equation (4) can be rewritten as follows:

$$V_{s^{\text{'}}}\left(\omega \right){=j\omega V}_s\left(\omega \right)={\sum }_{k=1}^K{{\displaystyle \sum }}_{n=0}^{N-1}\left({\alpha }_n{+\mathit{jp}}_k{\beta }_n\right)V_{s_k}^{t^n}\left(\omega \right),$$

(14)

where $V_{s_k}^{t^n}\left(\omega \right)$ denotes the Fourier transform of $t^n{s}_k\left(t\right)$. It can be found that the overdetermined system of linear equations cannot be constructed by Equation (14) because $p_k$ and $V_{s_k}^{t^n}\left(\omega \right)$ are unobtainable, and only $V_s^{t^n}\left(\omega \right)={{\displaystyle \sum }}_{k=1}^K{V}_{s_k}^{t^n}\left(\omega \right)$ can be calculated.

However, when still using $\mathit{y}$ and $\mathit{A}$ defined in Equation (6) to estimate $x$ according to Equation (7), an interesting result can be obtained as follows:

$$\mathrm{Im}\left(\widehat{\mathit{x}}\right)\approx P{\left[{\beta }_0{,\beta }_1{,\dots ,\beta }_{N-1}\right]}^{\mathrm{T}}\in {\mathrm{C}}^{N\times 1}$$

(15)

with

$$P=\frac{{{\displaystyle \sum }}_{k=1}^K{p}_k{a}_k^2}{{{\displaystyle \sum }}_{k=1}^K{a}_k^2},$$

(16)

where $P$ is a positive real number depending on $p_k$ and $a_k$, which means that we can obtain a special IF that is a multiple of the IF of the fundamental harmonic component.

It can be verified as follows:

Considering a multi-component nonstationary signal defined in Equation (11); then, we can obtain

$$\{\begin{array}{l}\mathit{y}={{\displaystyle \sum }}_{k=1}^K{\mathit{y}}_k\hfill \\ \begin{array}{l}\mathit{A}={{\displaystyle \sum }}_{k=1}^K{\mathit{A}}_k\hfill \\ {\mathit{y}}_k{=\mathit{A}}_k{\mathit{x}}_k\hfill \end{array}\hfill \end{array},$$

(17)

with

$$\{\begin{array}{l}{\mathit{y}}_k{=j\omega V}_{s_k}\left(\omega \right)\in C^{M\times 1}\hfill \\ \begin{array}{l}{\mathit{A}}_k=\left[V_{s_k}^{t^0}\left(\omega \right){,V}_{s_k}^{t^1}\left(\omega \right){,\dots ,V}_{s_k}^{t^{N-1}}\left(\omega \right)\right]\in C^{M\times N}\hfill \\ {\mathit{x}}_k={\left[{\alpha }_0{+\mathit{jp}}_k{\beta }_0{,\alpha }_1{+jp}_k{\beta }_1{,\dots ,\alpha }_{N-1}{+jp}_k{\beta }_{N-1}\right]}^{\mathrm{T}}\in C^{N\times 1}\hfill \end{array}\hfill \end{array}$$

(18)

where ${\mathit{y}}_k$, ${\mathit{A}}_k$, and ${\mathit{x}}_k$ corresponding to $s_k\left(t\right)$ are constructed according to Equation (6), respectively.

Then, ${\mathit{A}}^{\mathrm{H}}\mathit{A}$ and ${\mathit{A}}^{\mathrm{H}}\mathit{y}$ can be calculated as follows:

$$\begin{array}{c}{\mathit{A}}^{\mathrm{H}}\mathit{A}={\left({{\displaystyle \sum }}_{k_1=1}^K{\mathit{A}}_{k_1}\right)}^{\mathrm{H}}\left({{\displaystyle \sum }}_{k_2=1}^K{\mathit{A}}_{k_2}\right)={{\displaystyle \sum }}_{k_2=1}^K{{\displaystyle \sum }}_{k_1=1}^K{\mathit{A}}_{k_1}^{\mathrm{H}}{\mathit{A}}_{k_2}\\ {\mathit{A}}^{\mathrm{H}}\mathit{y}={\left({{\displaystyle \sum }}_{k_1=1}^K{\mathit{A}}_{k_1}\right)}^{\mathrm{H}}\left({{\displaystyle \sum }}_{k_2=1}^K{\mathit{y}}_{k_2}\right)={{\displaystyle \sum }}_{k_2=1}^K{{\displaystyle \sum }}_{k_1=1}^K{\mathit{A}}_{k_1}^{\mathrm{H}}{\mathit{y}}_{k_2}\end{array}$$

(19)

with

$$\begin{array}{l}{\mathit{A}}_{k_1}^{\mathrm{H}}{\mathit{A}}_{k_2}=\left[\begin{array}{ccc}{{\displaystyle \sum }}_{\omega }{V}_{s_{k_1}^{*}}^{t^0}\left(\omega \right)V_{s_{k_2}}^{t^0}\left(\omega \right)& \cdots & {{\displaystyle \sum }}_{\omega }{V}_{s_{k_1}^{*}}^{t^0}\left(\omega \right)V_{s_{k_2}}^{t^{N-1}}\left(\omega \right)\\ ⋮& \ddots & ⋮\\ {{\displaystyle \sum }}_{\omega }{V}_{s_{k_1}^{*}}^{t^{N-1}}\left(\omega \right)V_{s_{k_2}}^{t^0}\left(\omega \right)& \cdots & {{\displaystyle \sum }}_{\omega }{V}_{s_{k_1}^{*}}^{t^{N-1}}\left(\omega \right)V_{s_{k_2}}^{t^{N-1}}\left(\omega \right)\end{array}\right].\\ {\mathit{A}}_{k_1}^{\mathrm{H}}{\mathit{y}}_{k_2}=\left[\begin{array}{c}{{\displaystyle \sum }}_{\omega }{V}_{s_{k_1}^{*}}^{t^0}\left(\omega \right){j\omega V}_{s_{k_2}}\left(\omega \right)\\ ⋮\\ {{\displaystyle \sum }}_{\omega }{V}_{s_{k_1}^{*}}^{t^{N-1}}\left(\omega \right){j\omega V}_{s_{k_2}}\left(\omega \right)\end{array}\right]\end{array}$$

(20)

It can be found that the element of ${\mathit{A}}_{k_1}^{\mathrm{H}}{\mathit{A}}_{k_2}$ is the power spectral density (PSD) of $t^m{s}_{k_1}\left(t\right)$ and $t^n{s}_{k_2}\left(t\right)$, where $m=n$ on the primary diagonal and $V_{s_{k_1}^{*}}^{t^m}\left(\omega \right)$ denotes the conjugate of $V_{s_{k_1}}^{t^m}\left(\omega \right)$. Then, according to Parseval’s theorem, the PSDs can be calculated as follows：

$$\begin{array}{l}{{\displaystyle \sum }}_{\omega }{V}_{s_{k_1}^{*}}^{t^m}\left(\omega \right)V_{s_{k_2}}^{t^n}\left(\omega \right)={{\displaystyle \sum }}_t{t}^m{s}_{k_1}^{*}\left(t\right)t^n{s}_{k_2}\left(t\right){=a}_{k_1}{a}_{k_2}{{\displaystyle \sum }}_t{t}^{m+n}{A}_0^2\left(t\right)e^{j\left(p_{k_2}{-p}_{k_1}\right){\phi }_0\left(t\right)}\hfill \\ {{\displaystyle \sum }}_{\omega }{V}_{s_{k_1}^{*}}^{t^m}\left(\omega \right){j\omega V}_{s_{k_2}}\left(\omega \right)={{\displaystyle \sum }}_t{t}^m{s}_{k_1}^{*}\left(t\right)s_{k_2}^{\text{'}}\left(t\right){=a}_{k_1}{a}_{k_2}{{\displaystyle \sum }}_t{t}^m{A}_0^2\left(t\right)I_{k_2}^{\text{'}}\left(t\right)e^{j\left(p_{k_2}{-p}_{k_1}\right){\phi }_0\left(t\right)}\hfill \end{array}.$$

(21)

When $k_1\ne k_2$, Equation (21) represents the cross-correlation of two different components, which is usually small enough and can be ignored. Thus, ${\mathit{A}}^{\mathrm{H}}\mathit{A}$ and ${\mathit{A}}^{\mathrm{H}}\mathit{y}$ in Equation (19) can be rewritten as follows:

$$\begin{array}{l}{\mathit{A}}^{\mathrm{H}}\mathit{A}={{\displaystyle \sum }}_{k=1}^K{\mathit{A}}_k^{\mathrm{H}}{\mathit{A}}_k\approx {{\displaystyle \sum }}_{k=1}^K{a}_k^2\left[\begin{array}{ccc}{{\displaystyle \sum }}_t{t}^0{A}_0^2\left(t\right)& \cdots & {{\displaystyle \sum }}_t{t}^{N-1}{A}_0^2\left(t\right)\\ ⋮& \ddots & ⋮\\ {{\displaystyle \sum }}_t{t}^{N-1}{A}_0^2\left(t\right)& \cdots & {{\displaystyle \sum }}_t{t}^{2N-2}{A}_0^2\left(t\right)\end{array}\right]={{\displaystyle \sum }}_{k=1}^K{a}_k^2{\mathit{A}}_1^{\mathrm{H}}{\mathit{A}}_1,\\ {\mathit{A}}^{\mathrm{H}}\mathit{y}={{\displaystyle \sum }}_{k=1}^K{\mathit{A}}_k^{\mathrm{H}}{\mathit{y}}_k\approx {{\displaystyle \sum }}_{k=1}^K{a}_k^2\left[\begin{array}{c}{{\displaystyle \sum }}_t{t}^0{I}_k^{\text{'}}\left(t\right)A_0^2\left(t\right)\\ ⋮\\ {{\displaystyle \sum }}_t{t}^{N-1}{I}_k^{\text{'}}\left(t\right)A_0^2\left(t\right)\end{array}\right]={{\displaystyle \sum }}_{k=1}^K{a}_k^2\left[\mathrm{Re}\left({\mathit{A}}_1^{\mathrm{H}}{\mathit{y}}_1\right){+\mathit{jp}}_k\mathrm{Im}\left({\mathit{A}}_1^{\mathrm{H}}{\mathit{y}}_1\right)\right]\end{array}$$

(22)

where $\mathrm{Re}(·)$ and $\mathrm{Im}(·)$ denote taking the real and imaginary part, respectively.

Then, according to Equation (7), we can have the following:

$$\widehat{\mathit{x}}={\left({\mathit{A}}^{\mathrm{H}}\mathit{A}\right)}^{-1}{\mathit{A}}^{\mathrm{H}}\mathit{y}\approx {\left({\mathit{A}}_1^{\mathrm{H}}{\mathit{A}}_1\right)}^{-1}\mathrm{Re}\left({\mathit{A}}_1^{\mathrm{H}}{\mathit{y}}_1\right)+j\left[{\left({\mathit{A}}_1^{\mathrm{H}}{\mathit{A}}_1\right)}^{-1}\mathrm{Im}\left({\mathit{A}}_1^{\mathrm{H}}{\mathit{y}}_1\right)\right]\frac{{{\displaystyle \sum }}_{k=1}^K{p}_k{a}_k^2}{{{\displaystyle \sum }}_{k=1}^K{a}_k^2}.$$

(23)

Finally, we can obtain Equations (15) and (16).

3. Simulation and Results

In this part, to show and verify the performance of our proposed method, let us first consider a mono-component signal with the nonlinear law of the IA and IF, which is defined as follows during 1 s.

$$s\left(t\right)=A\left(t\right)e^{j2\pi \phi \left(t\right)}$$

(24)

with

$$\begin{array}{c}A\left(t\right)=\exp \left[0.2\sin \left(6\pi t\right)\right]\\ \phi \left(t\right){=32t+200t}^2{-36t}^3{-50t}^4.\end{array}$$

(25)

Both the original signal and noisy signal contaminated by Gaussian white noise with zero mean are shown in Figure 1, where the signal-to-noise rate (SNR) is 0 dB and the sampling frequency is 2048 Hz. Then, the IFs estimated by the PIFEM in the different iterations are shown in Figure 2, where the parameter N is set as 8 and the threshold $\epsilon$ is set as 0.1 for the PIFEM. Note that the MAPE defined in Equation (10) can be adapted to calculate the relative error between the real and estimated IF by using the real and estimated IF to replace ${\widehat{\mathit{IF}}}_i$ and ${\widehat{\mathit{IF}}}_{i-1}$, respectively.

It can be found that, in Figure 2a, compared with the real IF, the IF estimated by the first iteration of the PIFEM deviates significantly from the real IF because of the strong noise but still captures the trend of the real IF. After the second iteration shown in Figure 2b, the difference between the estimated IF and real IF has narrowed considerably, where the MAPE is quickly reduced to 2.44% because the band-pass filtering operation effectively filters out most of the noise. After the third iteration shown in Figure 2c, the estimated IF has been very close to the real IF, where the MAPE is 0.45%. Finally, as shown in Figure 2d, the break condition of the PIFEM algorithm is triggered after the fifth iteration, where the MAPE is 0.18%. It indicates that the proposed iterative procedure can accurately estimate the nonlinear IF for the mono-component nonstationary signal with noise.

In addition, the MAPEs under different SNRs are shown in Figure 3 to test the procedure’s reliability and practical limits, where each MAPE is the average of 100 estimations. It can be found that the estimation accuracy decreases correspondingly along with reducing of SNR, and the estimation accuracy has a rapid decline when the SNR is less than 0 dB. Thus, the PIFEM is not recommended for application in extremely noisy cases.

The TFRs generated by different TFA methods are shown in Figure 4, where the IF is estimated by finding the maximum of TFR at every time instant. It can be found that the IF estimated by the STFT in Figure 4b is not accurate enough because of the noise, where the MAPE is 1.48%. The energy concentration of the TFR generated by the improved multisynchrosqueezing (IMSST) [19] in Figure 4c is significantly higher than the STFT in Figure 4a, where the IMSST is an improved SST-based TFA postprocessing method, but the IF estimated by the IMSST is still inaccurate compared with the real IF in Figure 4d, where the MAPE is 1.55%. The PCT [22] shown in Figure 4e cannot enhance the energy concentration of the TFR but can obtain a more accurate IF estimate, where the MAPE is 0.73% in Figure 4f. This indicates that the performance of the PIFEM is better than these different TFA methods. Note that these TFA methods are all based on the STFT that needs to calculate the FFT at each time instant, whose computational complexities are higher than the PIFEM.

Then, the PIFEM is also compared with the other two current state-of-art procedures for the IF estimation, i.e., the first conditional spectral moment of the time-frequency distribution of the signal [25] and the derivative of the phase of the analytic form of the signal obtained by the Hilbert transform, which is shown in Figure 5. It can be found that, compared with the PIFEM, the IF estimated by the first time-frequency moment is inaccurate under the strong noise. The derivative of the phase based on the Hilbert transform is extremely sensitive to noise, and the estimated IF is mightily oscillating when the SNR is equal to just 20 dB in Figure 5b.

Then, let us consider a multi-component signal defined in Equation (11), where K = 5, $a_k=1/k$, $p_k=k$, $A_0\left(t\right)$, and ${\phi }_0\left(t\right)$ are defined in Equation (25).

For the signal without noise, the IF estimated by the PIFEM is shown in Figure 6a, where the coefficient P in Equation (16) should be 1.56 according to the $a_k$ and $p_k$ defined above, and the parameter N is set as 8 for the PIFEM. It can be found that the IF estimated by PIFEM exactly approximates the 1.56 multiples of the real IF of the fundamental harmonic component, where the MAPE is 0.53%, which verifies the effectiveness of our proposed method. For the signal with noise, the iteration procedure of the PIFEM cannot be used because the frequency demodulation operation cannot change all the harmonic components into the same narrow band signal and it is difficult to design the band-pass filter for de-noising. However, as shown in Figure 6b, when the intensity of the Gaussian white noise is not very strong (SNR = 10 dB), the difference between the IF estimated by the PIFEM and the real IF is slight.

Then, by resampling the signal in the angular domain according to the estimated IF, we can effectively separate all harmonic components. As shown in Figure 7a, before resampling the signal, the amplitude spectrums of all five components are overlapped with each other. However, in Figure 7b, the amplitude spectrums of all five components are well separated after resampling in the angular domain, and then it is easy to use a band-pass filter to pick out every component, where the passband range of the filter is from 75 Hz to 225 Hz. Thus, we can use the iteration procedure to estimate the IF of every separated component. As shown in Figure 8, the estimated IF of the separated fundamental harmonic component is very close to the real IF, where the MAPE is 0.13%. This indicates that the PIFEM can be used for the multi-component nonstationary signal constructed by multiple harmonics.

4. Experiment and Results

In this part, a vibration signal collected from a rotor test rig under the time-varying rotating speed condition is considered. The rotor test rig is shown in Figure 9. The rotor is driven by a DC motor with a rated current of 1.95 A and a maximum output power of 148 W, and the motor is controlled by a speed controller (DH5600) that changes the 220 V AC power into the PWM signal for motor speed control. The rotor with a diameter of 10 mm is supported by two bearing brackets, connected with the motor shaft by a coupling, and contains two mass disks with a weight of 500 g and a diameter of 75 mm. A vibration signal is acquired by an eddy current sensor (5E102) fixed on the sensor bracket; transmitted to a dynamic data acquisition instrument (DH5922) for amplification and filter with a sampling frequency of 2000 Hz and a sampling duration of 10 s under the speed-up and speed-down process; and finally passed to the computer for analysis and storage, which is shown in Figure 10a and its amplitude spectrum is shown in Figure 10b. The eddy current sensor [26] uses the principle of eddy current formation to sense displacement. It can be used to measure shaft displacement in rotating machinery and has been around for many years, offering manufacturers high-linearity, high-speed measurements, and high resolution.

A change in the rotating speed of the rotor causes a change in the amplitude and frequency of the rotor vibration. The acquired vibration signal is consequently a multi-component nonstationary signal with multiple AM-FM harmonic components. However, it can be found that in Figure 10b, different harmonic components are not easy to identify because of the blurry spectrum caused by the nonlinear frequency modulation.

Then, the PIFEM (N = 10) is adopted to estimate the IF of the acquired multi-component vibration signal, which is shown as a red dashed line in Figure 11a compared with the TFR generated by the STFT. It can be found that the estimated IF is close to the IF of the fundamental harmonic component but not perfectly matched. Subsequently, the acquired multi-component signal is resampled according to the estimated IF, and the amplitude spectrum of the resampled signal is shown in Figure 11b, where the passband range of the filter is from 10 Hz to 35 Hz for picking out the fundamental harmonic component. It can be found that all harmonic components are clearer to identify and separate in Figure 11b compared with Figure 10b, but the weak high-order harmonics are still a little blurry.

Thus, the iteration procedure of the PIFEM (N = 10, $\mathsf{\epsilon }$ = 0.1) can be used to estimate the IF of the separated fundamental harmonic component, which denotes the instantaneous rotating speed of the rotor system. As shown in Figure 12a, the estimated IF of the fundamental harmonic is closer to the real IF compared with Figure 11a after the third iteration. The amplitude spectrum of the resampled signal is shown in Figure 12b. It can be found that all harmonic components in Figure 12b are discernible compared with Figure 10b, and the weak high-order harmonics are clearer compared with Figure 11b, which means the estimated IF is more accurate.

This indicates that the proposed method can be used for the IF estimation of the vibration signal from the varying-speed rotor system.

5. Discussion

In this section, we mainly focus on the unsolved shortcomings and future development of the PIFEM.

Firstly, a sensitivity analysis of the effects of the parameters N and $\epsilon$, i.e., the order of the polynomial and the error threshold, is shown in Figure 13, where the signal is defined in Equations (24) and (25) and the SNR is set as 10 dB. In Figure 13a, the error threshold $\epsilon$ is set as 0.1; it can be found that the estimation error increases as the order of the polynomial increases because the high-order coefficient is more sensitive to small changes in the signal than the low-order coefficient. However, the choice of the order also depends on the nonlinearity degree of the IF. Stronger nonlinearity means a higher order is required. In Figure 13b, the order of the polynomial is set as 8; it can be found that the estimation error is insensitive to changes in the threshold. However, a smaller threshold means more iterations are required. It is thus unnecessary to set the error threshold as too small.

Then, in Figure 14, the proposed method is evaluated in the case of abrupt frequency changes by a piecewise signal defined as follows:

$$s_{pw}\left(t\right)=\{\begin{array}{cc}s\left(t\right),& 0\le t<0.5\\ s\left(t\right)e^{j100\pi t},& 0.5\le t<1\end{array}$$

(26)

where the $s\left(t\right)$ is defined in Equations (24) and (25).

It can be found that in Figure 14, the PIFEM cannot accurately estimate the abrupt frequency changes because the polynomial series expansion can only approximate the continuous function. It is a worthwhile further research effort work to extend the PIFEM to satisfy a much wider range of signals.

Finally, it is worth noting that not only the IF but also other important time-varying features are utilized for condition monitoring such as the closely-related IA and IP. In addition, recently, instantaneous spectral entropy has been proposed as well, with applications for both the condition monitoring of wind turbines [27] and the structural health monitoring of multi-storey frame structures [28]. Hence, it will be further research to combine the PIFEM and these time-varying features for fault diagnosis and condition monitoring in the future.

6. Conclusions

The vibration signal from the rotating machinery under time-varying speed is usually a multi-component nonstationary signal constructed by AM-FM harmonic components, which makes it difficult to analyze the signal in the frequency domain. To efficiently and accurately estimate the nonlinear IF law of the nonstationary signal, a novel parameterized TFA method was proposed, known as the PIFEM. The proposed method has high computational efficiency because it can obtain the estimation of IF over a period of time simultaneously, which just needs to carry out several FFT operations and once least squares estimation. In addition, an iteration procedure was developed to enhance the noise robustness for the mono-component nonstationary signal. The numerical and experimental validations show that the proposed method is suitable to analyze the vibration signal from the rotating machinery, such as the IF estimate of the rotor system.