# Intrinsic Mode Chirp Multicomponent Decomposition with Kernel Sparse Learning for Overlapped Nonstationary Signals Involving Big Data

Authors: Haixin Sun; Yongchun Miao; Jie Qi
Journal: Complexity
 (2018)

---

## Abstract

We focus on the decomposition problem for nonstationary multicomponent signals involving Big Data. We propose the kernel sparse learning (KSL), developed for the T-F reassignment algorithm by the path penalty function, to decompose the instantaneous frequencies (IFs) ridges of the overlapped multicomponent from a time-frequency representation (TFR). The main objective of KSL is to minimize the error of the prediction process while minimizing the amount of training samples used and thus to cut the costs interrelated with the training sample collection. The IFs first extraction is decided using the framework of the intrinsic mode polynomial chirp transform (IMPCT), which obtains a brief local orthogonal TFR of signals. Then, the IFs curves of the multicomponent signal can be easily reconstructed by the T-F reassignment. After the IFs are extracted, component decomposition is performed through KSL. Finally, the performance of the method is compared when applied to several simulated micro-Doppler signals, which shows its effectiveness in various applications.


---

## Body

## 1. Introduction

In practical applications such as oceanic investigation [1], radar [2], biomedical application [3], and mechanical fault diagnosis [4], we need to represent and process nonstationary signals. While Big Data can be explicitly regarded as a good fortune, great challenges also appear with extensive datasets. In general, Big Data introduce challenges in which properties, such as time domain, space-time domain, computational complexity, and information energy, are mingled in complicated ways with data resources, as shown in Figure 1. Particularly, a multicomponent signal (MCS) is a superposition of monocomponent signals whose may overlap (or cross) in time-frequency (T-F) domains. The fundamental problem in processing signals is to extract the useful information stored in the amplitude-modulated (AM) and/or the frequency-modulated (FM) signals. The time-frequency analysis (TFA) is a central tool to capture the time-varying features of multicomponent chirp signals. To analyze further characteristics of these signals, we are more eager to decompose them into monocomponent signals.Figure 1
Signal processing challenges and applications for Big Data.


### 1.1. Related Work

In the past few decades, many signal decomposition approaches have been developed to capture the accurate monocomponent signals. Commonly, the typical methods can be summarized as time-frequency distributions [5], the empirical mode decomposition [6] and its multivariate extensions [7], the local mean decomposition [8], the reassignment method [9, 10] such as the synchrosqueezing transform (SST) [9], the sparsification methods [11, 12], and so on.However, most of the decomposition methods are signal-dependent and presume that the signal components meet strict divisible conditions in the T-F domain [9], and thus they are affected by additive noise and close or overlapped instantaneous frequencies (IFs) of signal components. Recently, the intrinsic chirp component decomposition (ICCD) method proposed in [13] and used in [14] revealed the incapable of separating crossed or overlapped chirp signals by using a joint-estimation scheme. Essentially, the ICCD based on the short-time Fourier transform (STFT) can accurately reconstruct overlapped components. Unfortunately, according to the Heisenberg uncertainty principle, the STFT, as the liner TFA, cannot achieve an arbitrarily high T-F resolution. The overlapped components are very close in the T-F plane and the mask with low T-F resolution may cover all components and the unwanted cross terms, which limits the applications of the ICCD.To analyze overlapped signal components, TFA is an important issue which deserves to be mentioned for the extraction of feature information from the components. Currently, the extraction accuracy relies on the energy concentration of TFR generated by TFA methods. There exist a variety of types of TFA methods, and they can be normally divided into two categories: the parametric TFA (PTFA) methods and the nonparametric TFA (NPTFA) methods. PTFA methods, such as polynomial [1, 15], spline-kernelled chirplet transform (SCT) [16], and sinusoidal models [17], often involve the high-dimensional search of the IFs, which is very time consuming. Moreover, the predesigned parametric models may be only suitable for special applications. In this paper, we are more concerned about the NPTFA methods which are known for their more adaptability in practical applications. The traditional NPTFA methods by short-time Fourier transform (STFT), continuous wavelet transform (CWT), and Wigner-Ville distribution (WVD) are subjected to poor energy concentration. Recently, a class of TFA techniques, referred to as T-F reassignment method, proposes to improve the readability of TFR clearly. The synchrosqueezing transform (SST), belonging to this class, is a postprocessing technique based on the CWT to obtain better localized TFR of nonstationary signals [18, 19]. The SST provides a higher resolution that moderates the limitations of linear projection based traditional NPTFA methods, such as STFT, CWT, and SWT [19, 20]. Moreover, the SST reconstructs the energies of these transforms, so that the resulting energies of coefficients are condensed around the IF curves of the components. However, the interference of serious noise on TFA has attracted a major attention, since it brings about a dominant error source when it appears in the contaminated signal. Errors on account of the serious noise induces false maxima (maxima other than the autoterms) in the T-F plane. In [21], an ant colony optimization (ACO) algorithm and TFR by WVD is adopted to estimate IFs of the components under serious noise. This method is applicable to represent overlapped multicomponent signals. Nevertheless, we find that these algorithms mentioned above cannot further decompose the nonstationary overlapped components, which may be often seen in various signal processing of practical applications.


### 1.2. Proposed Method

Notice that despite all those efforts decomposing overlapped MCSs is still a challenging task. In this paper, we develop the intrinsic mode polynomial chirp transform (IMPCT) which takes into account IFs characteristics and a novel KSL method is proposed in order to obtain a robust decomposition overlapped components in high noise. For a nonstationary signal, instantaneous amplitudes (IAs) and IFs are two of the foremost properties. In [22], the IAs of the real signal are approximated by a third-order or more orders polynomial. To represent complex chirp signals, the IFs and IAs are both modeled as discrete chirp series in the IMPCT, which enables energy concentration and refines the most important coefficients in high noise. The coefficients for this model are distributed according to the time-varying band-pass filter in the T-F plane and not projected from either side. In addition, the distributed properties of the IMPCT method and carefully designed series make the representation more robust to swift IFs variations.For the practical implementation of a decomposition method, a very interesting solution is the kernel sparse learning (KSL) approach, which has not yet been developed for the multicomponent decomposition. While we propose the KSL, as an active learning, similar to kernel dictionary learning, which has been successfully applied to classification problems [23–25], this model works differently. In fact, the purpose of the algorithm is to maximize the accuracy to minimize the amount of training samples. For such aim, smart strategies are adopted to pick out the most prominent training samples. An advantage of our approach is that it separates closer crossing components by KSL in parallel.The remainder of this paper is organized as follows. In Section2, the theory of IMPCT is described in detail. The IFs estimation by T-F reassignment by the path penalty function and the decomposition model by kernel sparse learning are given in Section 3. Section 4 shows simulation results of decomposing some real and micro-Doppler signals are provided, and conclusions are followed in Section 5.


## 1.1. Related Work

In the past few decades, many signal decomposition approaches have been developed to capture the accurate monocomponent signals. Commonly, the typical methods can be summarized as time-frequency distributions [5], the empirical mode decomposition [6] and its multivariate extensions [7], the local mean decomposition [8], the reassignment method [9, 10] such as the synchrosqueezing transform (SST) [9], the sparsification methods [11, 12], and so on.However, most of the decomposition methods are signal-dependent and presume that the signal components meet strict divisible conditions in the T-F domain [9], and thus they are affected by additive noise and close or overlapped instantaneous frequencies (IFs) of signal components. Recently, the intrinsic chirp component decomposition (ICCD) method proposed in [13] and used in [14] revealed the incapable of separating crossed or overlapped chirp signals by using a joint-estimation scheme. Essentially, the ICCD based on the short-time Fourier transform (STFT) can accurately reconstruct overlapped components. Unfortunately, according to the Heisenberg uncertainty principle, the STFT, as the liner TFA, cannot achieve an arbitrarily high T-F resolution. The overlapped components are very close in the T-F plane and the mask with low T-F resolution may cover all components and the unwanted cross terms, which limits the applications of the ICCD.To analyze overlapped signal components, TFA is an important issue which deserves to be mentioned for the extraction of feature information from the components. Currently, the extraction accuracy relies on the energy concentration of TFR generated by TFA methods. There exist a variety of types of TFA methods, and they can be normally divided into two categories: the parametric TFA (PTFA) methods and the nonparametric TFA (NPTFA) methods. PTFA methods, such as polynomial [1, 15], spline-kernelled chirplet transform (SCT) [16], and sinusoidal models [17], often involve the high-dimensional search of the IFs, which is very time consuming. Moreover, the predesigned parametric models may be only suitable for special applications. In this paper, we are more concerned about the NPTFA methods which are known for their more adaptability in practical applications. The traditional NPTFA methods by short-time Fourier transform (STFT), continuous wavelet transform (CWT), and Wigner-Ville distribution (WVD) are subjected to poor energy concentration. Recently, a class of TFA techniques, referred to as T-F reassignment method, proposes to improve the readability of TFR clearly. The synchrosqueezing transform (SST), belonging to this class, is a postprocessing technique based on the CWT to obtain better localized TFR of nonstationary signals [18, 19]. The SST provides a higher resolution that moderates the limitations of linear projection based traditional NPTFA methods, such as STFT, CWT, and SWT [19, 20]. Moreover, the SST reconstructs the energies of these transforms, so that the resulting energies of coefficients are condensed around the IF curves of the components. However, the interference of serious noise on TFA has attracted a major attention, since it brings about a dominant error source when it appears in the contaminated signal. Errors on account of the serious noise induces false maxima (maxima other than the autoterms) in the T-F plane. In [21], an ant colony optimization (ACO) algorithm and TFR by WVD is adopted to estimate IFs of the components under serious noise. This method is applicable to represent overlapped multicomponent signals. Nevertheless, we find that these algorithms mentioned above cannot further decompose the nonstationary overlapped components, which may be often seen in various signal processing of practical applications.


## 1.2. Proposed Method

Notice that despite all those efforts decomposing overlapped MCSs is still a challenging task. In this paper, we develop the intrinsic mode polynomial chirp transform (IMPCT) which takes into account IFs characteristics and a novel KSL method is proposed in order to obtain a robust decomposition overlapped components in high noise. For a nonstationary signal, instantaneous amplitudes (IAs) and IFs are two of the foremost properties. In [22], the IAs of the real signal are approximated by a third-order or more orders polynomial. To represent complex chirp signals, the IFs and IAs are both modeled as discrete chirp series in the IMPCT, which enables energy concentration and refines the most important coefficients in high noise. The coefficients for this model are distributed according to the time-varying band-pass filter in the T-F plane and not projected from either side. In addition, the distributed properties of the IMPCT method and carefully designed series make the representation more robust to swift IFs variations.For the practical implementation of a decomposition method, a very interesting solution is the kernel sparse learning (KSL) approach, which has not yet been developed for the multicomponent decomposition. While we propose the KSL, as an active learning, similar to kernel dictionary learning, which has been successfully applied to classification problems [23–25], this model works differently. In fact, the purpose of the algorithm is to maximize the accuracy to minimize the amount of training samples. For such aim, smart strategies are adopted to pick out the most prominent training samples. An advantage of our approach is that it separates closer crossing components by KSL in parallel.The remainder of this paper is organized as follows. In Section2, the theory of IMPCT is described in detail. The IFs estimation by T-F reassignment by the path penalty function and the decomposition model by kernel sparse learning are given in Section 3. Section 4 shows simulation results of decomposing some real and micro-Doppler signals are provided, and conclusions are followed in Section 5.


## 2. IMPCT for TFR

### 2.1. The Motivation of IMPCT

First of all, a numerical example of a noisy environment is applied to demonstrate the motivation of our method proposed. A simulated signal is considered as follows:(1)st=s1t+s2twith (2)s1t=1+0.5cos⁡2πt×cos⁡2π0.2+532t-474t2+369t3,s2t=1+0.5cos⁡2πt×cos⁡2π0.8+50t+550t2-300t3,which is sampled at a sampling frequency of 2000 Hz and the signal is polluted by a noise with SNR=-2dB. The T-F representations shown in Figures 2(a) and 2(b) are generated by SCT [16] and SWT [19], respectively.Figure 2
TFs representations by (a) SWT, (b) SCT, (c) Iteration=0 IMPCT, and (d) Iteration=50 IMPCT.
(a)
SWT (b)
SCT (c)
Iteration=0 IMPCT (d)
Iteration=50 IMPCTThe SWT, the continuous wavelet transform- (CWT-) based SST, sharpens the time-frequency representation of a multicomponent signal but suffers from a lower T-F resolution, as shown in Figure2(a). The T-F coefficients of the IFs trajectories are obtained by SWT only in the frequency domains. Thus, it is no wonder that the SWT cannot accurately estimate time-varying IFs of an overlapped signal as well. In Figure 2(b), the SCT is able to generate a TFR with a poor energy concentration for the contaminated signal with nonlinearly time-varying IFs. If the chirp rate of SCT deviates from the true IFs, the T-F representation will smear heavily.An intuitive solution stems from a combination of high-energy concentrations that partially result from the different modulation frequencies. The new T-F representations of the IMPCT will have better performance than the SCT, as shown in Figures\ 2(c) and 2(d) with 50 iterations.Therefore, it is necessary to develop alternative approaches to deal with the challenging cases that there exist components overlapping in the T-F domain. We provide the following definition.Definition 1.
The two componentsxk(t) and xm(t) in model (3) are said to cross (or overlap) in the T-F domain. If there exist a time instant t^, a bound d≥0 and a step h shift up and down such that their IFs satisfy(3)IFkt^=IFmt^,IFkτ^±h=IFmτ^∓h,IFk′t^-IFm′t^≥d.Equation (3) illustrates a strict definition of these overlapped components for our algorithm to identify the case that some components are tangent to each other (i.e.,IFk′(t^)=IFm′(t^)). Figure 3 shows the strict overlapping and tangent of the two components. After shifting up and down, if there is a point such that IFk(τ^)±h=IFm(τ^)±h,t^≠τ^ at other intersecting time τ^, the components are strictly overlapping. It is worth noticing that the IFs of the arbitrary two components have different curvature at the intersecting time t^.Figure 3
(a) Two components are strictly overlapping. (b) Two components are tangent. Blue line, IF trajectory. Red line, IF trajectory after shift up. Purple line, IF trajectory after shift down.
(a) (b)A signalst=rmejφmt=Atcos⁡ft can be considered an intrinsic mode chirp function, with accuracy ε and A(t) and f(t) meeting the following conditions [14]:(4)A∈C1R∩L∞R,f∈C2Rinft∈RAt>0,inft∈Rft<∞supt∈RAt<∞,supt∈Rft<∞,supt∈Rf′t<∞A′t,f′t≤εft,∀t∈R.Here, to overcome these limitations, a robust technology architecture for overlapped nonstationary signals is proposed, as shown in Figure4. Initially, the input signal is passed by the N-point IMPCT stage. The stage brings about the IMPCT coefficients shown as r1,r2,⋯,rm. The inverse IMPCT operation for a signal sequence of length λ1,λ2,⋯,λm is performed. These characterized polynomial coefficients are concatenated and rearranged to the IFs estimation for each of the input nonstationary signals. Then, to decompose overlapped components in practical applications, we need to deal with a very significant but challenging issue of how to estimate IFs of overlapped components for the IMPCT. The second problem is how to determine the criterion which can separate the signals under micro-Doppler effect and noise interference. In the following sections, the solution for these two problems will be provided about detailed theoretical analysis.Figure 4
Block diagram of the proposed approach.


### 2.2. IMPCT Representation

The Fourier series analysis brings to light the sinusoidal property, so it may be utilized for the reconstruction of some signals with linear phase [14]. However, when processing nonlinear polynomial phase signals (NPPS) often involved in the practical application, a signal representations can be obtained by using the IMPCT. Namely, the IMPCT can be adopted to demodulating components of interest, reducing them to the sparse sinusoidal components. In this section, the Fourier series is stretched into the complex polynomial Fourier signal model as the IMPCT. Consider a multicomponent polynomial phase signal contaminated by a white Gaussian noise as(5)xt=∑m=1Mrmejφmt+εt,where rm is the amplitude coefficient function and φm(t)=α1mt+⋯+αnmtn/n! indicates the phase of the m-th component; ε(t) shows the noise. In general, the IMPCT of x(t) can be formulated as follows: (6)Xf1,⋯,fn=∑m=1Mrm×∫-∞∞ejα1m-f1t+⋯+αnm-fn/n!tndt.For all points in the n-dimensional frequency domain that correspond to the positions of signal components characterized by f1=α1m,...,fn=αnm,m=1,...,M, we can acquire (7)Xf1,⋯,fn=∑m=1M2πrmδf1-α1m,⋯,fn-αnm.That is to say, X(f1,⋯,fn)→∞ at the positions of signal components: f1=α1m,⋯,fn=αnm.Thus, X(f1,⋯,fn) can be deemed to a sparse M-component representation. In addition, for f1≠α1m,⋯,fn≠αnm,m=1,⋯,M has limited value, negligible when compared to (7).Let us inspect the signal vectorx with elements X(k), superimposed of a sum of M NPPS: (8)Xk=∑m=1MXmk+εk=∑m=1Mrmej2π/Kka1m+k2a2m/2+⋯+knanm/n!,where the polynomial coefficients are supposed to be bounded integers and K is the total number of discrete samples. To offer a sparse representation of nonlinear signals, the discrete form of the IMPCT is adopted, which may be expressed as (9)Xλ1,⋯,λn=∑k=0K-1∑m=1Mrmej2π/Kka1m+k2a2m/2+⋯+knanm/n!×e-j2π/Kk2λ2/2+⋯+knλn/n!e-j2π/Kkλ1.When the set of IMPCT parameters (λ1,⋯,λn) is selected to obtain the adaptive coefficients of the NPPS (a1i,a2i,⋯,ani), thus we can obtain (10)λ1,⋯,λn=a1m,a2m,⋯,anmthe m-th signal component is demodulated and the sinusoid rmej2π/K(ma1m) is paramount in the IMPCT spectrum. Under this circumstances, the energy of the spectrum is extremely concentrated at λ=a1m. On the contrary, while (λ1,⋯,λn)≠(a1m,a2m,⋯,anm) for all m∈1,⋯,M the spectrum consisted of some dispersed phase. Under high noise environment, we suppose that the minimum spectral component having amplitude coefficient Krmin can be also primary in the spectrum, meaning that it is above the aggregate value of all other components (in the worst case above their accumulation): K|rmin|>∑m=1M‍rm, where rmin is the minimum of all rm, for m=1,⋯,M. Therefore, a search is performed for M sets of parameters λ1m,⋯,λnm∣m=1,⋯,M so that (λ1,⋯,λn)=(a1m,a2m,⋯,anm) holds, and χλ1m,⋯,λnm is strongly concentrated. All in all, multicomponent set of parameters are calculated to obtain the optimal concentrated vector χ(λ1,⋯,λn) for all m=1,⋯,M. The maximized component will lie in the frequency λ1=a1m.The IMPCT representation explains the motivation for obtaining an initial IFs estimate in the algorithm. The IFs estimate will be based on the value ofχ(λ1,⋯,λn). In the next section, the T-F reassignment algorithm by the path penalty function is proposed in order to more accurately estimate IFs with the serious noise cases.


## 2.1. The Motivation of IMPCT

First of all, a numerical example of a noisy environment is applied to demonstrate the motivation of our method proposed. A simulated signal is considered as follows:(1)st=s1t+s2twith (2)s1t=1+0.5cos⁡2πt×cos⁡2π0.2+532t-474t2+369t3,s2t=1+0.5cos⁡2πt×cos⁡2π0.8+50t+550t2-300t3,which is sampled at a sampling frequency of 2000 Hz and the signal is polluted by a noise with SNR=-2dB. The T-F representations shown in Figures 2(a) and 2(b) are generated by SCT [16] and SWT [19], respectively.Figure 2
TFs representations by (a) SWT, (b) SCT, (c) Iteration=0 IMPCT, and (d) Iteration=50 IMPCT.
(a)
SWT (b)
SCT (c)
Iteration=0 IMPCT (d)
Iteration=50 IMPCTThe SWT, the continuous wavelet transform- (CWT-) based SST, sharpens the time-frequency representation of a multicomponent signal but suffers from a lower T-F resolution, as shown in Figure2(a). The T-F coefficients of the IFs trajectories are obtained by SWT only in the frequency domains. Thus, it is no wonder that the SWT cannot accurately estimate time-varying IFs of an overlapped signal as well. In Figure 2(b), the SCT is able to generate a TFR with a poor energy concentration for the contaminated signal with nonlinearly time-varying IFs. If the chirp rate of SCT deviates from the true IFs, the T-F representation will smear heavily.An intuitive solution stems from a combination of high-energy concentrations that partially result from the different modulation frequencies. The new T-F representations of the IMPCT will have better performance than the SCT, as shown in Figures\ 2(c) and 2(d) with 50 iterations.Therefore, it is necessary to develop alternative approaches to deal with the challenging cases that there exist components overlapping in the T-F domain. We provide the following definition.Definition 1.
The two componentsxk(t) and xm(t) in model (3) are said to cross (or overlap) in the T-F domain. If there exist a time instant t^, a bound d≥0 and a step h shift up and down such that their IFs satisfy(3)IFkt^=IFmt^,IFkτ^±h=IFmτ^∓h,IFk′t^-IFm′t^≥d.Equation (3) illustrates a strict definition of these overlapped components for our algorithm to identify the case that some components are tangent to each other (i.e.,IFk′(t^)=IFm′(t^)). Figure 3 shows the strict overlapping and tangent of the two components. After shifting up and down, if there is a point such that IFk(τ^)±h=IFm(τ^)±h,t^≠τ^ at other intersecting time τ^, the components are strictly overlapping. It is worth noticing that the IFs of the arbitrary two components have different curvature at the intersecting time t^.Figure 3
(a) Two components are strictly overlapping. (b) Two components are tangent. Blue line, IF trajectory. Red line, IF trajectory after shift up. Purple line, IF trajectory after shift down.
(a) (b)A signalst=rmejφmt=Atcos⁡ft can be considered an intrinsic mode chirp function, with accuracy ε and A(t) and f(t) meeting the following conditions [14]:(4)A∈C1R∩L∞R,f∈C2Rinft∈RAt>0,inft∈Rft<∞supt∈RAt<∞,supt∈Rft<∞,supt∈Rf′t<∞A′t,f′t≤εft,∀t∈R.Here, to overcome these limitations, a robust technology architecture for overlapped nonstationary signals is proposed, as shown in Figure4. Initially, the input signal is passed by the N-point IMPCT stage. The stage brings about the IMPCT coefficients shown as r1,r2,⋯,rm. The inverse IMPCT operation for a signal sequence of length λ1,λ2,⋯,λm is performed. These characterized polynomial coefficients are concatenated and rearranged to the IFs estimation for each of the input nonstationary signals. Then, to decompose overlapped components in practical applications, we need to deal with a very significant but challenging issue of how to estimate IFs of overlapped components for the IMPCT. The second problem is how to determine the criterion which can separate the signals under micro-Doppler effect and noise interference. In the following sections, the solution for these two problems will be provided about detailed theoretical analysis.Figure 4
Block diagram of the proposed approach.


## 2.2. IMPCT Representation

The Fourier series analysis brings to light the sinusoidal property, so it may be utilized for the reconstruction of some signals with linear phase [14]. However, when processing nonlinear polynomial phase signals (NPPS) often involved in the practical application, a signal representations can be obtained by using the IMPCT. Namely, the IMPCT can be adopted to demodulating components of interest, reducing them to the sparse sinusoidal components. In this section, the Fourier series is stretched into the complex polynomial Fourier signal model as the IMPCT. Consider a multicomponent polynomial phase signal contaminated by a white Gaussian noise as(5)xt=∑m=1Mrmejφmt+εt,where rm is the amplitude coefficient function and φm(t)=α1mt+⋯+αnmtn/n! indicates the phase of the m-th component; ε(t) shows the noise. In general, the IMPCT of x(t) can be formulated as follows: (6)Xf1,⋯,fn=∑m=1Mrm×∫-∞∞ejα1m-f1t+⋯+αnm-fn/n!tndt.For all points in the n-dimensional frequency domain that correspond to the positions of signal components characterized by f1=α1m,...,fn=αnm,m=1,...,M, we can acquire (7)Xf1,⋯,fn=∑m=1M2πrmδf1-α1m,⋯,fn-αnm.That is to say, X(f1,⋯,fn)→∞ at the positions of signal components: f1=α1m,⋯,fn=αnm.Thus, X(f1,⋯,fn) can be deemed to a sparse M-component representation. In addition, for f1≠α1m,⋯,fn≠αnm,m=1,⋯,M has limited value, negligible when compared to (7).Let us inspect the signal vectorx with elements X(k), superimposed of a sum of M NPPS: (8)Xk=∑m=1MXmk+εk=∑m=1Mrmej2π/Kka1m+k2a2m/2+⋯+knanm/n!,where the polynomial coefficients are supposed to be bounded integers and K is the total number of discrete samples. To offer a sparse representation of nonlinear signals, the discrete form of the IMPCT is adopted, which may be expressed as (9)Xλ1,⋯,λn=∑k=0K-1∑m=1Mrmej2π/Kka1m+k2a2m/2+⋯+knanm/n!×e-j2π/Kk2λ2/2+⋯+knλn/n!e-j2π/Kkλ1.When the set of IMPCT parameters (λ1,⋯,λn) is selected to obtain the adaptive coefficients of the NPPS (a1i,a2i,⋯,ani), thus we can obtain (10)λ1,⋯,λn=a1m,a2m,⋯,anmthe m-th signal component is demodulated and the sinusoid rmej2π/K(ma1m) is paramount in the IMPCT spectrum. Under this circumstances, the energy of the spectrum is extremely concentrated at λ=a1m. On the contrary, while (λ1,⋯,λn)≠(a1m,a2m,⋯,anm) for all m∈1,⋯,M the spectrum consisted of some dispersed phase. Under high noise environment, we suppose that the minimum spectral component having amplitude coefficient Krmin can be also primary in the spectrum, meaning that it is above the aggregate value of all other components (in the worst case above their accumulation): K|rmin|>∑m=1M‍rm, where rmin is the minimum of all rm, for m=1,⋯,M. Therefore, a search is performed for M sets of parameters λ1m,⋯,λnm∣m=1,⋯,M so that (λ1,⋯,λn)=(a1m,a2m,⋯,anm) holds, and χλ1m,⋯,λnm is strongly concentrated. All in all, multicomponent set of parameters are calculated to obtain the optimal concentrated vector χ(λ1,⋯,λn) for all m=1,⋯,M. The maximized component will lie in the frequency λ1=a1m.The IMPCT representation explains the motivation for obtaining an initial IFs estimate in the algorithm. The IFs estimate will be based on the value ofχ(λ1,⋯,λn). In the next section, the T-F reassignment algorithm by the path penalty function is proposed in order to more accurately estimate IFs with the serious noise cases.


## 3. Multicomponent Decomposition

### 3.1. IFs Estimation by the T-F Reassignment

There are several methods to implement the IFs estimation. Here, we will adopt a reassignment by the path penalty function to further reduce search space. The set of the IFs estimation is detected in the selected points with the best path penalty function gain. The algorithm is performed recursively until we are no longer able to search out a new point that can reduce values of the path penalty function.In T-F plane, the frequenciesF and time sequence T can be contained. The set of paths between two ending points is FT, which aids in searching out the optimal path impossible. Fortunately, the algorithm can be implemented recursively as an example about the generalized Viterbi algorithm [26]. The optimal paths can be obtained as Pi(t;fj) which join the moment t1 and entire points to the moment ti, t∈t1,ti, for j∈1,M. It can be defined as(11)Pit;fj=argminλt∈Mijpλt;ti,ti,fj,j∈1,M,where the set Mij is obtained in (9) and contains all paths between the moment t1 and entire points to the moment ti. Simultaneously, a sum of the path penalty functions is revealed by pλ(t);ti,ti,fj for each ridge λ(t). In the Viterbi algorithm, a chirplet path (11) is known as the partial optimal path. Currently, the IFs estimation, within the time interval ti,ti+1, can be expressed as (12)ϖit=argminPit;fj,j∈1,MpPit;fj;ti,ti,fj.At the next moment ti+1, the partial optimal paths will be represented as the cascade of (11) with the points at the new moment Pi+1(t;fj)=Pi(t;fk),(ti+1;fj), j∈1,M, k∈1,M.By selecting paths with shorter chirplets, one can search out chirplets with an increasing functionR(ϖk,ϖj) of the distance |ϖj-ϖk| to further decrease search space. A new minimum value at the next moment is represented as (13)argmin⁡pPit;fk;ti,ti,fk+Rϖk,ϖj+Cχti+1,fjfor each ϖj. It is worth noting that C(χ(ti+1,fj)) is positive constant to capture solutions of the constrained partial optimal path with different values of length. To illustrate algorithm, an example about two overlapping signals is shown in Figure 5. In the considered case, an optimal path (IFs estimation) connects points (2,1) and (1,1).Figure 5
Illustration of the IFs estimation in the T-F plane.To form a best chirplet path using the path penalty functions, two key constraints in the connections stage from each point should be satisfied as(14)ϖj-ϖk≤Δϖ,λj-λk≤Δλwhere λj,λk stand for the estimated polynomial chirp phase in (9), respectively. The parameters Δϖ and Δλ are determined by the boundaries of the frequency modulation rates of the considered signal. In Figure 3, it can be easily seen that two best chirplet paths are synthesized by connecting the partial optimal path to points (2,1) and (1,1), as well as (1,2) and (2,2), respectively. Note that the new IFs estimation can perform an update of the previous value.Furthermore, we only need to solve the 1D optimization problem in partial optimal path to estimate IFs ridges. Thus, our method can be efficiently used in various real applications and a numerical overlapped signal mentioned in Figure2(a) is analyzed in the high noise. The relative result of the IFs estimation is shown in Figure 6. Each mode occupies difference boundary region of the T-F domain as shown in Figure 6(a), and the corresponding energy is concentrated around the IFs ridges in Figure 6(b).Figure 6
The relative IFs estimation result of a numerical signal. (a) The spectrogram of IMPCT; (b) the IFs estimation result.
(a) (b)In [27], the grid search method is adopted to solving the optimization problem. Especially assuming that the chirp phase search space is λj,λk, the set of the discrete polynomial chirp phase can adjust the grid resolution. Considering the training set in the kernel sparse learning, our method searches out proper path from a dictionary which is reconstructed by possible partial optimal paths in the set (13). Each subset is generated from a different starting point of the sample, implemented independently to link the each other and used to update a different pool of optimal paths. Therefore, a parallel optimization calculation can be adopted in our learning method.


### 3.2. Decomposition by Kernel Sparse Learning

One can rewrite equation (5) (15)s=∑m=1Mrmxm+η,where s=s(t0)⋯s(tN-1)T; η=n(t0)⋯n(tN-1)T; xm=xmaTxmbTTwith (16)xma=am0,⋯amL,a¯m0,⋯a¯mL,(17)xmb=bm0,⋯bmL,b¯m0,⋯b¯mL,where am0,⋯amL,a¯m0,⋯a¯mL and bm0,⋯bmL,b¯m0,⋯b¯mL indicate complex polynomial Fourier coefficients which should be identified. For each signal component x(t), its inversion formula corresponding to (7) is written as (18)xt1x∗t2=1N+1∑k=-N/2N/2Xit1+t22,k×ej2π/N+1kt1-t2,i=1,2,⋯,M,where x(t1) shows a column vector whose elements are the signal values and x∗(t2) indicates the Hermitian transpose of x(t1). The signal x(t) is supposed to be time limited within |n|≤N/2. We have obtained the IMPCT Xi(t,k) of this i-th component, according to (8). Through accumulating the above relationship for i=1,2,⋯,M, the following is available: (19)∑i=1Mxt1x∗t2=1N+1∑k=-N/2N/2∑i=1MXit1+t22,k×ej2π/N+1kt1-t2.Considering formula (9), for the signals meeting the proposed model, this relation can be rewritten as (20)∑i=1Mxt1x∗t2=1N+1∑k=-N/2N/2Xit1+t22,k×ej2π/N+1kt1-t2.By indicating (21)Rt1,t2=1N+1∑k=-N/2N/2Xit1+t22,k×ej2π/N+1kt1-t2and the decomposition of the kernel matrix R with the elements R(t1,t2), we obtain (22)R=∑i=1Nλiκinκi∗n,where κi(n) will be equal to the i-th signal components.To recover the expected components, we need to deal with nonlinear equations in (15). Since this inverse problem is commonly ill-posed, the kernel sparse learning is employed to process real data. In traditional kernel dictionary learning the dictionary is expressed as a linear combination of a nonlinear representation of the data in [23]. This is defined as(23)κX=κXΨ︸dictionaryY,where Y is defined as a series of coefficients and κX is nonlinear function. Ψ is the sparse weight matrix combining fixed basis in order to engender the dictionary atoms.In kernel sparse learning, the fixed basis is replaced by a nonlinear combination of the input. In (22), a kernel could only be used while solving the inverse problem. This allows expressing (22) as (24)ΨκXT︸transformκX=Y,where the kernelized data matrix is formulated instead of the data matrix. One can give the kernel upfront KX,X=κXTκX, referred to as K. We design the learning as (25)minΨ,YΨK-YF2+λΨF2-log⁡detΨ+μY0For kernel sparse learning, one can notice that the gradients for different terms and updating the coefficients in (25) are easy to calculate. This is given by (26)Y⟵ΨK≥μ⊙ΨK.When the amount of samples is much bigger than the dimension of the data, storing and processing the kernel matrix are great difficulty in terms of memory. Consequently, a more efficient scheme needs to be implemented. Using the eigenvalue decomposition, the kernel matrix can be denoted as(27)K=UΛUT,where U shows the eigenvector matrix and Λ indicates the diagonal matrix of eigenvalues sorted in descending order. Alternating minimization of (25) brings about (28)Ψ⟵minΨΨUΛ-YF2+λΨF2-log⁡detΨ,(29)Y⟵minYΨUΛ-YF2+μUY′0.Solving (28) is straightforwardly updated along with the update of IMPCT, identically; UΛ plays the role of data. In this study, we denote how (29) can be efficiently addressed by adopting the separation variable method.An agent variableP=UY′ is introduced, which can ideally ensure that the agents and variables are equivalent in each iteration. In reality, the equality is only enforced at convergence. Therefore, under the slight relaxation of constraints, the augmented Lagrangian method [28] can be applied to improve scalability and availability of KSL. One can rewrite (29) as (30)Y⟵minYΨUΛ-YF2+μP0+εP-UΛF2.In the above description, the equivalence between the agents and the variables can be guaranteed by the hyperparameter μ. By regulating the level of the value, the strength of the constraint will be controlled to enforce the equivalence at convergence. Initially, the algorithm sets a low value of μ and gradually enhances the value of the hyperparameter to efficiently deal with the optimization problem.In the high noise environment, one may choose more complex methods to make sure the number of components is accurate, such as by the split Bregman approach [29, 30]. To further relax the constraint, a Bregman variable is used in the following mode: (31)Y⟵minY,PΨUΛ-YF2+μP0+εP-UY′-ΨF2.A relatively large Ψ is conducive to make up for IFs estimation errors. One starts a high value of the Bregman variable Ψ but it is automatically updated in each iteration. Therefore, an effective approach is to automatically adapt the equality between the agents and the variables to the noise level such as the methods in [29].In order to achieve parallel computing, through adopting the alternating directions method of multipliers, (31) can be decomposed into the following two easier subproblems. (32)Y⟵minYΨUΛ-YF2+εP-UY′-ΨF2,(33)Y⟵minPμP0+εP-UY-ΨF2.It can be easily seen that the first subproblem (32) is an easy least squares problem forming a convergence in a mode of the pseudo inverse problem. In the subproblem (33), a hard threshold needs to be taken to enforce the equality between the agents and the variables at convergence. The model is efficient to deal with two subproblems in parallel. In addition, the KSL has the ability of select optimal training samples in order to minimize the error of the decomposition process while minimizing the number of learning samples to label and thus to cut down the costs interrelated with the training sample collection.To test the result of the training samples, we need to generate the features for the sample signalXt. Simultaneously, we obtain a learned dictionary as shown in Figure 7, where signals are merged only to provide a more valuable time-frequency representation from the IMPCT. It is worth noting that, in all the following simulation experiments, the initial training samples can be obtained by the IMPCT based on different feature modalities, such as noise, micro-Doppler effect, and scattering. Feature modalities quantity and the number of the generated samples corresponding to each feature modality are determined according to the complexity of the signal to be decomposed.Figure 7
The dictionary of a test processed the IMPCT by kernel sparse learning.The question is how to determine the dictionary size hyperparameterμ. Note that although the same initial value μ=1/N4 is adopted for all proposed numerical results, a more appropriate choice of this parameter may bring about better estimation results. The choice is an empirical balance between the MSE acquired on the training features, the sparsity of the kernel dictionary, and the interpret-ability of the final dictionary (decision criteria that rely on the specific application can also be considered).Having estimated the IFs of the signal (SNR=-2dB) in the previous section, the decomposition of signal with severely overlapped components is executed by KSL as shown in Figure 8. It indicates that all the components are successfully separated and the SNR is obviously improved. The detected IFs ridges from the TFR of the signal are generated in Figure 8(a). Figures 8(b) and 8(c) denote the corresponding individual components. The improvement of the superimposed of all components (SNR=16.75dB) is shown in Figure 8(d).Figure 8
Component separation for the simulated signal by the KSL.
(a) (b) (c) (d)


## 3.1. IFs Estimation by the T-F Reassignment

There are several methods to implement the IFs estimation. Here, we will adopt a reassignment by the path penalty function to further reduce search space. The set of the IFs estimation is detected in the selected points with the best path penalty function gain. The algorithm is performed recursively until we are no longer able to search out a new point that can reduce values of the path penalty function.In T-F plane, the frequenciesF and time sequence T can be contained. The set of paths between two ending points is FT, which aids in searching out the optimal path impossible. Fortunately, the algorithm can be implemented recursively as an example about the generalized Viterbi algorithm [26]. The optimal paths can be obtained as Pi(t;fj) which join the moment t1 and entire points to the moment ti, t∈t1,ti, for j∈1,M. It can be defined as(11)Pit;fj=argminλt∈Mijpλt;ti,ti,fj,j∈1,M,where the set Mij is obtained in (9) and contains all paths between the moment t1 and entire points to the moment ti. Simultaneously, a sum of the path penalty functions is revealed by pλ(t);ti,ti,fj for each ridge λ(t). In the Viterbi algorithm, a chirplet path (11) is known as the partial optimal path. Currently, the IFs estimation, within the time interval ti,ti+1, can be expressed as (12)ϖit=argminPit;fj,j∈1,MpPit;fj;ti,ti,fj.At the next moment ti+1, the partial optimal paths will be represented as the cascade of (11) with the points at the new moment Pi+1(t;fj)=Pi(t;fk),(ti+1;fj), j∈1,M, k∈1,M.By selecting paths with shorter chirplets, one can search out chirplets with an increasing functionR(ϖk,ϖj) of the distance |ϖj-ϖk| to further decrease search space. A new minimum value at the next moment is represented as (13)argmin⁡pPit;fk;ti,ti,fk+Rϖk,ϖj+Cχti+1,fjfor each ϖj. It is worth noting that C(χ(ti+1,fj)) is positive constant to capture solutions of the constrained partial optimal path with different values of length. To illustrate algorithm, an example about two overlapping signals is shown in Figure 5. In the considered case, an optimal path (IFs estimation) connects points (2,1) and (1,1).Figure 5
Illustration of the IFs estimation in the T-F plane.To form a best chirplet path using the path penalty functions, two key constraints in the connections stage from each point should be satisfied as(14)ϖj-ϖk≤Δϖ,λj-λk≤Δλwhere λj,λk stand for the estimated polynomial chirp phase in (9), respectively. The parameters Δϖ and Δλ are determined by the boundaries of the frequency modulation rates of the considered signal. In Figure 3, it can be easily seen that two best chirplet paths are synthesized by connecting the partial optimal path to points (2,1) and (1,1), as well as (1,2) and (2,2), respectively. Note that the new IFs estimation can perform an update of the previous value.Furthermore, we only need to solve the 1D optimization problem in partial optimal path to estimate IFs ridges. Thus, our method can be efficiently used in various real applications and a numerical overlapped signal mentioned in Figure2(a) is analyzed in the high noise. The relative result of the IFs estimation is shown in Figure 6. Each mode occupies difference boundary region of the T-F domain as shown in Figure 6(a), and the corresponding energy is concentrated around the IFs ridges in Figure 6(b).Figure 6
The relative IFs estimation result of a numerical signal. (a) The spectrogram of IMPCT; (b) the IFs estimation result.
(a) (b)In [27], the grid search method is adopted to solving the optimization problem. Especially assuming that the chirp phase search space is λj,λk, the set of the discrete polynomial chirp phase can adjust the grid resolution. Considering the training set in the kernel sparse learning, our method searches out proper path from a dictionary which is reconstructed by possible partial optimal paths in the set (13). Each subset is generated from a different starting point of the sample, implemented independently to link the each other and used to update a different pool of optimal paths. Therefore, a parallel optimization calculation can be adopted in our learning method.


## 3.2. Decomposition by Kernel Sparse Learning

One can rewrite equation (5) (15)s=∑m=1Mrmxm+η,where s=s(t0)⋯s(tN-1)T; η=n(t0)⋯n(tN-1)T; xm=xmaTxmbTTwith (16)xma=am0,⋯amL,a¯m0,⋯a¯mL,(17)xmb=bm0,⋯bmL,b¯m0,⋯b¯mL,where am0,⋯amL,a¯m0,⋯a¯mL and bm0,⋯bmL,b¯m0,⋯b¯mL indicate complex polynomial Fourier coefficients which should be identified. For each signal component x(t), its inversion formula corresponding to (7) is written as (18)xt1x∗t2=1N+1∑k=-N/2N/2Xit1+t22,k×ej2π/N+1kt1-t2,i=1,2,⋯,M,where x(t1) shows a column vector whose elements are the signal values and x∗(t2) indicates the Hermitian transpose of x(t1). The signal x(t) is supposed to be time limited within |n|≤N/2. We have obtained the IMPCT Xi(t,k) of this i-th component, according to (8). Through accumulating the above relationship for i=1,2,⋯,M, the following is available: (19)∑i=1Mxt1x∗t2=1N+1∑k=-N/2N/2∑i=1MXit1+t22,k×ej2π/N+1kt1-t2.Considering formula (9), for the signals meeting the proposed model, this relation can be rewritten as (20)∑i=1Mxt1x∗t2=1N+1∑k=-N/2N/2Xit1+t22,k×ej2π/N+1kt1-t2.By indicating (21)Rt1,t2=1N+1∑k=-N/2N/2Xit1+t22,k×ej2π/N+1kt1-t2and the decomposition of the kernel matrix R with the elements R(t1,t2), we obtain (22)R=∑i=1Nλiκinκi∗n,where κi(n) will be equal to the i-th signal components.To recover the expected components, we need to deal with nonlinear equations in (15). Since this inverse problem is commonly ill-posed, the kernel sparse learning is employed to process real data. In traditional kernel dictionary learning the dictionary is expressed as a linear combination of a nonlinear representation of the data in [23]. This is defined as(23)κX=κXΨ︸dictionaryY,where Y is defined as a series of coefficients and κX is nonlinear function. Ψ is the sparse weight matrix combining fixed basis in order to engender the dictionary atoms.In kernel sparse learning, the fixed basis is replaced by a nonlinear combination of the input. In (22), a kernel could only be used while solving the inverse problem. This allows expressing (22) as (24)ΨκXT︸transformκX=Y,where the kernelized data matrix is formulated instead of the data matrix. One can give the kernel upfront KX,X=κXTκX, referred to as K. We design the learning as (25)minΨ,YΨK-YF2+λΨF2-log⁡detΨ+μY0For kernel sparse learning, one can notice that the gradients for different terms and updating the coefficients in (25) are easy to calculate. This is given by (26)Y⟵ΨK≥μ⊙ΨK.When the amount of samples is much bigger than the dimension of the data, storing and processing the kernel matrix are great difficulty in terms of memory. Consequently, a more efficient scheme needs to be implemented. Using the eigenvalue decomposition, the kernel matrix can be denoted as(27)K=UΛUT,where U shows the eigenvector matrix and Λ indicates the diagonal matrix of eigenvalues sorted in descending order. Alternating minimization of (25) brings about (28)Ψ⟵minΨΨUΛ-YF2+λΨF2-log⁡detΨ,(29)Y⟵minYΨUΛ-YF2+μUY′0.Solving (28) is straightforwardly updated along with the update of IMPCT, identically; UΛ plays the role of data. In this study, we denote how (29) can be efficiently addressed by adopting the separation variable method.An agent variableP=UY′ is introduced, which can ideally ensure that the agents and variables are equivalent in each iteration. In reality, the equality is only enforced at convergence. Therefore, under the slight relaxation of constraints, the augmented Lagrangian method [28] can be applied to improve scalability and availability of KSL. One can rewrite (29) as (30)Y⟵minYΨUΛ-YF2+μP0+εP-UΛF2.In the above description, the equivalence between the agents and the variables can be guaranteed by the hyperparameter μ. By regulating the level of the value, the strength of the constraint will be controlled to enforce the equivalence at convergence. Initially, the algorithm sets a low value of μ and gradually enhances the value of the hyperparameter to efficiently deal with the optimization problem.In the high noise environment, one may choose more complex methods to make sure the number of components is accurate, such as by the split Bregman approach [29, 30]. To further relax the constraint, a Bregman variable is used in the following mode: (31)Y⟵minY,PΨUΛ-YF2+μP0+εP-UY′-ΨF2.A relatively large Ψ is conducive to make up for IFs estimation errors. One starts a high value of the Bregman variable Ψ but it is automatically updated in each iteration. Therefore, an effective approach is to automatically adapt the equality between the agents and the variables to the noise level such as the methods in [29].In order to achieve parallel computing, through adopting the alternating directions method of multipliers, (31) can be decomposed into the following two easier subproblems. (32)Y⟵minYΨUΛ-YF2+εP-UY′-ΨF2,(33)Y⟵minPμP0+εP-UY-ΨF2.It can be easily seen that the first subproblem (32) is an easy least squares problem forming a convergence in a mode of the pseudo inverse problem. In the subproblem (33), a hard threshold needs to be taken to enforce the equality between the agents and the variables at convergence. The model is efficient to deal with two subproblems in parallel. In addition, the KSL has the ability of select optimal training samples in order to minimize the error of the decomposition process while minimizing the number of learning samples to label and thus to cut down the costs interrelated with the training sample collection.To test the result of the training samples, we need to generate the features for the sample signalXt. Simultaneously, we obtain a learned dictionary as shown in Figure 7, where signals are merged only to provide a more valuable time-frequency representation from the IMPCT. It is worth noting that, in all the following simulation experiments, the initial training samples can be obtained by the IMPCT based on different feature modalities, such as noise, micro-Doppler effect, and scattering. Feature modalities quantity and the number of the generated samples corresponding to each feature modality are determined according to the complexity of the signal to be decomposed.Figure 7
The dictionary of a test processed the IMPCT by kernel sparse learning.The question is how to determine the dictionary size hyperparameterμ. Note that although the same initial value μ=1/N4 is adopted for all proposed numerical results, a more appropriate choice of this parameter may bring about better estimation results. The choice is an empirical balance between the MSE acquired on the training features, the sparsity of the kernel dictionary, and the interpret-ability of the final dictionary (decision criteria that rely on the specific application can also be considered).Having estimated the IFs of the signal (SNR=-2dB) in the previous section, the decomposition of signal with severely overlapped components is executed by KSL as shown in Figure 8. It indicates that all the components are successfully separated and the SNR is obviously improved. The detected IFs ridges from the TFR of the signal are generated in Figure 8(a). Figures 8(b) and 8(c) denote the corresponding individual components. The improvement of the superimposed of all components (SNR=16.75dB) is shown in Figure 8(d).Figure 8
Component separation for the simulated signal by the KSL.
(a) (b) (c) (d)


## 4. Validation

It is known that micro-Doppler effect has been diffusely applied to target detection and recognition in sonar systems [31, 32]. In general, such signals from scattering points of the transmitted signal (or from the time delay and Doppler shift of the received echo) may seriously overlap in the T-F domain. To accurately characterize the TFR of micro-Doppler signals, it is key to decompose overlapped signal components, first and foremost.In this section, the presented approach will be used to decompose both simulated micro-Doppler signals to verify the performance of the proposed method. The maximum number of iterations is set to 80 and we adopt a more suitableμ=2/N4. To quantify the noise level of the inspected signal or the accuracy of IFs estimation and signal reconstruction, the SNR (unit: dB) of the signal is defined as (34)SNR=10log10⁡ωt22ωt-ω~t22where ω(t) indicates original signal or true IFs; ω~(t) is reconstructed signal or estimation IFs.In all simulated examples, a zero-mean Gaussian white noiseε(t) with variance σ2 and with i.i.d. real and imaginary parts is artificially added to the resulting micro-Doppler signal to analyze the effect of complicated AWGN.Example 2.
Consider an FM signal(35)st=sin⁡2π25t+10sin⁡t+sin⁡30+2π25t-10sin⁡t,induced by some typical moving targets [2] (such as rotating, tumbling, and coning targets); the signal is contaminated with noise with SNR=-5dB. Herein, the sampling frequency is 100 Hz and the time duration of the signal is set to 15 s. The IMPCT of the signal, as the starting point of the presented algorithm, is displayed in Figure 9(a).
The experimental results are observed after every 10 iterations, until the 80th iteration, are shown in Figures9(b)–9(i). It can be easily seen that after the 10th iteration a large amount of elements has concentrated the spectrum energy at many T-F points. However, due to the nonlinear polynomial phase effect in (6), the spectrum energy has already started to concentrate at the autoterm, as shown in Figure 9(b). With the iteration, the stronger the spectral energy is enhanced, in contrast, the weaker the noise interference is. The optimal TFR with a superior energy concentration from the 80th iteration is shown in Figure 9(i) due to the energy positioning update mechanism (7) and a sparse representation (8) avoiding the noise interference, used in (9).
The detected IFs ridges from the TFR of the signal are revealed in Figures10(a)–10(c). It reveals the results of IFs ridges with 0th iteration, 10th iteration, and 80th iteration, respectively. Before the 10th iteration, the T-F reassignment detector by the path penalty function cannot track precise IFs for the micro-Doppler signal with highly overlapped components. On the other hand, the proposed T-F reassignment can effectively distinguish these crossed components by considering the partial optimal path of the ridges. To abate the noise interference, one can further smooth the estimated IFs by nonlinear spline adaptive fitting [33] with a polynomial Fourier model, as illustrated in Figure 10(c).
Having estimated the IFs ridges, the separation of the multicomponent signal is executed by the KSL, as shown in Figures10(d) and 10(e). It displays that all the components are successfully discomposed and the SNR is availably improved.Figure 9
The dictionary of a test processed the IMPCT by kernel sparse learning.
(a) (b) (c) (d) (e) (f) (g) (h) (i)Figure 10
Estimated IFs and component separation for the micro-Doppler signal in Figure9 (SNR=10.16dB).
(a) (b) (c) (d) (e)Example 3.
Here, we consider again a poly-harmonic FM signal, expressed as(36)st=s1t+s2t+s3twith (37)s1t=cos⁡2π5sin⁡π4t+6t+1.2t2s2t=cos⁡2π0.2+45t-1.2t2s3t=cos⁡2π20sin⁡π4t+20t+0.2t2,where a sampling frequency is set to 100 Hz and the its time duration is 15 s. The signal is contaminated by high white noise and the SNR is 0 dB. The TFRs acquired by the presented IMPCT are compared with the results acquired by SWT and SCT, as shown in Figure 11. In addition, IFs ridges of the noise signal are estimated by the proposed T-F reassignment (see Figures 11(c), 11(f), and 11(i)).
The TFRs of nonnoisy and noise signal obtained by SWT, shown in Figures11(a) and 11(b), respectively, can barely show the inherent T-F pattern of the signal. Specifically, the TFRs only reveals the clear IFs trajectory for the autoterms of signal due to an easy approximation by a linear function. However, the TFRs of the cross terms are too blur to accurately estimate the IFs ridges of the overlapping components. In Figures 11(d) and 11(e), it can be noticed that the TFRs obtained by SCT scatter the energy around the IFs due to their rough frequency resolution induced by the high noise, and it is inaccurate to extract the IFs. In Figures 11(g) and 11(h), it is revealed that, due to the proposed optimal concentrated vector, the IMPCT surpasses SWT and SCT as it clearly shows the true T-F pattern of the signal. The IMPCT generates the TFRs with an outstanding concentration, based on which the IFs can be extracted accurately (see Figure 11(i)). However, in Figures 11(c) and 11(f), it can be seen that partial IFs ridges of cross terms corresponding to some points are truncated in disorder. The simulation results display that both methods suffer from high noise and cross terms. Obviously, one cannot accurately decompose the overlapped components by these two methods.
For all the following tests, the MSE versus variable SNR will be calculated (from -10 to 10 dB below, with a 1 dB step). The MSE is expressed as(38)MSE=Eξ^n-ξTn2,where ξ^(n) indicates the estimated IFs and ξT(n) the true IFs index.
We consider the IFs estimation from the TFRs of IMPCT, SWT, and SCT, using the T-F reassignment based on the path penalty function. The obtained MSEs are displayed in Figure12. Signals with fast IFs variations have been chosen to determine good performance of the presented method in such scenarios. Because of the distributed property of the optimization, the IMPCT is robust to characterize the TFRs, i.e., to the nonlinearity IFs ridges. Analyzing the results in Figure 12, we demonstrate that, through slightly promoting the IFs variations or adding the noise, the better improvement level of IMPCT over SWT and SCT is obtained.Figure 11
TFRs of the signal obtained SWT, SCT, and IMPCT, respectively, and estimated IFs by the proposed T-F reassignment.
(a) (b) (c) (d) (e) (f) (g) (h) (i)Figure 12
IFs estimation MSEs for signals in Experiment 2, using the proposed T-F reassignment.


## 5. Conclusion

In this paper, we presented a novel approach which combines the KSL algorithm with the T-F reassignment to decompose seriously overlapped components. The IMPCT is first proposed for the TFRs of a multicomponent nonstationary signal in high noise and can adapt to some practical applications. Comparing with traditional the decomposition approaches, the KSL based on the T-F reassignment can decompose signals with overlapped components. The modified T-F reassignment algorithm by the best path penalty function estimates the main IFs ridges of all the components and then reconstructs the obtained IFs ridges according to the partial optimal path. To estimate performance of the presented method, we analyze the numerical experiments including a real-world signal, which clearly illustrates that it outperforms the IMPCT based on the T-F reassignment and KSL in high noise environments. In addition, it outperforms the state-of-the-art TFRs methods when polynomial phase signal is contaminated by a white Gaussian noise. It should be emphasized that, in essence, the KSL algorithm involves postprocessing of ridges to optimize the IFs estimation of components contaminated by noise or induced by moving targets. Further, the applications in practical scenarios will be the part of our research development.



---

*Source: 8426790-2018-07-26.xml*