A Time–Frequency Correlation Analysis Method of Time Series Decomposition Derived from Synchrosqueezed S Transform

Traditional correlation analysis is analyzed separately in the time domain or the frequency domain, which cannot reflect the time-varying and frequency-varying characteristics of non-stationary signals. Therefore, a time–frequency (TF) correlation analysis method of time series decomposition (TD) derived from synchrosqueezed S transform (SSST) is proposed in this paper. First, the two-dimensional time–frequency matrices of the signals is obtained by synchrosqueezed S transform. Second, time series decomposition is used to transform the matrices into the two-dimensional time–time matrices. Third, a correlation analysis of the local time characteristics is carried out, thus attaining the time–frequency correlation between the signals. Finally, the proposed method is validated by stationary and non-stationary signals simulation and is compared with the traditional correlation analysis method. The simulation results show that the traditional method can obtain the overall correlation between the signals but cannot reflect the local time and frequency correlations. In particular, the correlations of non-stationary signals cannot be accurately identified. The proposed method not only obtains the overall correlations between the signals, but can also accurately identifies the correlations between non-stationary signals, thus showing the time-varying and frequency-varying correlation characteristics. The proposed method is applied to the acoustic signal processing of an engine–gearbox test bench. The results show that the proposed method can effectively identify the time–frequency correlation between the signals.


Introduction
When analyzing multi-input and multi-output vibration systems, most of the inputs are thought to be independent and to have no influence on each other. However, due to the physical and other connections in the system, mutual influence is inevitable, and a correlation inevitably exists [1]. In near-field measurements, the main signals measured can usually ignore the influence of others. However, when there are different sources with the same frequencies, such as the multiples of the fundamental frequency, correlation analysis becomes necessary. Moreover, blind source separation (BSS) [2] methods, such as independent component analysis (ICA) [3] and fast independent component analysis, have attracted more and more attention in the field of correlation analysis, especially in the cases of medium and strong correlations [4]. Thus, in order to separate vibration and noise sources more accurately, it is necessary to perform correlation analysis.
In the time domain, correlation analysis mainly uses an auto-correlation function and a cross-correlation function to denoise and filter signals and to describe the correlation delay of random By using the inverse Fourier transform, the SSST can be transformed from the TF domain into the time-time domain, and a time series decomposition (TD) [18,19] has been derived to achieve this.
In this paper, a time series decomposition algorithm derived from the SSST is proposed, and TF correlation is analyzed in the TF plane, thus obtaining the correlations varying with time and frequency. In Section 2, we derive formulas for the SSST and the TD, respectively. Then, the main steps of the time-frequency correlation analysis are presented. In Sections 3 and 4, we apply traditional and the SSST-TD correlation analysis to both synthetic and experimental data and compare the corresponding results. We present our conclusions in Section 5.

S Transform (ST)
Before presenting with the details of the synchrosqueezed S transform, the S transform is first introduced. The ST is a signal processing method proposed by Stockwell [20]. For x(t), the ST is defined as follows: The basic wavelet function in the ST is defined as follows: where x(t) represents a signal, ST x represents the ST of x(t), f denotes frequency, t denotes time, b denotes the time shift, and i is the imaginary unit.
The ST is invertible. The inverse transform corresponding to the ST is as follows: The ST is similar to the STFT but with a Gaussian window whose width scales inversely and whose height scales linearly with frequency. The ST can obtain a varying TF resolution of the signals and inherits the advantages of both the STFT and the WT.

Synchrosqueezed S Transform (SSST)
The ST of x(t) is defined in Equation (1). Thus, the following is true: According to the Parseval theorem, Equation (1) can also be expressed as follows: where ξ is the angular frequency, X(ξ) denotes the FT result of the signal x(t), and ψ(ξ) denotes the complex conjugate of the FT result of ψ(t). In particular, if ξ < 0, X(ξ) = 0. For a purely harmonic wave x(t) = A cos(2π f 0 t), its FT is presented as follows: where A and f 0 denote, respectively, the amplitude and frequency.
Then, we can obtain the instantaneous frequency of the signal x(t) as follows: For the harmonic wave x(t), the instantaneous frequency can be obtained by the following: that is, the instantaneous frequency of the standard sinusoidal signal is the frequency of the signal itself. Given the candidate instantaneous frequency f ( f , b) at each location ( f , b) on the spectrogram, the TF spectrogram now becomes the time-instantaneous frequency spectrogram, i.e., locations ( f , b) are mapped to ( f ( f , b),b) [21]. The next step is to reassign (squeeze) the value of the time-instantaneous frequency spectrogram.
Note that the ST x ( f , b) is only computed at discrete values f k , and frequencies on the "squeezed" TF spectrogram are f l . Hence, we squeeze the TF coefficients ST x ( f , b) together via the same instantaneous frequency f ( f , b). By superimposing the time-frequency spectra values among the on the frequency point f l , we then obtain a sharpened TF representation [22,23], thus getting the SSST. For signal x(t), the SSST [24] is defined as follows: where f k is the discrete frequency obtained by the ST; ∆ f k = f k − f k−1 ; f l denotes the center frequency of the SSST; ∆ f l denotes the bandwidth; and ∆ f l = f l − f l−1 . Equation (10) shows that the SSST has higher accuracy in its TF decomposition ability. Similar to the ST, the signal can be reconstructed from the SSST x [25]. The signal can be reconstructed as follows: where RE means the real part and Details of the derivation can be found in the Appendix A.

Time Series Decomposition Algorithm (TD)
In order to transform the multi-correlation analysis from the one-dimensional time domain to the TF domain, it is necessary to firstly transform the signal into a two-dimensional time domain. A time series decomposition algorithm (TD) can be derived from the SSST to achieve this. If the SSST is transformed by the inverse FT, a time series decomposition is obtained. For the TF spectrum attained by the SSST, if the time b in the TF domain is fixed but the frequency f is variable, then the TD is as follows: where t' denotes the time attained by the inverse FT on the frequency f.
Similar to the SSST, the TD is reversible: As the elements of the SSST are complex valued numbers, the elements of the TD are also complex valued. The TD can transform an arbitrary time series into a set of time-limited, localized two-dimensional time series composed of each component, which can highlight the local feature of the signal. The TD decomposes the signal into N-components, and the sum amplitude of these components is equal to that of the original signal. Therefore, in the process of decomposition, the amplitude of each decomposed component will be reduced compared with the original, but each component can highlight the local characteristics.

Correlation Analysis
In traditional correlation analysis, the ordinary coherence is the basis of the correlation analysis. The partial coherence is the ordinary coherence between the conditional signals, where the conditional signals mean the casual effects of the others removed in a linear least square sense and the virtual coherence is the ordinary coherence between a signal and a principle component. For x and y, the magnitude-squared coherence function (MSCF) [26] is defined as follows: where G xx ( f ) and G yy ( f ) denote the power spectral density and G xy ( f ) denotes the cross spectral density of x and y. The MSCF is a function of the frequency with values between 0 and 1. These values indicate how well the two signals correspond to each other at different frequencies. The greater the values, the stronger the correlation. The traditional correlation method is inaccurate for non-stationary signals. Thus, a time-frequency correlation method is proposed to process such non-stationary signals. As the TD is a two-dimensional time series matrix, if we fix time in the b-axis in the TF domain and compute time in the t'-axis using Equation (15), then we can get a local correlation between the signals at each time b and frequency f, and that is the synchrosqueezed S transform-time series decomposition time-frequency (SSST-TD TF) correlation, an extension to the traditional correlation analysis.

Steps of the SSST-TD TF Correlation Analysis
To obtain the local correlation between the signals, a synchrosqueezed S transform-time series decomposition time-frequency correlation method is proposed. The main steps are as follows: (1) For the time domain signal x and y, obtain the time-frequency spectrum by the SSST.
(2) Perform time series decomposition on the time-frequency spectrum. That is, the inverse Fourier transform is applied to the frequency-axis of the SSST to obtain the two-dimensional complex valued time series matrix.

Simulation of Cosine Signals
To verify the SSST-TD TF correlation method, a simulation was conducted. Because stationary signals, particularly the cosine ones, are a special case of the signals, analyzing stationary signals can show the advantage and stability of the proposed method. Thus, we took two cosine signals. Signal x was cosine with frequency 30 Hz, and signal y 1 was the sum of 30 Hz and 100 Hz. The time resolution was 0.1 ms, and the sampling points were 4096. The traditional correlation coefficient (directly computed in the time domain) between x and y 1 is shown in Figure 1.

192
Due to the existence of a window function and the time resolution, a correlation coefficient has Due to the existence of a window function and the time resolution, a correlation coefficient has a certain bandwidth. Here, the correlation coefficient was 1 near 30 Hz, which means the two signals were correlated. Moreover, it was below 0.2 near 100 Hz, which is considered uncorrelated. At other frequencies, the correlation coefficient fluctuated greatly.
Then, we applied the ST and the SSST to the signals x and y 1 . The magnitude of the TF spectrum is shown in Figure 2. Both the ST and the SSST obtained signal distributions in the TF domain. Moreover, TF representations gave an intuitive display of the signals. Comparing Figure 2A-D, we can conclude that the TF energy bandwidth of the ST was much wider than that of the SSST. The existence of the false spectral energies is likely to cause mistakes in interpreting the results of the spectral decomposition, especially in the wide energy bandwidth around 100 Hz. Using the SSST, the TF resolution was highly improved, and the energy concentrated where the actual frequencies were located.     Next, we decomposed the TF spectrum of the ST and the SSST by the TD, that is, the ST-TD and the SSST-TD, and carried out a correlation analysis. The correlation coefficient is shown in Figure 3.
Obviously, near 30 Hz, both the ST-TD and the SSST-TD detected the frequency correlation between the signals. However, the bandwidth of the ST-TD correlation was wider than that of the SSST-TD, even worse than the traditional method. Near 100 Hz, the result of the ST was inaccurate due to the low TF resolution; that is, the ST was suitable for acquiring the TF distribution of signals but not suitable for calculating the correlation coefficient. Thus, the SSST-TD TF correlation analysis not only obtained the same result as the traditional method, but also obtained the TF correlation distribution. Figure 4 shows three SSST-TD TF correlation coefficients at different times. It is obvious that for the cosine signals, the overall trend of the SSST-TD TF correlation coefficient is consistent with that of the traditional method with small differences at the different frequencies. The main difference is that the SSST-TD TF correlation coefficient changes with time.

230
Considering that the SSST-TD TF correlation changes with time, the phase change must be 231 considered. In order to analyze the influence of the phase, the phase of is advanced by π and π/2: 232 = cos(30 × 2π ) + cos(100 × 2π + 2 ) = cos(30 × 2π ) + cos(100 × 2π + ) According to the previous steps, the three traditional correlation coefficient curves are shown in

230
Considering that the SSST-TD TF correlation changes with time, the phase change must be 231 considered. In order to analyze the influence of the phase, the phase of is advanced by π and π/2: 232 = cos(30 × 2π ) + cos(100 × 2π + 2 ) = cos(30 × 2π ) + cos(100 × 2π + ) According to the previous steps, the three traditional correlation coefficient curves are shown in

Simulation of Different Phase Signals
Considering that the SSST-TD TF correlation changes with time, the phase change must be considered. In order to analyze the influence of the phase, the phase of y 1 is advanced by π and π/2: According to the previous steps, the three traditional correlation coefficient curves are shown in Figure 5. Under the influence of the phase change, the traditional correlation coefficient was consistent near 30 Hz and 100 Hz, with great differences at other frequencies.
= cos(30 × 2π ) + cos(100 × 2π + 2 ) = cos(30 × 2π ) + cos(100 × 2π + ) According to the previous steps, the three traditional correlation coefficient curves are shown in Using the SSST-TD TF correlation analysis, the correlation coefficients of x and y 2 and x and y 3 are shown in Figure 6. It can be seen that the overall correlation coefficient trend was consistent with the original phase with small differences at different times. Taking 63 Hz and 100 Hz as an example, the time domain curves of the correlation coefficients, after removing the two ends, are shown in Figure 7. It can be seen that at 100 Hz, the correlation coefficient did not follow the phase change, remaining unchanged, while at 63Hz, the correlation coefficient varied with the phase change. Therefore, compared with traditional correlation analysis, the SSST-TD TF correlation not only reflected the correlation changing with time, but also retained the correlation changing with the phases. At a high frequency, the phase changes quickly. However, when the frequency is low, the phase change will play an important role and cannot be ignored.

Simulation of Signals With Noise
Next, we added Gaussian noise to x and y 1 with a mean value of 0 and a standard variance of 1/3. The traditional correlation coefficient is shown in Figure 8A. It is clear that under noisy conditions, traditional correlation can recognize the correlation between x and y 1 near 30 Hz and the uncorrelation near 100 Hz. Due to the existence of noise, the correlation coefficients changed at other frequencies.

263
The correlation coefficient obtained by SSST-TD is shown in Figure 8B. The correlation coefficient obtained by SSST-TD is shown in Figure 8B. The overall correlation coefficient trend was consistent with the noiseless one. The correlation near 30 Hz and uncorrelation near 100 Hz can be accurately identified in the case of noise. Figure 9 shows the correlation coefficients at different times. Except at 30 Hz and 100 Hz, the correlation coefficients at the other frequencies fluctuated under the influence of noise. As the proposed method is a time-frequency extension to the ordinary coherence, the method inherits the characteristics of traditional correlation; that is, under the influence of noise, both the traditional correlation analysis and the proposed analysis are influenced.

281
The time domain signals are shown in Figure 10. The traditional correlation coefficient is shown 282 in Figure 11. Obviously, it is impossible to distinguish the correlation between and . For the
x 0 4 y Figure 9. The correlation coefficient between x and y 1 at different times.

Simulation of Non-Stationary Signals
For stationary signals with or without noise, both the traditional correlation analysis and the SSST-TD TF correlation analysis can identify the correlation between the signals. However, the SSST-TD TF correlation method can better identify the correlation at different times. Practical signals are often non-stationary, so the following signals were designed to analyze the advantages of the proposed method. They consisted of cosine, cosine frequency modulation (FM), and shock. We also added Gaussian noise to both with a mean value of 0 and a standard variance of 1. The time resolution was 1 ms, and the sampling points were 4096.
The time domain signals are shown in Figure 10. The traditional correlation coefficient is shown in Figure 11. Obviously, it is impossible to distinguish the correlation between x and y. For the frequency 30 Hz, it occurred at the second and third stage in x and only the third stage in y, so the correlation coefficients was expected to change and not be constant. Thus, the correlation coefficient was not accurate.
TD TF correlation method can better identify the correlation at different times. Practical signals are 277 often non-stationary, so the following signals were designed to analyze the advantages of the 278 proposed method. They consisted of cosine, cosine frequency modulation (FM), and shock. We also 279 added Gaussian noise to both with a mean value of 0 and a standard variance of 1. The time resolution 280 was 1 ms, and the sampling points were 4096.

281
The time domain signals are shown in Figure 10. The traditional correlation coefficient is shown 282 in Figure 11. Obviously, it is impossible to distinguish the correlation between and . For the 283 frequency 30 Hz, it occurred at the second and third stage in and only the third stage in , so the 284 correlation coefficients was expected to change and not be constant. Thus, the correlation coefficient 285 was not accurate.

281
The time domain signals are shown in Figure 10. The traditional correlation coefficient is shown 282 in Figure 11. Obviously, it is impossible to distinguish the correlation between and . For the 283 frequency 30 Hz, it occurred at the second and third stage in and only the third stage in , so the 284 correlation coefficients was expected to change and not be constant. Thus, the correlation coefficient 285 was not accurate.

290
The SSST of and is shown in Figure 12. It is clear that the TF distribution of and can 291 be attained by the SSST, and the energy is gathered at the actual frequency. Moreover, the SSST can 292 Figure 11. The traditional correlation coefficient between x and y.
The SSST of x and y is shown in Figure 12. It is clear that the TF distribution of x and y can be attained by the SSST, and the energy is gathered at the actual frequency. Moreover, the SSST can obtain the time-varying distribution more accurately. The correlation coefficient obtained by the SSST-TD TF correlation analysis is shown in Figure 13A.
Comparing the SSST spectrum with the correlation coefficient, the signals were correlated in the first and third stages. The 30 Hz uncorrelation was also obtained in the second stage. The correlation coefficient of the frequency modulation was obtained as well, as can be seen in the local, enlarged detail in Figure 13B. The outline was clear, and this cannot be attained by traditional correlation analysis. For the non-stationary signals, traditional correlation analysis cannot obtain the accurate correlation. Similarly, the correlation coefficients were affected by noise. The SSST-TD TF correlation analysis not only obtained the correlation between the signals, but also the correlation as it changed with time and frequency. Thus, by analyzing the correlation in the TF plane, the correlations between signals are more accurate and specific.  Figure 13A.

294
Comparing the SSST spectrum with the correlation coefficient, the signals were correlated in the 295 first and third stages. The 30 Hz uncorrelation was also obtained in the second stage. The correlation 296 coefficient of the frequency modulation was obtained as well, as can be seen in the local, enlarged 297 detail in Figure 13B.     Through a time series decomposition algorithm derived from the SSST, the SSST-TD TF correlation was carried out in the TF plane, and the local correlation between the signals was attained. Compared with the traditional correlation analysis method using sinusoidal signals, sinusoidal signals with noise, and non-stationary signals, the proposed method was validated by stationary and non-stationary signals simulation. The simulation results show that the traditional method can obtain the overall correlation between the signals but cannot reflect the local time and frequency correlations, especially the correlations between the non-stationary signals. The proposed method not only obtained the overall correlations between the signals, but also accurately identified the correlations between the non-stationary signals, revealing more specific and essential correlations between the signals. Moreover, by comparing the ST and the SSST, it can be concluded that a more energy-concentrated TF result denotes a better TF location and a better characterization of the time-varying feature.

Experimental Equipment
In order to verify the SSST-TD TF correlation analysis, experiments were carried out on an engine-gearbox test bench. The test bench was composed of a four-cylinder engine, a gearbox, a brake, a base, etc., as shown in Figure 14. This test bench was designed for identifying vibration and noise sources and the fault diagnosis of gearboxes, bearings, and shafts driven by engines or motors. Thus, the noise signals were computed.

324
In order to verify the SSST-TD TF correlation analysis, experiments were carried out on an 325 engine-gearbox test bench. The test bench was composed of a four-cylinder engine, a gearbox, a 326 brake, a base, etc., as shown in Figure 14. This test bench was designed for identifying vibration and 327 noise sources and the fault diagnosis of gearboxes, bearings, and shafts driven by engines or motors.

328
Thus, the noise signals were computed.

332
The gearbox was a primary planetary gearbox, and the details are shown in Table 1. The sun 333 gear was connected to the engine with the ring gear fixed. The operating condition was 2600 rpm.

334
The sensors were two sound pressure sensors. The sampling frequency was 5.12 kHz, and the 335 sampling points were 4096.  The gearbox was a primary planetary gearbox, and the details are shown in Table 1. The sun gear was connected to the engine with the ring gear fixed. The operating condition was 2600 rpm.

336
The sensors were two sound pressure sensors. The sampling frequency was 5.12 kHz, and the sampling points were 4096.

Experimental Signal Analysis
To reduce the influence between the sources, the signals were measured near-field. The measured near-field sound pressures of the gearbox x and the engine y are shown in Figure 15. The FT of x and y are shown in Figure 16, and the traditional correlation coefficient was calculated as shown in Figure 17. The most common signals of the rotating machinery were the fundamental frequency and the multiples of the fundamental frequency. As can be seen in Figure 16

332
The gearbox was a primary planetary gearbox, and the details are shown in Table 1. The sun 333 gear was connected to the engine with the ring gear fixed. The operating condition was 2600 rpm.

334
The sensors were two sound pressure sensors. The sampling frequency was 5.12 kHz, and the 335 sampling points were 4096.

338
To reduce the influence between the sources, the signals were measured near-field. The 339 measured near-field sound pressures of the gearbox and the engine are shown in Figure 15. The

340
FT of and are shown in Figure 16, and the traditional correlation coefficient was calculated as 341 shown in Figure 17. The most common signals of the rotating machinery were the fundamental 342 frequency and the multiples of the fundamental frequency. As can be seen in Figure 16, the amplitude     Then, we applied the ST and the SSST to x and y, and the TF spectra obtained are shown in Figure 18. As can be seen, the ST energy below 500 Hz was mainly gathered at the frequencies around 86 Hz and 173.65 Hz, which corresponds well with multiple rotating frequencies. For the SSST, the bandwidth around 86 Hz and 173.65 Hz was much smaller, showing the advantages of the SSST. The energy distributions on the TF spectrum were concentrated around the real instantaneous frequency of the signal. Then, we carried out the SSST-TD TF correlation analysis. The correlation coefficient is shown in Figure 19.   continuous. The result shows consistency with that of the traditional method, as shown in Figure  Figure 19. The SSST-TD TF correlation coefficient between x and y.
Taking 86.25 Hz and 173.75 Hz as an example. From the ST and the SSST, it can be seen that the signals were basically continuous at 86. 25 Hz, so the correlation coefficient was also expected to be continuous. The result shows consistency with that of the traditional method, as shown in Figure 20A. However, there were many discontinuities and fluctuations near 173.65 Hz, which means the correlation coefficient was not as continuous as that at 86. 25 Hz, as shown in Figure 20B. The other main frequencies can also be analyzed similarly. signals were basically continuous at 86. 25 Hz, so the correlation coefficient was also expected to be 373 continuous. The result shows consistency with that of the traditional method, as shown in Figure   374 20A. However, there were many discontinuities and fluctuations near 173.65 Hz, which means the 375 correlation coefficient was not as continuous as that at 86. 25 Hz, as shown in Figure 20B. The other 376 main frequencies can also be analyzed similarly.  As previously analyzed, the correlation coefficient varied with the phase. The ST/SSST time frequency spectrum should be considered to judge the correlation between the signals. Therefore, the traditional method is the overall correlation coefficient, which cannot reflect the time-varying characteristics, while the SSST-TD TF correlation method is able to obtain the time-varying correlation and obtain a better local correlation between the signals.

Conclusions
A more energy-concentrated TF result denotes a better TF location and better characterization of time-varying features. In this paper, a time series decomposition algorithm derived from the SSST was proposed, and TF correlation was analyzed in the TF plane. The local correlation between the signals in the TF domain were attained, and the correlations varying with time and frequency were able to reveal more specific and essential correlations between the signals. By comparing the traditional correlation analysis with the proposed method using sinusoidal signals, sinusoidal signals with noise, and non-stationary signals, the superiority of the proposed method is that the correlation can be obtained more intuitively and accurately. Similar to the traditional method, the SSST-TD TF correlation method is also affected by noise. The influence of noise must be further studied. As an extension to the traditional correlation analysis, the proposed method not only obtains the correlation between the signals in accordance with the traditional correlation analysis method, but also obtains the time-frequency correlation distribution and can be used for non-stationary signals. Thus, the following major conclusions can be drawn. 1 For stationary signals, both the traditional correlation and the SSST-TD TF correlation methods can identify the correlations between signals, and the SSST-TD TF correlation method can also obtain the correlation as it changes with time and frequency. 2 Noise and the phase will affect the correlation coefficient. The traditional correlation method cannot recognize the influence of the phase, but the SSST-TD TF correlation method can get the effect of the phase of the cosine signals. 3 For time-varying signals, the correlation between the signals cannot be recognized by the traditional correlation method. The SSST-TD TF correlation method can obtain the correlation of non-stationary signals and can also reflect the time-varying and frequency-varying characteristics between the signals. Acknowledgments: The authors wish to thank E. Brevdo for his MATLAB Synchrosqueezing Toolbox, and three anonymous reviewers for their constructive comments and suggestions.

Conflicts of Interest:
No conflicts of interest exist in the submission of this manuscript, and the manuscript has been approved by all the authors for publication. The work described herein is original research that has not been published previously and is not under consideration for publication elsewhere, in whole or in part. All the authors listed have approved the manuscript. The founding sponsors had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

Derivation of the Formula of Inverse SSST
The definition of the SSST is as follows: According to Parseval theorem, the ST of the signal can be expressed as follows: The following is true: We found the following: According to Equation (5) and using the continuous form of the right side of Equation (A4), we found the following: −ψ(ξ)ξ −1 dξ and assuming that the signal is real-valued, the signal can then be reconstructed by the following: In the piecewise constant approximation corresponding to the binning in f, the following is true: