Synchroextracting transform: The theory analysis and comparisons with the synchrosqueezing transform

Signal Processing 166 (2020) 107243

Contents lists available at ScienceDirect

Signal Processing journal homepage: www.elsevier.com/locate/sigpro

Synchroextracting transform: The theory analysis and comparisons with the synchrosqueezing transform Zhen Li a, Jinghuai Gao a,∗, Hui Li a, Zhuosheng Zhang b, Naihao Liu a, Xiangxiang Zhu b a

School of Information and Communications Engineering; National Engineering Laboratory for Offshore Oil Exploration, Xi’an Jiaotong University, Xi’an, Shaanxi 710049, China b School of Mathematics and Statistics, Xi’an Jiaotong University, Xi’an, Shaanxi 710049, China

a r t i c l e

i n f o

Article history: Received 16 December 2018 Revised 25 May 2019 Accepted 29 July 2019 Available online 30 July 2019 Keywords: Synchroextracting transform Fourier-based synchrosqueezing transform Theoretical analysis Synchroextracting operator Fixed frequencies

a b s t r a c t In this paper, we consider the theoretical analysis and new view of a time-frequency analysis method termed synchroextracting transform (SET), which is inspired by the Fourier-based synchrosqueezing transform (FSST) and the ideal time-frequency analysis theory. We first review the corresponding definitions of the FSST and SET. Then, we provide the theoretical analysis of the SET and study some properties with mathematical proofs. We show that the SET mainly benefits from the synchroextracting operator (SEO) that can only retain the TF information most related to time-varying features of the signal. We also give a new view that SEO is equivalent to the fixed frequencies of the frequency estimation operator of the FSST. Based on the SEO, SET greatly improves the concentration of the time-frequency representation and perfectly localizes the time-frequency features of chirp signals. Finally, numerical simulations are employed to illustrate the performance comparisons of the SET and FSST. © 2019 Elsevier B.V. All rights reserved.

1. Introduction The time-frequency analysis (TFA) technology plays a significant role in signal processing and interpretation. The TFA method can map a 1D signal into a time-frequency representation (TFR) that provides the information of time-varying features. Various TFA methods have been designed for well analysis of non-stationary signals recently. The classical linear TFA approaches, such as the short-time Fourier transform (STFT) [1] and the wavelet transform (WT) [2], have been extensively used [3]. However, these traditional methods are restricted by the Heisenberg uncertainty principle [4], namely, the optimal time and frequency resolutions cannot be achieved simultaneously. Consequently, the classical linear TFA methods usually generate a “blurred” TFR, failing to characterize TF features of non-stationary signals accurately. To improve the quality of the TFR, especially for multicomponent signals, some efforts have been made in the past decades such as the amplitude modulated-frequency modulated (AM-FM) method [5–7], parametric TFA method [8,9], and demodulated TFA method [10,11]. Among those advanced techniques, the reassignment method [12,13], empirical mode decomposition coupled with the Hilbert spectrum

∗

Corresponding author. E-mail address: [email protected] (J. Gao).

https://doi.org/10.1016/j.sigpro.2019.107243 0165-1684/© 2019 Elsevier B.V. All rights reserved.

[14], and synchrosqueezing transform [15,16] are the three most meaningful methods. The reassignment method (RM) was applied to the spectrogram analysis in the1970’s [12], and a new formulation and theoretical analysis of it were studied [13]. It is found that the RM is an effective approach to improve the concentration of TFRs obtained using the traditional TFA methods, e.g., the linear transforms [17,18] and bilinear transforms [13,19]. The RM has been widely applied in many practicalapplications [20–22]. However, it reassigns the TF coefficients along the time and frequency directions, leading to the inability to reconstruct the signal [23]. Empirical mode decomposition, proposed by Huang et al. [14], is a powerful and data-driven signal analysis technique to decompose the signal into a series of intrinsic mode functions (IMFs). A TFR can be achieved by combining EMD with the Hilbert transform to compute instantaneous frequencies, which is sometimes called the Hilbert-Huang transform (HHT) [14]. The HHT has shown its functional capability [24–26]. However, some features encumber its direct application: mode mixing, sensitivity to noise, aliasing, and end-point artefacts [27]. Furthermore, the obtained IMFs have an increasing frequency localization at the expense of a decreasing time localization [28]. The solutions such as ensemble empirical mode decomposition (EEMD) [29] and complete ensemble empirical mode decomposition (CEEMD) [30], address these issues in practice, but introduce new challenges to our mathematical understanding. Although attempts at a mathematical understanding of

2

Z. Li, J. Gao and H. Li et al. / Signal Processing 166 (2020) 107243

EMD and its results have been mostly exploratory [27,31,32], the properties of EMD and its variants are not fully understood. In the mid-1990s, Mase and Daubechies proposed a novel method, which was termed synchrosqueezing transform (SST), to analysis auditory signals [15]. A thorough theoretical studies of the SST was provided by Daubechies et al. [16], which gives new insights in understanding the principle of EMD, and more importantly, the SST can characterize components with time-varying spectrum more clearly. The purpose and manner of the SST are very similar to that of the RM family. Combining the sparsity of the RM with the reversibility of WT, the SST can concentrate the TFR as well as reconstruct the signal. Due to such properties, the SST has attracted a lot of interest and is widely studied in theory. One interesting aspect of the SST is the extension to another TFA method. The SST was extended to the frame work of STFT [33,34]. Liu et al. proposed the synchrosqueezing three parameter wavelet transform [35,36]. Two variants of SST based on the S transform and generalized S transform were introduced, respectively [37,38]. Daubechies et al. further proposed ConceFT, for the concentration of frequency and time, via STFT-based or WT-based SST [39]. For the signal f (t ) = A(t )e2iπ φ (t ) , all these advanced developments are designed to achieve the ideal TFR (ITFR) [40,41], which is in the form

IT F R(η, t ) = f (t )δ (η − φ  (t )) = A(t )e2iπφ (t ) δ (η − φ  (t )),

(1)

where δ is the Dirac distribution, and A(t) and φ (t) respectively stand for the time-varying amplitude and phase of f(t). φ  (t), the derivative of φ (t), is the instantaneous frequency (IF). Obviously, the ITFR is the one producing the Dirac pulse at the IF φ  (t) of the signal; elsewhere the value is zero [40]. The synchrosqueezing methods mentioned above are all based on the frequency estimation operator, which provides an unbiased estimation of the IF for the weakly FM signal [16]. Therefore, the SST-based methods can obtain an intensely concentrated TFR that is toward the ITFR. However, SST is not suitable for the strongly FM signal [42,43]. To address this issue, several advanced approaches have been proposed. One important scheme is to modify the frequency estimation operator of the original SST, such as the secondorder synchrosqueezing transform based-STFT or based-WT [44– 47], and the synchro-compensating chirplet transform [48]. Herein, we focus on one novel TFA method called synchroextracting transform (SET) [49], which combines the Fourier-based synchrosqueezing transform (FSST) with the ITFR theory. The SET can improve the TF concentration only using the frequency estimation operator of the original SST. However, the theoretical analysis of the SET was not provided in [49]. This paper mainly aims to provide the theoretical studies of the SET. We first recall some related definitions. Then the approximation theory is given and several properties of the SET are derived with mathematical proofs. In the theoretical analysis, the Gaus2 − σt 2

sian window function with generalized form g(t ) = σ1 e is employed. It is well known that the SST is not suitable for the strongly FM signal. In contrast, because the introducing of the synchroextracting operator helps to perfectly localize the chirp, the TF concentration in the SET is improved greatly and the SET is also able to provide a more concentrated TFR than the SST for the strongly FM signal. In the simulation, we show that compared with the SST, the SET achieves a better signal concentration for signals with slow varying IF, especially for that with fast varying IF, as well as for higher-order polynomial phase signals and other non-stationary signals. This paper is organized as follows. In Section 2, we review the reassignment method and the SST method based-STFT. Section 3 introduces the definition of the SET, and provides the corresponding theoretical analysis. In Section 4, we show several properties of the SET. Some numerical experiments are employed

to illustrate the performance of the SET in Section 5. The discussion and conclusion are drawn in Sections 6 and 7, respectively. 2. Definitions 2.1. Short-time fourier transform and reassigned method For a given signal f(t) ∈ L1 (R), the Fourier transform is defined by

fˆ(η ) =



+∞

f (t )e−2iπ ηt dt ,

−∞

(2)

which provides the frequency information of f(t) on the whole time domain. In order to characterize the time-localized frequency descriptors, one introduces the short-time Fourier transform (STFT) that is obtained by performing the Fourier transform within a sliding window g(t) ∈ L2 (R) according to:

S (η , t ) = =



+∞

f (ξ )g∗ (ξ − t )e−2iπ ηξ dξ

−∞  +∞

f (ξ )g(ξ − t )e−2iπ ηξ dξ

−∞

= M (η, t )e2iπ ϕ (η,t ) ,

(3) g∗ (t )

where g(t) is a finite-length, real-valued (g(t ) = and even window function, which has unit norm. M(η, t) is the magnitude of S(η, t), and ϕ (η, t) is the phase. By shifting the analysed signal f(t), instead of the window function, we can calculate the modified STFT by



Sgf (η, t ) =

+∞



−∞

 =

+∞

−∞

f (ξ  + t )g∗ (ξ  )e−2iπ ηξ dξ 

(4)

f (ξ )g(ξ − t )e−2iπ η (ξ −t ) dξ

(5)

= Mgf (η, t )e2iπ f (η,t ) . g

(6)

We should note that the STFT is invertible in L2 (R), the reconstruction formulae is

f (ξ ) =



=

+∞



+∞

−∞

−∞





+∞

−∞

Sgf (η, t )g(ξ − t )e2iπ η (ξ −t ) dt dη

+∞

−∞



Mgf (η, t )g(ξ − t )e

2iπ gf (η,t )+η (ξ −t )

(7) 

dtdη.

(8)

When the time variation of M(η, t) is slower than the phase variation, those points (η, t), providing the maximum reconstruction contribution, satisfy the phase stationary condition [42]

∂t { gf (η, t ) + η (ξ − t )} = 0, ∂η { gf (η, t ) + η (ξ − t )} = 0. (9) Based on those two formulas, the definitions of the so-called reassignment operators for the STFT frame are defined as follows:

ωˆ f (η, t ) = ∂t gf (η, t ), tˆf (η, t ) = t − ∂η gf (η, t ).

(10)

According to (5), (6), and (10), the frequency reassignment operator and the time reassignment operator can be described as:



 ∂t Sgf (η, t ) ωˆ f (η, t ) = ∂t (η, t ) =  , 2iπ Sgf (η, t )   ∂η Sgf (η, t ) tˆf (η, t ) = t − ∂η gf (η, t ) = t −  . 2iπ Sgf (η, t ) g f

(11)

(12)

Then, the reassignment method (RM) improves the concentration of the time-frequency representation (TFR) by reallocatg ing the energy distribution S f (η, t ) according to the map (η, t ) →

(ωˆ f (η, t ), tˆf (η, t )). The definition of the RM based on STFT is

Z. Li, J. Gao and H. Li et al. / Signal Processing 166 (2020) 107243

RM f (ω, τ ) =



+∞



−∞

+∞

0

|Sgf (η, t )|2 δ (ω − ωˆ f (η, t ))

δ (τ − tˆf (η, t ))dηdt,

(13)

where δ ( · ) denotes the Dirac distribution. It is well-known that the RM reassigns the information along the time and the frequency directions in the TF plane. Indeed, the RM can obtain a concentrated TFR toward to the ideal one. However, it does not allow for signal reconstruction. 2.2. STFT-based synchrosqueezing transform

fined in (13). Then, the expression of the FSST is:

1 g( 0 )



+∞ −∞

Sgf (η, t )δ (ω − ω ˆ f (η, t ))dη.

(14)

f (t ) =

+∞ −∞

T f ( ω , t )d ω .

(15)

Then, we recall the theoretical results of the FSST for signals with slowly time-varying amplitude and instantaneous frequency (IF) [34,42]. The following definitions define those signals: Definition 1. (Intrinsic mode type function). Let > 0, for a given continuous mapping f: R → C, f ∈ L∞ (R), it will be named intrinsicmode-type (IMT), if the amplitude A(t) and phase φ (t) of the signal f (t ) = A(t )e2iπ φ (t ) satisfy the following conditions:

A(t ) ∈ C 1 (R ) ∩ L∞ (R ), φ (t ) ∈ C 2 (R ), supt∈R φ  (t ) < ∞, φ  (t ) > 0, ∀t ∈ R, |A (t )| ≤ |φ  (t )|, |φ  (t )| ≤ , ∀t ∈ R.

f (t ) =

fk (t ) =

k=1

N 

(3) Moreover, for all k ∈ {1, . . . , N}, there exists a constant C such that

  

  λ,ε˜  lim   C ε , ∀t ∈ R. T ( ω , t ) d ω − f ( t ) k λ→0 |ω−φ  (t )|ε˜ f 

(18)

Theorem 1 provides a strong approximation result. The TF distribution of the FSST contains non-zero coefficients localized around the curves (φk (t ), t ). Furthermore, the mode can be reconstructed via the inverse transform of FSST from a small frequency band around the IF curve φk (t ). 3. The synchroextracting transform and approximation results To further improve the quality of the TFR, many different techniques have been carried out. Among these advanced methods, we focus on the synchroextracting transform (SET), which combines the reversibility of the STFT and the sparsity of the ideal TFR (ITFR) [49]. To our knowledge, no theoretical analysis is available for this new method. We aim to bridge the gap in this section.

The SET can be performed via three successive steps. The first g step is to calculate the STFT S f (η, t ). The second step is to calculate the estimation ω ˆ f (η, t ) of the IF as follows:



ωˆ f (η, t ) =

Re ∞,



∂t Sgf (η,t ) , 2iπ Sgf (η,t )

|Sgf (η, t )| > γ , |Sgf (η, t )| ≤ γ ,

(19)

where γ > 0 is to eliminate the unstable phenomenon or the influence of the noises. Here, γ is chosen as 10−8 in the noise free numerical experiment. In a  noisy environment, the adaptive threshold γ is determined by γ = 2log2 N · σ , where N is the signal length g g and σ = median(|S f (η, t ) − median(S f (η, t ))| )/0.6745 [50].

T e f (η, t ) = Sgf (η, t )δ (η − ω ˆ f (η, t )),

φk − φk −1 ≥ 2d.

δ (η − ωˆ f (η, t )) =

 Definition 3. Let function h(t) ∈ C∞ (R) and h(t )dt = 1, and set the threshold ε˜ and the accuracy λ. The FSST of f(t) ∈ Bε,d with ε˜ and λ can be defined by:



 ω − ωˆ f (η, t ) 1 1 Sgf (η, t ) h dη. g(0 ) |Sgf (η,t )|>ε˜ λ λ

(16)

When δ and ε˜ approach to zero, we obtain the usual formula in (14). 1

Theorem 1. Consider f ∈ Bε,d and set ε˜ = ε 3 . Let g ∈ S(R), the Schwartz class, be satisfy supp(gˆ) ⊂ [−d, d]. Further, if is small enough, then the following holds:

1, η = ω ˆ f ( η , t ), 0, otherwise,

(21)

which derives the following capability of extracting:



T e f (η, t )) =

Let the B ,d denotes the set of all SIMT.

(20)

where δ (η − ω ˆ f (η, t )) is synchroextracting operator (SEO) and can be interpreted as:



Ak (t )e2iπφk (t ) ,

k=1

where all the fk (t) are IMT. And further, for all k ∈ {2, . . . , N} and all t, those IMTs are separated with separation d ∈ (0, 1) i.e.

T fλ,ε˜ (ω, t ) =

(17)

The final step is energy extraction:

Definition 2. (Superposition of well-separated intrinsic mode components). The signal f is called to be a superposition of wellseparated IMT (SIMT), if there exists a finite N such that N 

|ωˆ f (η, t ) − φk (t )|  ε˜.

3.1. The definition of SET

And it enables to retrieve the signal by considering:



(1) |S f (η, t )| > ε˜ only when, there exists some k ∈ {1, . . . , N} such that (η, t) ∈ Zk := {(η, t ), s.t . |η − φk (t )| < d}. (2) For each k ∈ {1, . . . , N} and all (η, t) ∈ Zk such that |S f (η, t )| > ε˜, we have

k

The synchrosqueezing transform (SST) was originally introduced in the continuous wavelet transform (CWT) context for audio signal analysis [15]. In contrast to empirical mode decomposition (EMD) [14], the SST provides a good approximation to the ideal TFR and allows for the reconstruction of the weakly modulated signal [16]. The aim of the SST is twofold: (a) provide an energyconcentrated TFR of multicomponent signals, (b) enable to separate and recover each mode. As an extension of the SST, the Fourierbased SST (FSST) was proposed in [34]. Starting from the STFT calculated by (5), the FSST reassigns the g TF coefficients of S f (η, t ) from (η, t) to (ω ˆ f (η, t ), t ), which is de-

T f (ω , t ) =

3

Sgf (η, t ), η = ω ˆ f ( η , t ), 0, otherwise.

(22)

According to the definitions of the FSST and SET, i.e. the Eqs. (14) and (20), the only difference between them is the operation as illustrated in Fig. 1(a). The FSST (Fig. 1(b)) reassigns the STFT coefficients around the IF, while the SET (Fig. 1(c)) only extracts those in the curve η = ω ˆ f (η, t ), existing as an exact estimation of the IF. It is visually found that the TFR of the SET is more concentrated than that of the FSST. 3.2. The existence of the fixed points According to Eq. (22), the SET only extracts coefficients of the STFT in the trace η = ω ˆ f (η, t ), i.e., the fixed points of ω ˆ f (η, t ) with

4

Z. Li, J. Gao and H. Li et al. / Signal Processing 166 (2020) 107243

∂ g ∂ S (η , t ) = ∂η f ∂η 

=



+∞

−∞ +∞

−∞

= −2iπ = −2iπ

f (ξ ) · g(ξ − t ) · e−2iπ η (ξ −t ) dξ





f (ξ ) · g(ξ − t ) · −2iπ (ξ − t )e−2iπ η (ξ −t ) dξ 

+∞

−∞ Stg f

f (ξ ) · (ξ − t ) · g(ξ − t ) · e−2iπ η (ξ −t ) dξ

( η , t ),

(28)

where t g(t ) = t · g(t ) is the time-weighted analysis window, and tg S f (η, t ) is the short-time Fourier transform computed based this time-weighted analysis window. To arrive at an expression for the partial derivative of spectral phase with respect to frequency, we take the partial derivag tive of S f (η, t ). Applying the product rule of differential calculus to Eq. (25), we have

Fig. 1. The operators of processing the coefficient. (a) The modulus of the STFT, the blue line is the IF and red (respectively green) arrows symbolizing how the operator is performed in FSST (respectively SET). (b) The TF result obtained by FSST. (c) The same as (b) but by the SET. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

∂ Mgf (η, t ) 2iπ g (η,t ) g ∂ g ∂ g f S f (η , t ) = M f (η, t )e2iπ f (η,t ) = ·e ∂η ∂η ∂η ∂ gf (η, t ) 2iπ g (η,t ) f +Mgf (η, t ) · 2iπ e ∂ω ∂ (ln Mgf (η, t )) g ∂ gf (η, t ) g = · S f (η, t ) + 2iπ · S f ( η , t ). ∂η ∂η (29) According to Eqs. (28) and (29), multiplying by 1/S f (η, t ), we g

can obtain: respect to the frequency η. Therefore, it is necessary to prove the existence of those fixed frequency points. Theorem 2. For a given signal f(t), considering a Gaussian window function g(t), it holds

∂ Mgf (η, t ) = 0. ∂η

η = ωˆ f (η, t ) ⇔

Stg (η , t ) f Sgf



+∞

(24)

g

(25)

Taking the partial derivative of Eq. (24), we have:



=



−∞

+∞

−∞



+

+∞

−∞

f (ξ ) · g(ξ − t ) · e−2iπη (ξ −t ) dξ

d g(t ) dt



 ∂ g(ξ − t ) · e−iη (ξ −t ) dξ ∂t   ∂ −2iπη(ξ −t ) f ( ξ ) · g( ξ − t ) · e dξ ∂t

3.3. The approximation results

(26)

is the time-derivative analysis window, and

is the short-time Fourier transform computed based this g

= 0), it holds 

Sgf (η, t ) ∂t Sgf (η, t ) =− + η. g 2iπ S f (η, t ) 2iπ Sgf (η, t )

g

tion but the derivative must exist, which is shown in Fig. 2. Without loss of generality, we employ the Gaussian window function for the theoretical analysis and numerical simulation in this paper.

time-derivative analysis window. Multiplying by 1/S f (η, t ) (when g S f (η , t )

g

{S f (η, t )/S f (η, t )} = 0, that is owned by the real window functg



g S f (η , t )

Remark 1. This theorem means that the maximums of the STFT modulus, i.e. the ridge points, satisfy η = ω ˆ f (η, t ). In the g

= −Sgf (η, t ) + 2iπ ηSgf (η, t ), where g (t ) =

(31)

Eq. (31), we use the property, i.e. {S f (η, t )/S f (η, t )} = 0 ⇔

f (ξ ) ·

+∞

(30)



= Mgf (η, t )e2iπ f (η,t ) .

∂ g ∂ S (η , t ) = ∂t f ∂t

∂ (ln Mgf (η, t )) ∂ gf (η, t ) − . ∂η ∂η

  ∂t Sgf (η, t ) η = ωˆ f (η, t ) ⇔ η =  2iπ Sgf (η, t ) g  tg  S f (η , t ) S f (η , t ) ⇔ =0⇔ =0 Sgf (η, t ) Sgf (η, t ) ∂ (ln Mgf (η, t )) ∂ Mgf (η, t ) ⇔ =0⇔ = 0. ∂η ∂η

(23)

f (ξ ) · g(ξ − t ) · e−2iπη (ξ −t ) dξ

−∞

1 2iπ

Finally, combining Eq. (27) with (30), we have

Proof. We rewrite the expression of the STFT:

Sgf (η, t ) =

(η , t )

=−

(27)

In the similar way, we can also calculate the derivative of Eq. (24) with respect to frequency, and then obtain:

1

Theorem 3. Consider f ∈ Bε,d and set ε˜ = ε 3 . Let g ∈ S(R), the Schwartz class, satisfy supp(gˆ) ⊂ [−d, d]. Further, if is small enough, then the following holds: g (1) |S f (η, t )| > ε˜ only when, there exist some k ∈ {1, . . . , N} such that (η, t) ∈ Zk := {(η, t ), s.t . |η − φk (t )| < d}. (2) For each k ∈ {1, . . . , N} and all (η, t) ∈ Zk |Sgf (η, t )| > ε˜, we have

|ωˆ f (η, t ) − φk (t )| ≤ ε˜.

such that

(32)

(3) Moreover, for all k ∈ {1, . . . , N}, there exists a constant C such that for all t ∈ R,

 

1 gˆ(0 )





 |η−φ  (t )|<ε˜ T e f (η, t )dη − f k (t ) ≤ C .

(33)

Z. Li, J. Gao and H. Li et al. / Signal Processing 166 (2020) 107243

80

80

80 IF

40 20

0

0.2

0.4 0.6 Time (s) (d)

0.8

1

IF

60

Frequency (Hz)

60

Frequency (Hz)

Frequency (Hz)

IF

0

5

40 20 0

0

0.2

0.4 0.6 Time (s) (e)

0.8

1

60 40 20 0

0

0.2

0.4 0.6 Time (s) (f)

0.8

1

Fig. 2. Illustration of the TFR with different window functions for the signal defined by f (t ) = sin(2π (50t )). The TFR of the STFT with (a) Blackman window, (b) Hanning window, (c) Gaussian window. The TFR of the SET with (d) Blackman window, (e) Hanning window, (f) Gaussian window.

Proof. This theorem provides a strong approximation result of the class SE ,d . Note that the first two conclusions are the same as that of the FSST and has been proved in [46]. Herein we only give the detailed version about the third one: (3)We first give the following equality:



|

η − φ  (t )| < ε˜T e f (η, t )dη 

= Sgf (η, t )η=φ  (t ) k  +∞ = ( f (ξ ) · g(ξ − t ) · e−2iπη(ξ −t ) )dξ |η=φk (t ) −∞  +∞  = f (y + t ) · g(y ) · e−2iπφk (t )y dy, −∞

where y = ξ − t . Then, we have:

   1     T e ( η , t ) d η − f ( t ) f k  gˆ(0 ) |η−φ  (t )|<ε˜     + ∞  1  1 −2iπφk (t )y  = f k ( y + t ) · g( y ) · e dy − gˆ(0 ) · fk (t ) gˆ(0 ) −∞ gˆ(0 )  1  +∞  = fk (y + t ) · g(y ) · e−2iπφk (t )y dy gˆ(0 )  −∞   +∞  2iπφk (t ) − g(y ) · Ak (t ) · e dy −∞    1  +∞ 2iπφk (y+t )−2iπφk (t )y ≤ g(y )[Ak (y + t ) − Ak (t )]e dy gˆ(0 )  −∞  1  +∞  + g(y )[Ak (t ) · e2iπφk (y+t )−2iπφk (t )y gˆ(0 )  −∞  −Ak (t ) · e2iπφk (t ) ]dy  +∞   1 |g(y )||Ak (y + t ) − Ak (t )||e2iπφk (y+t )−2iπφk (t )y |dy ≤ gˆ(0 ) −∞  +∞   1 |g(y )||Ak (t )||e2iπ (φk (y+t )−φk (t )y) − e2iπ φk (t ) |dy + gˆ(0 ) −∞

≤

 +∞ 1 |g(y )| · |Ak (y + t ) − Ak (t )|dy gˆ(0 ) −∞  +∞  2π |g(y )| · |Ak (t )||φk (t ) + φk (t )y + gˆ(0 ) −∞



+φk (ξt ) 

 y2 − φk (t )y − φk (t )|dy 2

+∞

A (t )π |g(y )||y|dy + k gˆ(0 ) −∞ gˆ(0 )

= (U1 + π U2 Ak (t )). gˆ(0 )

≤



+∞ −∞

|g(y )||y|2 dy

−1

Let (t ) = maxk {U1 + π U2 Ak (t ))}. If ˜ ≤ (t )∞2 , for all t ∈ R and k ∈ {1, . . . , N}, it holds:

   1     T e ( η , t ) d η − f ( t ) k  gˆ(0 ) |η−φ  (t )|<ε˜ f 

≤ (U1 + π U2 Ak (t )) gˆ(0 ) 1

(t ) gˆ(0 ) 1 ≤

˜ . gˆ(0 ) ≤

 Remark 2. According to the Eqs. (18) and (37), to decompose a mono-component signal, both of the FSST and SET need the IF trajectories. However, apart from the calculation of IF trajectories, the FSST has to know the integration regions additionally. In contrast, the reconstruction of the SET is more convenient and straightforward. 4. The properties of the SET The SET shows some according to Eqs. (5) and (20). The corresponding proofs are presented in Appendix.

6

Z. Li, J. Gao and H. Li et al. / Signal Processing 166 (2020) 107243

(1) Linearity: As is well-known, the STFT is a linear transform. Though the frequency estimation operator defined in (20) is nonlinear, the SET still holds the linearity because of the extracting operator. (2) Time and frequency shifts: Considering f0 (t ) = f (t − t0 )e2iπ η0 t , we can obtain the following expressions: 2iπη0 t

S f 0 ( η , t ) = S f ( η − η0 , t − t 0 ) e , ωˆ f0 (η, t ) = ωˆ f (η − η0 , t − t0 ) + η0 , T e f0 (η, t ) = T e f (η − η0 , t − t0 )e2iπη0 t . This property denotes that the frequency estimation operator will be shifted and add amount of movement accordingly when the signal is shifted in the time or frequency direction. Moreover, the amplitude of the SET will be shifted in time or frequency direction, and the phase also generates a corresponding frequency shift. (3) Complex conjugate: If f0 (t ) = f (t ), we have S f0 (η, t ) = S f (−η, t ) and:

  |T e f (η, t )| − |IT F R(η, t )| ≤ ,

(37)

where is a small enough constant. Combining the localization with the norm, we can give a conclusion that the result of the SET towards to that of the ITFR for weakly amplitude-modulated and frequency-modulated (AM-FM) signals.

5. Numerical analysis of the behaviour of the SET In this part, we take several examples to verify the performance of the SET. Meanwhile, we also make a comparison between with STFT and FSST methods. For all the examples, the STFT is config1

ured with a Gaussian window function g(t ) = σ − 2 e

ωˆ f0 (η, t ) = −ωˆ f (−η, t ), and T e f0 (η, t ) = T e f (−η, t ). (4) Perfect localization on chirp signal: For a pure linear chirp signal f (t ) = Ae2iπ φ (t ) , where φ (t) is a second-order polynomial and A > 0, we can obtain

ωˆ f (η, t ) = b + ct +

approaching to IF φ  (t). Consequently, the localization of the SET is perfectly toward to that of the ITFR. We also pay attention to the norm of the TFR, which can characterize the TF structure of the signal. We have

c2 1/(bπ )2 + c2

(η − (b + ct )).

(34)

Then, it holds

η = ωˆ f (η, t ) ⇔ η = b + ct

(35)

which is exactly the IF of the chirp. Thus

T e f (η, t ) = Sgf (η, t )δ (η − ω ˆ f (η, t )) = Sgf (b + ct, t ).

(36)

It means that the SET of a chirp signal only appears in the IF of the signal, that is η = b + ct. (5)Perfectly toward to ITFR: According to Eqs. (1) and (20), we know that the localizations of the SET and ITFR are determined by the frequency estimation ω ˆ f (η, t )) and IF φ  (t), respectively. Furthermore, δ (η − ω ˆ f (η, t )) describes the ridge in the TF domain,

2

− π t2 σ

.

5.1. Quality of representation To quantify the enhancement of the SET on the concentration of TFR, we introduce the normalized energy that is associated with the first coefficients with the largest amplitude, which means: the faster it towards 1, the more concentrated the TFR [44]. In Fig. 3(a) and (b), we represent the results of normalized energies for two test-signals defined in (38) and (39). As expected, the results of the SET are almost perfect, since the energy are only contained in 2 coefficients. Comparatively, the energy of the FSST grows slightly slow. In order to investigate the performance of the SET in the presence of noise, we carry out the same experiments but for the signals added with Gaussian white noise (noise level is 0 dB). The images, as shown in Fig. 3(c) and (d), exhibit a slower increase of the normalized energy because of the diffusion corresponding to noise. However, we can still check that the concentration of the SET is still better than that of the STFT and FSST.

Fig. 3. Normalized energy as a function of the number of the sorted coefficients for (a) signal f1 (t), (b) signal f2 (t), (c) signal f1 (t) with 0 dB noise, (d) signal f2 (t) with 0 dB noise.

Z. Li, J. Gao and H. Li et al. / Signal Processing 166 (2020) 107243

7

Fig. 4. EMD corresponding to the STFT, FSST, or SET for (a) signal f11 (t), (b) signal f22 (t).

Fig. 5. Illustration of FSST and SET for an exponential chirp. The corresponding TFRs given by the (a) FSST, (b) SET.

Fig. 6. Illustration of FSST and SET for a hyperbolic chirp. The corresponding TFRs provided by the (a) FSST, (b) SET.



f11 (t )

f12 (t )

    f1 (t ) = sin(2π (32t + 10sin(t ))) + sin(2π (44t + 10sin(t ))), (38) 



f21 (t )



f22 (t )

    f2 (t ) = sin(2π (250t + 50(t ) )) + sin(2π (130t + 100(t )2 )) . 3

(39) In order to further evaluate the quality of TFR, we employ the means of the Earth mover’s distance (EMD) [39], which measures the dissimilarity between the obtained TFR and the ideal one. The smaller EMD denotes the more concentrated TFR and less noise fluctuations. Herein, we just take components f11 (t) and f22 (t) as the examples. The corresponding results are shown in Fig. 4. Whether for the weakly or strongly FM mode, the SET result achieves the minimum in each noise level, which means that the SET has better ability to improve the TF concentration in comparison with the other two methods.

5.2. Limiting modulations Taking the extracting operator into account, the SET is suitable for more types of signals. For verifying the validity of this infer, we first consider three possible strong frequency modulations, i.e. a phase φ (t) s.t. ε ≤ φ  (t). The first one is an exponential chirp with phase φ  (t ) = 10e3t , the second one is a hyperbolic chirp with phase φ  (t ) = −50 ∗ log(1.02 − t ), and the last one is a multi-component signal containing strongly non-linear sinusoidal frequency modulations defined in (40). The corresponding TFRs generated by using the FSST and SET are displayed in Figs. 5– 7, respectively. As expected, FSST provides a sharp representation where frequencies change slowly. However, the energy is spread out where frequencies vary quickly. We can check that the results of the SET are always with better concentration in the whole TF plane. Herein, we particularly test on a “crossover signal” in (41), which contains two components with intersecting IF trajectories. Fig. 8 displays the TFRs provided by using the FSST and SET. As

8

Z. Li, J. Gao and H. Li et al. / Signal Processing 166 (2020) 107243

Fig. 7. Illustration of FSST and SET for the multi-component signal f3 (t). The corresponding TFR obtained by the (a) FSST, (b) SET.

Fig. 8. Illustration of FSST and SET for f4 (t). The corresponding TFRs generated by the (a) FSST, (b) SET.

shown in the illustration, the SET-based image is more concentrated and has less cross-terms. f31 (t )





f32 (t )

    f3 (t ) = sin(2π (330t + 16cos3π t )) + sin(2π (190t + 9cos3π t )) 

f33 (t )





+ sin(80π t ),

f4 (t ) = cos(300π t ) + cos(2π (80(t ) + 150(t )2 +

(40) 3

π

cos(4π t )). (41)

5.3. Choice of the window parameter σ The choice of the window parameter σ is an important step for each time-frequency representation. To study the impact of the window on the representation, we introduce the Renyi entropy [49,51], which can measures the concentration of the TFR. The smaller the Renyi entropy is, the more concentrated the representation. As shown in Fig. 9, we illustrate the evolution of Renyi entropy with respect to σ for the signal f1 (t). For a given σ , the SET always provides a minimum value, which denotes the most concentrated TFR. The study of the approximation theorem of FSST suggests that the window should be taken narrow enough, to minimize the estimation error; while it must be wide enough to satisfy the separation condition between the different modes. This phenomenon is illustrated in Fig. 10 using three different window sizes (σ 1 < σ 2 < σ 3 ). It can be observed that the TFR achieved by the FSST contains cross-terms when the σ 3 is much small, and the TFR will be diffused dramatically when the σ 1 is much large. While, the corresponding result of the SET is always more concentrated and with less cross-terms or diffusion that of the FSST. That is to say that the SET is less sensitive to window parameter σ .

Fig. 9. Evolutions of Renyi entropies with respect to σ for the noisy free f1 (t) calculated by the STFT, FSST, and SET, respectively.

5.4. Perfect localization towards to ITFR As mentioned in Section 4, the SET is perfectly toward to the ITFR for weakly FM signals. To demonstrate this property, we employ mode f11 (t) with constant amplitude and signal f5 (t ) = e−0.5t f11 (t ) with time-varying amplitude. Their corresponding TFRs and modules of the SET and ITFR are presented in Figs. 11 and 12, respectively. It clearly shows that both the localization and magnitude of the SET are towards that of the ITFR for both f11 (t) and f5 (t). Herein, we use mode f21 (t), f22 (t) and f32 (t) to further verify the ability of the TF localization of the SET. Fig. 13 displays the corresponding IF and estimation calculated by using the SEO η = ω ˆ f (η, t ). As expected, Fig. 13(a) clearly shows that the SET achieves a perfect TF localization for the chirp f21 (t). As shown in Fig. 13(b) and (c), the SET still provides satisfactory results for the highly FM modes f22 (t) and f32 (t).

Z. Li, J. Gao and H. Li et al. / Signal Processing 166 (2020) 107243

9

Fig. 10. The TFRs of f1 (t) ontained by the FSST (top) and SET (bottom) with different window parameters, from left to right: σ 1, σ 2, σ 3.

1.1

Amplitude

ITFR SET

1

0.9

0

1

2 Time (s) (c)

3

4

Fig. 11. Comparisons of the localization and magnitude for f11 (t). The TFR provided by the (a) ITFR, (b) SET. (c) The corresponding magnitude of the ITFR and SET, (d) the relative error of the magnitudes in (c).

5.5. Illustration on real signals To better show the effectiveness of the SET, we employ two real-life signals. The first one is a portion of a speech signal, which is the same in [44]. Fig. 14 shows the corresponding TF distributions obtained by using the three different methods, as well as a zoom on a modulated portion. As shown in Fig. 14(d), the energy of STFT-based result is diffused seriously. For the FSST, the TF energy is well concentrated only where frequencies change slowly

(see Fig. 14(e)). As expected, Fig. 14(f) from the SET shows a result with high concentration along the ridges in the whole TF domain. Then, we apply the SET algorithm to characterize the whistles of killer whales [52], whose sampling time and sampling frequency are 3 seconds and 11025 Hz, respectively. The sound is very important for killer whales to echolocate and communicate. Thus, it is significant to study the TF feature of the sound signal, which is shown in Fig. 15(a). The corresponding TFRs calculated by

10

Z. Li, J. Gao and H. Li et al. / Signal Processing 166 (2020) 107243

1.2 ITFR SET

Amplitude

1 0.8 0.6 0.4 0.2 0

0

1

2 Time (s) (c)

3

4

Fig. 12. Comparisons of the localization and magnitude for f5 (t). The TFR provided by the (a) ITFR, (b) SET. (c) The corresponding magnitude of the ITFR and SET, (d) the relative error of the magnitudes in (c).

500

500

500

IF

IF

200 100

400 Frequency (Hz)

300

0

IF

400 Frequency (Hz)

Frequency (Hz)

400

300 200 100

0

0.2

0.4 0.6 Time (s) (a)

0.8

1

0

300 200 100

0

0.2

0.4 0.6 Time (s) (b)

0.8

0

1

0

0.2

0.4 0.6 Time (s) (c)

0.8

1

Fig. 13. The localization of the SET for highly FM signal (a) f21 (t), (b) f22 (t), (c) f32 (t).

using the STFT, FSST, and SET are respectively shown in Fig. 15(b)– (d). From the comparisons, we can clearly observe that this sound signal contains three components and is blurred by the background noise. Moreover, the SET provides a TFR with better energy concentration than the other two methods. The more concentrated TFR denotes the better ability of the TF localization and the better characterization of time-varying features. This indicates the potential of the SET in the application for whale research.

Then, the STFT will be easily obtained



Sgf

(η , t ) = e

2iπ φ (t )

σ1

1  (i/2π )n (n) ∂ n − σ2 π 2 (η−φ (t ))2 A (t ) n e 1−iσ2 πφ (t ) . n! ∂η

For the signal f (t ) = the assumed condition |A (t)| < ε |φ  (t)| is indispensable in the FSST and SET work. Now, we consider a more strongly amplitude-modulated (AM) signal, and the first-order expansion of the amplitude and second-order expansion of the phase are

A(τ + t ) = A(t ) + A (t )τ ,

(42)

φ (τ + t ) = φ (t ) + φ  (t )τ + φ  (t )τ 2 .

(43)

(44)

n=0

We will have the IF estimation

6. Discussion A(t )e2iπ φ (t ) ,

σ2 π

1 − iσ2 π φ  (t )

σ π 2 (η−φ  (t ))2 − 21−iσ πφ  (t ) 2

⎫ ⎪ ⎪ ⎬

A(t )e 1 (i/2π )n  ∂ n − σ2 π 2 (η−φ (t ))2 ⎪ ⎪ ⎪ A(n ) (t ) n e 1−iσ2 πφ (t ) ⎪ ⎩ ⎭ n! ∂η n=0 ⎧ g ⎫ ⎨ ∂ S f (η,t ) ⎬

ωˆ f (η, t ) = φ  (t ) + 

−φ  (t )

⎧ ⎪ ⎪ ⎨

∂η

,

⎩ Sgf (η, t ) ⎭

(45)

which indicates that the φ  (t) no longer satisfies η = ω ˆ f ( η , t ). Moreover, it also generates a relatively large error between ωˆ f (η, t ) and φ  (t). Thus, the SET will provide an imprecise TFR for

Z. Li, J. Gao and H. Li et al. / Signal Processing 166 (2020) 107243

11

Fig. 14. TFRs of the speech signal. The results are calculated by the (a) STFT, (b) FSST, (c) SET. The corresponding zooms on a strongly modulated part are shown in (d), (e), and (f), respectively.

Fig. 15. Illustrations of the whistles of killer whales. (a) The sound wave. The TFRs achieved by the (b) STFT, (c) FSST, and (d) SET.

the highly AM signal. Fig. 16(a) shows an example, whose definition is in (46) containing high amplitude modulation in time 3 − 4s. The corresponding TFR from the SET is displayed in Fig. 16(b) and is indeed not exact.

f6 (t ) = e0.6(t−2) sin(2π (25t + 10sin(1.5t ))).

lows:

Sgf

 = f (t )σ1

(46)

As proved in Theorem 2, the SET can well reconstruct the weakly AM-FM signal, i.e., |A (t)|, |φ  (t)| < ε |φ  (t)|. However, under the assumption (42) and (43), the STFT in the IF φ  (t) is given as fol-



(φ  (t ), t ) = e2iπ φ (t ) σ1

σ2 π

1 − iσ2 π φ  (t )

σ2 π

1 − iσ2 π φ  (t )

A(t )

.

(47)

√ Due to σ1 σ2 π = gˆ(0 ), we can rewrite (46) as:



Sgf

(φ  (t ), t ) = f (t )gˆ(0 )

1 , 1 − iσ2 π φ  (t )

(48)

12

Z. Li, J. Gao and H. Li et al. / Signal Processing 166 (2020) 107243

40

Frequency (Hz)

IF 35 30 25 20

3

3.2

3.4 3.6 Time (s) (b)

3.8

4

Fig. 16. The negative example 1. (a) The signal f6 (t) defined in (46), (b) the time-frequency representation from the SET in time 3 − 4s.

35

Frequency (Hz)

IF 30 25 20 15

0

1

2 Time (s) (b)

3

4

Fig. 17. The negative example 2. (a) The signal f7 (t)defined in (50), (b) the time-frequency representation from the SET.

which derives the following reconstruction formula

f (t ) =

Sgf

(φ  (t ), t )  gˆ(0 )

1 − iσ2 π φ  (t ).

suitable for more categories of signals and less sensitivity to the choice of the window length.

(49) Declaration of Competing Interest

According to Eq. (49), the SET is no longer suitable for the reconstruction of the signal containing highly AM modes. Furthermore, as shown in Fig. 17, the SET also cannot deal with the signal containing strong frequency modulations, i.e., the φ(t) is not negligible. The signal shown in Fig. 17(a) is defined in (50). It is clearly seen in Fig. 17(b) that the TFR of the SET is deviated from the IF. For the more exact spectrum and signal reconstruction, we should design better TF analysis method that considers higher order informations of the amplitude and phase of the signal. A novel TF analysis method termed high-order synchrosqueezing transform was proposed recently. It can obtain a more concentrated TFR for strongly AM-FM signals. In the future work, we can combine the idea of the SET with this high-order SST to better improve the quality of TFR for the general signal.

f7 (t ) = sin(2π (25t + 0.5sin(8t ))).

(50)

7. Conclusion We revisited the synchroextracting transformation (SET) in this work. The innovation of the SET is the synchroextracting operator (SEO). We provided the approximation theorem of the SET for weakly modulated signals. Several properties of the SET were studied and proved theoretically. We also gave a new view of the SEO that it is equivalent to fixed frequencies of the synchrosqueezing operator. Therefore, the time-frequency concentration and sparsity are greatly improved by the SET. Compared with the FSST, the SET algorithm is a more powerful tool for the time-frequency analysis of non-stationary signals, even for higher frequency modulated signals. In short, we proved that the SET can obtain a more concentrated TFR and is more robust against noise. As a result, it is

We declare that we do not have any commercial or associative interest that represents a conflict of interest in connection with the work submitted. Acknowledgment This work is supported in part by the Major Research Plan of the National Natural Science Foundation of China (91730306), in part by the Major National Science and Technology Projects (2016ZX05024-001-007 and 2017ZX05069), in part by the National Key R & D Program of the Ministry of Science and Technology of China (2018YFC0603501), in part by the National Postdoctoral Program for Innovative Talents under Grant BX201900279, and in part by the Fundamental Research Funds for the Central Universities under Grant xjh012019030. The authors appreciate the editor and reviewers for their valuable comments, which improve the presentation of this work. Appendix A. Properties of the SET Some properties of the SET presented in Section 4 are proved as the following. A1. Linearity Assume that f (t ) = α f1 (t ) + β f2 (t ), and f1 (t), f2 (t) ∈ Bε,d . From the linearity of STFT, it is known that

Sgf (η, t ) = α Sgf (η, t ) + β Sgf (η, t ). 1

2

(A.1)

Z. Li, J. Gao and H. Li et al. / Signal Processing 166 (2020) 107243

Then, substituting (A.1) into (20), we obtain the reassignment operator



 ∂t (α Sgf1 (η, t ) + β Sgf2 (η, t )) ωˆ f (η, t ) =  . 2iπ (α Sgf (η, t ) + β Sgf (η, t )) 1 2

ωˆ f (η, t ) =

   ωˆ f1 (η, t ), i f η ∈ φ1 (t ) − d, φ1 (t ) + d; ωˆ f2 (η, t ), i f η ∈ φ2 (t ) − d, φ2 (t ) + d .

(A.3)

Finally, we will obtain

=

α

Sgf 1 Sgf 1

(η , t ) + β

Sgf 2

α

=

α T e f 1 ( η , t ) + β T e f 2 ( η , t ).

(A.4)

Let us consider f0 (t ) = f (t − t0 )e2iπ η0 t . Then, we have 0

= = =

−∞  +∞ −∞  +∞ −∞ Sgf

f (τ ) · g(τ − (t − t0 )) · e−2iπ (η−η0 )(τ −(t−t0 ) dτ e2iπ η0 t 2iπη0 t

.

(A.5)

(A.6)

( η − η0 , t − t 0 ) δ ( ( η − η0 ) # −ω ˆ f (η − η0 , t − t0 ) e2iπη0 t 2iπη0 t

.

(A.7)

(A.8)

According to Eq. (5), the STFT of f(t) is then written as

Sgf (η, t )  +∞ = f (ξ ) · g(ξ − t ) · e−2iπω (ξ −t ) dξ f (τ + t ) · g(τ ) · e−2iπητ dτ 2

1   2 −τ f (t )e2iπ [φ (t )τ + 2 φ (t )τ ] · σ1 e σ2 · e−2iπητ dτ

=

σ1 · f (t )

=

σ1 · f (t )



+∞

e

 −( σ1 −iπφ  (t ))(τ + i(2π1 (η−φ (t ))) )2 2

2( σ −iπφ (t )) 2

−∞



chirp

signals:

(a)

Due

to

(A.13)

Because |ω ˆ (η, t ) − φ  (t )| ≤ ˜ , and the function gˆ is continuous

and gˆ(0 ) = 1, it holds

|T e f ( η , t )| ≈ | f (t )gˆ(η − φ  (t ))δ (η − ω ˆ f (η, t ))|

= | f (t )gˆ(ω ˆ f (η, t ) − φ  (t ))|

σ2 π 2 (η−φ  (t ))2 σ2 π e 1−iσ2 πφ (t ) ,  1 − iσ2 π φ (t )

(A.14)

dτ · e

2  − (2π1(η−φ (t)))

4( σ −iπφ (t )) 2

(A.9)

(A.15)

Thus

||T e f (η, t )| − |IT F R(η, t )|| ≤ .

. For such signal, wherever t and τ , it holds

−∞

linear

|IT F R(η, t )| = | f (t )δ (η − φ  (t ))| = | f (t )|.

1   2 f (τ + t ) = f (t )e2iπ [φ (t )τ + 2 φ (t )τ ] .

=

the

Furthermore

Let us consider a linear chirp signal f (t ) = Ae2iπ φ (t ) with φ (t ) = a + bt + 12 ct 2 and a Gaussian window function g(t ) =

−∞  +∞

3. For

(φ  (t ))2 < 1, FSST can make the TFR more concentrated 1/(σ2 π )2 +(φ  (t ))2  2 than STFT. (b) Since 1/(σ (πφ)2(+t ())φ  (t ))2 is increasing with increasing 2 φ  (t), the smaller the φ  (t) is, the more concentrated the FSST

≈ | f (t )gˆ(0 )|

A3. Perfect localization on chirp signal

=

(A.11)

(A.12)

= | f (t )|.

Sgf

= T e f ( η − η0 , t − t 0 ) e

−∞ +∞

(φ  (t ))2 (η − φ  (t )), 1/(σ2 π )2 + (φ  (t ))2

= |Sgf (η, t )δ (η − ω ˆ f (η, t ))|

0



ωˆ f (η, t ) = φ  (t ) +

Sgf (η, t ) ≈ f (t )gˆ(η − φ  (t )).

T e f0 (η, t ) = Sgf (η, t )δ (η − ω ˆ f0 (η, t ))

2

Then we achieve the estimation of the IF

f (ξ − t0 )e2iπη0 ξ · g(ξ − t ) · e−2iπη (ξ −t ) dξ

The result of SET will be derived as follows:

−σ

(A.10)

According to [33,36], we have the following approximate expression

ωˆ f0 (η, t ) = ωˆ f (η − η0 , t − t0 ) + η0 .

σ1 e

(η − φ  (t ))Sgf (η, t ).

A4. Perfectly toward to ITFR

Then, we obtain the IF estimation

t2

1/(σ2 π )2 + (φ  (t ))

f0 (ξ ) · g(ξ − t ) · e−2iπη (ξ −t ) dξ

( η − η0 , t − t 0 ) e

=

2( σ12 + iπ φ  (t ))(φ  (t ))2

representation.

A2. Time and frequency shifts

Sgf (η, t ) =

+

Remark

(η, t )δ (η − ωˆ f1 (η, t )) + β Sgf2 (η, t )δ (η − ωˆ f2 (η, t ))

+∞

∂t Sgf (η, t ) = φ  (t )Sgf (η, t )

ωˆ f (η, t ) = η ⇔ η = φ  (t ) = b + ct.

# (η, t ) δ (η − ωˆ f (η, t ))

=



g

which derives

T e f (η, t ) = Sgf (η, t )δ (η − ω ˆ f (η, t ))

"

with τ = ξ − t. Then, we calculate the derivative of S f (η, t )

(A.2)

Since these two signals satisfy that |φ1 (t ) − φ2 (t )| > 2d, we have

13

(A.16)

Remark 4. In the numerical simulations, we generically assume that



δ (t ) =

1, t = 0; 0, t = 0.

(A.17)

References [1] D. Gabor, Theory of communication, J. Inst. Electr. Eng. 93 (1946) 429–457. [2] E.P. Wigner, On the quantum correction for the rmodynamic equilibrium, Phy. Rev. 40 (1946) 749–759. [3] P. Flandrin, Time-Frequency/Time-Scale Analysis, 10, Academic Press, 1978, pp. 64–76. [4] D. Gabor, Theory of communication, J. Inst. Electr. Eng. 93 (26) (1946) 429–457. [5] F. Gianfelici, G. Biagetti, P. Crippa, C. Turchetti, AM-FM decomposition of speech signals: an asymptotically exact approach based on the iterated hilbert transform, in: Proceedings of the IEEE Workshop on Statistical Signal Processing Proceedings, 2005, pp. 333–337. [6] F. Gianfelici, G. Biagetti, P. Crippa, C. Turchetti, Multicomponent AM-FM representations: an asymptotically exact approach, IEEE Trans. Audio Speech Lang. Process. 15 (3) (2007) 823–837. [7] S. Pei, K. Chang, The mystery curve: a signal processing point of view, IEEE Signal Process. Mag. 34 (6) (2017) 158–163. [8] X. Li, G. Bi, S. Stankovi, A.M. Zoubir, Local polynomial Fourier transform: a review on recent developments and applications, Signal Process. 91 (6) (2011) 1370–1393.

14

Z. Li, J. Gao and H. Li et al. / Signal Processing 166 (2020) 107243

[9] Y. Yang, Z. Peng, X. Dong, W. Zhang, G. Meng, General parameterized time-frequency transform, IEEE Trans. Signal Process. 62 (11) (2014) 2751–2764. [10] Z. Yu, Y. Sun, W. Jin, A novel generalized demodulation approach for multi– component signals, Signal Process. 118 (2016) 188–202. [11] S. Wang, X. Chen, G. Cai, B. Chen, X. Li, Z. He, Matching demodulation transform and synchrosqueezing in time-frequency analysis, IEEE Trans. Signal Process. 62 (1) (2014) 69–84. [12] K. Kodera, R. Gendrin, C. Villedary, Analysis of time-varying signals with small BT values, IEEE Trans. ASSP 26 (2016) 1009. [13] F. Auger, P. Flandrin, Improving the readability of time-frequency and time-scale representations by the reassignment method, IEEE Trans. Signal Process. 43 (5) (1995) 1068–1089. [14] N.E. Huang, Z. Shen, S.R. Long, M.C. Wu, H.H. Shih, Q. Zheng, N.C. Yen, C.C. Tung, H.H. Liu, The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis, Proc. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci. 454 (1971) (1998) 903–995. [15] I. Daubechies, S. Maes, A nonlinear squeezing of the continuous wavelet transform based on auditory nerve models, in: A. Aldroubi, M. Unser (Eds.), Wavelets in Medicine and Biology, CRC Press, 1996, pp. 527–546. [16] I. Daubechies, J.F. Lu, H.T. Wu, Synchrosqueezed wavelet transforms: an empirical mode decomposition-like tool, Appl. Comput. Harmon. Anal. 30 (2) (2011) 243–261. [17] O. Rioul, P. Flandrin, Time-scale energy distributions: a general class extending wavelet transforms, IEEE Trans. Signal Process. 40 (7) (1992) 1746–1757. [18] I. Djurovic, L.J. Stankovic, Time frequency representation based on the reassigned S-method, Signal Process. 77 (1) (1999) 115–120. [19] X. Li, G. Bi, The reassigned local polynomial periodogram and its properties, Signal Process. 89 (2) (2009) 206–217. [20] K. Fitz, L.Haken, On the use of time frequency reassignment in additive sound modeling, J. Audio Eng. Soc. 50 (11) (2002) 879–893. [21] S.S. Ram, H. Ling, Analysis of microDopplers from human gait using reassigned jointtime frequency transform, Electron. Lett. 43 (23) (2007) 282–287. [22] X. Wu, T. Liu, Spectral decomposition of seismic data with reassigned smoothed pseudo Wigner Ville distribution, J. Appl. Geophys. 68 (3) (2009) 386–393. [23] F. Auger, P. Flandrin, Y.T. Lin, S. McLaughlin, S. Meignen, T. Oberlin, H.T. Wu, Time-frequency reassignment and synchrosqueezing: an overview, IEEE Signal Process. Mag. 30 (6) (2013) 32–41. [24] S.T. Quek, P.S. Tua, Q. Wang, Detecting anomalies in beams and plate based on the Hilbert-Huang transform of real signals, Smart Mater. Struct. 12 (3) (2003) 447–460. [25] B. Liu, S.D. Riemenschneider, Y.S. Xu, Gearbox fault diagnosis using empirical mode decomposition and Hilbert spectrum, Mech. Syst. Signal Process. 20 (3) (2006) 718–734. [26] B.M. Battista, C.C. Knapp, T. McGee, V. Gobel, Application of the empirical mode decomposition and Hilbert-Huang transform to seismic reflection data, Geophysics 72 (2) (2007) H29–H37. [27] D.P. Mandic, R.U. Rehman, Z. Wu, N.E. Huang, Empirical mode decomposition-based time-frequency analysis of multivariate signals: the power of adaptive data analysis, IEEE Signal Process. Mag. 30 (6) (2013) 74–86. [28] J.B. Tary, R.H. Herrera, J. Han, V.D. Mirko, Spectral estimation What is new? What is next? Rev. Geophys. 52 (4) (2014) 723–749. [29] Z. Wu, N. Huang, Ensemble empirical mode decomposition: a noise-assisted data analysis method, Adv. Adapt. Data Anal. 1 (1) (2009) 1–41. [30] M. Torres, M. Colominas, G. Schlotthauer, P. Flandrin, A complete ensemble empirical mode decomposition with adaptive noise, in: Proceedings of the IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 2011, pp. 4144–4147.

[31] G. Rilling, P. Flandrin, One or two frequencies? The empirical mode decomposition answers, IEEE Trans. Signal Process. 56 (1) (2008) 85–95. [32] M. Feldman, Analytical basics of the EMD: two harmonics decomposition, Mech. Syst. Signal Process. 23 (7) (2009) 2059–2071. [33] G. Thakur, H.T. Wu, Synchrosqueezing-based recovery of instantaneous frequency from nonuniform samples, SIAM J. Math. Anal. 43 (5) (2012) 2078–2095. [34] T. Oberlin, S. Meignen, V. Perrier, The Fourier-based synchrosqueezing transform, Proc. 39th Int. Conf. Acoust., Speech, Signal Process. (ICASSP), 2014, pp. 315–319. [35] N. Liu, J. Gao, Q. Wang, The extraction of instantaneous frequency from seismic data via synchrosqueezing three parameter wavelet transform, in: Proceedings of the SEG Technical Program Expanded Abstracts, 2015, pp. 2801– 2805. [36] N. Liu, J. Gao, Z. Zhang, X. Jiang, Q. Lv, High resolution characterization of geological structures using synchrosqueezing transform, Interpretation 5 (1) (2017) T75–T85. [37] Z. Huang, J. Zhang, Z. Zou, Synchrosqueezing S-transform and its application in seismic spectral decomposition, IEEE Trans. Geosci. Remote Sens. 54 (2) (2016) 817–825. [38] Q. Wang, J. Gao, N. Liu, X. Jiang, High-resolution seismic time-frequency analysis using the synchrosqueezing generalized S-transform, IEEE Geosci. Remote Sens. Lett. PP (99) (2018) 1–5. [39] I. Daubechies, Y. Wang, H.T. Wu, ConceFT: concentration of frequency and time via a multitapered synchrosqueezed transform, Philos. Trans. R. Soc. 374 (2065) (2016) 20150193. [40] L.J. Stankovi, S. Stankovi, An analysis of instantaneous frequency representation using time-frequency distributions-generalized Wigner distribution, IEEE Trans. Signal Process. 43 (2) (1995) 549–552. [41] L.J. Stankovi, I. Djurovi, S. Stankovi, M. Simeunovi, S. Djukanovi, M. Dakovi, Instantaneous frequency in time-frequency analysis: enhanced concepts and performance of estimation algorithms, Digit. Signal Process. 35 (2014) 1–13. [42] G. Thakur, H.T. Wu, Synchrosqueezing-based recovery of instantaneous frequency from nonuniform samples, SIAM J. Math. Anal. 43 (5) (2012) 2078–2095. [43] C. Li, M. Liang, A generalized synchrosqueezing transform for enhancing signal time-frequency representation, Signal Process. 92 (9) (2012) 2264–2274. [44] T. Oberlin, S. Meignen, V. Perrier, Second-order synchrosqueezing transform or invertible reassignment? Towards ideal time-frequency representations, IEEE Trans. Signal Process. 63 (5) (2015) 1335–1344. [45] T. Oberlin, S. Meignen, The second-order wavelet synchrosqueezing transform, in: Proceedings of the IEEE ICASSP, 2017. [46] R. Behera, S. Meignen, T. Oberlin, Theoretical analysis of the second-order synchrosqueezing transform, Appl. Comput. Harmon. Anal. 34 (2016) 1009. [47] S. Wang, X. Chen, I. Selesnick, Y. Guo, C. Tong, X. Zhang, Matching synchrosqueezing transform: a useful tool for characterizing signals with fast varying instantaneous frequency and application to machine fault diagnosis, Mech. Syst. Signal Process. 100 (2018) 242–288. [48] Y. Miao, H. Sun, J. Qi, Synchro-compensating chirplet Transform, IEEE Signal Process. Lett. 25 (2018) 1413–1417. [49] G. Yu, M. Yu, C. Xu, Synchroextracting transform, IEEE Trans. Ind. Electron. 64 (10) (2017) 8042–8054. [50] R. Herrera, J. Han, M. van der Baan, Applications of the synchrosqueezing transform in seismic time-frequency analysis, Geophysics 79 (3) (2014) V55–V64. [51] L. Stankovi, A measure of some time frequency distributions concentration, Appl. Comput. Harmon. Anal. 81 (3) (2001) 621–631. [52] Discovery of Sound in the Sea [Online]Available:http://www.dosits.org/.