Novel Method for Vibration Sensor-Based Instantaneous Defect Frequency Estimation for Rolling Bearings Under Non-Stationary Conditions

It is proposed a novel instantaneous frequency estimation technology, multi-generalized demodulation transform, for non-stationary signals, whose true time variations of instantaneous frequencies are unknown and difficult to extract from the time-frequency representation due to essentially noisy environment. Theoretical bases of the novel instantaneous frequency estimation technology are created. The main innovations are summarized as: (a) novel instantaneous frequency estimation technology, multi-generalized demodulation transform, is proposed, (b) novel instantaneous frequency estimation results, obtained by simulation, for four types of amplitude and frequency modulated non-stationary single and multicomponent signals under strong background noise (signal to noise ratio is −5 dB), and (c) novel experimental instantaneous frequency estimation results for defect frequency of rolling bearings for multiple defect frequency harmonics, using the proposed technology in non-stationary conditions and in conditions of different levels of noise interference, including a strong noise interference. Quantitative instantaneous frequency estimation errors are employed to evaluate performance of the proposed IF estimation technology. Simulation and experimental estimation results show high effectiveness of the proposed estimation technology.


Introduction
Non-stationary signals occur in many fields, such as radar system, sonar, speech and mechanical systems [1]. Effective nonstationary signal processing methods are important for industrial applications. In recent years, scholars have developed some effective methods to characterize time-varying signal properties.
Time-frequency analysis (TFA) [2] is an effective method to display time-varying instantaneous frequencies (IFs) and time-varying instantaneous amplitudes of nonstationary signal. The traditional TFA methods include linear TFA such as the short-time Fourier transform (STFT), the wavelet transform (WT), and the bilinear/quadratic TFA such as the Wigner Ville distribution (WVD). However, these methods are restricted by the Heisenberg uncertainty principle and cross terms. Time frequency representation (TFR) with insufficient time-frequency resolutions are not providing accurate estimations of nonstationary signal, which essentially restrict their engineering applications. For improving the time-frequency resolutions, the time-frequency reassignment methods are developed, e.g., the synchrosqueezing transform (SST), demodulated SST and synchro-extracting transform Cheng [32] improved the GDT, and it can be used for multi-component signal decomposition. The iterative GDT is developed to transform multi-frequency components [33]. In the iterative GDT, the true IF curve functions are required, and the multi-frequency curves are transformed into the constant frequency curves. In recent years, the GDT is widely applied in industrial applications, such as fault frequencies detection of bearings and planetary gearboxes under nonstationary conditions [34]. Cheng [35] proposed an envelope order spectrum-based GDT to detect gear fault frequencies. In the proposed method, time-varying fault frequencies are transformed into lines, and are also decomposed into some mono-component. However, the demodulation operator is difficult to obtain without a speed sensor. For example, if the signal component with IF, which is used to calculate the demodulation operator, is polluted by noise, it therefore cannot be extracted from the TFR of the original signal.
Based on the above analysis, the IF estimation technique is important for machineries' fault frequencies detection under nonstationary conditions. The tachometer-free-based methods have attracted scholars' attention. Feng proposed [36] an iterative GDT-based energy separation method; the method can be used to process rotational machinery signals with knowledge of true IF dependency. In this method, the GDT is used to separate an interested harmonic, and the harmonic frequency can be extracted. The GDT is applied only one time to the same specific harmonic, and the number of harmonics in a signal defines the number of GDT iterations. These methods are facing the following challenge. In addition to the time-varying IFs, nonstationary signals normally exhibit time-varying instantaneous amplitudes. As a result, frequency components with relatively low amplitudes are easily contaminated by noise, and they therefore cannot be accurately extracted from the TFA by the GDT. The above GDT methods do not consider IFs of components with relatively low amplitudes that are difficult to accurately estimate.
Another challenge is as follows. The GDT transforms a time-varying frequency curve into a constant frequency one, which can be quantitatively characterized by a peak in the spectrum. The transformed frequency curve has a better time-frequency concentration. However, without knowledge of a true IF curve function, the GDT does not produce IF estimations.
In order to address the above-mentioned challenges, the novel IF estimation technology, named multi-generalized demodulation transform (MGDT), is proposed here. The main innovations of the paper are as follows:

1.
A novel IF estimation technology is proposed for IF estimation. The proposed novel technology, MGDT, is based on the GDT and employs the novel proposition of multiple adaptive iterations of the GDT for a specific signal component to concentrate blurry TFR energy in a step-wise manner. A theoretical basis of the novel IF estimation technology is created. A novel termination criterion, that could control iterations of the MGDT method, is also proposed. By iteratively and adaptively applying multiple GDT operations to the same specific signal component, the energy of the TFR result should be concentrated in a step-wise manner. The technology does not require a priori knowledge of true IF time dependency in contrast to other advanced IF estimation technologies. This feature is important for industrial applications of the proposed technology. The main difference between the proposed MGDT and previously proposed transforms, e.g., [36], is that the MGDT is performing multiple adaptive GDT iterations for the same specific signal component and the previously proposed transforms are performing just one (initial) GDT iteration for the same specific signal component. Other differences are highlighted below. The important feature of the novel MGDT is that every new iteration of the GDT is adaptive, i.e., it is based on previous iteration results. Our extensive literature search clearly shows that nobody, in worldwide terms, has previously proposed the main important novelty of the MGDT: an adaptive iterative application of the GDT for the same signal component of a multi-component signal.

2.
Novel IF estimation results, obtained by simulation, for four types of AM-FM nonstationary signals under a strong background noise (signal to noise ratio is -5dB): (i) single-component signal with nonstationary amplitude and nonstationary, nonlinear (quadratic) variation of the instantaneous frequency; (ii) multi-component signal with components that have nonstationary amplitudes and nonstationary linear instantaneous frequencies with different chirp rates; (iii) multi-component signal with components that have nonstationary amplitudes and nonstationary, nonlinear (quadratic) instantaneous frequencies with different chirp rates and frequency accelerations; and (iv) multi-component signal with components that have nonstationary amplitudes and nonstationary, nonlinear (cubic) instantaneous frequencies with different chirp rates, frequency accelerations, and the first derivatives of frequency accelerations. 3.
Novel experimental IF estimation results for defect frequency of rolling bearings for multiple defect frequency harmonics, using the proposed technology in non-stationary conditions and in conditions of six different levels of a noise interference, including of two levels of a strong noise interference.
The objectives of the research are as follows: 1.
To develop, investigate, and create the theoretical basis of the novel IF estimation technology.

2.
To perform a novel validation of the proposed technology by simulation for four types of AM-FM nonstationary signals.

3.
To perform comparisons of the proposed IF estimation technology with four traditional TFAs technologies by simulation.

4.
To perform a novel experimental validation of the proposed technology for IF estimation for multiple harmonics of bearing defect frequency in non-stationary conditions and in conditions of different levels of noise interference, including a strong noise interference.
The rest of the paper is structured as follows. Section 2 presents the developed theoretical bases of the IF estimation technology, including brief review of the GDT, and details of the proposed technology. Section 3 presents novel validation of the proposed method, using four types of the simulated signals, and comparisons of the proposed IF estimation technology with four traditional TFAs technologies by simulation. Section 4 presents a novel experimental validation of the proposed technology for IF estimation for multiple harmonics of a bearing defect frequency in nonstationary conditions and in conditions of different levels of noise interference, including a strong noise interference. Finally, the conclusions are drawn in Section 5.

Theoretical Bases of the Proposed IF Estimation Technology
In this section, we, firstly, briefly review the GDT technology, and then the novel IF estimation technique, the MGDT, is presented in details.

The GDT Technology
The nonstationary signals, that are AM-FM functions, can be defined as: where A(t) > 0 is the instantaneous amplitude, f (s) > 0 is the IF of the nonstationary signal. ∅ denotes the initial phase. In general, the measured signal consists of multiple frequency components, and thus can be expressed by where K is the number of the frequency components.
Sensors 2020, 20, 5201 6 of 30 The analytic signal of signal (1) should be calculated as: where H(·) represents the Hilbert transform. Then, the demodulation operator and the modulation operator can be expressed by Equations (4) and (5), respectively.
The demodulation operator is used for demodulation transform and the modulation operator is used for inverse demodulation transform. The demodulated signal can be calculated by multiplying Let f d (t) = f (t), the demodulated signal will become a purely AM signal centered around the frequency f c . In the GDT technology, the frequency f c is equal to the starting frequency of frequency function f (t).
Similarly, the analytic signal x A (t) can be calculated by multiplying

Novel IF Estimation Technique
The main concept of the IF estimation technique, the MGDT, is to obtain demodulation operator Φ − (t) by multiple iterations and step-wise refinement, and it is introduced, based on the mono-component analytic signal, shown in Equation (3).
Before the demodulation, due to the modulation influence, the spectrum of the signal x A (t) widely spreads around its IF f (t), as shown in Figure 1. Based on the GDT, with the known demodulation operator Φ − (t), the time-varying IF curve can be transformed into a constant frequency f c labelled as a blue line in Figure 1. As a result, its bandwidth is narrow with high energy concentration. Hence, the spectrum of the demodulated signal is more concentrated, than the original spectrum.
Based on relatively higher energy concentration of the demodulated signal, the multiple iteration and stepwise refinement strategies are proposed as follows: (a) roughly estimate the IFs of the original signal from its TFR; (b) calculate the demodulation operator, based on the estimated IF function and roughly demodulate the raw signal; as a result, the bandwidth of the IF curve becomes narrow and time-frequency concentration of the demodulated signal is enhanced and (c) obtain the high-accuracy IFs estimation by the peak searching algorithm. The implementation procedure of the MGDT is presented in Figure 2.
iteration and stepwise refinement strategies are proposed as follows: (a) roughly estimate the IFs of the original signal from its TFR; (b) calculate the demodulation operator, based on the estimated IF function and roughly demodulate the raw signal; as a result, the bandwidth of the IF curve becomes narrow and time-frequency concentration of the demodulated signal is enhanced and (c) obtain the high-accuracy IFs estimation by the peak searching algorithm. The implementation procedure of the MGDT is presented in Figure 2. The specifics of the implementation are described in detail as follows 1.
IFs and demodulation operator extraction The rough IFs are firstly extracted from the TFR of the raw signal. The TFR is calculated, based on the STFT. The STFT of the signal x A (t) ∈ L 2 (R) with the window function g(t) ∈ L 2 (R) can be defined as The time-frequency resolution of the TFR can be adjusted by the window width, and a wide window results in better frequency resolution, but poor time resolution. An appropriate window width L should be determined for obtaining an accurate result.
Then, peak searching algorithm is employed to extract the IFs from the TFR. For improving the robustness against a noise, local maxima detection strategy is developed. The peak searching algorithm is defined as The rough IFs are firstly extracted from the TFR of the raw signal. The TFR is calculated, based on the STFT. The STFT of the signal ( ) ∈ (ℝ) with the window function ( ) ∈ (ℝ) can be defined as The time-frequency resolution of the TFR can be adjusted by the window width, and a wide window results in better frequency resolution, but poor time resolution. An appropriate window width L should be determined for obtaining an accurate result.
Then, peak searching algorithm is employed to extract the IFs from the TFR. For improving the robustness against a noise, local maxima detection strategy is developed. The peak searching algorithm is defined as where arg max represents maximize function; ∆ is the local search range, ( , ) is position of peak in the whole time-frequency amplitudes; ( , ) is the output result, i = 1, 2, ⋯, I. I is determined by the window width L and the signal length. For fitting the discrete IFs, there are two function models. The first function model is the polynomial function model, which can be defined as Based on the Weierstrass approximation theorem [37], every continuous function, defined on a closed interval can be uniformly approximated, as closely as need, by a polynomial function. This model is widely used for IF estimation of the chirp signals [38].
The second function model is the Fourier series model, which is defined as The Fourier series is an expansion of a periodic function in terms of a sum of sines and cosines, and it can break up and recombine an arbitrary periodic function with whatever needed or with a practical accuracy [39]. This model also can be used for IF estimation of the chirp signals [40].
In this paper, the polynomial function model is used to fit discrete IFs. Hence, the first IF curve function is ( ) = ∑ . The first demodulation and modulation operators can be calculated as where arg max represents maximize function; ∆p is the local search range, IF(t s , f s ) is position of peak in the whole time-frequency amplitudes; IF(t i , f i ) is the output result, i = 1, 2, · · · , I. I is determined by the window width L and the signal length. For fitting the discrete IFs, there are two function models. The first function model is the polynomial function model, which can be defined as Based on the Weierstrass approximation theorem [37], every continuous function, defined on a closed interval can be uniformly approximated, as closely as need, by a polynomial function. This model is widely used for IF estimation of the chirp signals [38].
The second function model is the Fourier series model, which is defined as The Fourier series is an expansion of a periodic function in terms of a sum of sines and cosines, and it can break up and recombine an arbitrary periodic function with whatever needed or with a practical accuracy [39]. This model also can be used for IF estimation of the chirp signals [40]. In this paper, the polynomial function model is used to fit discrete IFs. Hence, the first IF curve function is f 1 (t) = K k=0 a d1 k t k . The first demodulation and modulation operators can be calculated as 2.
Demodulation and IF accurate calculation The first demodulation operator Φ − 1 (t) is used to demodulate the analytic signal x A (t). The first demodulated signal x d1 A (t) is calculated as The first IF curve function f 1 (t) is not identical to the actual IF curve function f (t) of the analytic signal. Hence, the spectrum of the demodulated signal x d1 A (t) is not centered around the carrier frequency a d1 0 , but has narrower frequency band than the original spectrum. The second demodulation operator Φ − 2 (t) is calculated for further demodulation, where Φ − 2 (t) is obtained by the second IF function f 2 (t) K k=0 a d2 k t k , which is extracted from the TFR S 2 (u, ω) of the demodulated signal x d1 A (t). With the iteration procedure, the actual IF function f m (t) is calculated as Furthermore, the actual demodulation and modulation operators can be calculated as Hence, the MGDT result of the analytic signal x A (t) can be shown as

Termination criterion and quantitative evaluation
The iteration can stop if an obvious peak can be captured in the spectrum of the demodulated signal and the biggest peak frequency PF m is approximately equal to the starting frequency a dm 0 of the IF function f m (t). Namely, PF m − a dm 0 /a dm 0 ≤ δ. δ is the termination threshold. In this paper, δ = 0.01. The spectrum of the MGDT signal x dm A (t) is obtained, using the Fourier transform. The quantitative indicators of the IF estimation error and transform error are defined to evaluate performance of the MGDT.
The IF estimation error FE is defined as Sensors 2020, 20, 5201 9 of 30 The transform error TE is defined as where PF is peak frequency in the spectrum of the demodulated signal with theoretical demodulation operator; PF m denotes the biggest peak frequency in the spectrum of the demodulated signal x dm A (t). The sensitivity of the MGDT is determined by both the initialization parameters in Figure 2 and noise level. The noise level is the most important parameter. If the signal to noise ratio (SNR) is relatively high, the IFs can be calculated under few iterations. On the contrary, more iterations should be considered under a low SNR.

Single-Component Non-Stationary Signal Analysis with Nonlinear (Quadratic) FM
For evaluating the performance of the MGDT, a simulated signal, consisting of one frequency component and the Gaussian white noise n(t) is constructed as The IF function of the simulated signal is ( ) = −108 + 454 + 260; the sampling frequency is 2000 Hz and the time duration is 2s. It is noted, that this function is not used in the proposed technology; this function is just used for IF error estimation. The SNR is −5dB, which is calculated as where and denote power and standard deviation of the noiseless signal, respectively; and are power and standard deviation of the noise, respectively. Figure 3 shows the simulated signal with a noise, its IF curve and its instantaneous amplitude curve. The instantaneous amplitude of the simulated signal is time-varying, as shown in Figure 3b, i.e. the instantaneous amplitude is relatively low in the initial part of the signal and is relatively high in the final part of the signal. The analytic signal obtained by using the Hilbert transform, is presented in Figure 4a. TFR of the analytic signal, generated by the STFT, and its local zooms are shown in Figure 4b-d, respectively. From the local zooms it can be seen, that the time-frequency ridge in the initial part is much blurry, than that in the final part. The main reason is that the noise pollutes the

Single-Component Non-Stationary Signal Analysis with Nonlinear (Quadratic) FM
For evaluating the performance of the MGDT, a simulated signal, consisting of one frequency component and the Gaussian white noise n(t) is constructed as The IF function of the simulated signal is f (t) = −108t 2 + 454t + 260; the sampling frequency is 2000 Hz and the time duration is 2 s. It is noted, that this function is not used in the proposed technology; this function is just used for IF error estimation. The SNR is −5 dB, which is calculated as where p x and A x denote power and standard deviation of the noiseless signal, respectively; p n and A x are power and standard deviation of the noise, respectively. Figure 3 shows the simulated signal with a noise, its IF curve and its instantaneous amplitude curve. The instantaneous amplitude of the simulated signal is time-varying, as shown in Figure 3b, i.e., the instantaneous amplitude is relatively low in the initial part of the signal and is relatively high in the final part of the signal. The analytic signal obtained by using the Hilbert transform, is presented in Figure 4a. TFR of the analytic signal, generated by the STFT, and its local zooms are shown in Figure 4b-d, respectively. From the local zooms it can be seen, that the time-frequency ridge in the initial part is much blurry, than that in the final part. The main reason is that the noise pollutes the weak part of the signal.   The IF estimation technology, the MGDT, is employed to the envelope signal for calculating the IF. Figure 5 shows the first iteration results. Here, the length of the Gaussian window is L = 128, 50% overlap, the local search range is ∆ = 186Hz. The extracted IFs, fitted IF curve ( ( ) = 2.96 + 170.64 + 422.73) and the actual IF curve are presented in Figure 5a. It can be found, that due to the inaccurate extraction of IFs in the initial part of the signal, the fitted IF curve has a clear error, which is 62.6%, calculated by Equation (20). Based on Equation (12)     The IF estimation technology, the MGDT, is employed to the envelope signal for calculating the IF. Figure 5 shows the first iteration results. Here, the length of the Gaussian window is L = 128, 50% overlap, the local search range is ∆ = 186Hz. The extracted IFs, fitted IF curve ( ( ) = 2.96 + 170.64 + 422.73) and the actual IF curve are presented in Figure 5a. It can be found, that due to the inaccurate extraction of IFs in the initial part of the signal, the fitted IF curve has a clear error, which is 62.6%, calculated by Equation (20). Based on Equation (12) Figure 5b, and its TFR is presented in Figure 5c. The spectrum of the demodulated signal,  The IF estimation technology, the MGDT, is employed to the envelope signal for calculating the IF. Figure 5 shows the first iteration results. Here, the length of the Gaussian window is L = 128, 50% overlap, and the local search range is ∆p = 186 Hz. The extracted IFs, fitted IF curve ( f 1 (t) = 2.96t 2 + 170.64t + 422.73) and the actual IF curve are presented in Figure 5a. It can be found, that due to the inaccurate extraction of IFs in the initial part of the signal, the fitted IF curve has a clear error, which is 62.6%, calculated by Equation (20). Based on Equation (12), 73t is calculated. The first demodulated signal is shown in Figure 5b, and its TFR is presented in Figure 5c. The spectrum of the demodulated signal, generated by the Fourier transform, is shown in Figure 5d, from which it can be found, that there is no clear no-smeared peak. Figure 5 also shows the traditional GDT results. It can be found that the traditional GDT is not performing well for this case, i.e., some instantaneous amplitudes are very low and are polluted by a noise.
Sensors 2020, 20, x 11 of 31 generated by the Fourier transform, is shown in Figure 5d, from which it can be found, that there is no clear no-smeared peak. Figure 5 also shows the traditional GDT results. It can be found, that the traditional GDT is not performing well for this case, i.e. some instantaneous amplitudes are very low and are polluted by a noise.  Figure 6b. Its TFR and spectrum are shown in Figure 6c,d, respectively. There is a constant timefrequency ridge in the TFR. A clear non-smeared peak, whose frequency is 262.2 Hz, can be seen in the spectrum. Based on the termination criterion, i.e. |262.2 − 262.4|/262.4 ≤ 0.01, the iteration process is stopped. Let continue to the second iteration. The parameters of window function L and peak search range ∆p are the same, as the parameters in the first iteration procedure. The new IFs are extracted from the TFR of the first demodulated signal and they are plotted in Figure 6a. The fitted IF curve ( f 2 (t) = −109.1t 2 + 279.1t + 262.4) is calculated using Equation (10), and it is presented as blue line in Figure 6a. It could be seen a close coincidence between the extracted IFs and the fitted IF curve. Then, the new demodulation operator is calculated With the new demodulation operator, the second demodulated signal is obtained, as shown in Figure 6b. Its TFR and spectrum are shown in Figure 6c,d, respectively. There is a constant time-frequency ridge in the TFR. A clear non-smeared peak, whose frequency is 262.2 Hz, can be seen in the spectrum. Based on the termination criterion, i.e., |262.2 − 262.4|/262.4 ≤ 0.01, the iteration process is stopped.
The actual IF function f 2 (t) = −106.1t 2 + 449.7t + 262.4 can be calculated based on Equation (15). Figure 7a shows that the calculated IF curve (blue line) coincides with the actual IF curve (red line), and the error is 0.85%. Then, the actual demodulation operator 4t is calculated based on Equation (16). With the actual demodulation operator, the demodulated signal x d A (t) is calculated, as presented in Figure 7b. Its TFR and the spectrum are presented in Figure 7c,d, respectively. As shown in the TFR and the spectrum, the time-varying time-frequency ridge is transformed into a line, which can be quantitatively characterized. Its TFR and the spectrum are presented in Figure 7c,d, respectively. As shown in the TFR and the spectrum, the time-varying time-frequency ridge is transformed into a line, which can be quantitatively characterized.   Its TFR and the spectrum are presented in Figure 7c,d, respectively. As shown in the TFR and the spectrum, the time-varying time-frequency ridge is transformed into a line, which can be quantitatively characterized.  In addition, the frequency component can also be separated from noise, using the inverse MGDT algorithm for avoiding noise interference. The actual modulation operator Φ + 2 (t) can be determined, based on Equation (17). The TFR of the demodulated signal x d A (t) with one constant time-frequency ridge has well-localized zones of concentration, as shown in Figure 7c. Thus, the band-pass filter is employed to the demodulated signal with the central frequency 262.2 Hz. With the actual modulation operator Φ + 2 (t), the modulation signal x A f (t) of the filtered component , which is presented in Figure 8a. Figure 8b is the TFR of the modulation signal. Its local zoom is shown in Figure 8c, from which it can be found, that there are no interference components.
In addition, the frequency component can also be separated from noise, using the inverse MGDT algorithm for avoiding noise interference. The actual modulation operator ( ) can be determined, based on Equation (17). The TFR of the demodulated signal ( ) with one constant time-frequency ridge has well-localized zones of concentration, as shown in Figure 7c. Thus, the band-pass filter is employed to the demodulated signal with the central frequency 262.2 Hz. With the actual modulation operator ( ) , the modulation signal ( ) of the filtered component ( ) is calculated as ( ) = ( ) ( ), which is presented in Figure 8a. Figure 8b is the TFR of the modulation signal. Its local zoom is shown in Figure 8c, from which it can be found, that there are no interference components. For comparison of the MGDT with the traditional TFA methods, the GDT, WT, the SST and the SET are employed.
There are many different families of wavelet functions for various engineering applications. As the proposed MGDT is employing complex analytical signal, it requires a choice of complex wavelet type for comparison purpose [41]. In the present work, the complex Morlet mother wavelet function is selected due to the following main reasons: (i) the Morlet function is employing the Gaussian window and, therefore, is providing highly efficient balance between time and frequency resolutions [42] (ii) the Morlet wavelet is not sharp edge function; therefore, it minimizes the spectral leakage [43]; (iii) the Morlet wavelet allows selection of highly efficient balance between time-frequency resolutions by two parameters: the center frequency and the bandwidth parameter (iv) the wavelet transform with the Morlet function is computationally efficient due to its link to the fast Fourier transform. For the WT, the wavelet name is "cmor4-2", the central frequency is 2 Hz and bandwidth parameter is 4s 2 .
For the SST, the window length is 100. For the SET the window length is 120. The TFRs, calculated by the traditional methods, the WT, the SST, and the SET are shown in Figure 9a-c, respectively.
The IFs are extracted, using the peak search algorithm from the TFRs, and the extracted IFs are presented in Figure 9d. In addition, the actual IF and IF, estimated by the novel MGDT, are also presented in Figure 9d. All peak search parameters are the same, as applied to the proposed MGDT. From Figure 9d, it can be found, that the extracted IFs by GDT, the WT, the SST and the SET have relatively large errors in the initial duration of the signal. Their total estimation errors are 55.9%, 29.6%, 60.2%, and 113.8%, respectively. The IF estimation error of GDT, 55.9%, is obtained by use of For comparison of the MGDT with the traditional TFA methods, the GDT, WT, the SST and the SET are employed.
There are many different families of wavelet functions for various engineering applications. As the proposed MGDT is employing complex analytical signals, it requires a choice of complex wavelet types for comparison purposes [41]. In the present work, the complex Morlet mother wavelet function is selected due to the following main reasons: (i) the Morlet function employs the Gaussian window and, therefore, provides a highly efficient balance between time and frequency resolutions [42]; (ii) the Morlet wavelet is not a sharp edge function; therefore, it minimizes the spectral leakage [43]; (iii) the Morlet wavelet allows the selection of a highly efficient balance between time and frequency resolutions by two parameters: the center frequency and the bandwidth parameter; (iv) the wavelet transform with the Morlet function is computationally efficient due to its link to the fast Fourier transform. For the WT, the wavelet name is "cmor4-2", the central frequency is 2 Hz, and the bandwidth parameter is 4s 2 .
For the SST, the window length is 100. For the SET, the window length is 120. The TFRs, calculated by the traditional methods, the WT, the SST, and the SET are shown in Figure 9a-c, respectively.
The IFs are extracted, using the peak search algorithm from the TFRs, and the extracted IFs are presented in Figure 9d. In addition, the actual IF and IF, estimated by the novel MGDT, are also presented in Figure 9d. All peak search parameters are the same, as applied to the proposed MGDT. From Figure 9d, it can be found that the extracted IFs by the GDT, the WT, the SST, and the SET have relatively large errors in the initial duration of the signal. Their total estimation errors are 55.9%, 29.6%, 60.2%, and 113.8%, respectively. The IF estimation error of the GDT, 55.9%, is obtained by use of the traditional GDT just once for the raw signal. The IF estimation errors for different methods are listed in Table 1. It could be seen that the proposed method has a better ability to calculate IFs than the traditional TFA methods, for the simulated signal that features a strong nonlinear FM, a time-variable signal amplitude, and a very noisy environment, the signal to noise ratio is −5dB. It could be also seen from Figure 9d that the main contributions to the estimation errors of the WT, the SST, and the SET are in a very noisy area of a signal, characterized by a low signal amplitude, i.e., in the time slot of (0-0.6) s. the traditional GDT just one time for the raw signal. The IF estimation errors for different methods are listed in Table 1. It could be seen, that the proposed method has better ability to calculate IFs, than the traditional TFA methods, for the simulated signal, that is featuring of a strong nonlinear FM, a time variable signal amplitude and a very noisy environment, the signal to noise ratio is −5dB. It could be also seen from Figure 9d, that the main contributions to estimation errors of the WT, the SST and the SET are in a very noisy area of a signal, characterized by low signal amplitude: i.e. in the time slot of (0-0.6) s.
The IF functions of the three components of the simulated signal are ( ) = 90.8 + 260 , ( ) = 22.7 + 65, and ( ) = 4.54 + 13, respectively. It is noted, that these functions are not used in the proposed technology; these functions are just used for IF error estimations. The sampling frequency is 2000 Hz, signal time duration is 2s and the SNR is −5dB.

Multi-Component Nonstationary Signal with Linear FM
For further validation of the IF estimation technology, another multi-component signal with components that have time-varying amplitudes and linear FMs with different chirp rates is simulated on the background of a noise as The IF functions of the three components of the simulated signal are f a (t) = 90.8t + 260, f b (t) = 22.7t + 65, and f c (t) = 4.54t + 13, respectively. It is noted that these functions are not used in the proposed technology; these functions are just used for IF error estimations. The sampling frequency is 2000 Hz, signal time duration is 2s, and the SNR is −5 dB. Figure 10 shows the multi-component signal and the TFR of the envelope signal. From the TFR, it can be found that the time-frequency ridge in the initial part of the signal is much more blurry than that in the final part of the signal. Figure 10 shows the multi-component signal and the TFR of envelope signal. From the TFR, it can be found, that the time-frequency ridge in the initial part of the signal is much blurry, than that in the final part of the signal. It should be noted, that the amplitudes of components of the signal are different, as shown by Equation (23), and IF of the component with a higher amplitude is estimated more accurately.  The IF estimation technology, the MGDT, is employed to the envelope signal for calculating the IF with the highest amplitude. The extracted IF function f a (t) = 89.1t + 256.2 is obtained. Figure 11 shows that the extracted IF curve (blue line) coincides with the theoretical IF curve (red line), and the estimation error is 1.46%.

Multi-Component Non-Stationary Signal with Nonlinear(Quadratic) FM
Based on the estimated IF f a (t), the first component can be removed from the raw signal. Then, the MGDT is applied to the signal without the first component, and the IF of the second component is extracted as f b (t) = 22.2t + 66.8. The estimation error is 2.77%. Then, by removing the second component, the IF of the third component can be extracted as f c (t) = 4.65t + 13.58, whose estimation error is 4.46%. Hence, the proposed IF estimation technology could be also used for nonstationary multi-component signals with a time-varying amplitude and linear FM.
It should be noted that the amplitudes of components of the signal are different, as shown by Equation (23), and the IF of the component with a higher amplitude is estimated more accurately.
Sensors 2020, 20, x 15 of 31 Figure 10 shows the multi-component signal and the TFR of envelope signal. From the TFR, it can be found, that the time-frequency ridge in the initial part of the signal is much blurry, than that in the final part of the signal. It should be noted, that the amplitudes of components of the signal are different, as shown by Equation (23), and IF of the component with a higher amplitude is estimated more accurately.

Frequency(Hz)
IF curve Harmonic f a

Harmonic f b
Harmonic f c Calculated IF Figure 11. Calculated IF curves.

Multi-Component Nonstationary Signal with Nonlinear (Quadratic) FM
For further validation of the proposed IF estimation technology, another multi-component signal with components that have time-varying amplitudes and nonlinear (quadratic) FMs with different chirp rates and different frequency accelerations is simulated on the background of a noise as x(t) = 1 5 e 2t−4 sin 2π −7.2t 3 + 45.4t 2 + 260t +0.9 × 1 5 e 2t−4 sin 2π −14.4t 3 + 90.8t 2 + 520t + n(t)# The IF functions of the two components of the simulated signal are f a (t) = −21.6t 2 + 90.8t + 260 and f b (t) = −43.2t 2 + 181.6t + 520, respectively. It is noted that these functions are not used in the proposed technology; these functions are just used for IF error estimations. The sampling frequency is 2000 Hz and the time duration is 2s. The SNR is −5 dB. Figure 12 shows the multi-component signal and the TFR of the envelope signal. From the TFR, it can be found that the time-frequency ridge in the initial part of the signal is much more blurry than that in the final part of the signal. The IF functions of the two components of the simulated signal are ( ) = −21.6 + 90.8 + 260and ( ) = −43.2 + 181.6 + 520, respectively. It is noted, that these functions are not used in the proposed technology; these functions are just used for IF error estimations. The sampling frequency is 2000 Hz and the time duration is 2s. The SNR is −5dB. Figure 12 shows the multi-component signal and the TFR of envelope signal. From the TFR, it can be found that the time-frequency ridge in the initial part of the signal is much blurry than that in the final part of the signal.   The IF functions of the two components of the simulated signal are ( ) = −21.6 + 90.8 + 260and ( ) = −43.2 + 181.6 + 520, respectively. It is noted, that these functions are not used in the proposed technology; these functions are just used for IF error estimations. The sampling frequency is 2000 Hz and the time duration is 2s. The SNR is −5dB. Figure 12 shows the multi-component signal and the TFR of envelope signal. From the TFR, it can be found that the time-frequency ridge in the initial part of the signal is much blurry than that in the final part of the signal.

Multi-Component Nonstationary Signal with Nonlinear (Quadratic) FM
For further validation of the proposed IF estimation technology, another multi-component signal with components that have time-varying amplitudes and nonlinear (cubic) FMs with different chirp rates, different frequency accelerations, and different first derivatives of frequency acceleration is simulated on the background of a noise as x(t) = 1 5 e 2t−4 sin 2π t 4 − 7.2t 3 + 45.4t 2 + 260t +0.9 × 1 5 e 2t−4 sin 2π 2t 4 − 14.4t 3 + 90.8t 2 + 520t + n(t) The IF functions of the simulated signal are f a (t) = 4t 3 − 21.6t 2 + 90.8t + 260 and f b (t) = 8t 3 − 43.2t 2 + 181.6t + 520. It is noted that these functions are not used in the proposed technology; these functions are just used for IF error estimations. The sampling frequency is 2000 Hz and the time duration is 2 s. The SNR is −5 dB. Figure 14 shows the multi-component signal and the TFR of the envelope signal. From the TFR, it can be found that the time-frequency ridge in the initial part of the signal is much more blurry than that in the final part of the signal.
The IF functions of the simulated signal are ( ) = 4 − 21.6 + 90.8 + 260 and ( ) = 8 −43.2 + 181.6 + 520. It is noted, that these functions are not used in the proposed technology; these functions are just used for IF error estimations. The sampling frequency is 2000 Hz and the time duration is 2s. The SNR is −5dB. Figure 14 shows the multi-component signal and the TFR of envelope signal. From the TFR, it can be found that the time-frequency ridge in the initial part of the signal is much blurry than that in the final part of the signal. The IF estimation technology, the MGDT, is employed to the signal for estimating the IFs. Figure  15 shows that the estimated IF curves (blue lines) coincide with the theoretical IF curves (red lines) for both components. The estimation errors for the first component is 1.2% and for the second component is 1.6%.
Hence, it is validated, that the proposed IF estimation technology could be also used for multicomponent signals with time variable signal amplitude and nonlinear (cubic) frequency modulations with different chirp rates, different frequency accelerations and different first derivatives of frequency acceleration. The IF estimation technology, the MGDT, is employed to the signal for estimating the IFs. Figure 15 shows that the estimated IF curves (blue lines) coincide with the theoretical IF curves (red lines) for both components. The estimation error for the first component is 1.2% and for the second component it is 1.6%.
Hence, it is validated that the proposed IF estimation technology could be also used for multi-component signals with a time-variable signal amplitude and nonlinear (cubic) frequency modulations with different chirp rates, different frequency accelerations, and different first derivatives of frequency acceleration.
The IF functions of the simulated signal are ( ) = 4 − 21.6 + 90.8 + 260 and ( ) = 8 −43.2 + 181.6 + 520. It is noted, that these functions are not used in the proposed technology; these functions are just used for IF error estimations. The sampling frequency is 2000 Hz and the time duration is 2s. The SNR is −5dB. Figure 14 shows the multi-component signal and the TFR of envelope signal. From the TFR, it can be found that the time-frequency ridge in the initial part of the signal is much blurry than that in the final part of the signal. The IF estimation technology, the MGDT, is employed to the signal for estimating the IFs. Figure  15 shows that the estimated IF curves (blue lines) coincide with the theoretical IF curves (red lines) for both components. The estimation errors for the first component is 1.2% and for the second component is 1.6%.
Hence, it is validated, that the proposed IF estimation technology could be also used for multicomponent signals with time variable signal amplitude and nonlinear (cubic) frequency modulations with different chirp rates, different frequency accelerations and different first derivatives of frequency acceleration. The proposed method could be used for the effective online estimation of time-variable linear and nonlinear IFs. The method is designed to estimate IFs from vibration data automatically. Upon receiving, automatically, an IF estimation, related to the first iteration, a termination will be automatically used to determine if the method should continue to iterate or stop further iterations. If an estimation error is lower than a defined termination threshold, the method will stop working. Otherwise, the method will continue to iterate until an estimation error that is lower than a termination threshold is obtained.
The method is computationally fast and provides near real-time estimation results, e.g., the calculation time is 1.7 s for the IF estimation of 2-second vibration data, the frequency range of the STFT is 1000 Hz, the frequency resolution is 7.8 Hz, the overlapping for the STFT is 50%, the sampling frequency is 2000 Hz, and there are two iterations that are employed to achieve the required estimation accuracy.

Validation of the Proposed IF Estimation Technology via Experimental Trials
In this section, experimental validation of the proposed IF estimation technology for a nonstationary multi-component signal from the rolling element bearing is performed.
The rolling bearing is one of the most widely used mechanical parts in rotating machinery. Its faults are important reasons of shutdowns of machines. For reducing economic losses and avoiding disastrous accidents, novel IF estimations for rolling bearing has become an important topic [44][45][46][47].
The captured nonstationary bearing vibration signal from a faulty bearing is considered for experimental MGDT validation. The main purpose of the validation is to estimate the non-stationary instantaneous bearing defect frequencies for the fundamental harmonic and the higher harmonics, using the proposed IF estimation technology for different levels of noise interference at a very early stage of damage development. Figure 16 shows the structural sketch of the experimental set-up. The AC motor (Marathon, Beloit, WI, USA )(3 phase, Hz: 60/50, RPM: 3450/2850, voltage: 208-230/460) is used to drive the steel shaft, which is supported by two bearings (Rexnord, Milwaukee, WI, USA), including a faulty one.
The employed motor has a star connection. In delta connection, each winding will get phase to phase voltage. Alternatively, in star connection, there will be two windings in series across two phases. Hence, motor windings get a higher voltage in delta connection. In star connection, the stator phase voltage is about 58 % of the line voltage. The torque of the motor is directly proportional to the square of the applied voltage. In star connection, the torque of the motor is about 33 % of the maximum rated torque of the motor. The torque delivering capacity of the motor in the delta connection is more than the torque-delivering capacity in the star connection, i.e., torque (star connection) = 33 % of torque (delta connection). Motor vibration amplitudes are propositional to the motor torque. Therefore, it is easier to estimate IFs of narrowband motor components in the delta connection, compared with the star connection.
Two flywheels (SpectraQuest, Richmond, VA, USA) are installed on the shaft to apply load to the system. The RF is measured by a built-in tachometer (SpectraQuest, Richmond, VA, USA) (one pulse per revolution analogue transistor-transistor logic output for DAQ purposes) and the RF can be adjusted by the AC converter.
An accelerometer (PCB, Depew, NY, USA) (PCB 352A60) is mounted on the top surface of the faulty bearing housing by a tapped hole to capture the vibration signal. This miniature accelerometer achieves an extremely high frequency response from 5 Hz to 60 kHz. The sensitivity of the accelerometer is 10mV/g; electrical filter corner frequency is 45 kHz; electrical filter roll-off is 10 dB/decade; resonant frequency is ≥ 95 kHz; measurement range is up to 500 g; broadband resolution is 0.002 g rms; nonlinearity is ≤ 1%; transverse sensitivity is ≤ 5%; overload limit is 5000 g pk; temperature response is −54 to +121°C; base strain sensitivity is ≤ 0.05 g/µε; excitation voltage is 18 to 30 VDC; constant current excitation is 2 to 20 mA; output impedance is ≤ 100Ohm; output bias voltage is from 8 to 12 VDC; discharge time constant is from 0.02 to 0.06 s; spectral noise (10 Hz) is 160µg/ √ Hz ; spectral noise (100 Hz) is 40 µg/ √ Hz ; spectral noise (1 kHz) is 15 µg/ √ Hz ; spectral noise (10 kHz) is 10 µg/ √ Hz ; height is 21.6 mm; weight is 6 gm; sensing element is ceramic; size-hex is 3/8in; sensing geometry is shear; housing material is stainless steel; sealing is welded hermetic; electrical connector is 5-44 coaxial; electrical connection position is the top; mounting is integral stud; and mounting thread is 10-32 male. The accelerometer is shown in Figure 16. The signal collecting system is mainly consisted of a DAQ card (Model: NI-PCI6259, number of channels is 16 differential or 32 single-ended, ADC resolution is 16 bits, single channel maximum is 1.25 MS/s, and multichannel maximal sample rate is 1 MS/s) and the NI-DAQmx software.
The bearing type is ER-10K, number of balls is 8, ball diameter is 8.94 mm, pitch diameter is 33.5 mm, contact angle is 0, and its FCCs of the outer race, inner race, and ball are 3.05, 4.95, and 1.99, respectively. A defect in the form of a 0.1 mm diameter hole is considered in the outer raceway of the bearing; the relative damage size, related to the outer race circumference, is 0.2%. The bearing is shown in Figure 16. The signal collecting system is mainly consisted of a DAQ card (Model: NI-PCI6259, number of channels is 16 differential or 32 single ended, ADC resolution is 16 bits, single channel maximum is 1.25MS/s, multichannel maximal sample rate is 1 MS/s) and the NI-DAQmx software.
The bearing type is ER-10K, number of balls is 8, ball diameter is 8.94mm, pitch diameter is 33.5mm, contact angle is 0, and its FCCs of outer race, inner race, and ball are 3.05, 4.95, and 1.99, respectively. A defect in the form of 0.1 mm diameter hole is considered in the outer raceway of the bearing; the relative damage size, related to the outer race circumference, is 0.2%. The bearing is shown in Figure 16. First, raw bearing vibration signal is collected. The collected vibration signal from the bearing with outer race fault is presented in Figure 17a; the estimated from speed sensor RF function is ( ) = 2.9 + 19.9 Hz; the RF increases from 19.9 Hz to 31.5 Hz during 4s. The sampling frequency is 48,000 Hz and the time duration is 4s. It should be noted, that speed sensor data are not used for the proposed estimation technology, as the technology is not required a priori knowledge of a true time varying IF dependency. These data are used only for IF error estimations for the proposed technology.
For selection of the sampling frequency, the following frequencies should be taken into account: the fundamental outer race bearing defect frequency, that is varying from 60.7 Hz to 96.1 Hz and a frequency band, that is related to transient impulses, generated by bearing outer race local fault. Before the final selection of the sampling frequency, a preliminary assessment of a frequency band, that is related to transient impulses, generated by bearing local outer race fault, was performed by the spectral kurtosis [48][49][50], and it was found, that this band is (9-12) kHz. So, based on the estimated frequency band, it is clear, that the bearing outer race defect frequency is not critical in selection of the sampling frequency, and the estimated frequency band is critical in sampling frequency selection. In order to avoid aliasing, the sampling frequency is selected as four times of the maximum frequency of a frequency band, that is related to transient impulses, generated by bearing local faults, i.e. 48 kHz.
Firstly, the proposed technology is validated for IF estimation of the first harmonic of outer race defect frequency. The envelope signal, shown in Figure 17b, is estimated by the Hilbert transform. Figure 17c, d shows the TFR and the spectrum of the envelope signal, respectively. From the TFR, it can be found, that due to the noise, the weak IFCFs in the initial duration of the signal are polluted by noise. In addition, because of the time-varying operation condition, the IFCF curve is time-varying and cannot be quantitatively characterized by the classical Fourier spectrum, Figure 17d. First, the raw bearing vibration signal is collected. The collected vibration signal from the bearing with outer race fault is presented in Figure 17a; the estimated from speed sensor RF function is f (t) = 2.9t + 19.9 Hz; the RF increases from 19.9 to 31.5 Hz over 4 s. The sampling frequency is 48,000 Hz and the time duration is 4 s. It should be noted that speed sensor data are not used for the proposed estimation technology, as the technology does not require a priori knowledge of a true time-varying IF dependency. These data are used only for IF error estimations for the proposed technology.
For selection of the sampling frequency, the following frequencies should be taken into account: the fundamental outer race bearing defect frequency that varies from 60.7 to 96.1 Hz and a frequency band that is related to transient impulses, generated by a bearing outer race local fault. Before the final selection of the sampling frequency, a preliminary assessment of a frequency band that is related to transient impulses, generated by a bearing local outer race fault, was performed by spectral kurtosis [48][49][50], and it was found that this band is (9)(10)(11)(12) kHz. So, based on the estimated frequency band, it is clear that the bearing outer race defect frequency is not critical in the selection of the sampling frequency, and the estimated frequency band is critical in sampling frequency selection. In order to avoid aliasing, the sampling frequency is selected as four times the maximum frequency of a frequency band that is related to transient impulses, generated by bearing local faults, i.e., 48 kHz.
Firstly, the proposed technology is validated for IF estimation of the first harmonic of the outer race defect frequency. The envelope signal, shown in Figure 17b, is estimated by the Hilbert transform. Figure 17c,d show the TFR and the spectrum of the envelope signal, respectively. From the TFR, it can be found that due to the noise, the weak IFCFs in the initial duration of the signal are polluted by noise. In addition, because of the time-varying operation condition, the IFCF curve is time-varying and cannot be quantitatively characterized by the classical Fourier spectrum, as shown in Figure 17d The proposed IF estimation technology is employed to the envelope signal for obtaining the IFCF. Here, the length of the Gaussian window is L = 12000, 50% overlap, the local search range ∆ = 100Hz.
The main criterion for window length selection depends on an instantaneous frequency change, characterized by the chirp rate in the case of the linear frequency variation and by the chirp rate and the frequency acceleration in the case of quadratic frequency variation, and so on [51][52][53][54].
The employed criterion of window length selection is clearly explained below on the basis of the two main cases of combinations of instantaneous frequency change ( ) of components ( ) and frequency resolution, defined by a window length. Case 1. The instantaneous frequency change of IF of component ( ) during duration of the selected window is within frequency limits ( − /2 ; + /2 ), defined by the frequency resolution f at a particular central frequency of a frequency bin of the STFT. The considered case could be viewed as a quasi-stationary case within a window length and it is provided no additional errors, related to STFT estimation for non-stationary signal, because time variable IF is not moving from one frequency bin to another during the selected window length.
Case 2. The instantaneous frequency change of IF of component ( ) during duration of the selected window is more, than the frequency limits ( − /2; + /2), defined by the frequency resolution f at the particular central frequency of a frequency bin. This case is definitely not a quasi-stationary case and it is provided additional errors, related to STFT estimation for non-stationary signal, as time variable IF is moving from one frequency bin to another during the selected window length: i.e. IF will be, at least, in two different frequency bins (or even more, depending on speed of frequency change) for the selected duration of the window.
Based on the considered two cases, the clear criterion for window selection, that we employed, is as follows: the instantaneous frequency change of IF of component ( ) during duration of the selected window should be within the frequency limit ( − /2; + /2), defined by the frequency resolution with a "safety factor" within threshold range (1.5-3) for "safety factor", where a safety factor is defined as = / . The threshold range (1.5-3) is defined as follows: it is a danger to select a threshold around 1, as an IF still could move from one frequency bin to another during the selected window length. If the selected threshold range (1.5-3), then, it will be IF non-movement from one frequency bin to another during the selected window length. It is possible to have a safety factor more than 3; however, a The proposed IF estimation technology is employed to the envelope signal for obtaining the IFCF. Here, the length of the Gaussian window is L = 12000, 50% overlap, the local search range ∆p = 100 Hz.
The main criterion for window length selection depends on an instantaneous frequency change, characterized by the chirp rate in the case of linear frequency variation and by the chirp rate and the frequency acceleration in the case of quadratic frequency variation, and so on [51][52][53][54].
The employed criterion of window length selection is clearly explained below on the basis of the two main cases of combinations of the instantaneous frequency change f j (t) of the components j(t) and frequency resolution, defined by a window length. Case 1. The instantaneous frequency change f chj of IF of component f j (t) during the duration of the selected window is within the frequency limits ( f p − f /2; f p + f /2), defined by the frequency resolution f at a particular central frequency f p of a frequency bin of the STFT. The considered case could be viewed as a quasi-stationary case within a window length and no additional errors are provided that are related to STFT estimation for a nonstationary signal, because the time-variable IF is not moving from one frequency bin to another during the selected window length. Case 2. The instantaneous frequency change f chj of IF of component f j (t) during the duration of the selected window is more than the frequency limits ( f p − f /2; f p + f /2), defined by the frequency resolution f at the particular central frequency f p of a frequency bin.
This case is definitely not a quasi-stationary case and additional errors are provided that are related to STFT estimation for a nonstationary signal, as the time-variable IF is moving from one frequency bin to another during the selected window length, i.e., the IF will be, at least, in two different frequency bins (or even more, depending on the speed of the frequency change) for the selected duration of the window.
Based on the considered two cases, the clear criterion for window selection that we employed is as follows: the instantaneous frequency change f chj of IF of component f j (t) during the duration of the selected window should be within the frequency limit ( f p − f /2; f p + f /2), defined by the frequency resolution with a "safety factor" within the threshold range (1.5-3) for a "safety factor", where a safety factor is defined as S = f / f ch .
The threshold range (1.5-3) is defined as follows: it is a danger to select a threshold around 1, as an IF still could move from one frequency bin to another during the selected window length. If the selected threshold range is (1.5-3), then the IF will not move from one frequency bin to another during the selected window length. It is possible to have a safety factor of more than 3; however, a drawback of such a big safety factor is that the frequency resolution for this case becomes poorer and this will potentially compromise an IF estimation.
It is discussed here how this criterion was employed for the considered experimental trials. As the main effort is for estimation, firstly, the IF of the first harmonic of the bearing outer race defect frequency is considered a time variation of the IF of this harmonic, which is 3.05 × f (t), where f (t) = 2.9t + 19.9 f (t) is the time variation of the fundamental rotation frequency, and the coefficient of 3.05 is the FCC of the bearing outer race.
For the selected window length of 12000 (i.e., time duration of this window length is 0.25 s), the instantaneous frequency change f ch1 of the IF of the first harmonic of the bearing outer race defect frequency is f chj = (0.25s x 2.9 Hz/s x 3.05) = 2.21 Hz and the frequency resolution related to this window length is f = (1/0.25 Hz) = 4 Hz. So, for this selection of the window, the instantaneous frequency change (i.e., 2.21 Hz) during the duration of the window is within the frequency resolution (i.e., 4 Hz) and, therefore, it is possible to avoid additional errors related to STFT estimation for the first harmonic of the bearing outer race defect frequency. This avoidance is because the time-variable IF is not moving from one frequency bin to another during window length. The "safety factor" for this estimation is S = (4 Hz/2.21 Hz) = 1.81 and it is inside the specified threshold range (1.5-3).
It is possible to select a window length, e.g., of 6000 (i.e., the time duration of this window length is 0.125s). In this case, the instantaneous frequency change f ch1 of the IF of the first harmonic of the bearing outer race defect frequency is f chj = (0.125s x 2.9 Hz/s x 3.05) = 1.11 Hz and the frequency resolution related to this window length is f = (1/0.125 Hz) = 8 Hz. So, in this case, the instantaneous frequency change during the duration of the window (i.e., 1.11 Hz) is also within the frequency resolution (8 Hz) and it is also possible to avoid additional errors related to STFT estimation for the first harmonic of the bearing outer race defect frequency.
However, the "safety factor" for this estimation is too big at S = (8 Hz/1.11 Hz) = 7.21. A drawback of this big safety is that the frequency resolution for this case (8 Hz) is poorer, compared with window length of 12000 that, potentially, will compromise an IF estimation.
If a longer window length, e.g., 18000, is selected, then the instantaneous frequency change f chj of the IF of the first harmonic of the bearing outer race defect frequency is (0.3755s x 2.9 Hz/s x 3.05) = 3.32 Hz, and the frequency resolution related to this window length is f = (1/0.375 Hz) = 2.67 Hz. In this case, the instantaneous frequency change during the duration of the window (i.e. 3.32 Hz) is not within the frequency resolution (2.67 Hz) and it is not possible to avoid additional errors related to STFT estimation for the first harmonic of the bearing outer race defect frequency. This non-avoidance is because the time-variable IF is moving from one frequency bin to another during the window length. So, the selection of a longer window length of 18000 is not acceptable.
The same criterion was used for the selection of the window length for simulation trials. For the simulation of two harmonics with nonlinear frequency modulation, the selected window length is 128. For the selected window length of 128 (i.e., the time duration of this window length is 0.064 s), the instantaneous frequency change f ch1 of the IF of the first harmonic is 5.72 Hz and the frequency resolution related to this window length is f = (1/0.064 Hz) = 15.625 Hz. So, for this selection of the window, the instantaneous frequency change during the duration of the window (i.e., 5.72 Hz) is within the frequency resolution (i.e., 15.625 Hz) and it is possible to avoid additional errors related to STFT estimation for the first harmonic. The "safety factor" for this estimation S = (15.625 Hz /5.72 Hz) = 2.73 which is inside the specified threshold range (1.5-3).
So, in order to implement the described criterion for the proposed method, what need to be known, approximately, are the maximum chirp rates for signals with linear variation of the IF and the maximum chirp rates and frequency accelerations for quadratic variation of the IF, etc. Figure 18 shows the first iteration demodulation results. As shown in Figure 18a, the first extracted IFCF has a relatively big error; the error between the first fitted IFCF curve ( f 1 (t) = 7.6t 2 − 28.6t + 99) and the actual IFCF curve (3.05 × f (t)) is 63.1%. A clear peak, without smearing, cannot be captured in the spectrum, as shown in Figure 18d. The second iteration demodulation results are shown in Figure 19. From Figure 19a, it can be found that the extracted IFCF is coincident with the fitted IFCF curve ( f 2 (t) = −7.6t 2 + 37.1t + 58.8). There is a clear non-smeared peak in the spectrum, presented in Figure 19d, and the frequency of that peak is very close to the starting frequency of f 2 (t).
Based on f 1 (t), f 2 (t) and Equation (16), the estimated IFCF curve function is f 2 (t) = 8.5t + 58.8, and the IF estimation error between the estimated IF function f 2 (t) and the obtained one from the speed sensor function 3.05 × f (t) is 4.4%. Applying the termination threshold of 0.01, the third iteration is not needed.
Based on Equation (17), the demodulation operator is estimated as Φ − 2 (t) = exp −j2π t 0 8.5sds − 58.8t . With the obtained demodulation operator, the demodulated signal is estimated and is shown in Figure 20b. Figure 20c,d are the TFR and the spectrum of the demodulated signal. From the TFR, a constant time-frequency curve can be found, and its value is 58.8 Hz, which is a peak frequency in the spectrum, as shown in Figure 20d. There is a clear non-smeared peak in the spectrum, presented in Figure  19d, and the frequency of that peak is very close to the starting frequency of ( ). Based on ( ), ( ) and Equation (16), the estimated IFCF curve function is ( ) = 8.5 + 58.8, and the IF estimation error between the estimated IF function ( ) and the obtained from the speed sensor function 3.05 × ( ) is 4.4%. Applying the termination threshold of 0.01, the third iteration is not needed.
Based on Equation (17), the demodulation operator is estimated as ( ) = − 2 8.5 − 58.8 . With the obtained demodulation operator, the demodulated signal is estimated and is shown in Figure 20b. Figure 20c,d are TFR and the spectrum of the demodulated signal. From the TFR, it can be found a constant time-frequency curve; its value 58.8 Hz, which is a peak frequency in the spectrum, Figure 20d.     Then, the proposed technology is further validated for estimation of IFs of the second and the third harmonics of outer race defect frequency. Figure 21a,b is the TFR and the spectrum of the demodulated signal, related to the second harmonic. Figure 21c,d is the TFR and the spectrum of the demodulated signal, related to the third harmonic.
A constant time-frequency ridge and a peak can be detected in the TFR and the spectrum of demodulated signal, related to the second harmonic (Figure 21a,b). The peak frequency is 117.7 Hz, which is very close to the value 121. Then, the proposed technology is further validated for estimation of IFs of the second and third harmonics of the outer race defect frequency. Figure 21a,b is the TFR and the spectrum of the demodulated signal, related to the second harmonic. Figure 21c,d is the TFR and the spectrum of the demodulated signal, related to the third harmonic.
A constant time-frequency ridge and a peak can be detected in the TFR and the spectrum of demodulated signal, related to the second harmonic (Figure 21a,b). The peak frequency is 117.7 Hz, which is very close to the value 121.4 Hz, estimated by the speed sensor. Similarly, a constant time-frequency ridge and a peak can be detected in the TFR and the spectrum of the demodulated signal, related to the 3 third harmonic (Figure 21c,d). The peak frequency is 176.7 Hz, which is very close to the value 182.1 Hz, estimated by the speed sensor.  The IF estimation technology is further experimentally validated at different added noise levels. The Gaussian white noise is added to the experimental vibration signal and the SNRs are −4dB (strong background interference), 2dB (strong background interference), 8dB, 30dB, 40dB, and 50dB, respectively. The IF estimation technology is applied to these signals. All technology parameters are remaining the same. The obtained IFCF estimation errors under different noise levels are listed in Table 2. It can be found, that the results have no change under noise levels from 50dB to 8dB compared with the experimental signal without additional noise. The estimation error slightly increases to 5.8% for noise level 2dB. The estimation error further increases to 7.6%, for noise level -4dB. Hence, it is experimentally confirmed a high immunity level against noise, provided by the proposed technology in non-stationary conditions of very early stage of local bearing damage  The IF estimation technology is applied to these signals. All technology parameters remain the same. The obtained IFCF estimation errors under different noise levels are listed in Table 2. It can be found that the results have no change under noise levels from 50 to 8 dB compared with the experimental signal without additional noise. The estimation error slightly increases to 5.8% for noise level 2 dB. The estimation error further increases to 7.6%, for noise level −4 dB. Hence, it is experimentally confirmed that a high immunity level against noise is provided by the proposed technology in nonstationary conditions of the very early stage of local bearing damage development, and the relative damage size is 0.2%. It should be noted that the experimental raw signal has an initial noise level. As a result of the experimental investigation at different SNRs of the interference, it is obtained that a satisfactory estimation error of 5.8% is achieved at the additional noise level of 2 dB SNR. As experimental vibration data are already contaminated by noise, it is believed that a 5.8% error is achieved, at least, in the range of (0-2 dB) SNR. Thus, the (0-2) dB SNR range is a successful immunity level against interference achieved by the proposed method via experiments.
It is also shown that the −5 dB SNR is a successful immunity level against noise achieved by the proposed method via simulation, and the simulation error range is (1.2%-4.46%).
In the current experimental setup, the investigated faulty bearing does not slip. In the majority of bearing operation conditions under nominal loads, bearing slippage is insignificant and it is not needed to make any changes to the proposed method.
However, at high speeds and light loads or at no loads, bearing slippage may be essential. A minimum load must be applied to a bearing in order to avoid slippage. The value of this minimum load increases with speed.
If bearing slippage appears, the fault characteristic periods change and, therefore, the bearing defect frequencies are affected by the slippage. Bearing slippage needs to be taken into account for IF estimation of bearing defect frequencies. For bearing slippage conditions, a novel signal segmentation approach by slicing the bearing vibrations into small realizations (periods), considering simultaneously the relatively small number, e.g., 5-10, of these realizations (i.e., groups of realizations) for IF estimation/fault diagnosis purposes, and then tracking IF estimations/fault diagnosis results via multiple groups has been proposed, developed, investigated, and comprehensively experimentally validated [55][56][57][58][59][60]. This approach takes into account that bearing slippage is not a very fast developing process and, therefore, allows to consider simultaneously 5-10 realizations. It is possible to apply this novel approach for the proposed method for bearing slippage conditions by considering the STFT of a small number of realizations (i.e., groups), applying the method for multiple STFTs of groups, starting with groups characterized by relatively high SNR regions, and tracking the obtained estimation results via multiple groups in order to obtain the overall IF estimation for the whole nonstationary bearing vibration signal.
The proposed method is developed, firstly, as a generic IF estimation method and, secondly, with a specific view of application for IF estimation related to the bearing defect frequencies, for the purpose of bearing fault diagnosis. It is well known that all bearing defect frequencies, i.e., outer race defect frequencies, inner race defect frequencies, rolling element defect frequencies, and cage defect frequencies, are directly proportional to the fundamental rotational frequency and constants of proportionalities which are known as fault characteristic coefficients (FCCs).
Therefore, taking into account these direct proportionalities, it is known that for both stationary and nonstationary conditions of bearing operations and for a simultaneous appearance of multiple bearing local faults (e.g., outer race faults, inner race faults, etc.), the bearing defect frequencies of all harmonics are not intersecting. Even if a rotating shaft with two different rolling element bearings is considered, the bearing defect frequencies of all harmonics are also not intersecting.
Thus, the current version of the proposed IF estimation method is designed for non-intersecting IFs. Other known GDT approaches are also designed for non-intersecting IFs.
Taking into account the nature of the proposed method, the method does not have any limitations on the number of non-intersecting harmonics. It is shown above that the proposed method can successfully estimate multiple IFs for the following cases in a very noisy environment, where the SNR is −5 dB: • Therefore, the proposed technology has the potential ability to diagnose multiple bearing faults by IF estimations of bearing defect frequencies, related to local faults of outer race, inner race, cage, and roller elements. As the proposed technology has good immunity against interference, it allows to diagnose early stages of local bearing defects. In the performed experiment, the IF of the outer race defect is successfully estimated at the very early stage of outer race damage development, and the relative damage size is 0.2%.
As the proposed method does not require a priori knowledge about "true IF time dependency", it could be applied to two scenarios: (i) where bearing defect frequencies are known; and (ii) where bearing defect frequencies are unknown. The last case is very important for industrial applications.
It is possible to apply the proposed method for IF estimations also for induction motors, including inverter-fed induction motors. The proposed method could be used for IF estimations for harmonics related to issues other than bearing faults, such as mechanical faults in induction motors in stationary and nonstationary conditions. In most cases, such an application will be related to IF estimations for variable rotation speed harmonics. In nonstationary conditions of motor operations, these harmonics, normally, have time-variable amplitudes and time-variable nonlinear frequency modulations (pure linear frequency modulations do not often appear in the exploitation of induction motors). Therefore, the three important advantages of the proposed method: the ability of effective IF estimation in conditions of time-variable harmonic amplitude, the ability of effective IF estimation in conditions of time-variable nonlinear frequency modulations, and the ability of the estimation of IF of multiple harmonics, could be exploited for IF estimations for harmonics related to mechanical faults in induction motors.
The effects of inverter harmonics on motor vibration signatures are studied in detail in the literature. It is shown theoretically and experimentally that the vibration signatures, caused by the inverter, are related to the multiple harmonics of the grid supply frequency. It is a possibility that those harmonics will affect the reliability of IF estimations for induction motors as strong inverter-induced harmonics could dominate in vibration spectra and, therefore, will mask weaker harmonics of bearing defect frequencies. However, it is possible to eliminate the influence of inverter-induced harmonics from motor raw vibrations in the same way as gearbox mesh harmonics are removed from gearbox raw vibrations in order to obtain gearbox vibration residual signals, e.g., [48]. In these motor vibration residual signals, masking effects of inverter-induced harmonics will be essentially reduced and, thus, effective IF estimations of motor narrowband harmonics could be performed by the proposed technology.

Conclusions
The novel IF estimation technology, the multi-generalized demodulation transform (MGDT), is developed and investigated for nonstationary signals, whose nonstationary instantaneous frequencies are unknown and difficult to extract from time-frequency representations in an essentially noisy environment. The theoretical basis of the novel IF estimation technology is created. The main novel feature of the IF estimation technology is the proposition of multiple adaptive iterations of the GDT for the same specific signal component to concentrate blurry TFR energy in a step-wise manner. The important feature of the novel MGDT is that every new iteration of the GDT is adaptive, i.e., it is based on previous iteration results. Our extensive literature search clearly shows that nobody, in worldwide terms, has previously proposed the main important novelty of the MGDT: an adaptive iterative application of the GDT for the same signal component of a multi-component signal.
The proposed IF estimation technology allows to remove the GDT's restriction of knowledge (a priori or in-advance estimation from a speed sensor) of a true instantaneous frequency that is important for industrial applications of the proposed technology. With the proposed adaptive iterations of the generalized demodulation procedure for the same specific signal component, the energy concentration on the time-frequency ridge is enhanced and a specific signal component with the instantaneous frequency curve, that is polluted by noise, can be accurately estimated.
Novel validation of the proposed technology is successfully performed by the simulation for four types of amplitude-and frequency-modulated nonstationary signals under a strong background noise (signal to noise ratio is −5 dB): (i) single-component signal with a nonstationary amplitude and nonstationary, nonlinear (quadratic) time variation of the instantaneous frequency; (ii) multi-component signal with components that have nonstationary amplitudes and nonstationary, linear instantaneous frequencies with different chirp rates; (iii) multi-component signal with components that have nonstationary amplitudes and nonstationary, nonlinear (quadratic) instantaneous frequencies with different chirp rates and frequency accelerations; and (iv) multi-component signal with components that have nonstationary amplitudes and nonstationary, nonlinear (cubic) instantaneous frequencies with different chirp rates, frequency accelerations, and the first derivatives of frequency accelerations.
The achieved estimation errors for the simulated nonstationary signal in the presence of a strong interference (−5 dB signal to noise ratio) are in the following ranges: (i) (1.46-4.46)% for a multi-component signal with linear frequency modulations and nonstationary amplitudes of components, (ii) (1.2-1.4)% for a multi-component signal with nonlinear (quadratic) frequency modulation and nonstationary variation of amplitudes of components, (iii) (1.2-1.6)% for a multi-component signal with components that have nonstationary amplitudes and nonstationary, nonlinear (cubic) time variations of the instantaneous frequency with different chirp rates, frequency accelerations, and the first derivatives of frequency accelerations, and (iv) 0.85% for a single-component signal with a nonstationary amplitude and nonstationary, nonlinear (quadratic) time variation of the instantaneous frequency.
It is validated by experimental trails that the novel IF estimation technology successfully estimates nonstationary bearing defect frequencies for the first, second, and third harmonics of bearing defect frequencies (without a priori knowledge of true time-varying instantaneous frequency dependency) under time-varying rotation speeds and time-varying amplitudes of these harmonics in the presence of interference and at the very early stage of outer race damage development: the relative local damage size is 0.2%. It is shown that the estimation error for the instantaneous defect characteristic frequency estimation is 4.4%. The instantaneous bearing defect characteristic frequencies of the first, second, and third harmonics, that are polluted by noise in the time-frequency representation, are accurately estimated by the proposed technology.
The estimation errors for nonstationary bearing defect frequencies of the first, second, and third harmonics are experimentally evaluated for six difference levels of added Gaussian noise interference, characterized by a signal to noise ratio (SNR) from 50 to −4 dB, including for two levels of a strong noise interference. It is shown that the estimation errors are unchanged (i.e., staying at a level of 4.4%) under noise levels from 50 to 8 dB, compared with the case without additional noise. The estimation error slightly increases to 5.8% for noise level 2 dB (a strong noise interference) and it further increases to 7.6% for noise level −4 dB (a strong noise interference). These results experimentally confirm a high immunity level against noise delivered by the proposed technology. Based on the experimental validation, for an acceptable IF estimation accuracy, the SNR for the proposed technology should not be lower than (0-2 dB).
The IF estimation technology, the MGDT, is a promising novel concept that could be widely employed for sensor-based vibration instantaneous frequency estimation for nonstationary signals of electromechanical systems in nonstationary operations, in a noisy environment, and without knowledge (a priori or by in-advance estimation from a speed sensor) of true instantaneous frequency variations. The proposed IF estimation technology could be employed for rotating components/machineries such as rolling element bearings, gearboxes, compressors, and induction motors, and could be also effectively employed for other technologies, e.g., low-frequency acoustic non-destructive technologies, e.g., [61,62], motor current signature analysis for induction motors, e.g., [63], ultrasound fault detection technologies, etc.

Conflicts of Interest:
The authors declare no conflict of interest.