Mechanical varying non-stationary signal separation method based on instantaneous frequency estimation

Mechanical equipment often works on variable speed condition, its corresponding vibration signal presents multi-component, modulation coupling with fast time-varying instantaneous frequency (IF), how to effectively compute IF and realize fasting varying non-stationary signal decoupling separation plays an important role in mechanical system fault diagnosis. In this paper, a sparse representation method called multi-scale chirp sparse representation (MSCSR) is introduced to identify, extract, and trend IF for achieving a highly concentrated time-frequency energy. Simulation demonstrates that the proposed method performs better than traditional IF estimation method. Furthermore, an adaptive time-varying filter is constructed using the extracted instantaneous frequency to decouple non-stationary fast signal. Ultimately, by rapid instantaneous frequency fluctuation experiment, the effectiveness of proposed method for fast strong time-varying signal is validated, it can effectively extract rapid oscillation instantaneous frequency, and the error is less than 10%.


Introduction
Mechanical equipment under variable speed condition is almost everywhere, such as gearbox under raising speed, its fault vibration signal often reflects multicomponent coupling non-stationary property with fast IF. As a result, how to separate fast varying signal based on IF estimation has become a critical requirement for vibration-based mechanical equipment condition monitoring methods.
In general, the existing signal processing methods are majorly applicable to the steady speed working condition. 1 Under variable speed, the fault vibration signal will lose the periodic regularity, and present the frequency modulation, amplitude modulation property, leading to the invalidity of classical methods. 2 With the development of non-stationary signal processing technology, scholars have gradually carried out research on the mechanical fault vibration signal separation. 3,4 In general, they can be divided into four universal categories, order tracking, 5 order cyclostationary analysis 6 and timefrequency analysis, 7 tine-varying filter. 8 As for order tracking, 9 its main procedure is to transform non-stationary signal into stationary or cyclically stationary signal by equal interval sampling technology in angle domain, so as to establish a bridge between signal processing under variable speed condition and constant speed condition, so that the classical signal processing method can play its role again. However, it must be pointed out that the implementation of order tracking mostly relies on hardware device such as encoder or key phase sensor, nevertheless, in some critical condition, as sensor cannot be installed which brings challenges to order tracking implementation. Under this background, scholars propose tacholess order tracking (TLOT) method, 2 it can directly achieve order tracking on non-stationary by using time-frequency transformation without hardware, but it still needs to be further improved, mainly in terms of accurately tracking and estimating instantaneous frequency. As instantaneous frequency reflects the dynamic state of mechanical equipment, so improving the estimation accuracy of IF becomes the core of TLOT. Furthermore, considering the structure of most mechanical parts has spatial symmetry, its fault vibration signal presents essentially cyclostationary in the angle domain, by using order tracking and synchronous averaging, we can convert the non-stationary signal in time domain into cyclostationary signal in angle domain again and establish order cyclostationarity analysis 10 ultimately which is convenient for processing mechanical varying non-stationary signal, such as order spectrum, 11 order cepstrum, 12 envelope order spectrum, 13 order bispectrum, 14 and high-order spectrum. 15 As for time-frequency analysis, considering only from the time domain or frequency domain, we can not obtain the instantaneous time-frequency property which is the core of non-stationary signal processing, 16 the time-frequency analysis provides the joint distribution information in the time domain and frequency domain which can effectively give attention to both time resolution and frequency resolution, it is very suitable for estimating instantaneous frequency and separating mechanical varying non-stationary signal. Finally, as for time-varying filter, 17 its center frequency needs to match the instantaneous frequency and its parameters also should be set reasonably to achieve mechanical varying non-stationary signal separation.
In conclusion, instantaneous frequency estimation plays an important role in non-stationary signal decomposition. A series of methods have been introduced to achieve this goal, such as reassignment synchrosqueezing transform, 18 multi-scale chirplet path, 19 short-time Fourier transformation, 20 considering the application of MSCSR on nonlinear FM signal estimation, 21 parameter identification of linear time-varying SDOF system, 22 in this paper, we introduce it to estimate IF in mechanical varying non-stationary signal whose matched estimator is constructed according to fast varying IF, so as to improve the signal time-frequency representation energy aggregation. On the basis, by using the extracted IF acting as center frequency of the filter, we design a adaptive time-varying filter for decomposing mechanical varying non-stationary signal ultimately.
Finally, the process diagram of proposed methodology is presented in Figure 1. Firstly, the continuous wavelet transform is applied to obtain corresponding wavelet energy spectrogram, and then the number of sub-signal component in vibration signal is accurately judged according to the time-frequency distribution, which greatly improves the calculation efficiency. Secondly, MSCSR self-adaptive selects an appropriate dynamic time support region according to the IF varying law, tracks, and estimates the IF of sub-signal component which possesses largest energy. Finally, an adaptive time-varying filter is designed using estimated IF as the center frequency, and is applied for decomposing vibration signal into several independent components. A simulation and mechanical system high speed fluctuation experiment is to validate the proposed method, and compared with IF recognition result and signal decomposition result by using Hilbert-Huang transform (HHT), 23 continuous wavelet transform (CWT), 24 synchrosqueezing wavelet transform (SWT), 24 variational mode decomposition (VMD), 25 empirical mode decomposition (EMD) 26 respectively, the proposed method has a better effectiveness and accuracy.

An adaptive time-varying filter based on instantaneous frequency estimation
The key steps of proposed methodology as shown in Figure 1, are discussed in this section. The basic function library of MSCSR for tracking and estimating IF is shown as following.
In the above formulas, D represents base function library, h a;b;I t ð Þ is the multi-scale chirp base function, I = [kN2 2j ;(k + 1)N2 2j ] represents dynamic time support region, it describes time -varying characteristic of instantaneous frequency, j ¼ 0; 1; Á Á Á ; log 2 N À 1 ½ symbolizes analysis scale coefficient, N represents sampling length, it can reflect the description precision of time-varying characteristics. As for k, it is the sequence number of dynamic analysis section in j scale, and it varies from 0 to 2 j 21. K a,b,I is a normalization coefficient, a is symbolizes frequency offset coefficient, b is the frequency rake ratio, those parameters can well reflect the varying state of instantaneous frequency. Finally, 1 I t ð Þ is a rectangular window function. As for the sampling frequency f s , it must satisfy the following conditions according to Shannon sampling theory.
The multi-scale chirp basic function is applied to project fault response signal segment by segment, and calculating projection coefficient b I and corresponding basic function h a;b;I on each dynamic time support region I. If the fault signal is more similar to basic function h a;b;I , its projection coefficient b I will be larger, and the corresponding energy of basic function h a;b;I will be larger. As a result, fault response signal must be projected to multi-scale chirp basis function, a specific basis function will be obtained that is corresponding to maximum projection coefficient b I on each dynamic time support region I. Then, in accordance with maximum projection coefficient b I and its corresponding chirp base function, a dynamic linking algorithm called P must be designed to connect a series of basis functions, and there is no overlap between the basis functions each other on whole dynamic time support regions I.
At the same time, the corresponding projection coefficient set and basic function set are shown in following.
The steps of the connection algorithm P in MSCSSD can be summarized as follows.
Step 1: Initialization. i represents the serial number of sub-signal dynamic time support region, E d i represents the total energy of sub-signal before dynamic time support regions I i of sequence number i, l i represents previous serial number of dynamic time support region attached to i, E e i indicates energy of sub-signal corresponding to the projection coefficient on each dynamic time support region I i , initializing E d j ¼ 0, l i ¼ 0.
Step 2: For every dynamic time support region I i in the set I ¼ I 1 ; I 2 ; ::: f g , find all next dynamic time support regions adjacent to it to form a set I j È É , if There is the following formula, The dynamic linking algorithm called P guarantees extracted sub-signal connected by basic function set H over the entire dynamic time support region set I ¼ I 1 ; I 2 ; ::: f gis most similar to original fault response signal, and instantaneous frequency corresponding basic function h a;b;I can be computed f I i t ð Þ ¼ a m + 2b m t; t i eI i on the dynamic time support region I i , the instantaneous frequency f can be tracked and estimated ultimately by connecting the corresponding frequency curve set f I ¼ f I 1 ; f I 2 ; :::: f gover the whole dynamic time support region set I ¼ I 1 ; I 2 ; ::: f g in chronological order. As a result, the estimated algorithm of instantaneous frequency can be summarized as shown in Figure 2.
Considering adaptive time-varying filter is based on classical filter, the proposed adaptive time-varying filter based on MSCSR can make its center frequency consistent with estimated IF according to the mechanical varying non-stationary signal. As the signal often presents obvious amplitude modulation and frequency modulation phenomenon, we can make the estimated IF that it is corresponding to carrier frequency acting as center frequency and side frequency act as cut-off frequency respectively. Therefore, an adaptive filter H(s,t) design algorithm based on MSCSR is proposed in this paper, it can be obtained by transferring classical low-pass prototype filter H(s). As its central frequency and filter bandwidth can automatically vary according to the non-stationary property, as a result, it is very suitable for mechanical varying non-stationary signal separation and the design procedures is shown in Figure 3.
Step 1: As for signal s t ð Þ; s t ð Þ ¼ P l¼N l l¼1 s l t ð Þ, l symbolizes the serial number of sub-signal, l ¼ 1; 2; Á Á Á N l , N l is total amount of sub-signals. As for sub-signal with serial number l, its carrier frequency curve and side frequency are w l z t ð Þ, b l z t ð Þ; 0 ł t ł N À 1; which are obtained by MSCSR, N represents the signal length.
Step 2: Obtaining the spectrum S f ð Þ is corresponding to signal s t ð Þ.
Step 3: In order to extract sub-signal s l t ð Þ, the adaptive time-varying filter H t; f ð Þ is designed on Þas follows. w l z t 1 ð Þ acts as the center frequency, and n harmonic frequency of b l z t 1 ð Þ is half-bandwidth. Considering fast varying signal often presents modulation characteristic, in order to ensure extracted sub-signal without distortion, as Chebyshev II filter has fluctuation in the stop-band and a steep transition band, which is suitable for signal separation, this paper adopts it acting as prototype low-pass filter.
Step 4: The adaptive time-varying filter H t; f ð Þ is applied to obtain corresponding filtered sub-signal Step 5: The IFFT transformation is used to obtain corresponding sub-signal s l t 1 t ð Þ in time domain at the time t ¼ t 1 0 ł t 1 ł N À 1 ð Þ , repeating above process at any time to obtain filtered sub-signal s l t ð Þ.
Step 6: Repeating above procedures from step 1 to step 5, and achieving signal separation s t ð Þ completely.

Simulation analysis
In order to analyze the effectiveness of proposed method on varying non-stationary signal decomposition, a simulation signal x(t) containing three subsignals is presented in formula (9). Considering actual noise effect, adding 20% SNR level Gaussian white noise s t ð Þ to simulation signal, the sampling frequency is 512 Hz and sampling time is 1 s. Signal-to-noise ratio (SNR) is defined as in formula (10), A signal and A noise symbolize the RMS value of simulation signal and noise component respectively.
Firstly, we use the wavelet spectrum by using continuous wavelet transform to observe the number of subsignals as shown in Figure 4. It can be significantly seen that the simulation signal contains two sub-signals and their varying IFs conform to linear varying trend, and mainly concentrate in the range of [0 Hz, 40 Hz] and [0 Hz, 128 Hz].
Secondly, we apply the MSCSR to estimate IFs. We set the frequency variation range on MSCSR is from 10 to 400 Hz, f min = 30 Hz, f max = 300 Hz, search frequency slope is ranging from 0 to 200 Hz, search resolution is  1 Hz/s, and the extracted IFs f extract1 (t), f extract2 (t) are shown in Figure 5(b) and (d) respectively.
Time-frequency representation in Figure 5(a) and (c) present the varying law of instantaneous frequency in x 1 t ð Þ and x 2 t ð Þ, it can be seen that the IFs possess linear varying property. Compared with Figure 5(a) and (c), the estimated IF by MSCSR as shown in Figure 5(b) and (d) is consistent with the theoretical value, this result presents the estimated value by using the proposed method has a perfect precision.
For comparison, we also uses Hilbert-Huang Transform (HHT), Continuous Wavelet Transform (CWT), synchrosqueezing wavelet transform (SWT) to extract IF, as shown in Figure 6. Compared Figure 5 with Figure 6. It can be clearly seen that the tracking effect by using CWT in the beginning and end stages has certain deviation, and estimated value has a large fluctuation phenomenon. On the other hand, although SWT can achieve the tracking in general, the accuracy is not high. As for HHT, as IF presents a linear varying trend, it can not effectively achieve tracking, this phenomenon implies that the existing methods can not be effectively used to estimate IF under large fluctuation background, while the estimated IF based on MSCSR is more consistent with theoretical value, which indicates the proposed method in this paper adopts multi-scale chirp basic function whose parameters are adjustable, can dynamically and self-adaptive track and estimate fast varying IF. Furthermore, in order to quantify IF identification accuracy, we adopt root mean square (RMS) to act as precision index (Index of accuracy, IA). f i n ð Þ represents IF identification value, f t n ð Þ symbolizes IF theoretical value, N is number of points. Generally, if IA index has a smaller value, it indicates that the recognition value is connected with theoretical value excellently, that is, the identification accuracy is more higher. Table 1 shows corresponding IA values by using four methods, MSCSR, SWT, CWT, and HHT, IA1 and IA2 represent the recognition accuracy of sub-signals x 1 t ð Þ, x 2 t ð Þ respectively. As for simulation signal x t ð Þ, it can be seen from Table 1 that the proposed method is obviously superior to HHT, CWT, and SWT.
In order to further verify the anti-noise ability and accuracy of proposed method on instantaneous frequency estimation under strong noise, compared with SWT, CWT, and HHT, we add different SNR level Gaussian white noise s t ð Þ into simulation signal x(t) which is ranging from 210 to 20 dB, the IA index under different SNR is calculated and shown in Figure 7. As can be seen from the Figure 7, the corresponding IA value by using proposed method is less than compared methods which is all less than 5%, indicating that the proposed method has high accuracy and good noise resistance.
Through above simulation signal, the advantage of proposed method in IF tracking has been preliminarily validated, and it has potential application value in nonstationary signal separation based on filter technology. Considering that the actual mechanical equipment often run in large speed fluctuation condition, its IF often presents curvilinear varying trend. As a result, we simulate a mechanical varying non-stationary signal as shown in formula (12). It contains two sub-signals s 1 (t) and s 2 (t), they reflect modulation property and their instantaneous frequencies are overlapping as shown in formula (13) and (14), s t ð Þ represents gaussian white noise, sampling frequency is 1024 Hz, the corresponding time waveform is shown in Figure 8.
It can be seen that the simulation signal s t ð Þ has frequency-modulation and amplitude-modulation property, and the time-varying carrier frequencies for s 1 (t) and s 2 (t) are shown in formula (15) and (16). Besides, as for the IF corresponding to amplitude modulation component in s 1 (t) and s 2 (t) is shown in formula (17). As a result, the simulation signal s t ð Þ contains six instantaneous frequencies, Considering the simulation signal s t ð Þ represents remarkable non-stationary property, we firstly use short-time Fourier transform (STFT) spectrum to determine the number of sub-signal as shown in Figure  9, six time-varying instantaneous frequencies, especially two time-varying carrier frequencies can be highly legible, and the sub-signal's side-bands have similar varying trend, which sufficiently reflects the coupling appearance in simulation signal. Furthermore, from Figure 10, compared with side frequency, the energy of carrier frequency is more powerful, we will extract two carrier frequencies primarily.
Considering above issues, we apply the proposed method on simulation signal to identify instantaneous carrier frequency, the estimated frequency f extract1 (t), f extract2 (t) and correspondingly theoretical carrier frequency f 1 (t), f 2 (t) are shown in Figure 11. It can be seen that the proposed method in this paper also has good tracking accuracy for fast varying instantaneous frequency, as a result, other four side-band frequencies Additionally, a time-varying filter H 2 s; t ð Þ is designed based on extracted frequency for decomposing simulation signal s t ð Þ. According to formula (7) and (8), the estimated instantaneous carrier frequency is taken as central frequency and estimated side-band frequency is corresponding cutoff frequency, and the filter's timefrequency response is shown in Figure 12. It can be seen intuitively that the amplitude-frequency property of time-varying filter is very consistent with signal s t ð Þ. Considering amplitude-frequency property of simulation signal s t ð Þ in Figure 10, adopting time-varying filter H 2 s; t ð Þ decomposes this signal, it can well obtain a signal s p2 t ð Þ in Figure 13 corresponding to sub-signal s 2 t ð Þ according to the following formula (18). Compared with the theoretical simulation sub-signal s 2 t ð Þ, whose time-frequency spectrum presents in Figure 9, the spectrum of signal s p2 t ð Þ has some     distortion, as the adaptive time-varying filter exists amplitude loss, while the loss is not very serious from Figure 13.
Considering the adaptive decomposition method owns unique advantage without presetting basis function, we select empirical mode decomposition (EMD), variational mode decomposition (VMD) for comparison as shown in Figure 14. It can directly see that the number of sub-signals obtained by adaptive decomposition method is not equal to real number, and each sub-component presents modulation attribute, which indicates that the decomposition results reflect aliased phenomenon. In conclusion, we can know that when the signal has a fast varying IF, adopting signal decomposition based on presetting a series of basis functions will be more realistic.

Experimental application
In this section, we use the collected rolling bearing fault vibration signal under variable speed condition as shown in Figure 15 to validate the accuracy of the proposed method on instantaneous rotating frequency estimation. The rolling bearing fault test platform is composed of AC motor, motor speed controller, rotating shaft, loading system and experimental rolling bearing, NI CompactDAQ. It can carry out rolling bearing fault experiment under different rotating speeds. The test rolling bearing is 52732QT and we adopt electric discharge machining to cause outer ring damage. Here, we generate two non-stationary signals, whose instantaneous rotating frequencies are in accordance with fast linear varying trend and sinusoidal varying trend respectively.
Firstly, we collect the fault vibration signal with linearly varying frequency as shown in Figure 16 and the rotating speed increases from 900 to 1080 rad/min in 6 s. Furthermore, except the proposed method in this paper, we also use several different methods to analyze the signal and estimate its instantaneous frequency, the result is shown in Figure 17.
It can be obviously seen from Figure 17 that except HHT, the estimated values corresponding to other methods can well match with the theoretical result. Moreover, the proposed method performs best in estimation accuracy, this conclusion can also be proved by the IA index as shown in Table 2.
On the other hand, we analyze another fault vibration signal as shown in Figure 18 under a large speed fluctuation condition which presents sinusoidal varying instantaneous frequency, and the estimated IF is shown in Figure 19.
Although there is a certain deviation between the instantaneous frequency value identified by proposed method and theoretical value, it is basically consistent with theoretical value on the trend, and the recognition effect is obviously better than SWT and HHT. Although SWT can also identify the corresponding trend, the estimated instantaneous frequency curve still contains a lot of burrs. Besides, the IA value in Table 3 also can prove above conclusion. Furthermore, in order to evaluate the efficiency of proposed method in experiment, we set another evaluation index, which is the calculation time t, Table 4 shows the consumed time by each methods. We can see that the proposed method in this paper has certain advantage in computational efficiency, which indicates that the MSCSR can be applied to the field of rotating machinery online monitoring.
Considering that the rotating machinery works under the large speed fluctuation condition, its instantaneous rotating frequency often reflect significantly non-linear property, in order to solve this problem, we can divide the instantaneous frequency into a series of segments for converting the global non-linearity into local linearity, or we can use polynomial frequency factor instead of linear frequency factor in MSCSR whose formula is f t ð Þ ¼ a + 2bt.

Conclusion
An adaptive time-varying filter design method based on instantaneous frequency estimation is proposed and applied for achieving mechanical varying nonstationary signal separation. The proposed method is also validated by simulation signal, compared with other algorithms, the accuracy of proposed method is proved. Finally, an instantaneous rotating frequency estimation experiment based on rolling bearing fault vibration signal verifies application reliability. The main conclusion of this article is as follows.
(1) As the mechanical fault vibration signal presents significantly fast varying instantaneous     frequency and multi-component coupling modulation property, considering the multi-scale chirp basic function has the advantage of high performance computing, high accuracy, and strong noise resistance, the MSCSR using this basis function can effectively extract the instantaneous frequency that presents curvilinear varying law. Furthermore, an adaptive filter designed with the extracted instantaneous frequency has high calculation efficiency which solves the problem that the classical filter needs to select filter's parameters. Finally, the result proves that the proposed method can extract the modulation component and decomposing mechanical varying non-stationary signal. (2) Considering the instantaneous rotating frequency under large speed fluctuation condition shows fast varying nonlinear property, in order to use MSCSR to better estimate it, we propose two solutions which utilizes transformation idea and alternative idea respectively. As for the former, we can divide the globally nonlinear instantaneous frequency into several locally linear instantaneous frequencies, and then use MSCSR to achieve accurate estimation. As for the latter, we can replace the linear frequency factor f t ð Þ ¼ a + 2bt in the existing MSCSR with the polynomial frequency factor, which can better conform to the nonlinear varying frequency.

Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.