Enhanced Micro-Doppler Feature Extraction Using Adaptive Short-Time Kernel-Based Sparse Time-Frequency Distribution

Abstract

The extraction of the micro-Doppler (m-D) feature based on time-frequency distribution (TFD) is of great significance for target detection and identification. To improve the feature extraction performance, numerous TFDs have been developed, with the majority falling under Cohen’s class. Nevertheless, these TFDs basically face a trade-off between artifact suppression and energy concentration. The main reason is that each Cohen’s class TFD is constructed by applying the two-dimensional Fourier transform to a kerneled ambiguity function directly, while existing kernels generally attenuate artifacts at the expense of losing valuable information. In this paper, a TFD reconstruction method employing an adaptive short-time kernel (ASTK) is developed in the framework of sparse representation (SR) theory to overcome this trade-off and enhance the m-D feature. Firstly, the task of the optimal kernel is explained from the viewpoint of the instantaneous auto-correlation function (IAF). Secondly, based on the quasi-linear frequency modulation feature of most m-D signals during short-time periods, the distribution rule of the short-time IAF (STIAF) in the ambiguity plane is concluded. Guided by this rule, an ASTK that can effectively remove unwanted artifacts with the least information loss is designed. Finally, an SR-based reconstruction procedure is conducted on the kerneled STIAF to generate an artifact-free TFD with high energy concentration, which can effectively enhance the m-D feature. Experiments using both simulated and real-world m-D signals demonstrate the effectiveness of the proposed method.

1. Introduction

Micro-motions refer to the small-amplitude movements of radar targets or structures on the targets, in addition to the bulk motion. Micro-motions are commonly observed in the real world, such as the vibrations of a vehicle engine, the rotation of helicopter rotor blades, or the tumbling of a spacecraft in space [1]. From the perspective of radar, these micro-motions can cause unique frequency modulations in the radar echo, resulting in the generation of sidebands around the target’s Doppler frequency shift. This phenomenon is known as the micro-Doppler (m-D) effect.

In recent years, the value of the m-D feature as one of the most significant target signatures has been confirmed. On the one hand, the m-D feature enables the differentiation between targets and backgrounds, which can effectively aid target detection [2,3]. On the other hand, the m-D modulation parameters such as repetition period and maximum amplitude indicate the motion or physical properties of the target, which provides valuable additional information for precise target detection and identification [3,4,5]. Numerous algorithms have been developed to analyze and extract the m-D feature from radar measurements, among which the time-frequency distribution (TFD) has garnered significant attention [6]. The TFD is capable of visualizing the frequency content of an m-D signal over time, thus enabling a clear display of intricate time-varying m-D modulation information of the target. In [7], the Short-Time Fourier Transform (STFT) was utilized to reveal the m-D feature of rotationally symmetric ballistic targets with precession. In [8], radar measurements of diverse targets including humans, vehicles, dogs, and sheep were processed and compared using STFT. Marple proposed time-frequency sharpening techniques to extract the m-D characteristics of helicopter echo signals, which further supported the accurate classification of four-blade and two-blade helicopters [9,10]. In the analysis of human gait echoes, researchers pointed out that extracting the m-D feature from the low-resolution TFD generated by STFT is technically challenging since m-D modulation curves corresponding to different parts of the body overlap with each other. However, using the high-resolution TFD obtained by the reassignment method allows for a clear distinction of the m-D from nine different parts of the human body [11,12]. Therefore, it is crucial to employ a TFD with high energy concentration for the precise extraction of a target’s m-D feature.

At first glance, the Wigner-Ville distribution (WVD) may stand out as one of the most appealing options in m-D feature extraction because it is well-known that WVD has excellent TF concentration, particularly for mono-component linear frequency modulation (LFM) signals. However, as a quadratic TFD, WVD produces desired auto-terms (ATs) and unwanted artifacts simultaneously. To be more specific, the artifacts comprise inner and outer artifacts (commonly referred to as cross-terms), which are induced by the nonlinear frequency modulation (FM) characteristic and the multi-component nature of the signal, respectively [13,14]. Generally, the complex micro-motions of targets cause time-varying frequency modulations in the radar echo, and the observation scene rarely consists of a single target. As a result, the presence of artifacts in WVD is almost inevitable. These artifacts significantly obscure the target’s m-D information, thereby reducing the readability and practicality of WVD in m-D analysis [15,16].

Quantities of modified versions of WVD have been proposed to remove or at least attenuate artifacts, most of which can be categorized under Cohen’s class [17]. Based on the fact that ATs mainly cluster near the origin in the ambiguity plane, while artifacts tend to lie farther away from it [18,19], each member of Cohen’s class involves applying a low-pass kernel to the signal’s ambiguity function (AF) and then obtaining an artifact-suppressed TFD by taking the two-dimensional (2-D) Fourier transform (FT) of the weighted AF. Obviously, the performance of such a TFD relies heavily on the choice of the kernel.

So far, various kernels have been proposed according to the aforementioned distribution pattern of the AF [20,21,22,23,24,25,25]. However, existing kernels, including the well-known adaptive optimal kernel (AOK) [23], face a trade-off between AT preservation and artifact excision. The cause of this issue is that the words “near” and “farther” used to differentiate ATs and artifacts are inherently vague, and thus cannot provide precise guidance for kernel design. Moreover, applying a crude Fourier inversion to the kerneled AF without any constraint tends to reduce the energy concentration [26]. Consequently, the performance of Cohen’s class TFDs is fundamentally limited.

Inspired by the local quasi-LFM feature of most real-life signals [27,28,29], this paper proposes an effective TFD reconstruction approach to generate high-performance TFDs of radar echoes to enhance the target’s m-D characteristic. The main contributions are as follows.

• An explicit distribution rule of the signal’s short-time AF (STAF) is derived, providing a clear explanation of a phenomenon that was vaguely described in previous literature.
• The task of the optimal kernel is analyzed and clearly illustrated. Moreover, an adaptive short-time kernel (ASTK) is developed according to the concluded distribution rule, which achieves excellent artifact removal and accurate AT preservation simultaneously.
• A sparse representation (SR)-based reconstruction method for the instantaneous spectrum is proposed, which further facilitates the construction of an artifact-free TFD with high energy concentration to enhance the m-D characteristic.

This paper is organized as follows. In Section 2, the task of the optimal kernel is analyzed and the ASTK is derived based on the ambiguity domain distribution pattern of the m-D signal’s STIAF. Following this, the reconstruction process of the proposed high-performance TFD is detailed. Section 3 compares the proposed TFD with four other methods using simulated and real m-D signals to demonstrate its superiority. Finally, Section 4 presents conclusions and outlines future research directions.

2. Proposed M-D Feature Enhancement Method

2.1. Task of Optimal Kernel

In general, an m-D signal and its idealized TFD can be respectively represented as [30]

$$s\left(t\right)=\sum _{k=1}^K{s}_k\left(t\right)=\sum _{k=1}^K{a}_k\left(t\right)e^{j{\phi }_k\left(t\right)},$$

(1)

$$P_{\mathrm{ideal}}(t,\omega )=\sum _{k=1}^K{a}_k^2\left(t\right)\delta \left(\omega -{\phi }_k^{\prime }\left(t\right)\right),$$

(2)

where K is the number of signal components, $s_k\left(t\right)$ denotes the k-th component with time-varying amplitude $a_k\left(t\right)$, phase ${\phi }_k\left(t\right)$, $\delta (·)$ represents the Dirac delta function, and ${\phi }_k^{\prime }\left(t\right)$ stands for the first order derivative of ${\phi }_k\left(t\right)$. Conceptually, TFDs in Cohen’s class are generated. All integrals run from $-\infty$ to ∞ unless otherwise noted.

$$P(t,\omega )=\frac{1}{4{\pi }^2}\int \int A(\theta ,\tau )\mathsf{\Phi }(\theta ,\tau )e^{-j\theta t-j\tau \omega }d\theta d\tau ,$$

(3)

where $A(\theta ,\tau )$ is the AF of $s\left(t\right)$, $\mathsf{\Phi }(\theta ,\tau )$ denotes the TFD kernel, $\theta$ and $\tau$ signify the Doppler frequency and the correlation lag, respectively. More specifically, $A(\theta ,\tau )$ is given by

$$A(\theta ,\tau )=\int R(t,\tau )e^{j\theta t}dt,$$

(4)

where $R(t,\tau )$ stands for the instantaneous auto-correlation function (IAF) of $s\left(t\right)$, which is expressed as

$$R(t,\tau )=s(t+\frac{\tau }{2})s^{*}(t-\frac{\tau }{2}),$$

(5)

${(·)}^{*}$ represents the conjugate operator.

By comparing (2) and (3), it can be observed that if $A(\theta ,\tau )\mathsf{\Phi }(\theta ,\tau )$ can be replaced by $2\pi \int {\sum }_{k=1}^K{a}_k^2\left(t_1\right)e^{j{\phi }_k^{\prime }\left(t_1\right)\tau }{e}^{j\theta t_1}d{t}_1$, we have

$$\begin{array}{cc}\hfill P(t,\omega )& =\frac{1}{4{\pi }^2}\int \int A(\theta ,\tau )\mathsf{\Phi }(\theta ,\tau )e^{-j\theta t-j\tau \omega }d\theta d\tau \hfill \\ \hfill & =\frac{1}{4{\pi }^2}\int \int \left(2\pi \int \sum _{k=1}^K{a}_k^2\left(t_1\right)e^{j{\phi }_k^{\prime }\left(t_1\right)\tau }{e}^{j\theta t_1}d{t}_1\right)e^{-j\theta t-j\tau \omega }d\theta d\tau \hfill \\ \hfill & =\int \left(\int \left(\frac{1}{2\pi }\int \sum _{k=1}^K{a}_k^2\left(t_1\right)e^{j{\phi }_k^{\prime }\left(t_1\right)\tau }{e}^{j\theta t_1}d{t}_1\right)e^{-j\theta t}d\theta \right)e^{-j\tau \omega }d\tau \hfill \\ \hfill & =\int \sum _{k=1}^K{a}_k^2\left(t\right)e^{j{\phi }_k^{\prime }\left(t\right)\tau }{e}^{-j\tau \omega }d\tau \hfill \\ \hfill & =\sum _{k=1}^K{a}_k^2\left(t\right)\delta \left(\omega -{\phi }_k^{\prime }\left(t\right)\right)\hfill \end{array}$$

(6)

That is to say, Cohen’s class TFD can be transformed into the idealized TFD under the condition of (6). Therefore, the optimal kernel is supposed to be

$${\mathsf{\Phi }}_{\mathrm{opt}}(\theta ,\tau )=\frac{\int {\sum }_{k=1}^K{a}_k^2\left(t\right)e^{j{\phi }_k^{\prime }\left(t\right)\tau }{e}^{j\theta t}dt}{\int R(t,\tau )e^{j\theta t}dt},$$

(7)

where the amplitude of the kernel is neglected.

Expanding ${\phi }_k(t+\frac{\tau }{2})$ and ${\phi }_k(t-\frac{\tau }{2})$ into Taylor series at t, we have

$${\phi }_k(t+\frac{\tau }{2})={\phi }_k(t)+{{\phi }_k}^{\prime }(t)·\frac{\tau }{2}+\frac{{{\phi }_k}^{″}(t)}{2!}·{(\frac{\tau }{2})}^2+\cdots$$

(8)

$${\phi }_k(t-\frac{\tau }{2})={\phi }_k(t)+{{\phi }_k}^{\prime }(t)·(-\frac{\tau }{2})+\frac{{{\phi }_k}^{″}(t)}{2!}·{(-\frac{\tau }{2})}^2+\cdots$$

(9)

where ${\phi }_k^{″}\left(t\right)$ stands for the second order derivative of ${\phi }_k\left(t\right)$.

Substituting (1) as well as (8) and (9) into (5), $R(t,\tau )$ can be expanded as

$$R(t,\tau )=\sum _{k=1}^K{R}_{kk}(t,\tau )+\sum _{k=1}^K\sum _{l=1,l\ne k}^K{R}_{kl}(t,\tau )$$

(10)

where $R_{kk}(t,\tau )$ represents the IAF of $s_k\left(t\right)$ and $R_{kl}(t,\tau )$ represents the instantaneous cross-correlation function (ICF) between $s_k\left(t\right)$ and $s_l\left(t\right)$ as the l-th component of the m-D signal. $R_{kk}(t,\tau )$ and $R_{kl}(t,\tau )$ are respectively given by

$$\begin{array}{cc}\hfill R_{kk}(t,\tau )& =s_k(t+\frac{\tau }{2})s_k^{*}(t-\frac{\tau }{2})=a_k^2\left(t\right)e^{j\left[{\phi }_k(t+\frac{\tau }{2})-{\phi }_k(t-\frac{\tau }{2})\right]}\hfill \\ \hfill & \approx a_k^2\left(t\right)e^{j\left({{\phi }_k}^{\prime }\left(t\right)\tau +\underset{\text{high-order-derivative terms}}{\underbrace{{\displaystyle \sum _{m=1}^{\infty }}\frac{{{\phi }_k}^{(2m+1)}\left(t\right)}{(2m+1)!·2^{2m}}{\tau }^{(2m+1)}}}\right)}\hfill \end{array}$$

(11)

$$R_{kl}(t,\tau )=s_k(t+\frac{\tau }{2})s_l^{*}(t-\frac{\tau }{2})=a_k\left(t\right)a_l\left(t\right)e^{j\left[{\phi }_k(t+\frac{\tau }{2})-{\phi }_l(t-\frac{\tau }{2})\right]}$$

(12)

In (11) and (12), the amplitude of each component is assumed to vary slowly, and thus $a_k(t\pm \frac{\tau }{2})=a_k\left(t\right)$. ${\phi }_k^{(2m+1)}\left(t\right)$ stands for the $(2m+1)$ order derivatives of ${\phi }_k\left(t\right)$.

Referring to (7)–(12), the task of the optimal kernel can be more precisely interpreted. In essence, ${\mathsf{\Phi }}_{\mathrm{opt}}(\theta ,\tau )$ aims to eliminate the $R_{kl}(t,\tau )$ as well as the high-order-derivative terms within each $R_{kk}(t,\tau )$ in the ambiguity domain. This is because the former introduces outer artifacts while the latter contributes to the inner artifacts. However, designing such a kernel poses a technical challenge because it is still unclear how to accurately differentiate the terms that should be captured from the others that need to be excluded. The remainder of this section presents a practical scheme.

2.2. Distribution Rule of STIAF in Ambiguity Domain

In fact, a more explicit distribution rule for each term in $R(t,\tau )$ can be derived with the idea of local signal processing. A $t_0$-centered short-time m-D signal is represented as

$$s_{t_0}\left(t\right)=s\left(t\right)w(t-t_0),$$

(13)

where $w\left(t\right)$ is a symmetrical window function with length L, i.e., $w\left(t\right)\ne 0$ only for $t\in [-\frac{L}{2},\frac{L}{2}]$. In this paper, $w\left(t\right)$ is set to be a rectangular window.

With a proper value of L, $s_{t_0}\left(t\right)$ can be precisely approximated as the superposition of K LFM components, which is exactly the local quasi-LFM feature commonly observed in real-life signals [27,28,29]. In this case, if we expand each phase function into the corresponding Taylor series at $t_0$, all terms higher than the second can be omitted. That is,

$$s_{t_0}\left(t\right)\approx \sum _{k=1}^K{a}_k\left(t\right)\underset{\mathrm{LFM}\mathrm{signal}}{\underbrace{e^{j\left(x_k{t}^2+y_{k}t+z_k\right)}}},t\in [t_0-\frac{L}{2},t_0+\frac{L}{2}],$$

(14)

where $x_k={\phi }_k^{″}\left(t_0\right)/2$, $y_k=({\phi }_k^{\prime }\left(t_0\right)-{\phi }_k^{″}\left(t_0\right)t_0)$, $z_k={\phi }_k\left(t_0\right)-{\phi }_k^{\prime }\left(t_0\right)t_0+\frac{{\phi }_k^{″}\left(t_0\right)}{2}{t}_0^2$. Similarly, for the l-th component of $s_{t_0}\left(t\right)$, $x_l={\phi }_l^{″}\left(t_0\right)/2$, $y_l=({\phi }_l^{\prime }\left(t_0\right)-{\phi }_l^{″}\left(t_0\right)t_0)$, $z_l={\phi }_l\left(t_0\right)-{\phi }_l^{\prime }\left(t_0\right)t_0+\frac{{\phi }_l^{″}\left(t_0\right)}{2}{t}_0^2$.

The AF of $s_{t_0}\left(t\right)$, i.e., the short-time AF (STAF) of $s\left(t\right)$, is defined as

$$A_{t_0}(\theta ,\tau )=\int R_{t_0}(t,\tau )e^{j\theta t}dt,$$

(15)

where $R_{t_0}(t,\tau )$ stands for $s\left(t\right)$’s short-time IAF (STIAF) and is expressed as

$$\begin{array}{cc}\hfill R_{t_0}(t,\tau )& =s_{t_0}(t+\frac{\tau }{2})s_{t_0}^{*}(t-\frac{\tau }{2})\hfill \\ \hfill & =\sum _{k=1}^K{R}_{t_0,kk}(t,\tau )+\sum _{k=1}^K\sum _{l=1,l\ne k}^K{R}_{t_0,kl}(t,\tau ),\hfill \end{array}$$

(16)

$R_{t_0,kk}(t,\tau )$ and $R_{t_0,kl}(t,\tau )$ represent the time-localized version of $R_{kk}(t,\tau )$ and $R_{kl}(t,\tau )$, respectively. The range of t is $[t_0-\frac{L}{2},t_0+\frac{L}{2}]$, while $\tau$ is limited in $[-L,L]$.

Using (14), $R_{t_0,kk}(t,\tau )$ and $R_{t_0,kl}(t,\tau )$ can be obtained as

$$R_{t_0,kk}(t,\tau )\approx a_k^2\left(t\right){\mathrm{e}}^{j(2x_{k}t+y_k)\tau }={\mathrm{e}}^{j\left({\phi }_k^{\prime }\left(t_0\right)+{\phi }_k^{″}\left(t_0\right)(t-t_0)\right)\tau }$$

(17)

$$R_{t_0,kl}(t,\tau )\approx a_k\left(t\right)a_l\left(t\right){\mathrm{e}}^{j(x_{kl}{t}^2+y_{kl}t+z_{kl})}$$

(18)

where $x_{kl}=x_k-x_l$, $y_{kl}=(x_k+x_l)\tau +y_k+y_l$, $z_{kl}=(x_k-x_l)\frac{{\tau }^2}{4}+(y_k-y_l)\frac{\tau }{2}+z_k-z_l$.

For variable t, $R_{t_0,kk}(t,\tau )$ in (17) is a narrow-band single-frequency signal with constant frequency $2x_k\tau$ and initial phase $y_k\tau$, but $R_{t_0,kl}(t,\tau )$ in (18) is a wide-band LFM signal with chirp rate $2x_{kl}$, initial frequency $y_{kl}$, and initial phase $z_{kl}$. Since STAF is exactly the inverse FT of STIAF with respect to t, for each value of $\tau$, the energy of $R_{t_0,kk}(t,\tau )$ is highly concentrated in the ambiguity domain, whereas the energy of $R_{t_0,kl}(t,\tau )$ disperses across its instantaneous bandwidth. Considering all values of $\tau$, a more explicit distribution pattern in the ambiguity plane can be derived.

If the signal’s frequency exhibits approximately linear variation, the energy of auto-terms is distributed along lines passing through the origin, while the energy of artifacts is dispersed within a diamond-shaped region.

To illustrate the above conclusion, Figure 1 shows the STAF of a two-component LFM signal and the corresponding cross-section at $\tau =0.5$ s. The LFM signal is expressed as

$$s_1\left(t\right)=e^{j\left(-100\pi t+200\pi t^2\right)}+e^{j\left(100\pi t-200\pi t^2\right)}.$$

(19)

Through an inspection of the results shown in Figure 1, the correctness of the concluded distribution rule can be verified. It is worth noting that the prerequisite for applying this rule is the quasi-LFM feature of the short-time signal. Hence, the analysis window is supposed to break down the original m-D signal into a series of short-time segments with linearly varying frequencies. To achieve this, the adaptive window proposed in [27] is employed, which is generated through the minimization of the bandwidth of each short-time signal. The reason is that the more stationary the short-time signal, the slower its frequency varies, leading to a smaller bandwidth of the short-time signal. Therefore, it is reasonable to consider the bandwidth of each short-time signal as an indicator of LFM approximation precision. A smaller bandwidth indicates higher approximate precision.

2.3. ASTK Design

By comparing (11) and (17), it is apparent that the high-order-derivative terms within each $R_{kk}(t,\tau )$ are absent in its time-localized counterpart $R_{t_0,kk}(t,\tau )$. Therefore, the primary objective of the short-time kernel for $s_{t_0}\left(t\right)$ shifts toward the suppression of $R_{t_0,kl}(t,\tau )$. Based on the STIAF’s distribution rule concluded in Section 2.2, if we accumulate the energy of the STAF through lines passing through the origin as

$$E\left(\alpha \right)={\int }_0^{\infty }{A}_{t_0}\left(\rho cos\alpha ,\rho sin\alpha \right)d\rho ,\alpha \in [0,2\pi ),$$

(20)

energy peaks can be obtained in the directions that each $R_{t_0,kk}(t,\tau )$ is distributed along, while the energy of each $R_{t_0,kl}(t,\tau )$ cannot be effectively accumulated from any angle. In (20), $E\left(\alpha \right)$ denotes the accumulated energy along a line passing through the origin with an inclination angle of $\alpha$, and $\rho$ represents the polar radius. Prior to the energy accumulation, an interpolation of the STAF from the rectangular grid to the polar grid is necessary. A detailed description of the interpolation process can be found in [24].

Finally, the ASTK can be derived by assigning the value of 1 to directions that correspond to the energy peaks while maintaining 0 for all other locations. As a result, the shape of the kernel is dynamically determined through the energy accumulation procedure, allowing it to adaptively adjust to match various short-time signals. For example, an X-shaped kernel will be derived for the STAF shown in Figure 1. Using this kernel, each $R_{t_0,kl}(t,\tau )$ can be effectively removed while each $R_{t_0,kk}(t,\tau )$ is precisely preserved. As a result, we have

$$\begin{array}{cc}\hfill r_{t_0}(t,\tau )& =\underset{\theta \to t}{\mathcal{F}}\left\{A_{t_0}(\theta ,\tau ){\mathsf{\Phi }}_{t_0}(\theta ,\tau )\right\}\approx \sum _{k=1}^K{R}_{t_0,kk}(t,\tau )\hfill \\ \hfill & \approx \sum _{k=1}^K{a}_k^2\left(t\right){\mathrm{e}}^{j\left({\phi }_k^{\prime }\left(t_0\right)+{\phi }_k^{″}\left(t_0\right)(t-t_0)\right)\tau }\hfill \end{array}$$

(21)

where $r_{t_0}(t,\tau )$ stands for the kerneled STIAF, t is limited in $[t_0-\frac{L}{2},t_0+\frac{L}{2}]$, $\underset{\theta \to t}{\mathcal{F}}$ denotes the FT from Doppler frequency $\theta$ to time t.

2.4. SR-Based TFD Reconstruction

Traditionally, once $r_{t_0}(t,\tau )$ is determined, the instantaneous spectrum at time instant $t_0$ can be computed as [24]

$$P(t_0,\omega )=\underset{\tau \to \omega }{\mathcal{F}}\left\{r_{t_0}(t_0,\tau )\right\}=\underset{\tau \to \omega }{\mathcal{F}}\left\{\sum _{k=1}^K{a}_k^2\left(t_0\right){\mathrm{e}}^{j{\phi }_k^{\prime }\left(t_0\right)\tau }\right\}.$$

(22)

where $\underset{\tau \to \omega }{\mathcal{F}}$ denotes the FT from correlation lag $\tau$ to frequency $\omega$.

By continuously sliding the analysis window and computing the spectrum for each time instant, an artifact-suppressed TFD (called ASTK-FT TFD hereafter) can be generated by arranging the obtained spectrums chronologically. However, since the range of $\tau$ for each $r_{t_0}(t,\tau )$ is limited by the length of the short-time signal, the concentration of each instantaneous spectrum is usually unsatisfactory due to the uncertainty principle [22]. Hence, ASTK-FT TFD still inherits poor energy concentration.

Since K is usually a small number in practice, the signal’s inherent sparsity can be utilized to generate “sharper” spectrums [31,32]. By employing a $l_0$ norm sparsity constraint, the sparse spectrum at time instant $t_0$ can be computed as

$$\begin{array}{cc}\hfill \widehat{P}(t_0,\omega )=& \underset{P(t_0,\omega )}{argmin}{∥P(t_0,\omega )∥}_0\hfill \\ \hfill & \mathrm{s}.\mathrm{t}.{\underset{\tau \leftarrow \omega }{\mathcal{F}}}^{-1}\left\{P(t_0,\omega )\right\}=r_{t_0}(t_0,\tau )\hfill \end{array}$$

(23)

where ${∥·∥}_0$ signifies the $l_0$ norm, and ${\underset{\tau \leftarrow \omega }{\mathcal{F}}}^{-1}$ stands for the inverse FT from frequency $\omega$ to correlation lag $\tau$.

It is well-known that sparse representation problems based on the $l_0$ norm are NP-hard. In this paper, an effective sparse decomposition algorithm named SL0 is utilized to solve (23). The SL0 algorithm uses the smoothed $l_0$ norm to find sparse solutions and has been verified to perform even faster than $l_1$-norm-based algorithms [33]. The SL0 algorithm has been successfully used in time-frequency analysis and behaves well. Readers can refer to [34] for more details about the process of the SL0 algorithm for instantaneous spectrum reconstruction.

By combining the effective ASTK and SR-based spectrum reconstruction procedure, the desired high-performance TFD (named ASTK-SR TFD) can be derived. The proposed TFD differs from AOK TFD in three aspects. First, ASTK-SR TFD employs an adaptive window to decompose the signal into short-time segments, whereas the AOK TFD uses a fixed window. Secondly, ASTK-SR TFD devises an adaptive kernel to weigh the signal’s STAF, while the AOK TFD utilizes low-pass kernels with Gaussian radial cross sections. Lastly, ASTK-SR TFD utilizes SR techniques to construct the spectrum, whereas the AOK TFD relies on the Fourier transform. To conclude, the complete signal processing process of the proposed method is illustrated in Figure 2 and summarized as follows.

• Break the m-D signal to be analyzed into a series of short-time segments with the adaptive window.
• Compute the STAF of each short-time signal in both rectangular and polar coordinates.
• Accumulate the STAF energy along lines passing through the origin with different slopes to determine the shape of the ASTK.
• Apply the ASTK to the STAF to suppress unwanted artifacts.
• Take inverse FT to the kerneled STAF to obtain the artifact-free STIAF.
• Reconstruct the instantaneous spectrum by utilizing the sparsity of the STIAF slice in the corresponding Fourier dictionary.
• Obtain the high-performance TFD by arranging the spectrums of all time instants chronologically.

3. Experimental Results and Analysis

In this section, synthetic and real-life m-D signals are employed to demonstrate the superiority of the proposed ASTK-SR TFD in m-D feature enhancement. To provide a comprehensive comparison, WVD, which is one of the most well-known methods with excellent energy concentration, and RSPWVD, which is known to usually behave best in practical applications, are chosen as benchmarks in all comparisons. Moreover, since AOK TFD shares similar ideas and signal processing procedures with the proposed method, it is also included in the following experiments. Finally, the results of ASTK-FT TFD are also provided to demonstrate the benefit introduced by the SR technique.

3.1. Simulation

Tumbling refers to a type of complex micro-motion that consists of rotation, translation, and acceleration. Using the theoretical model deduced in [1], the radar return of a tumbling target composed of three ideal scatterers is computed. The simulation parameters are listed in

Table 1. Simulation Parameters.

Parameter	Value
Carrier frequency	10 GHz
Pulse repetition frequency (PRF)	3000 Hz
Target location	(1000 m, 5000 m, 5000 m)
Initial Euler angles	(20${}^{\circ }$, 20${}^{\circ }$, 20${}^{\circ }$)
Angular velocity	4$\pi$ rad/s
Initial velocity	5 m/s
Scatterer position (target local coordinate)	${\mathrm{P}}_1$ (0 m, 0 m, 0 m)
	${\mathrm{P}}_2$ (−0.5 m, 0.3 m, 0.4 m)
	${\mathrm{P}}_3$ (0.5 m, −0.3 m, −0.4 m)

.

The theoretical m-D feature of the simulated radar echo and five considered TFDs computed with a 0 dB signal-to-noise ratio (SNR) are depicted in Figure 3 for a visual comparison of their energy concentration and artifact suppression ability. For a more comprehensive comparison, the cross-sections of these TFDs at time = 0.5 s are also depicted in Figure 4. As the performance of AOK TFD relies on the kernel volume parameter [23] and RSPWVD is dependent on the choice of window lengths [30], experiments have been conducted with various parameters. The TFDs reported in Figure 3 have all been individually optimized to be the most visually appealing outcome for a fair and unbiased comparison.

As shown in Figure 3, WVD is seriously contaminated by artifacts and noise due to the lack of interference suppression ability, while the others significantly reduce artifacts. However, AOK TFD fails to provide a clear TF image due to the limited artifact suppression ability provided by the radially Gaussian kernel. Furthermore, the TF trajectories are heavily expanded as numerous AF data related to ATs have been deserted. By assigning the dispersed TF energy back to the center of gravity at each time instant, RSPWVD truly improves the energy concentration. On the other hand, it still encounters the same issue of residual artifacts as AOK TFD does. Moreover, RSPWVD is very sensitive to noise, and thus the obtained TF trajectories are slightly twisted. Due to the advanced capabilities of the proposed ASTK in attenuating artifacts and noise while accurately preserving ATs’ information, ASTK-FT TFD is able to distinctly reveal the target’s m-D characteristic. However, as evidenced by Figure 4, the TF trajectories are still widened due to the employment of the 2-D FT. By combining the strengths of SR and ASTK, the proposed ASTK-SR TFD effectively suppresses artifacts and maintains high TF concentration, which closely resembles the theoretical TFD and thus achieves the most accurate m-D feature expression among the methods under consideration.

To quantitatively measure the performance of different TFDs, the Pearson correlation coefficients between them and the theoretical TFD are computed as

$$\mathrm{Pcor}=\frac{\mathrm{vec}{\left(\widehat{\mathit{P}}\right)}^{\mathrm{T}}\mathrm{vec}\left({\mathit{P}}_{\mathrm{ideal}}\right)}{\parallel \mathrm{vec}\left(\widehat{\mathit{P}}\right){\parallel }_2{∥\mathrm{vec}\left({\mathit{P}}_{\mathrm{ideal}}\right)∥}_2}$$

(24)

where $\widehat{\mathit{P}}$ stands for the TFD generated by a certain method, $\mathrm{vec}\left(\widehat{\mathit{P}}\right)$ represents the vector form of $\widehat{\mathit{P}}$ with the mean value subtracted, ${(·)}^{\mathrm{T}}$ denotes the vector transpose operation, and ${\parallel ·\parallel }_2$ is the $l_2$ norm. The Pearson correlation can measure the similarity in shape between the computed TFDs and the theoretical TFD and a higher value means better shape similarity [32]. Conceptually, a value of 1 indicates an excellent match with the signal’s time-frequency feature, in other words, the most precise extraction of the target’s m-D feature.

The Pearson correlations between the results and the standard distribution are computed in

Table 2. Quantitative Quality Evaluation of Different TFDs.

TFD	WVD	AOK	RSPWVD	ASTK-FT	ASTK-SR
Pcor	0.1872	0.3209	0.3710	0.5191	0.6063

. From a closer inspection of the data, it is evident that the correlation coefficients for WVD, AOK TFD, and RSPWVD are all below 0.5. In contrast, notably higher values can be observed for ASTK-FT TFD and ASTK-SR TFD, which confirms the effectiveness of the ASTK and highlights the superiority of ASTK-SR TFD over the other methods.

Moreover, to explore the performance of these TFDs in noise tolerance, white Gaussian noise is added to the synthetic m-D signal to simulate different SNRs varying from −5 dB to 10 dB. At each SNR, 100 Monte-Carlo trials are conducted to obtain the average Pearson correlation coefficients of different methods. As shown in Figure 5, one can see that all TFDs experience a performance deterioration as the SNR decreases. However, the proposed ASTK-SR TFD and ASTK-FT TFD always have much higher correlation coefficients with the theoretical TFD. The two ASTK-based methods show better noise robustness because the noise is randomly distributed in the ambiguity domain, and thus can be effectively filtered by the well-designed ASTK. Taken together, the overall assessment is that the proposed ASTK-SR TFD behaves better in energy concentration, artifact suppression, precision, and noise robustness among all the considered TFDs.

3.2. Computational Complexity

In this subsection, the five considered TFDs are compared from a computational viewpoint.

Assume an N-point sampled m-D signal whose TFD is supposed to be computed over N frequency bins. As previously mentioned, the reconstruction of the proposed ASTK-SR TFD can be summarized into seven steps. From a computational perspective, the cost of steps one and seven can be omitted. For step two, assuming the average length of the short-time window is $\overline{L}$, the dimension of each STAF in rectangular coordinate is $\overline{L}\times \overline{L}$ and can be obtained by $\overline{L}$ FFT procedures of size $\overline{L}$. Therefore, the computation of N STAFs in rectangular coordinates requires $\mathcal{O}(N{\overline{L}}^{2}log\overline{L})$ operations. Bilinear interpolation is an effective method for calculating the STAF in polar coordinates. The proposed method sets the dimension of the STAF in polar coordinates to be $2\overline{L}\times 2\overline{L}$, which means that it takes $\mathcal{O}({\overline{L}}^2)$ operations to generate it from a rectangularly sampled STAF. Therefore, the cost of computing N STAFs in polar coordinates is $\mathcal{O}(N{\overline{L}}^2)$. To sum up, the total computational cost of step two is $\mathcal{O}(N{\overline{L}}^{2}log\overline{L})$. For step three, the energy accumulation in each direction is simply achieved by $2\overline{L}$ additions. Consequently, the cost of energy accumulation in $2\overline{L}$ angles is $\mathcal{O}({\overline{L}}^2)$, which means the complexity of the determination of N adaptive kernels is $\mathcal{O}(N{\overline{L}}^2)$. For step four, each rectangularly sampled STAF is multiplied by the corresponding adaptive kernel, requiring $\mathcal{O}(N{\overline{L}}^2)$ operations in total. The transformation from each STAF to STIAF involves $\overline{L}$ FFT procedures of size $\overline{L}$. Therefore, the computational cost of step five is $\mathcal{O}(N{\overline{L}}^{2}log\overline{L})$. It has been derived that the computational cost of the SL0 algorithm is $O\left(U^{2}N\right)$ for a $U\times N$ dictionary matrix [34]. As a result, the generation of N instantaneous spectrums in step six needs $\mathcal{O}(N^2{\overline{L}}^2)$ operations since the average size of the adaptive Fourier dictionaries is $\overline{L}\times N$. To be more intuitive, the computational cost of the main steps of the proposed method is listed in

Table 3. Computational Costs of the Main Steps of the Proposed Method.

Step	2	3	4	5	6
Computational cost	$\mathcal{O}(N{\overline{L}}^{2}log\overline{L})$	$\mathcal{O}(N{\overline{L}}^2)$	$\mathcal{O}(N{\overline{L}}^2)$	$\mathcal{O}(N{\overline{L}}^{2}log\overline{L})$	$\mathcal{O}(N^2{\overline{L}}^2)$

.

Summing all data in Table 3, the computational complexity of the reconstruction of ASTK-SR TFD can be derived as $\mathcal{O}(N^2{\overline{L}}^2)$. Since the ASTK-FT TFD differs from ASTK-SR TFD only by using FFT to construct the spectrum in step six and the N FFT procedures require $\mathcal{O}(N^{2}logN)$ operations, the computational cost of ASTK-FT TFD can be obtained as $\mathcal{O}(N{\overline{L}}^{2}log\overline{L})$.

After consulting relevant literature [23,30,34], the computational complexities of the five methods are listed in

Table 4. Computational Performance of Different Methods.

Method	WVD	AOK	RSPWVD	ASTK-FT	ASTK-SR
Computational cost	$\mathcal{O}(N^{2}logN)$	$\mathcal{O}(N{\overline{L}}^2)$	$\mathcal{O}(N^{2}logN)$	$\mathcal{O}(N{\overline{L}}^{2}log\overline{L})$	$\mathcal{O}(N^2{\overline{L}}^2)$
Execution time	0.025 s	14.343 s	0.507 s	18.422 s	37.516 s

. Furthermore, the execution time of each method for computing the TFD of the simulated m-D signal is also provided for verification. All computations are performed in a computer equipped with an AMD Ryzen 7 5800H CPU and 32 GB RAM.

From the data in Table 4, it is clear that the proposed ASTK-SR TFD requires a little more computational cost due to the numerous complex iterations required for solving the SR problem. In the future, we will conduct more in-depth studies of the SR theory and explore adjustments to existing sparse decomposition algorithms to enhance the efficiency of our method.

3.3. Real Data Test

The real micro-Doppler (m-D) signal of two turntable-driven rotating four-sided corner reflectors collected in an anechoic chamber is utilized to illustrate the potential of ASTK-SR TFD in practical applications. As shown in Figure 6, a radar with 220-GHz carrier frequency and 1000 Hz pulse repetition frequency is employed. Two corner reflectors are rotating with a speed of 11 r/min, and the rotation radii are set to be 0.24 m and 0.16 m, respectively.

The waveform of the m-D signal is displayed in Figure 7, where five kinds of TFDs are also presented. In Figure 8, the cross-sections of these TFDs at time = 3 s are displayed for comparison. It can be observed that WVD still suffers from serious interference, posing a significant challenge to the m-D feature extraction. Similar to the result of the previous simulation, AOK TFD effectively enhances the m-D characteristic but shows some residual artifacts and unsatisfactory energy concentration. As shown in Figure 8d, ASTK-FT TFD attenuates interference more thoroughly than AOK TFD due to the use of the ASTK. By further improving the energy concentration of ASTK-FT TFD with SR techniques, the ASTK-SR TFD realizes excellent energy concentration and artifact suppression, even as compared to that of RSPWVD, which is known to behave best in practice. In summary, the proposed method can yield a TF image of the highest quality, which is beneficial to accurate micro-motion parameter estimation and target identification.

It is worth noting that ASTK-SR TFD may also exhibit a loss of TF information when dealing with signals that contain components with varying energy levels. This occurs because the energy of certain weak components might be fainter than artifacts, thus being incorrectly omitted when designing the ASTK. Research to address this issue will be conducted in the near future.

4. Conclusions

In this paper, an effective TFD reconstruction method has been proposed to enhance the m-D feature of targets with micro-motions. Benefiting from the well-designed ASTK, unwanted artifacts and noise can be effectively suppressed. Then, a high-performance TFD is produced by utilizing SR techniques. Simulations demonstrate the superior performance and robustness of the proposed TFD when compared to four other TFDs, including the well-known AOK TFD and RSPWVD. Moreover, a real m-D signal processing application verifies the practicability of the proposed method. Our future work will focus on exploring more effective criteria for extracting weak components in the ambiguity domain and improving the computational efficiency of the proposed method.