A Chirplet Transform-based Mode Retrieval Method for Multicomponent Signals with Crossover Instantaneous Frequencies

This paper studies blind source nonstationary signal separation in which a nonstationary signal is represented as a superposition of Fourier-like oscillatory modes:

$$x(t)=A_{0}(t)+\sum_{k=1}^{K}x_{k}(t),\quad x_{k}(t)=A_{k}(t)\cos\big{(}2\pi\phi_{k}(t)\big{)},$$

(1)

with $A_{k}(t),\phi_{k}^{\prime}(t)>0$ and $A_{k}(t)$ varying slowly. Such a representation, called an adaptive harmonic model (AHM) representation, is important for extracting information, such as the underlying dynamics, hidden in the nonstationary signal, with the trend $A_{0}(t)$, instantaneous amplitudes (IAs) $A_{k}(t)$ and the instantaneous frequencies (IFs) $\phi^{\prime}_{k}(t)$ being used to describe the underlying dynamics.

In nature, many real-world phenomena that can be formulated as signals (or in terms of time series) are often affected by a number of factors and appear as time-overlapping multicomponent signals in the form of (1). A natural approach to understand and process such phenomena is to decompose, or even better, to separate the multicomponent signals into their basic building blocks $x_{k}(t)$ (called components, modes or sub-signals) for extracting the necessary features. Also, for radar, communications, and other applications, signals often appear in multicomponent modes. Since these signals are mainly nonstationary, meaning the amplitudes and/or phases of some or all components change with the time, there have been few effective rigorous methods available for decomposition of them.

The empirical mode decomposition (EMD) algorithm along with the Hilbert spectrum analysis (HSA) is a popular method to decompose and analyze nonstationary signals [1]. EMD works like a filter bank [2, 3] to decompose a nonstationary signal into a superposition of intrinsic mode functions (IMFs) and a trend, and then the IF of each IMF is calculated by HSA. There are many articles studying the properties of EMD and variants of EMD have been proposed to improve the performance, see e.g. [2]-[14]. In particular, the separation ability of EMD was discussed in [4], which shows that EMD cannot decompose two components when their frequencies are close to each other. The ensemble EMD (EEMD) was proposed to suppress noise interferences [5]. The original EMD was extended to multivariate signals in [6, 8, 3]. Alternative sifting algorithms and formulations for EMD were introduced in [7, 9, 11]. Similar to the EMD filter bank, the wavelet filter bank for signal decomposition was proposed in [12], called empirical wavelet transform. An EMD-like sifting process was recently proposed in [13] to extract signal components in the linear time-frequency (TF) plane one by one. EMD is an efficient data-driven approach and no basis of functions is required. A weakness of EMD or EEMD is that it can easily lead to mode mixtures or artifacts, namely undesirable or false components [15]. In addition, there is no mathematical theorem to guarantee the recovery of the components.

Time-frequency analysis is another method to separate multicomponent signals, which is widely used in engineering fields such as communication, radar and sonar as a powerful tool for analyzing signals [16]. Time-frequency signal analysis and synthesis using the eigenvalue decomposition method have been studied [17, 18, 19]. Recently the synchrosqueezing transform (SST) was developed in [20] to provide mathematical theorems to guarantee the component recovery of nonstationary multicomponent signals. The SST, which was first introduced in 1996, intended for speech signal separation [21], is based on the continuous wavelet transform (CWT). The short-time Fourier transform (STFT)-based SST was also proposed in [22, 23] and further studied in [24]. SST provides an alternative to the EMD method and its variants, and it overcomes some limitations of the EMD and EEMD schemes [25, 26]. SST has been used in machine fault diagnosis [27, 28], crystal image analysis [29, 30], welding crack acoustic emission signal analysis [31], medical data analysis [32, 33, 34], etc.

SST works well for sinusoidal signals, but not for broadband time-varying frequency signals. To provide sharp representations for signals with significantly frequency changes, two methods were proposed. One is the matching demodulation transform-based SST (or called the instantaneous frequency-embedded SST) proposed in [35, 36] which changes broadband signals to narrow-band signals (see also [37]). The instantaneous frequency-embedded SST proposed in [36] preserves the IFs of the original signals. The other method is the 2nd-order SST introduced in [38] and [39]. The higher-order FSST was presented in [40] and [41], which aims to handle signals containing more general types. Very recently an adaptive SST with a time-varying parameter were introduced in [42, 43]. They obtained the well-separated condition for multicomponent signals using the linear frequency modulation to approximate a nonstationary signal at any local time. In addition, theoretical analysis of adaptive SST was obtained in [44, 45]. SST with a time-varying window width has also been studied in [46, 47].

Recently, a direct time-frequency approach, called signal separation operator (SSO), was introduced in [48] for multicomponent signal separation. In SSO approach, the components are reconstructed simply by substituting the time-frequency ridge to SSO. SSO based on CWT was proposed in [49, 50]. More recently, [51] proposed SSO of “linear chirp-based model” with the phase functions of the signal components in the source signal model (1) approximated by quadratic polynomials at any local time. This model provides a more accurate component recovery formulae, with analysis established in [52, 50].

While SST and SSO are mathematically rigorous on IF estimation and component recovery, both of them require that the components $x_{k}(t)$ are well-separated in the time-frequency plane, namely IFs of $x_{k}(t)$ satisfy

$$\phi^{\prime}_{k}(t)-\phi^{\prime}_{k-1}(t)\geq 2\triangle,\;2\leq k\leq K,$$

(2)

for some $\triangle>0$. In many applications, multicomponent signals are overlapping in the time-frequency plane, that is the IFs of its components are crossover. For example, in radar signal processing, the micro-Doppler

We say two components $x_{k-1}(t)$ and $x_{k}(t)$ in (1) are overlapping in time-frequency plane when the IF curves of them are crossing at $t=t_{0}$, i.e., $\phi_{k-1}^{{}^{\prime}}(t_{0})=\phi_{k}^{{}^{\prime}}(t_{0})$. In this paper we consider multicomponent signals of the form (1) satisfying

$$|\phi^{\prime}_{k}(t)-\phi^{\prime}_{\ell}(t)|+\rho|\phi^{\prime\prime}_{k}(t)-\phi^{\prime\prime}_{\ell}(t)|\geq 2\triangle,$$

(3)

where $\rho\geq 0$ is the regularization parameter and $\triangle>0$ is called the separation resolution. When $\rho=0$, (3) is reduced to the well-separated condition (2) required for SST and SSO. Condition (3) allows the IFs of some components $x_{k}(t)$ to cross as long as their chirp rates $\phi_{k}^{\prime\prime}(t)$ are different near the time $t_{0}$ where IFs crossing occurs.

Recently, many studies have been devoted to estimations of crossing instantaneous frequencies. An algorithm based on linear least square fitting technique was introduced [55]. In [56], the authors employed a three-dimensional parameter space (time, frequency, chirp rate) to estimate instantaneous frequencies and chirp rates jointly of signal, which can make the crossing of IFs appear as separated in the time-frequency-chirprate domain. Furthermore, [56] proposed a three-dimensional ridge detection algorithm. [57] used time-frequency analysis and Radon transform for the unsupervised separation of frequency modulated signals in noisy environment. The chirp rates were also introduced in [58], where estimations of IFs were turned into the resolution of a two-dimensional linear system. Effective quasi maximum likelihood–random samples consensus algorithm was also proposed for IF estimation of overlapping signals in the time-frequency (TF) plane [59].

Moreover, [54] developed a compressive sensing approach to recover stationary narrowband signals contaminated by strong nonstationary signals. An algorithm for estimations of IFs of signals with intersecting signatures in the time-frequency domain is proposed in [60]. Using the estimated IF, the corresponding component is removed from the mixture signal based on de-chirping method [60]. [61] combined parameterized de-chirping and band-pass filter to obtain components of multi-component signal, which avoided dealing with time-frequency representation of the signal and worked well under heavy noise. [62] introduced a novel nonparametric algorithm called ridge path regrouping to extract the IFs of the overlapped components and recover the corresponding component by using the intrinsic chirp component decomposition method. The time-frequency distribution (with a cross-term) [63] was used to estimate IFs. In [64], the authors developed a general model to characterize multi-component chirp signals, where IFs and instantaneous amplitudes of the intrinsic chirp components were modeled as Fourier series. Similarly, a multivariate intrinsic chirp mode was defined in [65]. Then IFs can be estimated using the framework of the general parameterized time-frequency transform and the corresponding multivariate intrinsic chirp modes were reconstructed by solving multivariate linear equations through an extended least square method [65]. Variational nonlinear chirp mode decomposition was employed to separate close or even crossed modes in [66]. In [67], the blind source was decomposed into a series of nonlinear chirp components and the reconstructed nonlinear chirp components can be clustered into corresponding intrinsic chirp sources. The frequency-domain intrinsic component decomposition method was proposed to characterize nonlinear and non-monotonic nonlinear group delays in frequency domain by using different kernel functions and separate and reconstruct each mode simultaneously as a time-frequency filter [68]. An adaptive chirp mode pursuit was proposed to capture signal modes one by one in a recursive framework [69]. In [70], the authors transformed a two-dimensional sinusoidal pattern into a single point in a two-dimensional plane and reconstructed signals from a reduced set of observations (back-projections).

In this paper, using the idea of SSO, we propose a chirplet transform-based signal separation scheme (CT3S) to retrieve modes of multicomponent signals with crossover IFs. Our algorithm addresses the inverse transform and aims to retrieve modes of multicomponent signals with fast varying IFs, even crossing IFs, while the aforementioned papers use the chirplet transform to obtain IFs. Most importantly, we provide a mathematically rigorous theorem which guarantees the recovery of components of multicomponent signals satisfying (3) with CT3S. To our best knowledge, there is no existing literature that establishes mathematical theorems to guarantee the recovery of the components overlapping in the time-frequency plane. Furthermore, we use a matched-filter along the specific time-frequency lines respect to the chirp rate to make different components be further separated and more concentrated in the three dimensional space of the chirplet transform. Based on the approximation of source signals with linear chirps at any local time, we propose an innovative signal reconstruction algorithm, called group filter-matched CT3S (GFCT3S), which takes a group of components into consideration simultaneously. GFCT3S is more suitable for signals with crossing IFs. When the IFs curves of different components are not crossover, but fast-varying and close to each other, our reconstruction algorithm will decrease the recovery errors significantly.

The remainder of this paper is organized as follows. In Section 2, we first review the signal separation methods by SST and SSO. After that we state the problem of recovering sub-signals with crossing IFs. In Section 3, we state our theorem which guarantees the recovery of components for multicomponent signal satisfying (3) with CT3S. We also introduce the filtered chirplet transform and propose GFCT3S for separating multicomponent signals with fast-varying and crossing IFs. We present the numerical experiments in Section 4. A conclusion is then presented in Section 5.

The (modified) STFT of signal $x(t)\in L_{2}({\mathbb{R}})$ with a window function $g(t)\in L_{2}({\mathbb{R}})$ is defined by,

$$V_{x}(t,\eta)=\int_{{\mathbb{R}}}x(\tau)g(\tau-t)e^{-i2\pi\eta(\tau-t)}d\tau,$$

(4)

If $g(0)\not=0$, then the original signal $x(t)$ can be recovered back from its STFT:

$$x(t)=\frac{1}{g(0)}\int_{\mathbb{R}}V_{x}(t,\eta)d\eta.$$

(5)

For multicomponent signal $x(t)$ in (1) satisfying the separation condition (2), the sub-signal $x_{k}(t)$ can be recovered by

$$x_{k}(t)\approx\frac{1}{g(0)}\int_{|\eta-\phi_{k}^{\prime}(t)|<\Gamma_{1}}V_{x}(t,\eta)d\eta,$$

(6)

for some $\Gamma_{1}>0$.

To enhance the time-frequency resolution and concentration, the idea of SST is to reassign the frequency variable. As in [22], denote

$$\omega_{x}(t,\eta)=\frac{\frac{\partial}{\partial t}V_{x}(t,\eta)}{i2\pi V_{x}(t,\eta)}.$$

(7)

The quantity $\omega_{x}(t,\eta)$ is called the “phase transformation” [20]. The STFT-based SST is to reassign the frequency variable $\eta$ by transforming the STFT $V_{x}(t,\eta)$ of $x(t)$ to a quantity, denoted by $R_{x}(t,\eta)$, on the time-frequency plane:

$$R_{x}(t,\eta)=\int_{\{\xi:V_{x}(t,\xi)\not=0\}}V_{x}(t,\xi)\delta\left(\omega_{x}(t,\xi)-\eta\right)d\xi.$$

(8)

Observe that signal reconstructions with the STFT and SST depend on the IF estimation of $\phi_{k}^{\prime}(t)$ and a given threshold $\Gamma_{1}$, hence it is indirect. In contrast, signal separation operator (SSO) [48] extracts signal components via local frequencies directly. The SSO $T_{a,\delta}$, which is applied to signals $x$ in (1), is defined by

$$\left(T_{a,\delta}x\right)(t,\eta)=\frac{1}{\hbar_{a}}\sum\limits_{n\in{\mathbb{Z}}}x(t-n\delta)h\left(\frac{n}{a}\right)e^{i2\pi n\eta},$$

(9)

where $h$ is a compactly supported window function, $\eta\in[0,1]$ and $\delta,a>0$ are parameters, with $a$ so chosen that

$$\hbar_{a}=\sum\limits_{n\in{\mathbb{Z}}}h\left(\frac{n}{a}\right)>0.$$

(10)

For a multicomponent signal $x(t)$ defined in (1), satisfying certain conditions including the well-separation condition (2), the set

$$\left\{\eta\in[0,1]:|\left(T_{a,\delta}x\right)(t,\eta)|>\Gamma_{2}\right\}$$

can be expressed as a disjoint union of exactly $K$ non-empty sets $\Theta_{\ell},\ell=1,2,...,K$, corresponding to the $K$ components of $x(t)$. The sub-signal $x_{k}(t)$ can be reconstructed by,

$$\widehat{x}_{k}(t)=2\Re e\left\{\left(T_{a,\delta}x\right)(t,\widehat{\eta}_{k})\right\},$$

(11)

where

$$\widehat{\eta}_{k}(t)=\arg\mathop{\max}\limits_{\eta\in\Theta_{k}}|\left(T_{a,\delta}x\right)(t,\eta)|.$$

(12)

As mentioned in Section 1, to recover components $x_{k}(t)$ with SST or SSO, it is required all IFs of different components be separated from each other, namely they be far away from each other and non-crossing as shown in (2). In particular, there is no mathematical theorem to guarantee the recovery of the waveforms of components when their IFs are crossover with only one observation $x(t)$ available. This paper is to provide a method to retrieve modes of such nonstationary multicomponent signals and establish retrieve error bounds. Next let us consider an example to show the performances of EMD, SST, SSO and our proposed method GFCT3S in retrieving the modes of a signal with crossover IFs.

Let $f(t)$ be the two-component signal introduced in [20]:

$$f(t)=f_{1}(t)+f_{2}(t)=\cos\left(t^{2}+t+\cos(t)\right)+\cos(8t).$$

(13)

Here we let the sampling rate $F_{s}=20$Hz and we only analyze the truncation signal on $t\in[0,256/F_{s}]$, with $256$ discrete sampling points. The instantaneous frequencies of $f_{1}(t)$ and $f_{2}(t)$ are $\phi_{1}^{\prime}(t)=(2t+1-\sin(t))/(2\pi)$ and $\phi_{2}^{\prime}(t)=4/\pi$, respectively.

Figure 2 shows the recovery results of $f_{1}(t)$ and $f_{2}(t)$. Observe that compared to the STFT and the STFT-based SST, the 2nd-order SST of $f(t)$ represents this two-component signal with crossing IF curves much sharper and clearer. However, the existing methods including EMD [1], SST [22], the 2nd-order SST [38] and SSO [48] are unable to recover the waveforms $f_{1}(t)$ and $f_{2}(t)$ accurately, see the recovered $f_{1}$ and $f_{2}$ by these methods in Figure 2. Our proposed GFCT3S in Algorithm 1 provided in Section 3.2 can recover the two components accurately as shown in the 4th panels (from the left) in Row 3 and Row 4 respectively in Figure 2. Note that EMD, SST, 2nd-order SST and SSO all result in big recovery errors for either $f_{1}$ or $f_{2}$ around $t_{0}=3.38$ where IFs crossing occurs, while our method produces very small errors near $t_{0}$. The boundary effect is unavoidable for all methods since we only use the truncation signal for $0\leq t\leq 12.8$ and the boundary extension has not be considered in this example. In addition, in this example we simply use the same Gaussian window with constant variance $\sigma=1.6$ for STFT, SST, the 2nd-order SST, SSO and GFCT3S .

In this section we propose a chirplet transform-based signal separation scheme (CT3S) to retrieve modes. We provide the main theorem on component recovery analysis. In addition, we introduce filtered chirplet transform (CT) to make IFs crossover components further separated and more concentrated in the three-dimensional space of CT. Finally, to improve the performance of CT3S in mode retrieval, we present an algorithm, called group filter-matched CT3S or GFTC3S for short, based on filtered CT and the approximation of source signals with linear chirps at any local time.