Instantaneous Frequency-Embedded Synchrosqueezing Transform for Signal Separation

The synchrosqueezing transform (SST) and its variants have been developed recently as an alternative to the empirical mode decomposition scheme to model a non-stationary signal as a superposition of amplitude- and frequency-modulated Fourier-like oscillatory modes. In particular, SST performs very well in estimating instantaneous frequencies (IFs) and separating the components of non-stationary multicomponent signals with slowly changing frequencies. However its performance is not desirable for signals having fast-changing frequencies. Two approaches have been proposed for this issue. One is to use the 2nd-order or high-order SST, and the other is to apply the instantaneous frequency-embedded SST (IFE-SST). For the SST or high order SST approach, one single phase transformation is applied to estimate the IFs of all components of a signal, which may yield not very accurate results in IF estimation and component recovery. IFE-SST uses an estimation of the IF of a targeted component to produce accurate IF estimation. The phase transformation of IFE-SST is associated with the targeted component. Hence the IFE-SST has certain advantages over SST in IF estimation and signal separation. In this article, we provide theoretical study on the instantaneous frequency-embedded short-time Fourier transform (IFE-STFT) and the associated SST, called IFE-FSST. We establish reconstructing properties of IFE-STFT with integrals involving the frequency variable only and provide reconstruction formula for individual components. We also consider the 2nd-order IFE-FSST.


INTRODUCTION
Recently the continuous wavelet transform-based synchrosqueezed transform (WSST) was developed in [1] as an empirical mode decomposition (EMD)-like tool to model a non-stationary signal x(t) as with A k (t), φ ′ k (t) > 0, where A k (t) is called the instantaneous amplitudes and φ ′ k (t) the instantaneous frequencies (IFs). The representation (1) of non-stationary signals is important to extract information hidden in x(t). WSST not only sharpens the time-frequency representation of a signal, but also recovers the components of a multicomponent signal. The synchrosqueezing transform (SST) provides an alternative to the EMD method introduced in [2] and its variants considered in many articles such as [3][4][5][6][7][8][9][10][11][12], and it overcomes some limitations of the EMD and ensemble EMD schemes such as mode-mixing. Many works on SST have been carried out since the publication of the seminal article [1]. For example, the short-time Fourier transform (STFT)-based SST (FSST) [13][14][15], the 2nd-order SST [16][17][18], the higher-order FSST [19,20], a hybrid EMD-WSST [21], the WSST with vanishing moment wavelets [22], the multitapered SST [23], the synchrosqueezed wave packet transform [24] and the synchrosqueezed curvelet transform [25] were proposed. Furthermore, the adaptive SST with a window function having a changing parameter was proposed in [26][27][28][29][30][31]. SST has been successfully used in machine fault diagnosis [32,33], and medical data analysis applications [see [34] and references therein]. [35] proposed a direct time-frequency method (called SSO) based on the ridges of spectrogram for signal separation. This method has been extended recently to the linear chirp-based models [36,37] and the models based on the CWT scaleogram [38,39]. A hybrid EMD-SSO computational scheme was developed in [40].
If the IFs φ ′ k (t) of the components x k (t) of a nonstationary multicomponent signal change slowly or change slowly compared with φ k (t), then SST performs very well in estimating φ ′ k (t) and separating the components x k (t) from x(t). However its performance is not desirable for signals having fastchanging frequencies. The 2nd-order and high-order SSTs were proposed for this issue and they do improve the accuracy of IF estimation and component recovery. The problem with the 2ndorder and high-order SSTs is that, like the convectional SST, one single phase transformation is applied to estimate the IFs of all components of a signal, which may not yield desirable results in IF estimation or component recovery.
Another approach is to demodulate the original signal to change a wide-band component into a narrow-band component. Li and Liang [41] and Meignen et al. [42] demodulate the original signal into a pure carrier signal and apply WSST and the 2nd-order FSST to the demodulated signal, respectively. FSST based on another demodulation was proposed in [ , we derive in this article mathematically rigorous phase transformation for IFE-FSST. In addition, in this article we also consider the 2nd-order IFE-FSST and derive the associate phase transformation.
The rest of this article is organized as follows. In Section 2, we briefly review FSST and the 2nd-order FSST. After that, we consider in Section 3 the IFE-STFT and establish reconstructing properties of IFE-STFT with integrals involving the frequency variable only. In Section 4, we derive mathematically rigorous phase transformations for IFE-FSST and the 2nd-order IFE-FSST. In addition, we provide reconstruction formula for individual components. Implementations and IFE-FSST-based component recovery algorithms are discussed in Section 5. Some experimental results are also provided in Section 5.

SHORT-TIME FOURIER TRANSFORM-BASED SST
where g(t) is a window function with g(0) = 0. x(t) can be reconstructed from its STFT: x(t) can also be recovered back from its STFT with an integral involving only the frequency variable η: In addition, one can show that if g(t) and x(t) are realvalued, then Furthermore, one can verify that STFT can be written as The STFT V x (t, η) of a slowly growing x(t) is well-defined and the above formulas still hold if the window function g(t) has certain smoothness and certain decaying order as t → ∞, for example g(t) is in the Schwarz class S. In this article, unless otherwise stated, we always assume that a window function g(t) has certain smoothness and decaying properties and g(0) = 0, and assume that a signal x(t) is a slowly growing function.

FSST
The idea of FSST is to re-assign the frequency variable η of V x (t, η). First we look at the STFT of x(t) = Ae i2π ξ 0 t , where ξ 0 is a positive constant. With we can obtain the IF ξ 0 of x(t) by where throughout this article, ∂ t denotes the partial derivative with respect to variable t. For a general x(t), at (t, η) for which V x (t, η) = 0, a good candidate for the IF of x(t) is In the following, denote which is called the "phase transformation" [1], "instantaneous frequency information" [13], or the "reference IF function" in [21]. FSST is to re-assign the frequency variable η by transforming the STFT V x (t, η) of x(t) to a quantity, denoted by R λ,γ x (t, ξ ), on the time-frequency plane defined by R λ,γ where ξ is the frequency variable, h(t) a compactly supported function with certain smoothness and ∞ −∞ h(t)dt = 1, γ > 0 is the threshold for zero and λ > 0 is a dilation. As λ, γ → 0, FSST is rewritten as For simplicity of presentation, throughout this article SSTs will be expressed as (7). Due to (4), we have that the input signal x(t) can be recovered from its FSST by If in addition, g(t) and x(t) are real-valued, then by (5), For a multicomponent signal x(t) given by (1), when A k (t), φ k (t) satisfy certain conditions, each component x k (t) can be recovered from its FSST: for certain Ŵ > 0, where IF k (t) is an estimate to φ ′ k (t). See [13][14][15] for the details.

Second-Order FSST
The 2nd-order FSST was introduced in [16]. The main idea is to define a new phase transformation ω 2nd x such that when x(t) is a linear frequency modulation (LFM) signal (also called a linear chirp), then ω 2nd x is exactly the IF of x(t). We say x(t) is a LFM signal or a linear chirp if with phase function φ(t) = ct + 1 2 rt 2 , IF φ ′ (t) = c + rt and chirp rate φ ′′ (t) = r. In [16], the reassignment operators are used to derive ω 2nd x . Different phase transformation ω 2nd x for the 2nd-order SST can be derived without using the reassignment operators see [28,29].
Let g be a given window function. Denote Recall that V x (t, η) denotes the STFT of x(t) with g defined by (2). In this article, we let V g 1 x (t, η) denote the STFT of x(t) with g 1 (t), namely, the integral on the right-hand side of (2) with g(t) replaced by g 1 (t). Define where q 0 (t, η) := 1 Then one can show that ω 2nd is an LFM signal given by (11), see [19,28]. Thus, we may define ω 2nd x (t, η) in (13) as the phase transformation for the 2ndorder FSST. Very recently a simple phase transformation for the 2nd-order FSST was proposed in [18].

INSTANTANEOUS FREQUENCY-EMBEDDED STFT
IFE-FSST is based on the IFE-STFT, which is defined below.
In the above definition, we assume x(t) ∈ L 2 (R). The definition of IFE-STFT can be extended to slowly growing functions x(t) if g(t) has certain smoothness and certain decaying order as t → ∞. Li and Liang [41] proposed the modulation x(τ ) → x(τ ) = x(τ )e −i2π (ϕ(τ )−η 0 τ ) and applied WSST to the modulated signal x(τ ), while [42] applied the 2nd-order FSST to x(t). The modulation: for IFE-FSST and also used in [44] for IFE-WSST is different from that used in [41,42]. IFE-STFT and IFE-CWT with such a modulation not only have more concentrated time-frequency representation than the conventional STFT and CWT respectively, but also well keep the IF of the signal. The reader is referred to [43] and [44] for detailed discussions.
[43] provides a reconstruction formula with IFE-STFT for the whole signal x(t), which is similar to (3) and involves an integral with both the time and frequency variables.
[43] does not consider individual component recovery formula with IFE-FSST. In this article, we provide such a component recovery formula.
To this regard, in this section we establish a reconstruction formula with IFE-STFT like (4), which involves an integral with the frequency variable only. First we have the following property about the IFE-STFT.
Proof: We have where the last equality follows from (6).
The next theorem shows that x(t) can be recovered from its IFE-STFT with an integral involving η only. Theorem 1. Let x(t) be a function in L 2 (R). Then Proof: Let x(t) be the function defined by (16). From (15), we have Thus, Equation (17) holds.
If one is interested in V I x (t, η) with the positive frequency η > 0 only, then we have the following result on how to recover . If y(η) = 0, η ≤ B for some constant B, then for any η 0 ≥ −B, Proof: Let x(t) be the function defined by (16). Then Frontiers in Applied Mathematics and Statistics | www.frontiersin.org Thus, Equation (18) holds.
Next theorem shows that when the condition y(η) = 0, η ≤ B in Theorem 2 does not hold, the integral in the right-hand side of (18) can still approximate x(t) well if η 0 is large.

IFE-STFT BASED SYNCHROSQUEEZING TRANSFORM
In this section, we consider IFE-FSST, the synchrosqueezing transform based on IFE-STFT. First we show how to derive the phase transformation associated with (the 1st-order) IFE-FSST. After that we introduce the 2nd-order IFE-FSST.

IFE-FSST
To define IFE-FSST, first we need to define the corresponding phase transformation ω I x (a, b). Let us consider the case On the other hand, where V I,g ′ x (t, η) denotes the IFE-STFT of x(t) defined by (14) with ϕ(t) and the window function g ′ given by (12). Thus, if V I x (t, η) = 0, then Based on the above discussion, for a general signal x(t), we define the phase transformation ω I x (a, b) of the IFE-FSST of x(t) to be The IFE-FSST of a signal x(t) with ϕ and ξ 0 is defined by where ω I x (t, η) is the phase transformation defined by (21).
The IFE-FSST is called the demodulation transform-based SST in [43]. The corresponding phase transformation in [43] is different from our ω I x (t, η) defined in (21). By (18) in Theorem 1, we know the input signal x(t) can be recovered from its IFE-FSST as shown in the following: For x(t) ∈ L 2 (R), and if, in addition, the conditions in Theorem 2 hold, then For a multicomponent signal x(t) in the form (1), if R I x k (t, ξ ), 1 ≤ k ≤ K lie in different time-frequency zones, then following (18), we know x k (t) can be recovered from its IFE-FSST: for certain Ŵ 1 > 0, where IF k (t) is an estimate of φ ′ k (t). If x k (t) and g(t) are real-valued, then

2nd-Order IFE-FSST
In this subsection, we propose the 2nd-order IFE-FSST. The key point is, based on IFE-STFT, to define a phase transformation ω I,2nd x (t, η) which is the IF φ ′ (t) of x(t) when x(t) is a linear chirp given by (11). As above, for g 1 (t) = tg(t), we use V I,g 1 x (t, η) to denote the IFE-STFT of x(t) with the window function g 1 (t), namely, the integral on the right-hand side of (14) with g(t) replaced by g 1 (t). Next we define the phase transformation ω I,2nd x (t, η) for the 2nd-order IFE-FSST to be: where ω I x (t, η) is defined by (21), and Theorem 4. If x(t) is a linear chirp signal given by (11), then at (t, η) where V I x (t, η) = 0, ∂ η V I,g 1 x (t, η)/V I x (t, η) = 0, ω I,2nd x (t, η) defined by (26) is the IF of x(t), namely ω I,2nd x (t, η) = c + rt.
Proof: Here, we provide the proof of ω I,2nd x (t, η) = c+rt for more general linear chirp signals given by where p, q are real numbers. For the simplicity of presentation, we denote and thus, V I x (t, η) can simply be written as Observe that for x(t) given by (28), we have Thus, On the other hand, as shown above, V I x ′ (t, η) is equal to the quantity in (20). Therefore, Hence, at (t, η) on which V I x (t, η) = 0, we have Taking partial derivative ∂ η to the both sides of (29), we have , which leads to q i2π + r = Q 0 (t, η), for (t, η) with ∂ η V I,g 1 x (t, η)/V I x (t, η) = 0, where Q 0 (t, η) is defined by (27).
Since c + rt is real, taking the real parts of the quantities in the above equation, we have , which is ω I,2nd x (t, η). This completes the proof of Theorem 4.
With the phase transformation ω I,2nd x (t, η) in (26), we have the corresponding 2nd-order IFE-FSST of a signal x(t) with ϕ, ξ 0 and window function g defined by One has reconstruction formulas with R I,2 x (t, ξ ) similar to (22)-(25).

Calculating ω I
x (t, η) and ω I,2nd x (t, η) First we consider the IFE-FSST. We need to calculate ω I x (t, η). We will use (15) so that FFT can be applied to (discrete signals) x and xϕ ′ to calculate V I (t, η), V I xϕ ′ (t, η) and V I,g ′ x (t, η). V I xϕ ′ (t, η) can be obtained by (15) with x replaced by xϕ ′ . As long as After obtaining V I (t, η), V I xϕ ′ (t, η) and V I,g ′ x (t, η), we get ω I x (t, η) and then the IFE-FSST.
Frontiers in Applied Mathematics and Statistics | www.frontiersin.org

IFE-FSST Algorithms for IF Estimation and Component Recovery and Experiments
To apply IFE-STFT or IFE-FSST, first of all we need to choose ϕ(t) and ϕ ′ (t). For the purpose of estimating the IF φ ′ k (t) of the kth component x k (t) and/or recover x k (t) of a multicomponent signal x(t), we should choose ϕ(t) and ϕ ′ (t) close to φ k (t) (up to a constant) and φ ′ k (t) respectively. One way is to use the ridges of the STFT. More precisely, suppose {t n } 0≤n<N , {η j } 0≤j<J , {ξ m } 0≤m<M are the sampling points of t, η, ξ respectively for STFT V x (t, η), FSST R x (t, ξ ), and IFE-FSST R I x (t, ξ ). Let η j n ,k , 0 ≤ n < N be the STFT ridge corresponding to x k (t) given by for each n, 0 ≤ n < N, where for each n, G t n ,k is an interval containing φ ′ k (t n ) (with convention: φ 0 (t) ≡ 0) at the time instant t n , and G t n ,k , 0 ≤ k ≤ K form a disjoint union of η : |V x (t n , η)| > γ , namely for each t n , See more details on G t,k in [37]. x (t, η)| (right).
To recover a component by either FSST or IFE-FSST, we need an estimate IF k (t) for φ ′ k (t) so that (10) or (24)/(25) can be applied. One way is to use the ridges of FSST and IFE-FSST to approximate φ ′ k (t n ). More precisely, let ξ m n ,k , 0 ≤ n < N be the FSST ridge defined similarly to the STFT ridge in (34): Then Equation (10) becomes Similarly, Equation (24) implies that x k (t) can be recovery from (discrete) IFE-FSST: where m I n , 0 ≤ n < N are the indices for IFE-FSST ridge defined as (36) with R x (t n , ξ m ) replaced by R I x (t n , ξ m ): For real-valued x k (t) and g(t), the recovery formulas (37) and (38) are respectively Re and To summarize, we have the following algorithm to estimate IF φ ′ k (t) and recover x k (t) by IFE-FSST.  (1). To estimate φ ′ k (t) and recover x k (t) by IFE-FSST, do the following. Step 1. Obtain the STFT ridge η j n ,k , 0 ≤ n < N by (34) and ϕ ′ (t n ), ϕ(t n ), 0 ≤ n < N by (35).
Step 2. Calculate IFE-FSST with ϕ ′ , ϕ obtained in Step 1. The ridge ξ m I n ,k , 0 ≤ n < N defined by (39) is an estimate of φ ′ k (t) and x I,rec k (t) in (38) is an approximation to x k (t).
We can use Algorithm 1 to recover each component x k (t) one by one. We can also apply Algorithm 1 to the remainder x(t) − x I,rec k (t) to recover another component after x k (t) is recovered; and we can repeat this procedure. The procedure of this iterative method is described as follows. Step 1. Apply Algorithm 1 to obtain x I,rec 1 (t).
Step 2. Let y 1 = x − x I,rec 1 . Apply Algorithm 1 to y 1 to obtain x I,rec 2 (t).
Step 3. Let y 2 = x − x I,rec 1 − x I,rec 2 . Apply Algorithm 1 to y 2 to obtain x I,rec 3 (t). Repeat this process to obtain x I,rec 4 (t), · · · , and finally x I,rec K (t).  Step 4. Apply Algorithm 1 to x − K k=2 x I,rec k to recover x 1 (t). Let x I,rec 1 (t) be the recovered x 1 (t). Then Apply Algorithm 1 to x − x I,rec 1 (t) − K k=3 x I,rec k to recover x 2 (t). Let x I,rec 2 (t) be the recovered x 2 (t). Obtain x I,rec 3 (t) by applying Algorithm 1 to x− x I,rec 1 (t)− x I,rec 2 (t)− K k=4 x I,rec k . Repeat this process to obtain x I,rec 4 (t), · · · , and finally x I,rec K (t).
We can repeat the procedure in Step 4 of Algorithm 2. That is why we call Algorithm 2 an iterative algorithm.
Next we consider two examples. We let be the window function, where σ > 0. First we consider a mono-component signal where x(t) is uniformly sampled with sample points t n = n△t, 0 ≤ n < N = 128, △t = 1 128 . The IF of x(t) is φ ′ (t) = 16 + 32t and it is shown in the 1st row of Figure 1. The FSST and IFE-FSST of x(t) are provided in the 2nd row; and the 2nd-order FSST and IFE-FSST are shown in the 3rd row. In this example we let σ = 1 16 . As mentioned above, discrete ϕ ′ (t) and ϕ(t) defined by (35) are used to define IFE-STFT and the 2nd-order IFE-STFT. Obviously IFE-FSST provides a much sharper time-frequency representation of x(t) than FSST. Both the 2nd-order FSST and the 2nd-order IFE-FSST as well give sharp time-frequency representations of x(t).
For a mono-component signal x(t) as given by (43), since x(t) can be recovered from FSST or IFE-FSST as shown in (8) and (22) respectively, theoretically, either (40) or (41) gives high accurate approximation to x(t) as long as M 0 is large enough. We choose a small M 0 so that the recovery errors with it show how sharp the time-frequency representations with FSST and IFE-FSST are. Here and below we set M 0 = 8.
In Figure 2, we provide the recovery errors x rec (t n ) − x(t n ), x I,rec (t n ) − x(t n ) for x(t) by FSST and IFE-FSST, where x rec (t n ) and x I,rec (t n ) are given by (40) and (41) respectively with M 0 = 8. Here, we show the error on [0.125, 0.875) only to ignore the boundary effect. Obviously, IFE-FSST provides a much sharper time-frequency representation than FSST.
Next we consider a two-component signal given by x(t) = x 1 (t) + x 2 (t), x 1 (t) = cos 2π 32t + 10 π cos(2πt) , x 2 (t) = cos 2π 64t + 10 π cos(2πt) , where t ∈ [0, 2), and x(t) is uniformly sampled with sample points t n = n△t, 0 ≤ n < N = 512, △t = 1 256 . Thus, IFs of x 1 (t), x 2 (t) are φ ′ 1 (t) = 32−20 sin(2πt), φ ′ 2 (t) = 64−20 sin(2πt), which are shown on the top-left panel of Figure 3. In this example we let σ = 1 32 for the window function. To this two-component signal, we apply Algorithm 2 to obtain x I,rec 1 (t) and x I,rec 2 (t). In the 3rd row of Figure 3 we show the IFE-FSSTs of x I,rec 1 (t) and x I,rec 2 (t). The FSST of x(t) is provided in the top-right panel of Figure 3. Of course, we can also apply iterative method to FSST to recover components one by one. Namely, we apply FSST to obtain x rec 1 (t), then apply FSST to x(t) − x rec 1 (t) to obtain x rec 2 (t). After that we apply FSST to x(t) − x rec 2 (t) to obtain x rec 1 (t), and finally to obtain x rec 2 (t) by applying FSST to x(t) − x rec 1 (t). The 2nd row of Figure 3 shows the FSSTs of x rec 1 (t) and x rec 2 (t). Comparing the FSST of x in the top-right panel with the individual FSSTs in the 2nd row, we see there is not much improvement of the time-frequency representation of FSST of x after we apply the iterative component recovery procedure.
In the 4th row of Figure 3, we provide the recovery errors for x 1 (t), x 2 (t) by FSST and IFE-FSST. Here, we show the error on [0.125, 1.875). From Figure 3, we see IFE-FSST provides a much sharper time-frequency representation for x(t). We also consider FSST and IFE-FSST of two-component x(t) in the noisy environment and our experiments show that IFE-FSST provides a sharp time-frequency representation in the noisy environment. In addition, we consider the 2nd-order IFE-FSST for component recovery. It does not provide much improvement than IFE-FSST. This may be due to that the results from IFE-FSST are hard to improve. Due to that only 15 pictures are allowed to be included in a article in this journal, we do not present these results here.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author/s.

ETHICS STATEMENT
This article was approved by AFRL for public release on 03 Dec. 2021, Case Number: AFRL-2021-4285, Distribution unlimited.