Parametric Sparse Recovery and SFMFT Based M-D Parameter Estimation with the Translational Component

Abstract

The micro-Doppler effect (m-D effect) provides unique signatures for target discrimination and recognition. In this paper, we consider a solution to the m-D parameter estimation. This method mainly consists of two procedures, with the first being the radar returns decomposition to extract the m-D components in Bessel domain. Then the parameter estimation issue is transformed as a parametric sparse recovery solution. A parametric sparse dictionary, which depends on m-D frequencies, is constructed according to the inherent property of the m-D returns. Considering that the m-D frequency is unknown, the discretizing m-D frequency range for the parametric dictionary matrix is calculated by the sinusoidal frequency modulated Fourier transform (SFMFT). In this manner, the finer m-D frequency, initial phases, maximum Doppler amplitudes and scattering coefficients are obtained by solving the sparse solution of the m-D returns. The simulation results verify the effectiveness.

This work was supported by the National Natural Science Foundation of China under Grant 61631019 and 61471386.

1 Introduction

Micro-Doppler effects (m-D effects) depict targets’ refined movement features, and provide the important information for target recognition and discrimination [1]. The rotation, vibration of a target induces the time-varying Doppler modulation on the received signals, which are most sinusoidal frequency-modulated (SFM) signals [2]. Hence, the parameters of the returned signals, including the scattering coefficients, modulation amplitudes, modulation frequencies and initial phases, suggest significant targets’ m-D features.

Time-frequency analysis (TF analysis) provides locations of such non-stationary signals in the TF plane, which has been widely utilized in m-D parameter estimation. However, most classic TF analyses either suffer from cross-term interferences, or are hardly to obtain the high resolution of time and frequency simultaneously [3]. Recently, some parametric methods which decompose signals into some different domain, such as the Bessel domain, are proposed and verified effective for the signal estimation and separation [4, 5]. Moreover, the sinusoidal frequency modulated Fourier transform (SFMFT) decomposes the phase term of a signal in the sinusoidal domain, so that the modulation frequency can be estimated by its projections on different sinusoids [6].

However, based on the Shannon-Nyquist theory, the received signals suffer from aliasing when the sampling frequency is less than twice of the maximum Doppler frequency shift, which makes it fail to estimate the motion features. In recent years, the compressive sensing (CS) technique is introduced to characterize the TF domain for m-D signatures extraction, and it is proved to be a promising tool with high-resolution [7,8,9,10]. CS theory supposes the sparsity of a signal in some transform domain, and the signal will be recovered from the limited samples. As the received signals induced by m-D are a sum of frequency modulated signals, they are sparse in TF domain [11]. As the approach proposed in [12], m-D signals can be first demodulated into some sinusoids, and then be reconstructed by exploiting their sparse solution in the Fourier domain. In [13], to extract accurate parameters, the demodulation procedure is conducted on the basis of a CS reconstruction process, but the initial phase cannot be obtained. In [14], the problem of m-D parameter estimation is transformed as a parametric sparse solution. It takes different m-D frequencies as the variables of the parametric dictionary matrix, and the sparse signal is solved by the pruned orthogonal matching pursuit (POMP) algorithm. Similarly, m-D signatures are extracted by parametric sparse recovery and OMP algorithm in passive radar systems [15]. In these methods, the m-D frequency is first being discretized as a frequency series. However, the estimation results are very sensitive to the m-D frequency, and meanwhile, the prior knowledge of the m-D frequency is unknown. In this case, the inaccurate division of m-D frequency series will lead to wrong estimation results. On the other hand, the wrong entries of the m-D frequency series are meaningless for estimation, which will lead to large computation amount. Additionally, some translational components induced by the bulk of non-rigid target or the scatters on the rotation axis, will make the traditional sparse recovery-based methods invalid for m-D parameter estimation.

To solve these problems, an approach based on the SFMFT in conjunction with the parametric sparse recovery and the k-resolution Fourier-Bessel (k-FB) series expansion is proposed. The received signal is first decomposed into different orders of k-FB series to extract the m-D signals of interest, and meanwhile, the phase shift ambiguity of the extracted m-D signals is revised. Considering the m-D frequency is unknown, it is first estimated coarsely by the SFMFT to provide a probable discretized frequency range for the parametric matrix. In this manner, the m-D frequency, as well as the initial phases and the modulation amplitudes are all discretized into series. The sparse representation model is solved by the OMP algorithm, and the finer modulation frequency, initial phases, modulation amplitudes and the scattering coefficients are obtained.

The following parts of the paper are organized as follows. In Sect. 2, the background of m-D signals and some basic theory is introduced. In Sect. 3 the proposed approach is described in details, while in Sect. 4 the simulation experimental results are analyzed. Section 5 draws the conclusion.

2 Background

In this section, some background including the m-D signal model and the basic theory of CS are briefly reviewed.

2.1 M-D Signal Model

The m-D effect is known as the time-varying Doppler modulations induced by the micro-motion of a target or its parts. The basic mathematical description of the m-D effect induced by rotation is discussed in this subsection.

The radar return of a target can be modeled as the signal reflected by a set of point scatters on the target. In this case, we suppose that there exist $K$ point scatters on the target. Suppose that the target rotates with the angular frequency $\omega$. The modulations induced by a rotation target consist of two parts: the m-D modulation induced by scatters rotating around the rotation axis; and the Doppler modulation induced by scatters on the rotation axis or the bulk of a non-rigid target. The received signal induced by rotation after the translation compensation is presented as

$$ s\left( t \right) = \sum\limits_{i = 1}^{K} {a_{i} \exp \left\{ {\text{j} \frac{4\pi }{\lambda }\left[ {f_{\text{D}} t + d_{i} \sin \left( {\omega t + \theta_{i} } \right)} \right]} \right\}} ,t = t_{1} , \cdots ,t_{M} $$

(1)

where $M$ denotes the signal sequence length, $\lambda$ is the radar wavelength and $f_{D}$ is the Doppler modulation frequency. $d_{i}$ is the distance between the rotation axis and the scatter $i$ in the LOS direction. $\varphi_{i}$ and $a_{i}$ are the initial phase and the scattering coefficient of scatter $i$, respectively. The Doppler frequency of scatter $i$ can be given according to the derivative of the phase term in [1] as:

$$ f_{i} \left( t \right) = \frac{2}{\lambda }\left[ {f_{\text{D}} + d_{i} \omega \,\sin \left( {\omega t + \theta_{i} } \right)} \right] $$

(2)

Note the fact that the rotation frequencies of all the scatters on one target are the same, with the rotation radius, initial phases and scattering coefficients are different. From (2) we can see that the translational component is different from that of the m-D component. The Doppler frequency of the translational component is constant, while the frequency of the m-D component is sinusoidal.

2.2 The CS Theory

The CS theory utilizes a signal’s sparsity, so that the discrete sampling can be obtained under the sampling rate much less than the Nyquist sampling condition. In accordance with the CS theory, a signal can be recovered from the limited samples if it is sparse in some domain. The classic CS theory is described as

$$ \varvec{s} = \varvec{\varPsi g} $$

(3)

where $\varvec{s}$ is an $M \times 1$ observation vector, $\varvec{\varPsi}$ is an $M \times N$ matrix, and $\varvec{g}$ is an $N \times 1$ sparse solution to be solved. The observation model in (3) can be solved under the constraint conditions, and thus, the sparse signal can be recovered from the observation vector by the following solution

$$ \hat{\varvec{g}} = \arg \mathop {\hbox{min} }\limits_{\varvec{g}} \left\| \varvec{g} \right\|_{0} {\text{ s}} . {\text{t}} .\left\| {\varvec{s} - \varvec{\varPsi g}} \right\|_{2}^{2} \, < \,\varepsilon $$

(4)

where $\left\| \cdot \right\|_{0}$ and $\left\| \cdot \right\|_{2}$ denote the $L_{0}$ norm and the $L_{2}$ norm, respectively, and $\varepsilon$ is the permitted error.

3 Parameter Estimation Method

As the modulation form of the translational component is different from that of the m-D component, the translational component will destroy the inherent property of the m-D signals, so that makes the sparse recovery based methods invalid. Hence, it requires an approach for the extraction of the m-D components from the received signals. Then the SFMFT and the parametric sparse recovery solution are adopted. The proposed parameter estimation method is introduced in the following subsections in detail.

3.1 Extraction of M-D Components

The discrete form $s\left( n \right)$ of the received signal in (1) with $N$ samples is presented as

$$ s\left( n \right) = \sum\limits_{i = 1}^{K} {a_{i} \exp \left[ {\text{j} \;pha\left( n \right)} \right]} , \, n = 1, \cdots ,N $$

(5)

where $pha\left( n \right) = \frac{4\pi }{\lambda }\left[ {f_{\text{D}} n + d_{i} \cos \left( {\omega n + \theta_{i} } \right)} \right]$ denotes the phase term.

As the discrete signals are mainly subjected to the problem of the phase shift ambiguity [16], the phase shift ambiguity should be revised for the received signal. Since the phase measurements are only possible in $\left( { - \pi ,\pi } \right)$, the phase shift ambiguity occurs when the total phase shift is more than $2\pi$ when sampling, which can be revised by comparison of adjacent samples.

$$ phar(n) = \left\{ {\begin{array}{*{20}l} {pha(n) - 2\pi ,} \hfill & {pha(n) - pha(n - 1) > \pi } \hfill \\ {pha(n) + 2\pi ,} \hfill & {phar(n - 1) - phar(n) > \pi } \hfill \\ {pha(n),} \hfill & {\left| {phar(n) - phar(n - 1)} \right| < \pi } \hfill \\ \end{array} } \right. $$

(6)

where $i > 1$ and $phar\left( n \right)$ denotes the difference of the real phase shift and the phase shift measurement.

After the phase shift ambiguity is revised, the m-D returns are extracted on the basis of the Fourier-Bessel transform (FBT) and its related theories. A signal can be decomposed into a weighted sum of Bessel functions by FBT [4]. Similar to the Fourier transform and the Fourier series, FBT presents as the FB series within a finite interval. For a better resolution, a resolution metric $k$ is introduced into the traditional FB series recently, which is called the k-FB series [5]. As the k-FB series has a one-to-one relationship between the order of the series and the frequency of a signal, it is effective for different signals separation overlapped in the TF domain.

Since the received signal $s\left( n \right)$ is in the neighborhood of zero-frequency after the translation compensation, the first several orders of the series correspond to the translational component. Expanded the signal by k-FB series and eliminate the first several orders of the series, the m-D signals $s_{\text{m - D}} \left( n \right)$ extracted from the radar returns can be reconstructed as follows

$$ s_{\text{m - D}} \left( n \right) = \sum\limits_{{m = m_{0} }}^{M} {\sum\limits_{n = 1}^{N} {C_{m} } J_{0} \left( {\frac{{\lambda_{m} }}{kN}n} \right)} $$

(7)

where $J_{0} \left( \cdot \right)$ is the Bessel function of the zero order, $\lambda_{m}$ is the $m$ th positive root of the function $J_{0} \left( t \right) = 0$, $m_{0}$ is the maximum order of k-FB series corresponding to the translational component, and $M$ is the order corresponding to the pulse repetition frequency (PRF). The series of order $m$ in (7) can be calculated as

$$ C_{m} = \frac{{2{\text{PRF}}^{2} }}{{\left[ {NJ_{1} \left( {\lambda_{m} } \right)} \right]^{2} }}\sum\limits_{n = 1}^{N} {ns\left( n \right)J_{0} \left( {\frac{{\lambda_{m} }}{kN}n} \right)} $$

(8)

3.2 Parameter Estimation

Considering that the parameters of the extracted m-D signals are $\left\{ {\omega ,a,d,\theta } \right\}$, discretize the candidates $\left\{ {d,\theta } \right\}$ in a certain range as $d\,{ \in }\,\left\{ {d_{1} ,d_{2} , \cdots ,d_{P} } \right\}$ and $\theta \,{ \in }\,\left\{ {\theta_{1} ,\theta_{2} , \cdots ,\theta_{Q} } \right\}$, respectively. Then the received m-D signal can be presented as

$$ \varvec{s} =\varvec{\varPsi}\left( \omega \right)\varvec{g}\left( \omega \right) $$

(9)

where $\varvec{s} = \left[ {s\left( 1 \right),s\left( 2 \right), \cdots ,s\left( N \right)} \right]^{\text{T} }$, $\varvec{g}\left( {\hat{\omega }} \right) = \left[ {g\left( 1 \right),g\left( 2 \right), \cdots ,g\left( {PQ} \right)} \right]^{\text{T} }$ and $\varvec{\varPsi}\left( \omega \right)$ is a $M \times PQ$ dictionary matrix. Then the sparse solution $\varvec{g}\left( \omega \right)$ can be obtained by

$$ \varvec{g}\left( \omega \right) = \mathop {\arg \hbox{min} }\limits_{{\varvec{g}\left( {\hat{\omega }} \right)}} \left\| {\varvec{g}\left( \omega \right)} \right\|_{0} {\text{ s}} . {\text{t}} . { }\left\| {\varvec{s} -\varvec{\varPsi}\left( \omega \right)\varvec{g}\left( \omega \right)} \right\|_{2}^{2} \, \le \,\varepsilon $$

(10)

$\varvec{\varPsi}\left( \omega \right)$ is constructed by the m-D signal’s intrinsic property

$$ \varvec{\varPsi}\left( \omega \right) = \left[ {\varphi \left( {1,\omega } \right), \cdots ,\varphi \left( {n,\omega } \right), \cdots ,\varphi \left( {N,\omega } \right)} \right]^{\text{T} } $$

(11)

where

$$ \varphi \left( {n,\omega } \right) = \left[ {\varphi_{0} \left( {n,\omega ,d_{1} ,\theta_{1} } \right), \cdots ,\varphi_{0} \left( {n,\omega ,d_{1} ,\theta_{P} } \right),\varphi_{0} \left( {n,\omega ,d_{2} ,\theta_{1} } \right), \cdots ,\varphi_{0} \left( {n,\omega ,d_{Q} ,\theta_{P} } \right)} \right] $$

(12)

and

$$ \varphi_{0} \left( {n,\omega ,d_{p} ,\theta_{q} } \right) = \exp \left[ {j\frac{4\pi }{\lambda }d_{p} \cos \left( {\omega n + \theta_{q} } \right)} \right] $$

(13)

As the error of the sparse recovery solution is more sensitive to $\omega$ than $\left\{ {d,\theta } \right\}$ in (9), serious recovery errors will be resulted in inaccurate or wrong $\omega$ [14]. However, without the prior knowledge of $\omega$, the discrete values of $\omega$ will not be fine enough, so that inaccurate even wrong estimation results will be obtained. Meanwhile, to avoid the unnecessary computation amount of wrong $\omega$ in the sparse recovery process, the SFMFT is first utilized to obtain a coarse estimation of $\omega$, which provides the candidate $\omega$ for a correct division range. In this manner, more accurate results of $\left\{ {\omega ,a,d,\theta } \right\}$ by the sparse recovery technique will be obtained.

SFMFT decomposes the modulation frequency $\omega$ on different sinusoidal basis, and the spectrum is obtained by the projections, referred to as the coefficient $e_{i}$ on the different sinusoidal basis. The discrete SFMFT [6] of the signal $s_{\text{m - D}} \left( n \right)$ is presented as

$$ s_{\text{m - D}} \left( n \right) = \left\langle {\sum\nolimits_{i = 0}^{N - 1} {e_{i} \left\langle \times \right\rangle \,\exp \left[ {\text{j} \,\exp \left( {\text{j} \omega_{0} ni} \right)} \right]} } \right\rangle $$

(14)

where the coefficient $e_{i}$ indicates the modulation spectrum of the signal, and $\omega_{0} = {{2\pi } \mathord{\left/ {\vphantom {{2\pi } N}} \right. \kern-0pt} N}$ is the modulation frequency unit. Coefficient $e_{i}$ can be calculated as

$$ e_{i} = \left\langle {\left\langle {s_{\text{m - D}} \left( n \right),\exp \left[ {\text{j} \exp \left( {\text{j} \omega_{0} ni} \right)} \right]} \right\rangle } \right\rangle = \frac{1}{N}\sum\limits_{n = 0}^{N - 1} {\ln \left\{ {x\left( n \right)} \right\}} \cdot \left[ { - \text{j} \exp \left( { - \text{j} \omega_{0} ni} \right)} \right] $$

(15)

The frequency estimation of $\omega$ is obtained by $e_{i}$ of the maximum value. Suppose that the maximum coefficient $e_{i}$ locates at $i = i_{0}$. Then the estimation of the rotation angular frequency is $i_{0} \omega_{0}$, and the discrete range of the candidate $\omega$ is set as $\left[ {\left( {i_{0} - 1} \right)\omega_{0} ,\left( {i_{0} + 1} \right)\omega_{0} } \right]$.

One of the solutions of (10) is the OMP algorithm [17]. It makes the dictionary matrix orthogonal by the Schmidt orthogonalization method. And then the process of the signal decomposition is iterated on the over-complete orthogonal basis. In each round of iteration, the most matched atom is found as the sparse approximation of the signal. The OMP algorithm is widely used in sparse signal processing because of its advances in the decomposition efficiency. The OMP algorithm is selected to solve (10). For clarity, the m-D parameter estimation algorithm is stated as follows.

4 Simulation Experiments

Suppose that the PRF of the radar is 1024 Hz, and with 1024 samples being collected. Considering that a target rotates with the angular velocity $\omega = 38.96\;{\text{rad/s}}$ with three scatters on it. One scatter locates at the rotation axis and thus induces a translational component. The rotation radiuses of other two scatters are 4 mm and 44 mm, the initial phases are 2.09 rad and 4.18 rad, and the scattering coefficients are 0.7 and 1, respectively. The received signal of the target under SNR = 12 dB is presented in Fig. 1(a).

Fig. 1.

TF analysis of the received signals. (a) The original received signal, (b) The extracted m-D signal.

From Fig. 1(a) we can see that the scatter with small rotation radius is almost contaminated by the translational component. The TF analysis of extracted the m-D signal is shown in Fig. 1(b), where the two m-D components can be both recognized. The m-D signals are extracted when the order is selected as $m_{0} = 10$ and $k = 2$ in (7). After the m-D component extraction, the spectrum of the m-D frequency calculated by SFMFT is shown in Fig. 2.

Fig. 2.

The SFMFT spectrum of the m-D frequency.

It is shown in Fig. 2 that the extracted m-D components only contain a single modulation frequency at around $37.7\,{{\text{rad}} \mathord{\left/ {\vphantom {{\text{rad}} {\text{s}}}} \right. \kern-0pt} {\text{s}}}$. The candidate $\omega$ is divided into the series of 40 entries between $\left[ {31.42,43.98} \right]$, and the candidates $\left\{ {d,\theta } \right\}$ of the dictionary matrix are divided into $25 \times 30$. After the iteration, the estimation results of $\left\{ {d,\theta } \right\}$ and $a$ are shown in Fig. 3.

Fig. 3.

The estimation result of $\left\{ {d,\theta } \right\}$ and $a$. (a) The extracted m-D signal, (b) The original received signal.

As denoted in Fig. 3(a), the results of the two groups of the estimated parameters are $\left\{ {d_{1} ,\theta_{1} } \right\} = \left\{ {4\;{\text{mm}},2.094\;{\text{rad}}} \right\}$ and $\left\{ {d_{2} ,\theta_{2} } \right\} = \left\{ {44\;{\text{mm}},4.189\;{\text{rad}}} \right\}$. The estimated scattering coefficients are $\left\{ {a_{1} ,a_{2} } \right\} = \left\{ {0.7049,1.006} \right\}$, which reflect on the magnitude of the pinnacles. The estimated m-D angular frequency is $38.96\;{{\text{rad}} \mathord{\left/ {\vphantom {{\text{rad}} {\text{s}}}} \right. \kern-0pt} {\text{s}}}$. The estimation results all agree well with the true values. On contrast, in Fig. 3(b), there are some interference pinnacles in the estimation results of the original received signal. The interference pinnacles are most in the location where the modulation amplitude approaches to zero. It agrees with (1) that the scatter on the rotation axis whose rotation radius is $d = 0$.

5 Conclusion

A novel approach of the m-D parameter estimation is proposed. We consider the m-D received signals containing both the m-D components and the translational components, which are invalid for the traditional sparse recovery based methods. The m-D components are first extracted in the procedure of signal decomposition and reconstruction by the k-FB series expansion. Then for finer estimation results, the SFMFT is adopted to provide the discretized modulation frequency range for parametric sparse recovery. By these means, interference pinnacles exist in traditional sparse recovery based methods are eliminated, and the accurate estimation results are obtained simultaneously.