A Novel Fault Feature Recognition Method for Time-Varying Signals and Its Application to Planetary Gearbox Fault Diagnosis under Variable Speed Conditions

The existing time-frequency analysis (TFA) methods mainly highlight the time-frequency ridges of the interested components by optimizing the time-frequency plane to facilitate the extraction of the relevant components. Generalized demodulation (GD), order tracking (OT), and other methods are generally used in conjunction with the TFA methods to realize the transition from a time-varying signal to a stationary signal, and finally identify the fault feature through a time-frequency plane. Generally, it is necessary to clarify the accuracy of the estimated components such as the rotational frequency or the fault characteristic frequency (FCF) during the operation of the GD or OT methods. Unfortunately, it is not only difficult to extract and locate rotational frequency or FCF, but also complicated in the whole estimation process. In this paper, a simple yet readable method is proposed to reveal the fault feature of time-varying signals. First, the method only needs to extract an arbitrary instantaneous frequency (IF). This is different from the GD method which needs to estimate and locate all phase functions. Then, it converts all variable frequency curves into corresponding lines parallel to the frequency axis based on the extracted IF to determine the proportional relationship between the components. Finally, to further improve the readability of the final results, we reduce the dimension of the transformed time-frequency representation to generate a two-dimensional (2D) energy-frequency map with high resolution and the same proportion. Subsequently, the performance is validated by simulated and experimental data.


Introduction
With the increasingly serious energy crisis and environmental pollution, as a green renewable energy, wind energy has attracted more and more attention throughout the world. As the main transmission component of wind turbines, planetary gearboxes are prone to failure due to long-term impact and alternating loads [1,2]. Therefore, accurate identification of gearbox fault characteristics is the prerequisite for fault diagnosis and health management of mechanical equipment [3,4].
Under the condition of variable speed, the signal of planetary gearbox has typical aperiodicity and non-stationarity. Undoubtedly, directly obtaining the fault characteristic frequencies of non-stationary signals by using the traditional Fourier spectrum analysis such fast Fourier transform is almost 1.
Can the process of estimating and locating all phase functions in GD with only one instantaneous frequency be simplified and achieve similar smooth transformation? 2.
Can the process be further simplified? For example, can the instantaneous frequency that needs to be estimated be arbitrary? 3.
Can more intuitive and readable results on the basis of TFR be presented?
In order to realize this idea, we mainly improve the following algorithms: 1.
To obtain stable high resolution TFR, an improved multi-synchrosqueezing transform (IMSST) algorithm based on MSST is proposed. According to the intermediate frequency extracted in the first step, the time-varying full frequency can be directly converted into stable full frequency.

2.
To transform the instantaneous frequency ridges into a series of lines parallel to the frequency axis, we improve the instantaneous frequency estimation operator based on the MSST algorithm, so that we achieve the result which is similar to GD using only an arbitrary extracted IF. 3.
To present more intuitive and readable results, we propose a simple data dimension reduction method, which generates a more readable two-dimensional (2D) energy-frequency diagram.
This paper contributes to the following three major aspects:

1.
A three-step model is used to enhance readability of the final results.

2.
The proposed method is validated in multicomponent and planetary gearbox simulation signals in the form of increasing signal complexity.

3.
The proposed method is applied to diagnose the planetary gearbox fault under a time-varying condition and directly recognize the fault type from the 2D energy-frequency map without using any other method.

Theoretical Description
According to the previous ideas, the proposed method is mainly divided into three steps: firstly, extract an arbitrary instantaneous frequency; then, obtain steady and high resolution TFR from time-varying signals; finally, produce a 2D energy-frequency map. Therefore, in this section, we will introduce the relevant theories of these three steps in turn.

Extraction Algorithm of IF
At present, the peak searching algorithm is the most common algorithm in IF extraction. It obtains the IF of interested components by searching the peak of time-frequency ridge energy. However, due to the interference of modulation phenomenon and strong noise in measured signals, it is generally difficult to extract high accuracy time-frequency ridges directly.
In Ref. [33], G. Yu proposed a general linear chirplet transform (GLCT). He used a time-varying demodulation operator to produce a rotational effect on the time-frequency plane, so as to solve the problem of energy divergence caused by modulation. The GLCT constructs discrete time-varying demodulation operators by defining rotating degree α which is written as follows: where, T is sampling time and f s is sampling frequency. N denotes the number of discrete time-varying demodulation operator c. From Equation (1), the rotating degree α is equally divided into N parts. However, considering the difference of the instantaneous change trend of general time-varying signals, it is difficult to completely represent the change trend of the time-frequency curve only through the division of equal intervals.
In this paper, we abandon the above estimation form and reconstruct the time-varying demodulation operator directly by estimating the trend of the instantaneous frequency. In this way, the energy concentration of the instantaneous frequency curve is maintained to the greatest extent.
An IF extraction algorithm based on time-varying demodulation and the STFT (TDSTFT) is proposed, which combines the idea of time-varying demodulation and the peak searching algorithm. It repeatedly estimates the time-varying demodulation operator through the changing trend of IF to optimize the TFR, and then extracts ridge line using the peak searching algorithm. The following is a detailed theoretical introduction.
An AM-FM signal is defined as: where, A(t) and Φ(t) are their instantaneous amplitude and instantaneous phase, respectively. And ϕ(t) = Φ (t) is regarded as the ideal IF of the analytical signal. In a short time u, the IF of signal s(t) can be approximately regarded as a linear equation, which is expressed by the first-order Taylor expansion: Then, the corresponding STFT is expanded into the following form: where, g(t) is a window function. Here we choose the Gauss window, whose expression is written as: According to Equation (4), the signal is modulated and the energy of STFT diverges due to the presence of e iϕ(t)(u−t) 2 /2 . Therefore, we introduce the time-varying demodulation operator e i c(t)(u−t) 2 /2 to eliminate the influence of e iϕ(t)(u−t) 2 /2 , which makes the TFR energy more concentrated. After introducing the time-varying demodulation operator, the STFT is rewritten as: Because the c(t) is unknown, here, we adopt the approximation strategy and achieve the accuracy allowed by the algorithm through the change trend of the approximate IF step by step. Considering the uniformity of energy concentration, this paper directly uses the change trend of instantaneous frequency to express c(t).
The specific process of the whole extraction algorithm (Algorithm 1) is as follows: length(s(t)); P = peak − searching(TDSTFT); I(t) = poly f it(P); I(t) = poly f it(P); I (t) = di f f (I(t)); 3: Iteration: for j = 1 : n; The IF of a component can be obtained more accurately after the above steps. The following is an example to illustrate the precision of the IF extraction algorithm.
We use a numerical example to illustrate the effectiveness of this method and the simulated signal is s(t) = sin 2π 10t + 10t 2 − 19 9 t 3 + 9 50 t 4 − 1 200 t 5 + n(t), where n(t) is −2 dB while Gaussian noise and the sampling frequency is set to 200 Hz. Figure 1a,c shows the TFRs by STFT and TDSTFT after four iterations, respectively. The energy of TFR, as shown in Figure 1c, is more concentrated. Figure 1b,d shows the extracted IFs obtained from Figure 1a,c, respectively, using the peak searching algorithm (where, the blue line is the extracted IF and the red line is the real value). Through the comparative analysis, it is illustrated that the IF extracted from graph d in Figure 1 is obviously closer to the true value, which also proves the effectiveness of the proposed method.

The TFRs by IMSST
According to the expression of time-varying signals in Equation (2), a standard STFT is written as follows: where, g * (t) is the complex common of g(t).
Then, the SST based on the STFT is expressed as: where, γ is a very small threshold to prevent excessive error when V g f (t, η) = 0. And δ function is defined as:

The TFRs by IMSST
According to the expression of time-varying signals in Equation (2), a standard STFT is written as follows: where, ( ) g t * is the complex common of ( ) g t .
Then, the SST based on the STFT is expressed as: where, γ is a very small threshold to prevent excessive error when And δ function is defined as: In Equation (9), ( , ) f t ω η represents the IF estimation of signals on t and η . Since the phase after STFT is not affected by the window function ( ) g t , we use ( , ) g f V t η to estimate the IF: In Equation (9), ω f (t, η) represents the IF estimation of signals on t and η. Since the phase after STFT is not affected by the window function g(t), we use V g f (t, η) to estimate the IF: where, arg[A] denotes the angle of the complex A. Considering that ω f (t, η) is a complex number, we only take the real part as the estimation of IF.
Although SST cannot deal with strong time-varying signals well, TFR with high concentration can also be obtained by continuously utilizing SST concentrating energy. This is the core idea of MSST method.
In a short time, the phase Φ(τ) of the original signal s(τ) is expanded as: The purpose of (15) is to simplify the STFT coefficient V g f (t, η) and V g f (t, η) is rewritten as follows: Similarly, the IF estimation is expressed as: In order to facilitate the writing, here, we useω [3] f (t, η) instead ofω f (t,ω f (t,ω f (t, η))). Similarly, the estimation of IFω f (t,ω f (t,ω f (t,ω f (t, · · · )))) after the Nth iteration is written asω [N] f (t, η). According to the law, it is not difficult to deduce the estimation of IF after the Nth iteration. Therefore, MSST is expressed as: To illustrate the advantages of MSST in sharpening time-frequency domain, we take a simple time-varying signal, s(t) = 10 cos 2π(500t − 2e −2t sin(20πt)) as an example, and the sampling frequency of this signal is set as 4096 Hz. Figure 2 shows its TFRs obtained by SST, 2-SST, 4-SST and MSST (N = 6), respectively. And their resolution is quantified by Rényi entropy shown in Table 1. In general, the smaller the value of Renyi entropy [34], the higher the resolution obtained. It proved that MSST method has the smallest value and produces more energy-intensive TFRs as shown in Table 1. In order to obtain stable and high resolution TFR, we directly improve MSST by referring to the idea of demodulation, so that the final TFR is a series of lines parallel to the frequency axis.
Suppose we obtained an arbitrary instantaneous frequency ( ) IF f t using TDSTFT and then the proportional coefficient at any time can be expressed as:   In order to obtain stable and high resolution TFR, we directly improve MSST by referring to the idea of demodulation, so that the final TFR is a series of lines parallel to the frequency axis.
Suppose we obtained an arbitrary instantaneous frequency f IF (t) using TDSTFT and then the proportional coefficient at any time can be expressed as: According to the proportional coefficient K, the transformed IF is written as: According to (22), the IFs are redistributed on the time-frequency plane to obtain the target time-frequency representation.
Finally, the expression of IMSST is written as follows:

Two-Dimensional Energy-Frequency Map Obtained by IMSST
In the fault diagnosis of time-varying signals, scholars not only expect to retain as much original information as possible, but also hope to be able to produce familiar analysis graphs to directly judge, such as the Fourier spectrum, order spectrum, and so on. We can obtain high resolution time-frequency plane using IMSST. In order to further clarify the proportional relationship between signal components, it is necessary to eliminate the influence of time and transform the TFR into a 2D frequency-energy map.
The specific conversion is as follows: where, ITsre(ω j ) represents the amplitude (energy) at frequency ω j and M is the number of samples, i.e., the length of the original time series. By (25), the influence of t is eliminated and the corresponding 2D energy-frequency map is obtained. The flowchart of the proposed method is shown in Figure 3, which is also summarized as follows: 1.
Apply acceleration sensor to obtain the vibration signal s(t).

2.
Extract an arbitrary IF f IF (t) from s(t) using TDSTFT.

3.
Calculate the proportionality coefficient K according to f IF (t).

5.
Transforming instantaneous frequency curves to a series of lines which are parallel to the frequency axis according to the proportionality coefficient K using IMSST, and then obtaining reconstructed TFR ITs [N] (t, ω).

6.
Eliminate the influence of time t to obtain ITsre(ω j ) which is a result of dimensionality reduction comparing with ITs [N] (t, ω), and then derive the 2D energy-frequency map. 7.
Identify the fault pattern using the 2D energy-frequency map.
frequency axis according to the proportionality coefficient K using IMSST, and then obtaining reconstructed TFR [ ] ( , ) N ITs t ω .
6. Eliminate the influence of time t to obtain ( ) j ITsre ω which is a result of dimensionality reduction comparing with [ ] ( , ) N ITs t ω , and then derive the 2D energy-frequency map.
7. Identify the fault pattern using the 2D energy-frequency map.

Simulation Validation
In this section, we describe a multicomponent signal and a planetary gear simulation signal which are employed to illustrate the effectiveness of the proposed method.

Multicomponent Signal
To better illustrate the principle of this method, a simple time-varying multicomponent signal is established in this section, which consists of four sinusoidal functions with proportional frequencies.
where, f (t) is the IF and n(t) is the white Gaussian noise. In this signal, the sampling frequency is where, ( ) f t is the IF and ( ) n t is the white Gaussian noise. In this signal, the sampling frequency is 1000 Hz and the sampling number is 6000. Similarly, a SNR of −5 dB is added to form a synthetic signal, as shown in Figure 4a. The Fourier spectrum of the synthetic signal is shown in Figure 4b. For the influence of variable speed, the energy of Fourier spectrum is divergent and blurry, which means the useful information is hard to find. For comparison, the TFR results of the signal obtained using STFT, MSST, and IMSST are presented in Figure 4c-e. Due to the influence of the Heisenberg uncertainty principle [17], the TFR of STFT has the characteristics of a large frequency bandwidth and a low frequency resolution. Especially in noise background, the whole time-frequency plane is basically blurred, which corresponds to the result in Figure 4c. Figure 4d,e are TFRs obtained using MSST and IMSST, respectively. In comparison to Figure 4c, Figure 4d and e shows sharper time-frequency ridges. Moreover, Figure 4e transforms the time-varying instantaneous frequency ridge into a line parallel to the frequency axis on the basis of retaining the high resolution of Figure 4d, which makes the proportional relationship between the components clearer. Figure 4f  For further comparison, we draw the final results of the OT method and the proposed method in Figure 5. It should be noted that the result of the OT method is obtained through three steps: firstly, obtain an IF using the peak searching algorithm on the TFR by STFT. Secondly, apply the angle domain resampling technology to obtain the stable angle domain signal. Thirdly, draw the Fourier spectrum of the angle domain signal. As compared with the results of OT and the proposed method, we see that the latter not only highlights the feature components, but also reduces the influence of noise, which indicates that the proposed method produces clearer results.
Moreover, Figure 4e transforms the time-varying instantaneous frequency ridge into a line parallel to the frequency axis on the basis of retaining the high resolution of Figure 4d, which makes the proportional relationship between the components clearer. Figure 4f is a 2D energy-frequency map derived from the TFR of Figure 4e which shows the four proportional frequencies: 20.33 Hz, 40.67 Hz, 61 Hz, and 81.67 Hz, which correspond to the integral times of (0) 20.32 f = Hz basically.
There is no doubt that Figure 4f is more readable than Figure 4e. For further comparison, we draw the final results of the OT method and the proposed method in Figure 5. It should be noted that the result of the OT method is obtained through three steps: firstly, obtain an IF using the peak searching algorithm on the TFR by STFT. Secondly, apply the angle domain resampling technology to obtain the stable angle domain signal. Thirdly, draw the Fourier spectrum of the angle domain signal. As compared with the results of OT and the proposed method, we see that the latter not only highlights the feature components, but also reduces the influence of noise, which indicates that the proposed method produces clearer results.

The Simulated Planetary Gear Signal under Time-Varying Rotating Speed
To evaluate the effectiveness of the proposed method for analyzing planetary gear vibration signals, the simulated planetary gear signals are generated according to [35]:

The Simulated Planetary Gear Signal under Time-Varying Rotating Speed
To evaluate the effectiveness of the proposed method for analyzing planetary gear vibration signals, the simulated planetary gear signals are generated according to [35]: where, f (r) s (t) is the instantaneous rotational frequency of sun gear and f ss (t) is the instantaneous fault characteristic frequency (IFCF). f m (t) denotes the instantaneous meshing frequency of the simulated signal. A 1 and A 2 are expressed as modulation factors. θ 1 , θ 2 and θ 3 are initial phases. n(t) corresponds to the while Gaussian noise. All parameters are shown in Table 2.
Considering the variable speed, the instantaneous rotation frequency f  Figure 6a,b shows the time domain waveforms and the Fourier spectrum of the simulated signals, respectively. The TFR obtained using STFT is shown in Figure 6c, where we can roughly see that the energy is concentrated between 400 and 600 Hz. However, it is difficult to further distinguish the meshing frequency and its sideband frequency. Figure 6d is the extracted IF by TDSTFT, which is basically close to the true value (where the blue line is the extracted IF and the red line is the real IF).
Figure 6e,f are TFRs obtained using MSST and IMSST, respectively. As comparing with Figure 6d, they are clearer to find the meshing frequency and its sideband frequency. Figure 6a,b shows the time domain waveforms and the Fourier spectrum of the simulated signals, respectively. The TFR obtained using STFT is shown in Figure 6c, where we can roughly see that the energy is concentrated between 400 and 600 Hz. However, it is difficult to further distinguish the meshing frequency and its sideband frequency. Figure 6d is the extracted IF by TDSTFT, which is basically close to the true value (where the blue line is the extracted IF and the red line is the real IF). Figure 6e,f are TFRs obtained using MSST and IMSST, respectively. As comparing with Figure 6d, they are clearer to find the meshing frequency and its sideband frequency. To better illustrate the advantage of the method, we respectively generate the 2D energy-frequency maps by SST and 4-SST base on the same algorithmic theory in Sections 2.2 and 2.3, in Figure 7. Due to the aliasing of each component, it is difficult to capture useful peaks from the graph. This is mainly because two methods cannot provide enough clear TFRs, and therefore it results in blurred 2D energy-frequency maps. Figure 8 is the 2D energy-frequency map using IMSST (N = 20) of a planetary gear simulated signal, which clearly shows the meshing frequency and its corresponding sideband. In comparison to Figure 6b, the phenomenon of frequency divergence has been greatly improved. Furthermore, the proportional relationship among the components exhibits more clearly in the energy-frequency map as comparing with those in Figure 7. Moreover, the results perfectly realize the transformation from non-stationary signal to stationary signal, which is also convenient for fault diagnosis. corresponding sideband. In comparison to Figure 6b, the phenomenon of frequency divergence has been greatly improved. Furthermore, the proportional relationship among the components exhibits more clearly in the energy-frequency map as comparing with those in Figure 7. Moreover, the results perfectly realize the transformation from non-stationary signal to stationary signal, which is also convenient for fault diagnosis.   has been greatly improved. Furthermore, the proportional relationship among the components exhibits more clearly in the energy-frequency map as comparing with those in Figure 7. Moreover, the results perfectly realize the transformation from non-stationary signal to stationary signal, which is also convenient for fault diagnosis.   For further comparison, we also draw the results of the OT method and the proposed method in Figure 9 which shows that the proposed method is better for highlighting fault characteristic frequency and eliminating background noise.
To further testify the robustness of the proposed method, more simulated signals under different noise levels are analyzed. Figure 10 presents the TFRs using IMSST and the corresponding 2D energy-frequency maps of the simulated signals with SNR of 0 dB, −5 dB, and −10 dB, respectively. With the enhancement of noise, it becomes more and more difficult to directly recognize the interest components from TFR using IMSST. Correspondingly, the peak value of characteristic frequencies becomes less prominent in the 2D energy-frequency maps. Therefore, when the noise SNR of a simulated signal is greater than −5 dB, the proposed method performs well with respect to fault feature recognition. If the noise level is too high, the effectiveness of the proposed method is not ideal. It is worth noting that, when the signal noise is very strong, it can be combined with the existing singular spectrum analysis (SSA), convex optimization, time-varying filter, and other denoising methods to reduce the noise level of the signal. For further comparison, we also draw the results of the OT method and the proposed method in Figure 9 which shows that the proposed method is better for highlighting fault characteristic frequency and eliminating background noise. To further testify the robustness of the proposed method, more simulated signals under different noise levels are analyzed. Figure 10 presents the TFRs using IMSST and the corresponding 2D energy-frequency maps of the simulated signals with SNR of 0 dB, −5 dB, and −10 dB, respectively. With the enhancement of noise, it becomes more and more difficult to directly recognize the interest components from TFR using IMSST. Correspondingly, the peak value of characteristic frequencies becomes less prominent in the 2D energy-frequency maps. Therefore, when the noise SNR of a simulated signal is greater than −5dB, the proposed method performs well with respect to fault feature recognition. If the noise level is too high, the effectiveness of the proposed method is not ideal. It is worth noting that, when the signal noise is very strong, it can be combined with the existing singular spectrum analysis (SSA), convex optimization, time-varying filter, and other denoising methods to reduce the noise level of the signal.

Experimental Validation
As we all know, the simulated signal is usually a vibration model based on the meshing characteristics of gearbox which is a simplified AM-FM model. The real acquisition signal contains not only the coupling vibration information of key components such as gears, bearings, and so on, but also the other interference information caused by background noise and vibration of other

Experimental Validation
As we all know, the simulated signal is usually a vibration model based on the meshing characteristics of gearbox which is a simplified AM-FM model. The real acquisition signal contains not only the coupling vibration information of key components such as gears, bearings, and so on, but also the other interference information caused by background noise and vibration of other components of the mechanical equipment. It is more complex than the simulation signal. Therefore, it is necessary to verify the proposed method by experiment.

Experimental Rig
In this section, the proposed method is applied to detect the gear damage. The experiment is carried out on the test rig shown in Figure 11. The experimental rig is composed of a driving motor, planetary gear box, helical gearbox, magnetic powder brake, and other components. The internal structure of the test rig is shown in Figure 12. The test fault is located on the sun gear of planetary gearbox, as shown in Figure 13. A PCB accelerometer is mounted on the top of the planetary gearbox to measure the vibration signal.        In this experiment, the rotational speed of the motor decreases linearly from 3000 r/min with a sampling frequency of 2560 Hz and a sampling number of 7680. The applied torque by the magnetic powder brake is 0. The parameters of the planetary gearbox are shown in Table 3, and the corresponding fault characteristic frequencies are calculated with reference to [36], as shown in Table 4.   Figure 14a,b are the time domain waveform and the Fourier spectrum of the fault vibration signals of the sun gear, respectively. As shown in Figure 14a, the impact decreases with a decrease of rotational speed. Due to the influence of variable speed, it is still difficult to identify the fault type directly from the Fourier spectrum. Figure 14c shows the TFR using STFT, where we can roughly see the changing trend of the vibration signal. f IF (t) is the IF extracted using TDSTFT (see Figure 14d), which shows that the speed of the vibration signal basically decreases in a linear way, which is also in line with our Introduction.
According to the extracted IF, select f IF (0) as the molecule of the proportional coefficient K, and then we obtain K = f IF (0) f IF (t) . At this time, the variable speed signal will be converted to a constant frequency according to the proportional relationship. Then the rotation frequency of the sun gear in the energy-frequency map is f (r) s = 23.5 Hz, and the other frequencies are calculated according to Table 4. rectly from the Fourier spectrum. Figure 14c shows the TFR using STFT, where we can roughly s e changing trend of the vibration signal.
( ) IF f t is the IF extracted using TDSTFT (see Figure 14d hich shows that the speed of the vibration signal basically decreases in a linear way, which is al line with our Introduction. At this time, the variable speed signal will be converted to s − f c ), the above frequency can also be written as f m ± k f c + (n ± k/N) f ss . However, manufacturing error is inevitable in the actual gearbox, which means the three planetary gears meshing with the sun gear cannot be exactly identical, leading to the difference of the impact sequence between the fault teeth of the sun gear and three planetary gears. If these impacts are considered as three different sequences, the FCF of the sun gear can be written as f ss /N. Correspondingly, the peak value of the sideband appears at frequency f m ± k f c + ((n ± k)/N) f ss . Figure 15a is the final 2D energy-frequency map, where the rotational frequency f  Figure 15b is an enlarged view of Figure 13 near the frequency f m and shows that the corresponding peak values are related to one-third times the characteristic frequency of the sun gear fault characteristic frequency, which could be caused by the manufacturing or installation errors of the three planetary gears. However, this does not prevent us from judging that the fault of the planetary gearbox is the broken teeth of sun gear. frequency and shows that the corresponding peak values are related to one-third times the characteristic frequency of the sun gear fault characteristic frequency, which could be caused by the manufacturing or installation errors of the three planetary gears. However, this does not prevent us from judging that the fault of the planetary gearbox is the broken teeth of sun gear.

Conclusions
In this paper, a novel fault diagnosis method for time-varying signals was proposed, and the performance of this method has been verified using simulation signals and a planetary gearbox experiment. Its advantages are as follows: Firstly, this method does not need to estimate the instantaneous shaft frequency. It only needs to estimate an arbitrary instantaneous frequency with time-varying characteristics, such as gear meshing frequency, which is easier than extracting the instantaneous shaft frequency.
Secondly, this method can smooth the time-varying frequency ridges and achieve the effect similar to GD and OT. Moreover, this method omitted the steps of filtering and resampling, so that it can retain the complete information of the vibration signal.
Finally, this method can generate a 2D energy-frequency map, which is similar to the spectrum map. In addition, in comparison to direct analysis on time-frequency plane, this form is more intuitive and friendly.