Adaptive short-time Fourier transform and synchrosqueezing transform for non-stationary signal separation

Signal Processing 166 (2020) 107231

Contents lists available at ScienceDirect

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

Adaptive short-time Fourier transform and synchrosqueezing transform for non-stationary signal separationR Lin Li a, Haiyan Cai b, Hongxia Han a, Qingtang Jiang b,∗, Hongbing Ji a a b

School of Electronic Engineering, Xidian University, Xi’an 710071, PR China Department of Math & CS, University of Missouri-St. Louis, St. Louis, MO 63121, USA

a r t i c l e

i n f o

Article history: Received 6 March 2019 Revised 5 July 2019 Accepted 20 July 2019 Available online 23 July 2019 Keywords: Instantaneous frequency Adaptive short-time Fourier transform Adaptive synchrosqueezing transform Well-separated condition for multicomponent non-stationary signal Component recovery of non-stationary signal

a b s t r a c t The synchrosqueezing transform, a kind of reassignment method, aims to sharpen the time-frequency representation and to separate the components of a multicomponent non-stationary signal. In this paper, we consider the short-time Fourier transform (STFT) with a time-varying parameter, called the adaptive STFT. Based on the local approximation of linear frequency modulation mode, we analyze the wellseparated condition of non-stationary multicomponent signals using the adaptive STFT with the Gaussian window function. We propose the STFT-based synchrosqueezing transform (FSST) with a time-varying parameter, named the adaptive FSST, to enhance the time-frequency concentration and resolution of a multicomponent signal, and to separate its components more accurately. In addition, we also propose the 2nd-order adaptive FSST to further improve the adaptive FSST for the non-stationary signals with fast-varying frequencies. Furthermore, we present a localized optimization algorithm based on our wellseparated condition to estimate the time-varying parameter adaptively and automatically. Simulation results on synthetic signals and the bat echolocation signal are provided to demonstrate the effectiveness and robustness of the proposed method. © 2019 Elsevier B.V. All rights reserved.

1. Introduction To model a non-stationary signal as a superposition of locally band-limited, amplitude and frequency-modulated Fourier-like oscillatory modes:

x(t ) =

K 

xk (t ),

xk (t ) = Ak (t )ei2πφk (t ) ,

(1)

k=1

where Ak (t ), φk (t ) > 0, has been a very active research area over the past few years. Note that the number of component K may change with time t, but it should be constant for long enough time intervals. The representation of x(t) in (1) with Ak (t) and φk (t ) varying slowly or more slowly than φ k (t) is called an adaptive harmonic model (AHM) representation of x(t), where Ak (t) are called the instantaneous amplitudes and φk (t ) the instantaneous frequencies (IFs). To decompose x(t) as an AHM representation (1) is important to extract information, such as the underlying dynamics, hidden in x(t).

R This work was supported in part by the National Natural Science Foundation of China (Grant No. 61803294) and Simons Foundation (Grant No. 353185). ∗ Corresponding author.: E-mail address: [email protected] (Q. Jiang).

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

Time-frequency (TF) analysis is widely used in engineering fields such as communication, radar and sonar as a powerful tool for analyzing time-varying non-stationary signals [1]. Time-frequency analysis is especially useful for signals containing many oscillatory components with slowly time-varying amplitudes and instantaneous frequencies. The short-time Fourier transform (STFT), the continuous wavelet transform (CWT) and the WignerVille distribution are the most typical TF analysis, see details in [1– 6]. Other TF distributions of Cohen’s class include the exponential distribution [7], a smoothed pseudo Wigner distribution [8] and the complex-lag distribution [9]. In addition, the TF signal analysis and synthesis using the eigenvalue decomposition method has been studied [10,11]. In particular, an eigenvalue decompositionbased approach which enables the separation of non-stationary components with overlapped supports in the TF plane has been proposed in [12]. Recently a number of new TF analysis methods such as the Hilbert spectrum analysis with empirical mode decomposition (EMD) [13], the reassignment method [14] and synchrosqueezed wavelet transform (SST) [15] have also been proposed to obtain Ak (t) and φk (t ). EMD is a data-driven decomposition algorithm which separates the time series signal into a set of monocomponents, called intrinsic mode functions (IMFs) [13]. EMD has been studied by many

2

L. Li, H. Cai and H. Han et al. / Signal Processing 166 (2020) 107231

researchers and has been used in many applications, see e.g. [16– 23]. Because of the presence of widely disparate scales in a single IMF, or a similar scale residing in different IMF components, named as mode mixing [24], two close IMFs are hardly distinguished by EMD. The CWT-based synchrosqueezing transform (WSST), introduced in [15] and future studied in [25], is a special case of reassignment methods, which aims to sharpen the TF representation of the signal by allocating the CWT coefficient value to a different point in the TF plane. A variant of WSST, the STFT-based SST (FSST) was proposed in [26] and further studied in [27,28]. Both WSST and FSST have been proved to be robust to noise and small perturbations [29–31]. However for frequency-varying signals, the squeezing effect of SST is not desirable. In this regard, a 2nd-order SST was introduced in [32,33] and further studied in [34,35]. The 2ndorder SST improves the concentration of the TF representation well on perturbed linear chirps with Gaussian modulated amplitudes. The higher-order FSST was presented in [36], which aims to handle signals containing more general types. Other SST related methods include the generalized SST [37], a hybrid EMD-SST computational scheme [38], the synchrosqueezed wave packet transform [39], the S-transform-based SST [40], SST with vanishing moment wavelets [41], the multitapered SST [42] and the demodulation-transform based SST [43,44]. In addition, the synchrosqueezed curvelet transform for two-dimensional mode decomposition was introduced in [45], the signal separation operator which is related to FSST was proposed in [46] and the empirical signal separation algorithm was introduced in [47]. The statistical analysis of synchrosqueezed transforms has been studied in [48] and a new IF estimator within the framework of the signal’s phase derivative and the linear canonical transform was introduced in [49]. SST has been used in machine fault diagnosis [50,51], crystal image analysis [52,53], welding crack acoustic emission signal analysis [54], and medical data analysis [55–57]. Most of the FSST algorithms available in the literature are based on the short-time Fourier transform (STFT) with a fixed window, which means high time resolution and frequency resolution cannot be obtained simultaneously. For broadband signals, a wide window is suitable for the low-frequency parts. On the contrary, a narrow window is suitable for the high-frequency parts. To enhance the TF resolution and energy concentration, we propose in this paper the adaptive FSST based on the STFT with a time-varying window. More precisely, let Vx (t, η) be the (modified) STFT of x(t ) ∈ L2 (R ) with a window function h(t ) ∈ L2 (R ) defined by

Vx (t, η ) :=



∞ −∞

 =

∞ −∞

x(τ )h(τ − t )e−i2πη (τ −t ) dτ

(2)



∞ −∞

 =

∞ −∞

x(t + τ )h(τ )e−i2πητ dτ ,

(3)

x(τ )gσ (t ) (τ − t )e−i2πη (τ −t ) dτ x(t + τ )

 τ  e−i2πητ dτ , σ (t )

1 g σ (t )

(4)

(5)

where σ = σ (t ) is a positive function of t, and gσ (t) (τ ) is defined by

gσ (t ) (τ ) :=

xk (t + τ ) = Ak (t + τ )ei2π φk (t+τ ) ≈ Ak (t )e 

= xk (t )ei2π φk (t )τ

 τ  , σ (t )

1 g σ (t )

(6)

with g ∈ L2 (R ). The window width of gσ (t) (τ ) is σ (t) (up to a constant), depending on the time variable t. In this paper we consider

for



i2π φk (t )+φk (t )τ



τ ≈ 0,

(7)

then the STFT Vxk (t, η ) of xk (t) with a window function h(t) is

Vxk (t, η ) ≈ xk (t ) h(η − φk (t )), where h denotes the Fourier transform of h(t). Hence if supp( h) ⊆ [−, ] for some  > 0, then Vxk (t, η ) lies in the TF zone given by

Z k = {(t, η ) : |η − φk (t )| < , t ∈ R}. Therefore, if

φk (t ) − φk −1 (t ) ≥ 2, t ∈ R, 2 ≤ k ≤ K,

(8)

then Zk ∩ Z = ∅, k = , which means the components of x(t) are well separated in the TF plane. (8) is a required condition for the study of FSST in [27,28] and even for the study of the 2nd-order FSST in [32,35]. Here we call (8) the sinusoidal signal model-based well-separated condition for x(t). In this paper we use the linear frequency modulation (LFM) signal to approximate a non-stationary signal at any local time to x (t, η ). More precisely, study the TF zone of the adaptive STFT V k we assume that each xk (t) is well approximated by an LFM at any local time: for any t ∈ R, 1   2 xk (t + τ ) = Ak (t + τ )ei2π φk (t+τ ) ≈ Ak (t )ei2π (φk (t )+φk (t )τ + 2 φk (t )τ ) 1   2 = xk (t )ei2π (φk (t )τ + 2 φk (t )τ )

for τ ≈ 0,

(9)

where for a given t, the quantity in (9) as a function τ is called an LFM signal (or a linear chirp signal). Thus we have

x (t, η ) ≈ V k

where t and η are the time variable and the frequency variable respectively. In this paper we consider the STFT with a time-varying parameter σ (t) (called the adaptive STFT) defined by

x (t, η ) := V

x (t, η ) (called the adaptive FSST) and study the the FSST based on V choice of the time parameter σ (t) so that the adaptive FSST gives a better instantaneous frequency estimation of the component of a multicomponent signal, and provides more accurate component recovery. To recover/separate the components xk (t) of a multicomponent signal as given by (1) with the SST approach, Iatsenko et al. [30] indicates that if STFTs Vxk−1 (t, η ) and Vxk (t, η ) of two components xk−1 (t ) and xk (t) are mixed, then FSST cannot separate these two components xk−1 (t ) and xk (t) either, and hence it cannot recover/separate components accurately. Thus it is desirable that appropriate window width of the window function h(t) can be chosen (if possible) so that the STFTs of different components do not overlap. When xk (t) are sinusoidal signals Ak ei2π ck t for some constants Ak , ck > 0 or they are well approximated by sinusoidal functions at any local time in the sense that for any t ∈ R,



 τ  e−i2π ητ dτ . σ (t )

1 1   2 xk (t )ei2π (φk (t )τ + 2 φk (t )τ ) g σ (t ) R

(10) In this paper we will obtain the LFM model-based well-separated condition which guarantees that for different k, the quantities as functions of (t, η) on the right-hand side of (10) lie within nonoverlapping zones in the TF plane when g is the Gaussian window function. We will also discuss how to select the time-varying parameter σ (t) such that the corresponding adaptive FSST and 2ndorder adaptive FSST have sharp TF representation. In particular, we propose a localized optimization method based on our wellseparated condition to estimate the time-varying adaptive window width σ (t). The idea of using an optimal time-varying window parameter (window width) has been studied or considered extensively in the literature, see e.g. [58–63]. In particular, the authors in [63] introduced a method to select the time-varying window width for sharp SST representation by minimizing the Rényi entropy. In addition, after we completed our work, we were aware of the very

L. Li, H. Cai and H. Han et al. / Signal Processing 166 (2020) 107231

recent work [64] on the adaptive STFT-based SST with the window function containing time and frequency parameters. Our motivation is different from others in that we do not focus on the optimal parameter such that the corresponding STFT has the sharpest representation in the TF plane. Instead, we pursue the establishment of the LFM model-based well-separated condition for multicomponent signals based on the adaptive STFT and we propose how to select window width σ (t) such that the STFTs of the components lie in non-overlapping regions of the TF plane. The Rényi entropybased optimal parameter may give the overall sharp representation of STFT or FSST, but it does not guarantee all the components to be separated. The selected σ (t) proposed by us does not necessarily result in a sharp representation of the associated STFT. Instead, it is selected in such a way that the adaptive STFTs of the components are well separated, and the corresponding adaptive FSSTs have sharp representation and hence the components can be recovered more accurately. The remainder of this paper is organized as follows. We introduce the adaptive STFT and adaptive FSST with a time-varying parameter in Section 2, where we also introduce the 2nd-order adaptive FSST. We derive the optimal time-varying parameter for a monocomponent signal based on the LFM model in Section 3. In Section 4, we establish the LFM model-based well-separated condition for multicomponent signals. In Section 5 we propose a localized optimization method on the selection of window parameters based on our well-separated condition. Experimental results are provided in Section 6. Finally we give the conclusion in Section 7.

3

where σ > 0 is a parameter, g(t) is a positive function in L2 (R ) with g(0 ) = 0 and having certain decaying order as t → ∞. If

1 − t2 g(t ) = √ e 2, 2π

(14)

then gσ (t) is the Gaussian window function. The parameter σ is also called the window width in the time-domain of the window function gσ (t) since the time duration gσ of gσ is σ (up to a constant): gσ = σ g , where g is the time duration of g. The parameter σ affects the shape of gσ and hence, the representation of the STFT of a signal with gσ . As mentioned in Section 1, Iatsenko et al. [30] states that for a multicomponent signal as given by (1), if STFTs Vxk−1 (t, η ) and Vxk (t, η ) of two components xk−1 and xk are mixed, then FSST cannot separate these two components xk−1 (t ) and xk (t) either. Thus it is desirable that an appropriate σ can be chosen so that the STFTs of different components do not overlap. In this paper we introduce STFT with a time-varying parameter and then establish the separability condition of a multicomponent signal based on this type of STFT. x (t, η ) (called The STFT of x(t) with a time-varying parameter V the adaptive STFT) we consider is defined by (4). One can verify x (t, η ) can be written as that V

x (t, η ) = V =



∞

−∞  ∞ −∞

x(ξ ) gσ (t ) (η − ξ )ei2π t ξ dξ





x(ξ ) g σ (t )(η − ξ ) ei2π t ξ dξ .

(15)

2. STFT and FSST with a time-varying parameter

where for a signal x(t), its Fourier transform x(ξ ) is defined by

In this section we first provide a brief review of FSST, then we propose the adaptive FSST based on the STFT with a time-varying parameter.

x (ξ ) =

2.1. Short-time Fourier transform-based synchrosqueezed transform Recall that Vx (t, η) is the STFT of x(t) defined by (2), which can be extended to a slowly growing x(t) provided that the window function h(t) is in the Schwarz class S. The idea of FSST is to reassign the frequency variable. As in [26], for a signal x(t), at (t, η) for which Vx (t, η ) = 0, denote

ωx (t, η ) =

∂ ∂ t Vx (t, η ) .

(11)

2π iVx (t, η ) Aei2π ct ,

When x(t ) = where A, c are constants with c > 0, then ωx (t, η) is exactly c, the IF of x(t). The quantity ωx (t, η) is called the “phase transformation” [25]. FSST is to reassign the frequency variable η by transforming STFT Vx (t, η) of x(t) to a quantity, denoted by Rx (t, η), on the TF plane:

Rx (t, ξ ) :=



{ζ :Vx (t,ζ ) =0}

Vx (t, ζ )δ



 ωx (t, ζ ) − ξ dζ ,

∞ −∞

x(t )e−i2π ξ t dt .

x (t, η ) as shown in We can obtain that x(t) can be recovered from V the following theorem. x (t, η ) be the time-varying STFT of x(t ) ∈ L (R ) deTheorem 1. Let V 2 fined by (4). Suppose x, g ∈ L1 (R ). Then

x(t ) =

σ (t ) g( 0 )



∞ −∞

x (t, η )dη. V

(16)

If in addition g(t) is real-valued, then for a real-valued x(t), we have

x(t ) =

2σ (t ) Re g( 0 )



∞ 0



x (t, η )dη . V

(17)

The proof of Theorem 1 is presented in Appendix.

2.3. Adaptive FSST with a time-varying parameter

(12)

where ξ is the frequency variable. For a multicomponent signal x(t) given by (1), when Ak (t), φ k (t) satisfy certain conditions, each component xk (t) can be recovered from its FSST:

 1 xk (t ) ≈ Rx (t, ξ )dξ , h(0 ) |ξ −φk (t )|<



(13)

Next we introduce the synchrosqueezing transform (SST) associated with the adaptive STFT. First we need to define the phase ad p transformation ωx associated with the adaptive STFT. In the following we use g(τ ) replaced by τ g (τ ), namely,

τ g (τ ) (t, η ) := V x



∞ −∞

x(t + τ )

τ   τ  −i2π ητ g e dτ . σ 2 (t ) σ (t )

for certain > 0. For more mathematically precise definition of FSST and the conditions on Ak (t), φ k (t) for (13), see [26–28].

To define the phase transformation ωx s(t ) = Aei2π ct . From

2.2. Adaptive STFT with a time-varying parameter

s (t, η ) = V

We consider the window function given by

t  1 gσ (t ) = g , σ σ

ad p



∞

s(t + τ )gσ (t ) (τ )e−i2π ητ dτ

−∞  ∞

=A

−∞

ei2π c(t+τ )

, we first consider

 τ  e−i2π ητ dτ , σ (t )

1 g σ (t )

(18)

4

L. Li, H. Cai and H. Han et al. / Signal Processing 166 (2020) 107231

we have

One can use the following formula to recover the kth component xk (t) of a multicomponent signal from the adaptive FSST:

 τ  ∂ 1 Vs (t, η ) = A (i2π c )ei2π c(t+τ ) g e−i2πητ dτ ∂t σ (t ) σ (t ) −∞

  ∞ σ  (t )  τ  −i2πητ i2π c (t+τ ) +A e − g e dτ σ (t ) σ (t )2 −∞

  ∞  τ  σ  (t )τ 

∞

i2π c (t+τ )

g

xk (t ) =

e

−∞

−

σ (t )3

s (t, η ) i2 π V



∂ ∂ t Vs (t, η )



1 2 s(t ) = Aei2π φ (t ) = Aei2π (ct + 2 rt )

s (t, η ) i2π V

x

adaptive STFT defined by (4) with g(τ ) replaced by τ g (τ ). In the τ g(τ ) (t, η ) to denote the adaptive STFT defined following, we use V x by (4) with g(τ ) replaced by τ g(τ ). That is,

(19) x (t, η ) = 0, the quanHence, for a general x(t), at (t, η) for which V tity in the right-hand side of the above equation is a good candidate for the IF of x. This quantity is also called the phase transforad p mation, and we denote it by ωx (t, η ):

ωxadp,2nd (t, η ) =

ω

ad p x

⎧  ∂ V (t,η )  ⎪ Re i2∂πt Vx (t,η ) + ⎪ x ⎪ ⎪ ⎨ ⎪ ⎪  ⎪ ⎪ ⎩Re

∂  ∂ t Vx (t,η )

x (t,η ) i2π V





  Vτx g (τ ) (t,η ) σ  (t ) σ (t ) Re i2π  Vx (t,η )

 + σσ ((tt)) Re

 τ g (τ )



− Re

∂ if ∂η

Vx (t,η ) x (t,η ) i2π V



τ g(τ ) (t, η ) := V x

Vx (t,η ) P i2π  Vx (t,η ) 0

Vx (t,η ) x (t,η ) V

∂ , if ∂η



 τ g(τ )



{η∈R: Vx (t,η ) =0}

x (t, η )δ V

 adp  ωx (t, η ) − ξ dη,

(21)

where ξ is the frequency variable. The reconstruction formulas in (16) and (17) lead to that x(t) can be reconstructed from its adaptive FSST:

x(t ) =

σ (t ) g( 0 )



∞ −∞

Radp x (t, ξ )d ξ ;

(22)

⎧  ∂   τ g (τ )  V (t,η ) ⎪ ⎨Re i2∂πt Vxx (t,η ) − Re i2Vxπ Vx((tt,,ηη)) p0 (t, η ) , ωx2nd (t, η ) =  ∂  ⎪ ⎩Re i2∂πt VVx ((t,t,ηη)) , x

∂ if ∂η ∂ if ∂η

 

τ  τ  −i2π ητ g e dτ . σ (t ) σ (t ) 2

(26)

(27)

x (t, η ) = 0, = 0, V

where

P0 (t, η ) =

1

∂ ∂η

 τ g(τ )

Vx (t,η ) x (t,η ) V



∂  ∂∂t Vx (t, η )  ∂η Vx (t, η )

 σ  (t ) ∂  Vxτ g (τ ) (t, η )  + . x (t, η ) σ (t ) ∂η V

The FSST with a time-varying parameter (called the adaptive FSST of x(t)) is defined by

Radp x (t, ξ ) :=

−∞

x(t + τ )

 ( t, η ) ,

(20)



∞

x (t, η ) = 0; = 0 and V

Vx (t,η ) x (t,η ) V

      τ g (τ ) (t, η ) ∂t Vx (t, η ) σ  (t ) V x x (t, η ) = 0. (t, η ) = Re + Re , for V x (t, η ) x (t, η ) σ (t ) i2π V i2π V



For a signal x(t), we define the phase transformation for the 2nd-order adaptive FSST as

 τ g(τ )

 τ g(τ )

(25)

with phase function φ (t ) = ct + 12 rt 2 , IF φ  (t ) = c + rt and chirp rate φ  (t ) = r. τ g (τ ) (t, η ) to denote the Recall that in Section 2.3, we use V

τ g (τ )

 σ  (t ) V (t, η ) s (t, η ) = 0. Re s , for V s (t, η ) σ (t ) i2π V

+

(24)

The 2nd-order FSST was introduced in [32]. The main idea is to define a new phase transformation ωx2nd such that when x(t) is a linear frequency modulation (LFM) signal, then ωx2nd is exactly the IF of x(t). We say s(t) is an LFM signal or a linear chirp if

σ  (t ) σ  (t ) Vsτ g (τ ) (t, η ) − . i2π σ (t ) σ (t ) i2π Vs (t, η )



2.4. Second-order adaptive FSST

Therefore, the IF of s(t), which is c, can be obtained by

c = Re

|ξ −φk (t )|< 1

Radp x (t, ξ )d η

for some 1 > 0.

s (t, η ) = 0, we have Thus, if V ∂ ∂ t Vs (t, η ) = c −



−i2πητ

e dτ σ (t )   s (t, η ) − σ (t ) V s (t, η ) − σ (t ) V τ g (τ ) (t, η ). = i2 π c V σ (t ) σ (t ) s +A

2σ (t ) Re g( 0 )

(28)

Then we have the following theorem with its proof given in Appendix. Theorem 2. If x(t) is an LFM signal given by (25), then at (t, η)  τ g( τ )  (t,η ) ∂ Vx x (t, η ) = 0, ωad p,2nd (t, η ) defined by where ∂η = 0 and V x  Vx (t,η )

(27) is the IF of x(t), namely ωx

ad p,2nd

(t, η ) = c + rt.

Observe that when σ (t) ≡ σ is a constant ωxad p,2nd (t, η ) is reduced to ωx2nd (t, η ) given by Vxτ g(τ ) (t,η ) Vx (t,η ) Vxτ g(τ ) (t,η ) Vx (t,η )

 

function,

= 0, Vx (t, η ) = 0; (29)

= 0, Vx (t, η ) = 0,

where and if in addition g(t) is real-valued, then for real-valued x(t), we have

x(t ) =

2σ (t ) Re g( 0 )



0

∞



Radp x (t, ξ )d ξ .

(23)

p0 (t, η ) =

∂ ∂η



1 τ g( τ )

Vx (t,η ) Vx (t,η )



∂  ∂∂t Vx (t, η )  . ∂η Vx (t, η )

ωx2nd in (29) is one of the phase transformations considered in [36] for the conventional 2nd-order FSST.

L. Li, H. Cai and H. Han et al. / Signal Processing 166 (2020) 107231

With the phase transformation ωx (t, η ) in (27), we define the 2nd-order FSST with a time-varying parameter, called the 2ndorder adaptive FSST, of a signal x(t) as in (21): ad p,2nd

2nd Radp, x

(ξ , t ) :=

 {η∈R: Vx (t,η ) =0}

  x (t, η )δ ωadp,2nd (t, η ) − ξ dη, V x (30)

where ξ is the frequency variable. We also have the reconstruction formulas for x(t) and xk (t) similar to (22), (23) and (24) with ad p ad p,2nd Rx (ξ , t ) replaced by Rx (ξ , t ). Note that the conventional 2nd-order FSST is defined by

R2x nd (ξ , t ) :=



{η∈R: Vx (t,η ) =0}

  Vx (t, η )δ ωx2nd (t, η ) − ξ dη,

(31)

where one can use ωx2nd (t, η ) defined by (29) or choose one of several different ωx2nd in [32].

The parameter σ for the window function gσ affects the sharpness of the STFT of a signal. In this section, we study how the timex (t, η ) varying parameter σ (t) controls the representation of STFT V of a monocomponent signal x(t) and provide the parameter σ (t) with which STFT has the sharpest representation in the TF plane. In the next section, we will consider the following problem: under which condition (if any) for a multicomponent signal as given by x (t, η ), 1 ≤ k ≤ K (1), with a suitable choice of σ (t), the STFTs of V k are well separated. To study the sharpness of the STFT of a monocomponent signal or the separability of STFTs (including STFTs with a time-varying parameter) of different components xk of x(t), we need to consider the support zone of STFT Vxk (t, η ) in the TF plane, the region outside which Vxk (t, η ) ≈ 0. Since the support zone of Vxk (t, η ) is determined by the support of g outside which g(ξ ) ≈ 0, first of all, we need to define the “support” of g if g is not band-limited. More precisely, for a given threshold 0 <  < 1, if | g(ξ )|/ maxξ | g( ξ ) | <  for |ξ | ≥ ξ 0 , then we say g(ξ ) is “supported” in [−ξ0 , ξ0 ]. We use L g and g = 2ξ0 to denote the length of the “support” interval of we call it the duration of g. Note that ξ0 = ξ0, depends on  . For simplicity, here and below we drop the subscript  . Also in applications,  is quite small. In the remainder of this paper, we consider g given by (14) and thus gσ (t) is the Gaussian window function defined by

gσ ( τ ) =

τ2 1 e− 2σ 2 , √ σ 2π

(32)

gσ ( ξ ) = e

.

For g given by (14), | g( ξ ) | = where

α=

(33) 2 2 e−2π ξ

1  2 ln(1/ ). 2π

Vs (t, η ) = √

A 1 − i2 π σ 2 r

ei2π (ct +rt

2

/2 )

h



 η − (c + rt ) ,

(36)

where

h (ξ ) = e

−

2π 2 σ 2 1+(2π rσ 2 )2

(1+i2π σ 2 r )ξ 2

.

One can obtain (36) by applying the following formula (see [1]): for real α and β with α > 0



∞

e

−(α +iβ )t 2 +iωt

−∞

dt =



√

π

ω2

α + iβ

e− 4(α+iβ ) .

(37)

Observe that |h(ξ )| is a Gaussian function with duration



1 + ( 2π r σ 2 )2

σ2



= 2α

1

σ2

+ ( 2π r σ )2 .

Thus the ridge of Vs (t, η) concentrates around η = c + rt in the TF plane, and Vs (t, η) lies within the zone of TF plane of (t, η):

1 1 − L|h| ≤ c + rt − η ≤ L|h| , 2 2 or equivalently

 c + rt − α



1

σ2

+ (2π rσ )2 ≤ η ≤ c + rt + α

1

σ2

+ ( 2π r σ )2 . (38)

1

L|h| gains its minimum when σ 2 =

σ=

1 2π |r |

=



1 2π |φ  (t )|

( 2π r σ )2 ,

namely,

.

(39)

The choice of σ given in (39) results in the sharpest representation of Vs (t, η). For a monocomponent signal x(t ) = A(t )ei2π φ (t ) , if its STFT with gσ , which is also given by (refer to (3)),

Vx (t, η ) =



∞ −∞

A(t + τ )ei2π φ (t+τ ) gσ (τ )e−i2π ητ dτ

can be well approximated by

Vx (t, η ) ≈



∞ −∞

A(t )e



i2π φ (t )+φ  (t )τ + 12 φ  (t )τ 2



gσ (τ )e−i2π ητ dτ ,

then the choice of σ given

with its Fourier transform given by −2π 2 σ 2 ξ 2

Proposition 1. Let s(t) be an LFM given by (25). The STFT of s(t) with the Gaussian window function gσ (τ ) is given by

L|h| = 2α

3. Support zones of STFTs of LFM signals

5

<  if and only if |ξ | > α ,

(34)

Thus we regard that g is “supported” in [−α , α ]. Hence, gσ , given 2α by (33), is “supported” in [− σα , σα ], and L gσ = σ . For s(t ) = Aei2π ct , since its STFT with gσ is

σ=

1 2π |φ  (t )|

(40)

results in the sharpest representation of Vx (t, η). The choice of σ = σ (t ) in (40) coincides with the result derived in [1]. Observe that σ in (40) is the optimal parameter for the representation of a monocomponent signal. In the next section, we will use the obtained TF zone in (38) for the STFT of an LFM signal to study the well-separated condition for a multicomponent signal.

Vs (t, η ) = Aei2π tc gσ ( η − c ) , and gσ (η − c ) is “supported” in c − σα ≤ η ≤ c + σα , Vs (t, η) concentrates around η = c and lies within the zone (a strip) of the TF plane (t, η):



(t, η ) :

c−

 α α ≤η≤c+ , t ∈R . σ σ

(35)

Next we consider LFM signals with IF φ  (t ) = c + rt > 0. First we find the STFT of an LFM signal.

4. Separability of multicomponent signals and selection of time-varying parameter In this section, we will consider the problem that under which condition (if any), for a multicomponent signal as given by (1), x (t, η ), 1 ≤ k ≤ K of differwith a suitable choice of σ (t), STFTs V k ent components xk defined in (4) are well separated, and the associated adaptive FSST of x(t) has a sharp representation.

6

L. Li, H. Cai and H. Han et al. / Signal Processing 166 (2020) 107231

From (38), STFT Vxk (t, η ) of xk with Gaussian window function gσ lies within the zone of TF plane (t, η):

4.1. Sinusoidal signal model First we consider the sinusoidal signal model. Recall that the STFT of s(t ) = Aei2π ct with gσ (t) is supported in the zone of the TF plane given by (35). Suppose x(t) is a finite summation of sinusoidal signals:

x(t ) =

K 

Ak e i2π c k t ,

(41)



ck + rk t − α

α α ≤ ck − , σ σ

or equivalently

2α , ck − ck−1

for k = 2, 3, · · · , K.

More generally, for x(t) given by (1), suppose each xk (t) is well approximated by sinusoidal functions at any local time, that is x (t, η ) of x (t) can be well(7) holds. Then the time-varying STFT V k k approximated by

x (t, η ) ≈ V k



∞ −∞

xk (t )e

i2πφk (t )τ



gσ (t ) (τ )e

−i2πητ



for all t. Thus xk−1 (t ) and xk (t) are separable in the TF plane if



(t, η ) : φk (t ) −



 η − φk (t ) ≤ α}

2α

k−1

k

(t )

, t ∈ R, k = 2, 3, · · · , K.

σ1 (t ) = max

2≤k≤K

The condition that (45) holds for k = 2, · · · , K is the wellseparated condition for a multicomponent signal consisting of LFM signals. One of the main goals of this paper to obtain an explicit σ (t) such that Vxk (t, η ), 1 ≤ k ≤ K lie within non-overlapping TF zones. To this end, we replace the TF zone of Vxk (t, η ) in (44) by a larger zone for Vxk (t, η ) by using σ1 + 2π |rk |σ to replace



+ (2π rk σ )2 in (44):

1

σ2

ck + rk t − α

1

+ 2π |rk |σ

σ

ck−1 + rk−1t + α



≤ η ≤ ck + r k t + α

1 σ



+ 2π |rk |σ .



∞





≤ ck + rk t − α

1 σ



+ 2π |rk |σ .

Ak (t + τ )ei2π φk (t+τ ) gσ (t ) (τ )e−i2π ητ dτ ,

xk (t )

−

1 − i2π σ (t )2 φk (t )

e

1

σ (t )2

2π 2 +(2πφ  (t )σ (t ))2 k

(1+i2π σ 2 (t )φk (t ))(η−φk (t ))2

.

x (t, η ) lies within the zone of TF plane: Thus V k



φk (t ) −

α

1

σ (t )2

+ (2π φk (t )σ (t ))2 ≤ η



≤ φk (t ) + α

(43)

1

σ (t )2

+ (2π φk (t )σ (t ))2 ,

(47)

for t ∈ R, and the well-separable condition for x(t) is



φ

k−1

(t ) + α

In this subsection we will derive the well-separated condition based on the LFM model. More precisely, we consider x(t ) = K k=1 xk (t ), where each xk (t) is an LFM signal, namely,

1

σ (t )2

+ (2π φk−1 (t )σ (t ))2

≤ φk (t ) − α

1 2 xk (t ) = Ak ei2π (ck t+ 2 rk t )

with the phase function φk (t ) = ck t + 12 rk t 2 and φk −1 (t ) < φk (t ).

σ

+ 2π |rk−1 |σ

can be well-approximated by the quantity on the right-hand side of (10), which is (by applying Proposition 1 or (37))



4.2. Linear frequency modulation (LFM) model

1

More generally, for x(t) given by (1), suppose each xk is well approximated by an LFM at any local time, namely (9) holds. Then x (t, η ) of x with g , which is (refer to the time-varying STFT V k σ (t) k (5))

−∞

2α . φk (t ) − φk −1 (t )

+ ( 2π rk σ )2 .

(45)

(42)

(42) is the sinusoidal signal model-based well-separable condition for x(t) with the adaptive STFT. When σ (t) ≡ σ is a positive constant function, (42) is reduced to (8) with  = α /σ . We observe in our experiments that in general a big σ will result in low time-resolution and unreliable representation of the FSST of a signal x(t). Actually, the error bounds derived in [25– 28,35] imply that for a signal, its synchrosqueezed representation is sharper when the window width in the time domain of the window function gσ , which is σ (up to a constant), is smaller. Thus we should choose σ (t) as small as possible. Hence, we propose the sinusoidal signal model-based choice for σ , denoted by σ 1 (t), to be



1

σ2

the zone given by (46) is slightly larger than that given by (44). Clearly, xk−1 (t ) and xk (t) are separable in the TF plane if

α α  ≤ η ≤ φk (t ) + . σ (t ) σ (t )

α α ≤ φk (t ) − , t ∈ R, k = 2, · · · , K, σ (t ) σ (t )

φ  (t ) − φ 

+ (2π rk−1 σ ) ≤ ck + rk t − α

√    1 1 2 1 + 2π |rk |σ ≤ + ( 2π rk σ )2 ≤ + 2π |rk |σ , 2 σ σ σ2

or equivalently

σ (t ) ≥

σ2

2

Since

Thus, the components of x(t) will be well-separated in the TF plane (namely, Ok , 1 ≤ k ≤ K do not overlap) if

φk −1 (t ) +

ck−1 + rk−1 t + α



1

(46)

x (t, η ) lies within the zone of the TF plane (t, η): Hence V k

=

+ ( 2π rk σ )2 ,

dτ

= xk (t ) g σ (t )(η − φk (t )) . Ok = {(t, η ) : −α ≤ σ (t )

1

σ2

(44)



where Ak , ck are positive constants with 0 < ck < ck+1 . Since the STFT of the k-component of x(t) lies within the zone of the TF plane (t, η): ck − α /σ ≤ η ≤ ck + α /σ for any t, the components of x(t) will be well-separated in the TF plane if

σ≥

+ ( 2π rk σ ) ≤ η ≤ ck + rk t + α 2

σ2

k=1

ck−1 +



1



1

σ (t )2

+ (2π φk (t )σ (t ))2 , 2 ≤ k ≤ K, (48)

for t ∈ R.

L. Li, H. Cai and H. Han et al. / Signal Processing 166 (2020) 107231

x (t, η ) by a larger As above, we replace the TF zone (47) of V k zone given by

  1 + 2π |φk (t )|σ (t ) ≤ η ≤ φk (t ) σ (t )   1 +α + 2π |φk (t )|σ (t ) . σ (t )

φk (t ) − α

plane. As discussed above, since a smaller σ (t) gives a sharper synchrosqueezing representation, we should choose σ (t) as small as possible. Hence, we propose the LFM model-based choice for σ , denoted by σ 2 (t), to be

σ2 (t ) = max

Then the corresponding well-separable condition for x(t) is

7



bk (t ) +



4α bk (t )2 − 8α ak (t )



: 2≤k≤K ,

(55)

where ak (t) and bk (t) are defined by (51), and α is defined by (34). Next we show some experimental results. We consider a twocomponent LFM signal,

  1  (t )|σ (t ) ( t ) + α + 2 π | φ k−1 k−1 σ (t )   1 ≤ φk (t ) − α + 2π |φk (t )|σ (t ) , 2 ≤ k ≤ K, σ (t )

φ

y(t ) = y1 (t ) + y2 (t ) (49)

= cos(2π (12t + 25t 2 )) + cos(2π (34t + 32t 2 )),

t ∈ [0, 1], (56)

which is equivalent to

ak (t )σ (t )2 − bk (t )σ (t ) + 2α ≤ 0, 2 ≤ k ≤ K,

(50)

where

ak (t ) = 2π α (|φk−1 (t )| + |φk (t )| ), bk (t ) = φk (t ) − φk −1 (t ). (51) If

 φk (t ) − φk −1 (t ) 2   − 16π α 2 |φk (t )| + |φk−1 (t )| ≥ 0,

bk (t )2 − 8α ak (t ) =



then (49) (or (50)) is equivalent to

bk (t ) +



4α bk (t ) − 8α ak (t ) 2

≤ σ (t ) ≤

bk (t ) −

2 ≤ k ≤ K,



4α bk (t )2 − 8α ak (t )

,

(52)

for t ∈ R. Otherwise, if bk (t )2 − 8α ak (t ) < 0, then there is no suitable solution of the parameter σ for (49) or equivalently (50). In this case we say that components xk−1 (t ) and xk (t) of multicomponent signal x(t) cannot be separated in the TF plane. Note that when ak (t ) = 0, i.e. φk (t ) = φk−1 (t ) = 0, (52) is reduced to (42). In the next theorem, we summarize the LFM model-based wellseparated condition we have derived above.  Theorem 3. Let x(t ) = Kk=1 xk (t ), where each xk (t ) = Ak (t )ei2π φk (t ) x (t, η ) with g is an LFM signal or its adaptive STFT V σ (t) can be well k approximated by (10), and φk −1 (t ) < φk (t ). If

4α

 π |φk (t )| + |φk−1 (t )| ≤ φk (t ) − φk −1 (t ), k = 2, · · · , K,

√

and

(53)

 max

2≤k≤K

bk (t ) +

≤ min

2≤k≤K







4α bk (t )2 − 8α ak (t )

bk (t ) −



4α bk (t )2 − 8α ak (t )



,

(54)

for t ∈ R, then the components of x(t) are well-separable in TF plane x (t, η ), 1 ≤ k ≤ K with σ (t) chosen to satisfy in the sense that V k (52) lie in non-overlapping regions in the TF plane. We call (53)–(54) the LFM model-based well-separated condition for a multicomponent signal x(t). Observe that our LFM model-based well-separated condition (53) requires the boundedness of the 2nd-order derivatives φk (t ), while it seems the sinusoidal signal-based well-separated condition (8) or (42) does not have such a constraint. Actually the sinusoidal signal model assumption (7) requires φk (t ) be small. In addition, to make the recovery error in (13) small, φk (t ) must be very small (see [27,28,35] for the details about the recovery error estimates). Any σ (t) between the two quantities in the two sides of the inequality (54) can separate the components of x(t) in the TF

where the number of sampling points is 256, namely the sampling rate is 256 Hz. The IFs of y1 (t) and y2 (t) are φ  (t ) = 12 + 50t and φ2 (t ) = 34 + 64t, respectively. The top-left panel of Fig. 1 shows the instantaneous frequencies of y1 (t) and y2 (t). With σ 2 (t), both the proposed adaptive FSST defined by (21) and 2nd-order adaptive FSST defined by (30) can represent this signal sharply. Here and below, we choose  = 15 , and hence α which is defined by (34) and used in (55) is α ≈ 0.2855. Observe that the 2nd-order adaptive FSST further improves the TF energy concentration of the adaptive FSST. Here we also give the results of conventional FSST studied in [26–28], and conventional 2nd-order FSST defined in [32] with σ = 0.057. This σ is obtained by minimizing the Rényi entropy of the STFT (refer to the next section about the definition of Rényi entropy). Observe that the 2nd-order FSST is better than the FSST with the same σ . When σ = 0.057 the TF representation of the conventional 2nd-order FSST is not as sharp or clear as that of the 2nd-order adaptive FSST. 5. Selecting the time-varying parameter automatically Suppose x(t) given by (1) is separable, meaning (53) and (54) hold. If we know φk (t ) and φk (t ), then we can choose a σ (t) such as σ 2 (t) in (55) to satisfy (52) to define the adaptive STFT and adaptive FSST for sharp representations of xk (t) in the TF plane and for accurate recovery of xk (t). However in practice, we in general have no prior knowledge of φk (t ) and φk (t ). Hence, we need to have a method which provides suitable σ (t). In this section, we propose an algorithm to estimate σ (t) which is based on the wellseparated condition of (49). First for temporarily fixed t and σ , denote Vx,(t,σ ) (η ) = Vx (t, η, σ ), the STFT of x(t) with a time-varying parameter defined by (4). We extract the peaks (local maxima) of |Vx,(t,σ ) (η)| with certain height. More precisely, assuming γ 1 > 0 is a given threshold, we find local maximum points η1 , η2 , , ηm of |Vx,(t,σ ) (η)| at which |Vx,(t,σ ) (η)| attains local maxima with

|Vx,(t,σ ) (ηk )| > γ1 , k = 1, · · · , m. maxη |Vx,(t,σ ) (η )|

(57)

Note that m may depend on t and σ . We assume η1 < η2 <  < ηm . The threshold γ 1 is used to remove the local maxima with smaller amplitudes, which are regarded as noises and interferences. For each local maximum point ηk , we regard ηk is the local maximum of the adaptive STFT Vxk ,(t,σ ) (η ) of a potential component, denoted by xk (t) of x(t). To check whether xk is indeed a component of x(t) or not, we consider the support interval [lk , hk ] for Vxk ,(t,σ ) (η ) with |Vxk ,(t,σ ) (η )| > 0 for η ∈ [lk , hk ]. If there is no overlap among [lk , hk ], [lk−1 , hk−1 ], [lk+1 , hk+1 ], then we decide that xk (t) is indeed a component of x(t), where [lk−1 , hk−1 ], [lk+1 , hk+1 ] are the support intervals for xk−1 and xk+1 defined similarly. With

8

L. Li, H. Cai and H. Han et al. / Signal Processing 166 (2020) 107231

Fig. 1. Experimental results on the two-component LFM signal in (56): IFs (Top-left); adaptive FSST with time-varying parameter σ 2 (t) (Top-middle); 2nd-order adaptive FSST with time-varying parameter σ 2 (t) (Top-right); conventional FSST (Bottom-left) and conventional 2nd-order FSST (Bottom-right) with σ = 0.057.

our LFM model, if the estimated IF φk (t ) of xk (t) is ck + rk t, then by (49),

(58) and (59). Especially when rk = 0, recalling the support zone of a sinusoidal signal mode in (35), we have

hk = ck + α

(58)

hk = ck +

(59)

η, σ ) for fixed t and σ :

1 σ

lk = ck−1 − α



+ 2π | rk |σ ,

1 σ

This way we obtain the collection of support intervals for Vx (t,



+ 2π | rk−1 |σ .

Notice that ck = ηk . Thus we need to estimate the chirp rate rk of xk (t). To this end, we extract a small piece of curve in the TF plane passing through (t, ηk ) which corresponds to the local ridge on |Vxk ,(t,σ ) (η )|. More precisely, letting

tk1 = t −

1 Lg = t − 2π ασ , 2 σ

tk2 = t −

1 Lg = t + 2π ασ , 2 σ

define

d k (τ ) =

argmax

η: η is near ηk

|V(τ ,σ ) (η )|, τ ∈ [tk1 , tk2 ].

In the above we have used the fact that the duration of gσ (t) is (refer to (34))

Lgσ = 2σ



α α , l = ck − . σ k σ

2 ln(1/ ) = 4π σ α .

Note that d k (t ) = ηk and (t, ηk ) is a point lying on the curve in the TF of (τ , η) given by

L = {(τ , d k (τ )) : τ ∈ [tk1 , tk2 ]} = {(τ , η ) : η = d k (τ ), τ ∈ [tk1 , tk2 ]}. Most importantly, {|Vxk ,(τ ,σ ) (η )| : (τ , η ) ∈ L} is the local ridge on |Vxk ,(τ ,σ ) (η )| near (t, ηk ), and thus, it is also the local ridge on |Vx (t, η, σ )|. Observe that from the STFT of an LFM given by Proposition 1, the local ridge on |Vxk (t, η, σ )| occurs when φk (t ) = ck + rk t. Thus we use the linear function

dk ( τ ) = rk ( τ − t ) + ck , τ ∈ [tk1 , tk2 ] to fit d k (τ ). The obtained rk is the estimated chirp rate rk of xk (t). With this rk and ck = ηk as given above, we have hk , lk given in

s = {[l1 , h1 ], · · · , [lm , hm ]}.

(60)

If adjacent intervals of s do not overlap, namely,

hk ≤ lk+1 , for all k = 1, 2, · · · , m − 1

(61)

holds, then this σ is a right parameter to separate the components and such a σ is a good candidate which we consider to select. Otherwise, if a pair of adjacent intervals of s overlap, namely, (61) does not hold, then this σ is not the parameter we shall choose and we need to consider a different σ . In the above description of our idea for the algorithm, we start with a σ and (temporarily fixed) t, then we decide whether this σ is a good candidate to select or not based our proposed criterion: (61) holds or does not. The choice of the initial σ plays a critical role for the success of our algorithm due to that on one hand, as we have mentioned above, in general a smaller σ will result in a sharper representation of SST, and hence, we should find σ as small as possible such that (61) holds; and on the other hand, different σ with which (61) holds may result in different number of intervals m in (60) even for the same time instance t. To keep the number m (an estimation of the number of modes K for a given time t) unchanged when we search for different σ with a fixed t, the initial σ is required to provide a good estimate of the number of the components of a multicomponent signal x(t). To this end, in this paper we propose to use the Rényi entropy to determine the initial σ (t). The Rényi entropy is a commonly used measurement to evaluate the concentration of a TF representation such as STFT, SST, etc. of a signal of x(t), see [62,65,66]. Taking STFT Vx (t, η) of a signal

L. Li, H. Cai and H. Han et al. / Signal Processing 166 (2020) 107231

x(t) as an example, the Rényi entropy with Vx (t, η) is

Eζ (t ) :=

1 log2 1−

|Vx (b, η )|2 dηdb  , t+ζ  ∞ 2 V ( b, η ) d η d b | | x t−ζ 0 t−ζ

0

(62)

where  is usually greater than 2. In this paper we choose  = 2.5, a common value used in other papers, see for examples [62,63]. Parameter ζ > 0 determines the local duration to calculate the local Rényi entropy. We choose ζ = 4. One may choose some large ζ for non-stationary signals with slow-varying IFs. Note that the smaller the Rényi entropy, the better the TF resolution. So for a fixed time t, we can use (62) to find a σ (denoted as σ u (t)) with the best TF concentration of Vx (t, η, σ ), where Vx (t, η, σ ) is the regular STFT of x(t) defined by (2) with the window function h(τ ) = gσ (τ ) = 1 τ σ g( σ ) having a parameter σ . More precisely, replacing Vx (b, η) in (62) by Vx (b, η, σ ), we define the Rényi entropy Eζ (t, σ ) of Vx (t, η, σ ), and then, obtain

  σu (t ) = argmin Eζ (t, σ ) .

(63)

σ >0

We set σ u (t) as the upper bound of σ (t) for a fixed t. With these discussions, we propose an algorithm to estimate σ (t) as follows. Algorithm 1. (Separability parameter estimation) Let {σ j , j = 1, 2, · · · , n} be an uniform discretization of σ with σ 1 > σ 2 >  > σ n > 0 and sampling step σ = σ j−1 − σ j . The discrete sequence s(t), t = t1 , t2 , · · · , tN (or t = 0, 1, · · · , N − 1) is the signal to be analyzed. Step 1. Let t be one of t1 , t2 , , tN . Find σ u in (63) with σ ranging over {σ j , j = 1, 2, · · · , n}. Step 2. Let s be the set of the intervals given by (60) with σ = σu . Let z = σu . If (61) holds, then go to Step 3. Otherwise, go to Step 5. Step 3. Let σ = z − σ . If the number of intervals m in (60) with this new σ remains unchanged, σ ≥ σ n and (61) holds, then go to Step 4. Otherwise, go to Step 5. Step 4. Repeat Step 3 with z = σ . Step 5. Let C (t ) = z, and do Step 1 to Step 4 for different time t of t1 , t2 , , tN . Step 6. Smooth C(t) with a low-pass filter B(t):

σest (t ) = (C ∗ B )(t ).

(64)

We call σ est (t) the estimation of the separability time-varying parameter σ 2 (t) in (55). In Step 6, we use a low-pass filter B(t) to smooth C(t). This is because of the assumption of the continuity condition for Ak (t) and φ k (t). With the estimated σ est (t), we can define the adaptive STFT, the adaptive FSST and the 2nd-order adaptive FSST with a time-varying parameter σ (t ) = σest (t ). In [63], the Rényi entropy-based optimal time-varying window was proposed for the sharp representation of SST. More precisely, let Rx (b, ξ , σ ) and R2x nd (b, ξ , σ ) be the regular FSST and the regular 2nd-order FSST of x(t) (with the phase transformation ωx2nd (a, b) given in [32]) defined by (12) and (31) respectively with the window function h(t ) = gσ (t ) given by (32) containing σ > 0. Denote SST (t ) the Rényi entropies of Rx (b, ξ , σ ) and R2x nd (b, ξ , σ ) by E, ζ ,σ

SST 2 (t ) respectively, which are defined by (62) with V (b, ξ ) and E, x ζ ,σ to be replaced by Rx (b, ξ , σ ) and R2x nd (b, ξ , σ ) for certain fixed , ζ . The optimal time-varying parameter is obtained by minimizing SST (t ) and E SST 2 (t ): E, ζ ,σ ,ζ ,σ

σRe (t ) = argmin E,SSTζ , σ (t ), σRe2 (t ) = argmin E,SSTζ 2, σ (t ). σ >0

σ >0

the phase transformation ωx (t, ξ ) in (20) replaced by the regular phase transformation ωx (t, ξ ) defined by the formula (11) for the conventional FSST. Similarly, the 2nd-order time-varying-window FSST with σ = σRe2 (t ) in [63] is defined by (30) but with the phase ad p,2nd transformation ωx (t, ξ ) in (27) replaced by a regular phase transformation ωx2nd (t, ξ ) defined by a formula in [32] for the conventional 2nd-order FSST. With PT representing phase transformation, we call them the regular-PT adaptive FSST and the 2nd-order regular-PT adaptive FSST, respectively. We use the Rényi entropy-based adaptive FSST and our proposed adaptive FSST with σ = σest (t ) to process the twocomponent linear chirp signal in (56). The different time-varying parameters are shown in the top-left panel of Fig. 2, where σ 1 (t), σ 2 (t), σ u (t), σ est (t), σ Re (t) and σ Re2 (t) are defined by (43), (55), (63), (64) and (65), respectively. Here we let σ ∈ [0.001, 0.2] with σ = 0.001, namely σ1 = 0.2 in Algorithm 1. We set  = 2.5, ζ = 4 (sampling points, for discrete signal) and γ 1 in (57) to be 0.3. Note that we set the same values of , ζ , and γ 1 for all the following experiments. We use a simple rectangular window B = {1/5, 1/5, 1/5, 1/5, 1/5} as the low-pass filter. One can use some other filters, such as an FIR filter or a window of Gaussian or Hamming. Note that the length of the filter we use is 5, which is related to the parameter ζ = 4. And we use a constant  = 1/5 in (34), namely constant α in (59). The estimation σ est (t) by Algorithm 1 is very close to σ 2 (t) except for the start near t = 0. So the estimation algorithm is an efficient method to estimate the well-separation time-varying parameter σ 2 (t). Fig. 2 shows the proposed adaptive FSST and 2nd-order adaptive FSST with σ est (t). The proposed 2nd-order adaptive FSST gives energy concentration. In Fig. 2, we also provide the regular-PT adaptive FSST with σ = σRe (t ) and the 2nd-order regular-PT adaptive FSST with σ Re2 (t) as described above. The regular-PT adaptive FSST performs well in the TF energy concentration of this two-component signal. The Matlab routines for Algorithm 1, FSST, the adaptive FSST and regular-PT adaptive FSST can be downloaded at the website of one of the authors: www.math.umsl.edu/∼jiang. For most well-separated signals, Algorithm 1 results in a suitable σ est (t) with which the 2nd-order adaptive FSST is clear, sharp and concentrated. However, when IFs of different components are too close, then two adjacent components at Step 2 of Algorithm 1 may merge into one, which results in component mixing. In addition, the regular-PT adaptive FSST method is unable to separate such components either, see an experimental example in the next section. To tackle this problem, we propose to use a varying  or α in (34), which defines the bandwidth of g(ξ ) and hence determines the support zones of STFTs. Although a greater  may result in a larger recovery error, some components with extremely close IFs can be separated with a large  . Suppose  ∈ [ s ,  o ] for some 0 <  s <  o < 1. Our method is first we choose the maximum  = o for fixed t and σ first, and obtain the support intervals s in (60) satisfying (61). Then we decrease  step by step. This way the support intervals in (60) will increase gradually. We stop our procedure when  reaches the minimum value  s or the condition in (61) does not hold. The following is the revised algorithm to estimate σ (t). ad p

 t+ζ  ∞



9

(65)

With σ Re (t) and σ Re2 (t) obtained by (65), Sheu et al. [63] defines the time-varying-window FSST with σ = σRe (t ) by (21) but with

Algorithm 2. Let {σ j , j = 1, 2, · · · , n} be an uniform discretization of σ with σ 1 > σ 2 >  > σ n > 0 and sampling step σ = σ j−1 − σ j . Let { j , j = 1, 2, · · · , m} be an uniform discretization of  with  1 >  2 >  >  m > 0 and sampling step  =  j−1 −  j . The discrete sequence s(t), t = t1 , t2 , · · · , tN (or t = 0, 1, · · · , N − 1) is the signal to be analyzed. Step 1. Let t be one of t1 , t2 , , tN . Find σ u in (63) with σ ∈ {σ j , j = 1, 2, · · · , n}.

10

L. Li, H. Cai and H. Han et al. / Signal Processing 166 (2020) 107231

Fig. 2. Experimental results of different time-varying parameters for the two-component LFM signal in (56): various time-varying parameters (Top-left); adaptive FSST with σ est (t) (Top-middle) and 2nd-order adaptive FSST with σ est (t) (Top-right); regular-PT adaptive FSST with σ Re (t) (Bottom-left) and 2nd-order regular-PT adaptive FSST with σ Re2 (t) (Bottom-right).

Step 2. Let s be the set of the intervals given by (60) with σ = σu and  = 1 . Let z = σu . Step 3. If (61) holds and  >  m , update  with  −  , and repeat Step 3. Step 4. If  >  m , go to Step 7. Otherwise, if  = m , go to Step 5. Step 5. Let σ = z − σ . If the number of intervals m in (60) with this new σ remains unchanged, σ ≥ σ n and (61) holds, then go to Step 6. Otherwise, go to Step 7. Step 6. Let  = 1 , go to Step 3 with z = σ . Step 7. Let C (t ) = z, and do Step 1–Step 6 for different time t of t1 , t2 , , tN . Step 8. Smooth C(t) with a low-pass filter B(t):

σest2 (t ) = (C ∗ B )(t ). 6. Further experiments and results In this section, we provide more numerical examples to further illustrate the effectiveness and robustness of our method in the IF estimation and component recovery. 6.1. Experiments with a three-component synthetic signal The three-component signal we consider is given by

z(t ) = z1 (t ) + z2 (t ) + z3 (t ), where







t ∈ [1/2, 1],



z2 (t ) = cos 94π t + 13 cos(4π t − π /2 ) + 110π t , z3 (t ) = cos 194π t + 112π t

6.2. Application to bat echolocation signal

(66)

z1 (t ) = cos 118π (t − 1/2 ) + 100π (t − 1/2 )2 ,



for this experiment is 512 Hz, namely we have 512 discrete samples for z(t). The IFs of the three components are φ1 (t ) = 59 + 100(t − 1/2 ), φ2 (t ) = 47 − 26 sin(4π t − π /2 ) + 110t and φ3 (t ) = 97 + 112t, respectively. Fig. 3 shows the waveform and IFs of z(t). We calculate various time-varying parameters as those shown in Fig. 2 and σ est2 (t) as well. In this experiment, to obtain σ est2 (t) with Algorithm 2, we consider a time-varying  (t) with  ∈ [0.2, 0.8] and  = 0.01. For this three-component signal, we observe that the conventional FSST, regular-PT adaptive FSST and adaptive FSST with σ est (t) cannot separate the three components well due to that the frequencies of two components are close to each other, see Fig. 4, while the 2nd-order adaptive FSST with σ est2 (t) provides quite sharp and clear representations of the three components. In Fig. 4, for the conventional FSST and conventional 2nd-order FSST, we use σ = 0.04 which is obtained by minimizing the Rényi entropy of the STFT. We also consider FSSTs in noise environment. We add Gaussian noises to the original signal given in (66) with different signal-tonoise ratios (SNRs). Fig. 5 shows the conventional 2nd-order FSST and our proposed 2nd-order adaptive FSST under different noise levels. Observe that our method works well under noisy environment.

2

2



,

t ∈ [0, 1],

t ∈ [0, 3/4].

Note that the durations of the three components in (66) are different, that is K in (1) can be time-varying. The sampling rate

In order to further verify the reliability of the proposed algorithm, we test our method on a bat echolocation signal emitted by a large brown bat. There are 400 samples with the sampling period 7 μs (sampling rate Fs ≈ 142.86 kHz). For a given real-world signal, how to select an appropriate constant σ such that the resulting conventional SST or 2nd-order SST has a sharp representation is probably not very simple. Here we choose σ = 8 × 10−5 , which is close to the mean of σ est (t) obtained by Algorithm 1. Fig. 6 shows the TF representations of the echolocation signal: STFT, conven-

L. Li, H. Cai and H. Han et al. / Signal Processing 166 (2020) 107231

11

Fig. 3. Three-component signal in (66): its waveform (Left panel) and IFs of its components (Right panel).

Fig. 4. Experimental results on the three-component signal in (66). Top-row (from left to right): conventional FSST with constant σ = 0.04, regular-PT adaptive FSST with σ Re (t) and adaptive FSST with σ est2 (t); Bottom-row (from left to right): conventional 2nd-order FSST with constant σ = 0.04, 2nd-order regular-PT adaptive FSST with σ Re2 (t) and 2nd-order adaptive FSST with σ est2 (t).

tional 2nd-order FSST with σ = 8 × 10−5 and the 2nd-order adaptive FSST with time-varying parameter σ est (t). Unlike the threecomponent signal in Fig. 4, the four components in the bat signal are much well separated. Thus, both the conventional 2nd-order FSST and the 2nd-order adaptive FSST can separate well the components of the signal. In addition, they give sharp representations in the TF plane. Comparing with the conventional 2nd-order FSST, the 2nd-order adaptive FSST with σ est (t) gives a better representation for the fourth component (the highest frequency component) and the two ends of the signal. Furthermore, σ est (t) provides a hint how to select σ for the conventional 2nd-order FSST. 6.3. Signal separation Finally we consider the separation of a multicomponent signal: to recover/reconstruct its components. We use (13), (24) and similar formulas to recover the signal components for conventional FSST and adaptive FSST. Here we use the maximum values on the FSST plane to search for the IF ridges φk (t ) one by one, see details in [31]. Then integrate around the ridges with = 1 = 15 (discrete value, unitless). We use the relative “root mean square error”

(RMSE) to evaluate the separation performance, which is defined by





K 1  zk − zˆk 2 RMSE = , K z k 2

(67)

k=1

where zˆk is the reconstructed zk , K is the number of components. We also consider signal separation in noise environment. As before we add Gaussian noises to the original signal given in (66) with different signal-to-noise ratios (SNRs). Due to the page limitation of the paper, we just provide the pictures of the reconstructed components of the three-component signal z(t) in (66) by the 2nd-order adaptive FSST with σ est2 (t) under the noiseless environment, while we provide RMSEs of four different methods, all in Fig. 7. In the bottom-right panel of Fig. 7 for RMSEs, Method 1, 2, 3 and 4 denote the 2nd-order adaptive FSST with σ est2 (t), the regular-PT adaptive FSST with σ Re (t), the 2nd-order regular-PT adaptive FSST with σ Re2 (t) and the conventional 2nd-order FSST with constant σ = 0.01, respectively. This panel gives the RMSEs of these 4 methods when SNR varies from 0 dB to 20 dB. Under each SNR, we do Monte-Carlo experiment

12

L. Li, H. Cai and H. Han et al. / Signal Processing 166 (2020) 107231

Fig. 5. FSSTs of three-component signal in (66) under different noise levels. Top row (from left to right): Conventional 2nd-order FSSTs with constant σ = 0.01 under SNRs of 5 dB, 10 dB and 15 dB. Bottom row (from left to right): 2nd-order adaptive FSSTs under SNRs of 5 dB, 10 dB and 15 dB.

Fig. 6. Example of the bat echolocation signal. Top-left: waveform; Top-right: conventional STFT with σ = 8 × 10−5 ; Bottom-left: conventional 2nd-order FSST with σ = 8 × 10−5 ; Bottom-right: 2nd-order adaptive FSST with time-varying parameter σ est (t) obtained by our proposed Algorithm 1.

L. Li, H. Cai and H. Han et al. / Signal Processing 166 (2020) 107231

13

Fig. 7. Reconstruction results of the signal in (66). Top-left, top-right and bottom-left: the reconstructed z1 (t), z2 (t) and z3 (t) by the 2nd-order adaptive FSST. Bottom-right: RMSEs under different SNRs with various methods (Method 1: 2nd-order adaptive FSST; Method 2: regular-PT adaptive FSST; Method 3: 2nd-order regular-PT adaptive FSST; Method 4: conventional 2nd-order FSST).

for 50 runs. Obviously, the reconstruction error with the 2nd-order adaptive FSST is less than those with other methods. Observe that when the noise level is high, for example SNR=0dB, RMSEs for all methods are large. This is mainly due to the fact that in a high level noise environment, the IFs of the modes are hardly estimated by the ridge detection process with the local maxima in the TF plane.

7. Conclusion In this paper, we introduce the adaptive short-time Fourier transform (STFT) with a time-varying parameter and the adaptive STFT-based synchrosqueezing transform (called the adaptive FSST). We also introduce the 2nd-order adaptive FSST. We analyze the support zones of the STFTs of linear frequency modulation (LFM) signals with the Gaussian window function. We develop the wellseparated condition for non-stationary signals by using LFM signals to approximate non-stationary signals during at local time. We propose a method to select the time-varying parameter automatically. The experimental results on both synthetic and real data demonstrate that the adaptive FSST is efficient for the instantaneous frequency estimation, sharp representation in the TF and the separation of multicomponent non-stationary signals with fastvarying frequencies. We will study the theoretical analysis of the adaptive FSST in our future work. In addition, our further study will consider other types of time-varying window functions besides the Gaussian window function. In this paper we consider signals of components without crossover IF curves. In the future, we will consider how to recover components with crossover IF curves.

Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgments The authors would like to thank the anonymous reviewers for their valuable comments. In addition, the authors wish to thank Curtis Condon, Ken White, and Al Feng of the Beckman Institute of the University of Illinois for the bat data in Fig. 6 and for permission to use it in this paper. Appendix Proof of Theorem. 1. From (15), we have



 ∞ ∞ x (t, η )dη = V x(ζ ) gσ (t ) (η − ζ )ei2π t ζ dζ dη −∞ −∞ −∞  ∞  ∞ = x ( ζ ) e i2π t ζ gσ (t ) (η − ζ )dηdζ −∞ −∞  ∞  ∞  ∞ = x ( ζ ) e i2π t ζ gσ (t ) (η )dηdζ = gσ (t ) (η )ei2π ·0·η dη ∞

−∞ ∞



−∞

−∞

g( 0 ) × x(ζ )ei2π t ζ dζ = gσ (t ) (0 )x(t ) = x(t ), σ (t ) −∞ where exchanging the order of dη and dζ follows from the Fubini’s theorem. This shows (16).

14

L. Li, H. Cai and H. Han et al. / Signal Processing 166 (2020) 107231

To prove (17), note that for real-valued x(t), since gσ (t) (τ ) is x (t, −η ) = V x (t, η ). Hence, from (16), we real-valued, we have V have

g( 0 ) x(t ) = σ (t )

 

= 

∞ 0 ∞ 0

x (t, η )dη + V x (t, η )dη + V



0



−∞ ∞



∞

0

x (t, η )dη V

φ  (t ) = c + rt = Re +

x (t, −η )dη V



∂ ∂ t Vs (t, η )



s (t, η ) i2 π V  Vτ g (τ ) (t, η ) 

− Re

 Vτ g(τ ) (t, η ) s

s (t, η ) i2 π V

P0 (t, η )



σ  (t ) Re s . s (t, η ) σ (t ) i2 π V τ g(τ ) (t,η ) V

∂ ( x Thus for a signal x(t) given by (68), at (t, η) where ∂η x (t,η ) ) = V ad p, 2 nd  x (t, η ) = 0, ω 0 and V (t, η ) defined by (27) is φ (t ) = c + rt,

∞

x (t, η )dη + x (t, η )dη V V 0  ∞  x (t, η )dη . = 2Re V =

Hence,

0

x

the IF of x(t). This shows Theorem 2.



0



Thus (17) holds.

References

Proof of Theorem. 2. Here we will show that ωs rt for s(t) given by

ad p,2nd

(t, η ) = c +

q 2 1 2 s(t ) = A(t )ei2πφ (t ) = Ae pt + 2 t ei2π (ct + 2 rt )

(68)

where p, q are two  real constants.  From s (t ) = p + qt + i2π (c + rt ) s(t ) and

s (t, η ) = V



∞

s(t + τ )

−∞

we have

∂ V (t, η ) = ∂t s



1

σ (t )

g(

τ )e−i2πητ dτ , σ (t )

1 τ g( )e−i2πητ dτ σ (t ) σ (t )  ∞ σ  (t ) τ + s(t + τ )(− ) g( )e−i2πητ dτ 2 σ (t ) σ ( t ) −∞  ∞ σ  (t )τ  τ + s(t + τ )(− )g ( )e−i2πητ dτ σ (t ) σ (t )3 −∞  s (t , η ) + (q + i2π r ) = ( p + qt + i2π (c + rt ) V  ∞ 1 τ × τ s(t + τ ) g( )e−i2πητ dτ σ ( t ) σ (t ) −∞ σ  (t )  σ  (t ) τ g (τ ) − Vs (t, η ) − V (t, η ) σ (t ) σ (t ) s  σ  (t )  = p + qt + i2π (c + rt ) − V (t, η ) σ (t ) s  τ g(τ ) (t , η ) − σ (t ) V τ g (τ ) (t , η ) + (q + i2π r )σ (t )V s σ (t ) s ∞

−∞

s (t + τ )

s (t, η ) = 0, we have Thus, if V ∂ ∂ t Vs (t, η ) = p + qt −

s (t, η ) V

×

σ  (t ) + i2π (c + rt ) + (q + i2π r )σ (t ) σ (t )

τ g(τ ) (t, η ) V s − s (t, η ) V



σ  (t ) Vsτ g (τ ) (t, η ) . σ (t ) Vs (t, η )

(69)

∂ to both sides of (69), Taking partial derivative ∂η

∂  ∂∂t Vs (t, η )  ∂  Vsτ g(τ ) (t , η )  = (q + i2π r )σ (t ) s (t , η ) ∂η Vs (t, η ) ∂η V  σ  (t ) ∂  Vτ g (τ ) (t, η )  −

σ (t ) ∂η

 τ g(τ )

s

s (t, η ) V

(t,η )

.



∂ s Therefore, if in addition, ∂η = 0, then (q + i2π r )σ (t ) = s (t,η ) V P0 (t, η ), where P0 (t, η) is defined by (28). Back to (69), we have V

∂ ∂ t Vs (t, η ) = p + qt −

s (t, η ) V

τ g(τ ) (t , η ) σ  (t ) V + i2π (c + rt ) + P0 (t , η ) s s (t , η ) σ (t ) V 

−

σ  (t ) Vsτ g (τ ) (t, η ) . σ (t ) Vs (t, η )

[1] L. Cohen, Time-Frequency Analysis, Prentice Hall, New Jersey, 1995. [2] P. Flandrin, Time-frequency/time-scale analysis, Wavelet Analysis and its Applications, vol. 10, Academic Press Inc., San Diego, CA, 1999. [3] L. Stankovic´ , M. Dakovic´ , T. Thayaparan, Time-Frequency Signal Analysis with Applications, Artech House, Boston, 2013. [4] F. Hlawatsch, G.F. Boudreaux-Bartels, Linear and quadratic TF signal representations, IEEE Signal Proc. Mag. 9 (2) (1992) 21–67. [5] S. Mallat, A Wavelet Tour of Signal Processing, Academic press, 1999. [6] S. Meignen, T. Oberlin, P. Depalle, P. Flandrin, S. McLaughlin, Adaptive multimode signal reconstruction from time frequency representations, Phil. Trans. R. Soc. A 374 (2065) (2016). [7] H. Choi, W. Williams, Improved TF representation of multicomponent signals using exponential kernels, IEEE Trans. Acoust. Speech, vol. ASSP 37 (6) (1989) 862–871. [8] L. Stankovic´ , A method for TF signal analysis, IEEE Trans. Signal Proc. 42 (1) (1994) 225–229. [9] S. Stankovic´ , I. Orovic, C. Ioana, Effects of cauchy integral formula discretization on the precision of IF estimation: unified approach to complex-lag distribution and its l-form, IEEE Signal Proc. Lett. 16 (4) (2009) 307–310. [10] H. Hassanpour, M. Mesbah, B. Boashash, SVD-Based TF feature extraction for newborn EEG seizure, EURASIP J. Adv. Signal Proc. 16 (2004) 2544–2554. [11] L. Stankovic´ , T. Thayaparan, M. Dakovic´ , Signal decomposition by using the s-method with application to the analysis of HF radar signals in sea-clutter, IEEE Trans. Signal Proc. 54 (11) (2006) 4332–4342. [12] L. Stankovic´ , D. Mandic´ , M. Dakovic´ , M. Brajovic´ , Time-frequency decomposition of multivariate multicomponent signals, Signal Proc. 142 (2018) 468–479. [13] N.E. Huang, Z. Shen, S.R. Long, M.L. Wu, H.H. Shih, Q. Zheng, N.C. Yen, C.C. Tung, H.H. Liu, The empirical mode decomposition and hilbert spectrum for nonlinear and nonstationary time series analysis, Proc. R. Soc. Lond. A 454 (1971) (1998) 903–995. [14] F. Auger, P. Flandrin, Improving the readability of TF and TF representations by the reassignment method, IEEE Trans. Signal Proc. 43 (5) (1995) 1068–1089. [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] P. Flandrin, G. Rilling, P. Goncalves, Empirical mode decomposition as a filter bank, IEEE Signal Proc. Lett. 11 (2004) 112–114. [17] N.E. Huang, Z. Wu, A review on Hilbert–Huang transform: method and its applications to geophysical studies, Rev. Geophys. 46 (2) (2008). [18] G. Rilling, P. Flandrin, One or two frequencies? the empirical mode decomposition answers, IEEE Trans. Signal Proc. 56 (2008) 85–95. [19] Z. Wu, N.E. Huang, Ensemble empirical mode decomposition: a noise-assisted data analysis method, Adv. Adapt. Data Anal. 1 (1) (2009) 1–41. [20] L. Li, H. Ji, Signal feature extraction based on improved EMD method, Measurement 42 (2009) 796–803. [21] L. Lin, Y. Wang, H.M. Zhou, Iterative filtering as an alternative algorithm for empirical mode decomposition, Adv. Adapt. Data Anal. 1 (4) (2009) 543–560. [22] J.D. Zheng, H.Y. Pan, T. Liu, Q.Y. Liu, Extreme-point weighted mode decomposition, Signal Proc. 42 (2018) 366–374. [23] R.R. Sharma, R.B. Pachori, Improved eigenvalue decomposition-based approach for reducing cross-terms in wignerville distribution, Circuits Syst. Signal Proc. 37 (8) (2018) 3330–3350. [24] T. Oberlin, S. Meignen, V. Perrier, An alternative formulation for the empirical mode decomposition, IEEE Trans. Signal Proc. 60 (5) (2012) 2236–2246. [25] I. Daubechies, J. Lu, H.T. Wu, Synchrosqueezed wavelet transforms: an empirical mode decomposition-like tool, Appl. Comput. Harmon. Anal. 30 (2) (2011) 243–261. [26] G. Thakur, H.T. Wu, J. Math, Synchrosqueezing based recovery of instantaneous frequency from nonuniform samples, SIAM. Anal. 43 (5) (2011) 2078–2095. [27] H.T. Wu, Adaptive Analysis of Complex Data Sets, Princeton Univ., Princeton, NJ, 2012 Ph.d. dissertation. [28] T. Oberlin, S. Meignen, V. Perrier, The fourier-based synchrosqueezing transform, in: Proc. 39th Int. Conf. Acoust., Speech, Signal Proc. (ICASSP), 2014, pp. 315–319. [29] G. Thakur, E. Brevdo, N. Fucˇ kar, H.T. Wu, The synchrosqueezing algorithm for time-varying spectral analysis: robustness properties and new paleoclimate applications, Signal Proc. 93 (5) (2013) 1079–1094.

L. Li, H. Cai and H. Han et al. / Signal Processing 166 (2020) 107231 [30] D. Iatsenko, P.V.E. McClintock, A. Stefanovska, Linear and synchrosqueezed TF representations revisited: overview, standards of use, resolution, reconstruction, concentration, and algorithms, Digital Signal Proc. 42 (2015) 1–26. [31] S. Meignen, D.H. Pham, S. McLaughlin, On demodulation, ridge detection and synchrosqueezing for multicomponent signals, IEEE Trans. Signal Proc. 65 (8) (2017) 2093–2103. [32] T. Oberlin, S. Meignen, V. Perrier, Second-order synchrosqueezing transform or invertible reassignment? Towards ideal TF representations, IEEE Trans. Signal Proc. 63 (5) (2015) 1335–1344. [33] T. Oberlin, S. Meignen, The 2nd-order wavelet synchrosqueezing transform, in: 2017 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 2017. New Orleans, LA, USA [34] D. Fourer, F. Auger, K. Czarnecki, S. Meignen, P. Flandrin, Chirp rate and instantaneous frequency estimation: application to recursive vertical synchrosqueezing, IEEE Signal Process. Lett. 24 (11) (2017) 1724–1728. [35] R. Behera, S. Meignen, T. Oberlin, Theoretical analysis of the 2nd-order synchrosqueezing transform, Appl. Comput. Harmon. Anal. 45 (2) (2018) 379–404. [36] D.H. Pham, S. Meignen, High-order synchrosqueezing transform for multicomponent signals analysis - with an application to gravitational-wave signal, IEEE Trans. Signal Proc. 65 (12) (2017) 3168–3178. [37] C. Li, M. Liang, A generalized synchrosqueezing transform for enhancing signal TF representation, Signal Proc. 92 (9) (2012) 2264–2274. [38] C.K. Chui, M.D. van der Walt, J. Geomath, Signal analysis via instantaneous frequency estimation of signal components, Int’l. 6 (1) (2015) 1–42. [39] H.Z. Yang, Synchrosqueezed wave packet transforms and diffeomorphism based spectral analysis for 1d general mode decompositions, Appl. Comput. Harmon. Anal. 39 (1) (2015) 33–66. [40] Z.L. Huang, J.Z. Zhang, T.H. Zhao, Y.B. Sun, Synchrosqueezing s-transform and its application in seismic spectral decomposition, IEEE Trans. Geosci. Remote Sensing 54 (2) (2016) 817–825. [41] C.K. Chui, Y.T. Lin, H.T. Wu, Real-time dynamics acquisition from irregular samples - with application to anesthesia evaluation, Anal. Appl. 14 (4) (2016) 537–590. [42] I. Daubechies, Y. Wang, H.T. Wu, ConceFT: concentration of frequency and time via a multitapered synchrosqueezed transform, Phil. Trans. R. Soc. A, 374 (2065) (2016). [43] S. Wang, X. Chen, G. Cai, B. Chen, X. Li, Z. He, Matching demodulation transform and synchrosqueezing in TF analysis, IEEE Trans. Signal Proc. 62 (1) (2014) 69–84. [44] Q.T. Jiang, B.W. Suter, Instantaneous frequency estimation based on synchrosqueezing wavelet transform, Signal Proc. 138 (2017) 167–181. [45] H.Z. Yang, L.X. Ying, Synchrosqueezed curvelet transform for two-dimensional mode decomposition, SIAM J. Math Anal. 46 (3) (2014) 2052–2083. [46] C.K. Chui, H.N. Mhaskar, Signal decomposition and analysis via extraction of frequencies, Appl. Comput. Harmon. Anal. 40 (1) (2016) 97–136. [47] L. Li, H.Y. Cai, Q.T. Jiang, H.B. Ji, An empirical signal separation algorithm based on linear TF analysis, Mech. Syst. Signal Proc. 121 (2019) 791–809. [48] H.Z. Yang, Statistical analysis of synchrosqueezed transforms, Appl. Comput. Harmon. Anal. 45 (3) (2018) 526–550.

15

[49] Z.C. Zhang, T. Yu, M.K. Luo, K. Deng, Estimating instantaneous frequency based on phase derivative and linear canonical transform with optimised computational speed, IET Signal Proc. 12 (5) (2018) 574–580. [50] C. Li, M. Liang, Time frequency signal analysis for gearbox fault diagnosis using a generalized synchrosqueezing transform, Mech. Syst. Signal Proc. 26 (2012) 205–217. [51] S.B. Wang, X.F. Chen, I.W. Selesnick, Y.J. Guo, C.W. Tong, X.W. Zhang, Matching synchrosqueezing transform: a useful tool for characterizing signals with fast varying instantaneous frequency and application to machine fault diagnosis, Mech. Syst. Signal Proc. 100 (2018) 242–288. [52] H.Z. Yang, J.F. Lu, L.X. Ying, Crystal image analysis using 2d synchrosqueezed transforms, Multiscale Model. Simul. 13 (4) (2015) 1542–1572. [53] J.F. Lu, H.Z. Yang, J. Imaging, Phase-space sketching for crystal image analysis based on synchrosqueezed transforms, SIAM Sci. 11 (3) (2018) 1954–1978. [54] K. He, Q. Li, Q. Yang, J. Testing, Characteristic analysis of welding crack acoustic emission signals using synchrosqueezed wavelet transform and evaluation 46 (6) (2018) 2679–2691. [55] H.T. Wu, Y.H. Chan, Y.T. Lin, Y.H. Yeh, Using synchrosqueezing transform to discover breathing dynamics from ECG signals, Appl. Comput. Harmon. Anal. 36 (2) (2014) 354–459. [56] H.T. Wu, R. Talmon, Y.L. Lo, Assess sleep stage by modern signal processing techniques, IEEE Trans. Biomed. Eng. 62 (4) (2015) 1159–1168. [57] C.L. Herry, M. Frasch, A.J. Seely, H.T. Wu, Heart beat classification from single-lead ECG using the synchrosqueezing transform, Physiol. Measur. 38 (2) (2017). [58] D.L. Jones, R.G. Baraniuk, A simple scheme for adapting TF representations, IEEE Trans. Signal Proc. 42 (12) (1994) 3530–3535. [59] N. Czerwinski, D.L. Jones, Adaptive short-time fourier analysis, IEEE Signal Proc. Lett. 4 (2) (1997) 42–45. [60] V. Katkovnik, L. Stankovic´ , Instantaneous frequency estimation using the wigner distribution with varying and data-driven window length, IEEE Trans. Signal Proc. 46 (9) (1998) 2315–2325. [61] J.G. Zhong, Y. Huang, Time-frequency representation based on an adaptive short-time fourier transform, IEEE Trans. Signal Proc. 58 (10) (2010) 5118– 5128. [62] L. Stankovic´ , A measure of some TF distributions concentration, Signal Proc. 81 (3) (2001) 621–631. [63] Y.L. Sheu, L.Y. Hsu, P.T. Chou, H.T. Wu, Entropy-based time-varying window width selection for nonlinear-type TF analysis, Int. J. Data Sci. Anal. 3 (2017) 231–245. [64] A. Berrian, N. Saito, Adaptive synchrosqueezing based on a quilted short-time Fourier transform, 2017, arXiv:1707.03138v5. [65] R. Baraniuk, P. Flandrin, A. Janssen, O. Michel, Measuring TF information content using the Rényi entropies, IEEE Trans. Inform. Theory 47 (4) (2001) 1391–1409. [66] V. Sharma, A. Parey, Performance evaluation of decomposition methods to diagnose leakage in a reciprocating compressor under limited speed variation, Mech. Syst. Signal Proc. (2018), doi:10.1016/j.ymssp.2018.07.029.