A Novel Velocity Estimation and Compensation Approach for Moving Target in Stepped Frequency Radar Based on Fractional Fourier Transform

1. Introduction

Noncooperative target recognition has been a highly concerned issue in recent years, leading to modern radar performs target imaging and recognition tasks besides conventional tasks of detection and tracking [1,2,3]. High range resolution profile (HRRP) achieved by larger working bandwidth is considered as a promising technology for noncooperative target identification [4,5]. Within HRRP radar systems, targets are discerned as a constellation of multiple scattering centers, each occupying distinct range cells, a phenomenon characteristically described as range-spread targets. Such detailed HRRP signatures markedly enhance the radar system's capabilities in both recognition and tracking of targets [6,7,8].

HRRP radar can be realized by two approaches, instantaneous wideband radar and synthetic wideband radar [9]. Instantaneous wideband radar provides big bandwidth in each pulse. It is simple on theory but expensive in cost. synthetic wideband radar divides the big bandwidth into several narrow-band pulses, therefore, more adaptable to resource-constrained platform. One practical form of synthetic wideband is the stepped-frequency waveform (SFW) technique [10,11,12], which transmits a series of coherent pulse sequences whose carrier frequencies increasing or decreasing by step size. Compared to the instantaneous wideband radar, the SFW operates within narrow instantaneous bandwidth and works well under the condition of narrowband transmitter and receiver [13]. It is widely used in satellite radar, warhead radar and even automotive radar. Unfortunately, the SFW is quite sensitive to the radial motion of the target for its large time-bandwidth product [1,14]. Especially in the context of high-speed target detection, the range-doppler coupling of the SFW causes range cell shift and power divergence of the target echoes.

The range cell shift and power divergence derived from the additive linear phase term and quadratic phase term of the Doppler effect has been the most concerned issues of SFW radar. Numerous studies on the SFW have been presented in recent years. These studies mainly focus on waveform designed [4,15,16,17,18] and motion compensation [19,20,21,22,23,24]. In [4], a novel SFW with double pulse repetition interval (PRI) and a Doppler phase cancellation approach is designed. It does not need estimation of the velocity and is easy to be realized. Referring to the opinion in [4], methods by using multiple stepped-frequency pulse trains and the robust phase unwrapping theorem to estimate the range and velocity of a target were introduced in [7,25]. These methods can well solve the range cell shift of the SFW radar. However, they reduce the utilization of the radar signal for the fact two sets of sequent pulse trains are used to unwrap phase. In [26,27,28], random stepped frequency signal is used instead of the linear frequency stepped signal to solve the range-doppler coupling of the SFW. But these methods introduce the phase compensation problem which may degenerate the correlation peak of the signal [29,30]. Besides, all of these methods based on waveform design suffer from the pairing of scattering centers when there are more than two targets with different velocities [14,23].

On the other hand, motion compensation-based methods focus on estimating the velocity of the target with the quadratic phase term of the Doppler effect, and compensating the Doppler effect with the estimated result. In [31], an algorithm based on the maximum likelihood (ML) estimation is proposed to estimate the target velocity accurately. However, the planar search of the ML makes its computation complexity is large, causing it difficult to be implemented in real time. Velocity estimation based on cross-correlation Fast Fourier transform (FFT) and the least burst error function is invested in [32,33,34]. They use an estimation strategy of coarse to fine, which firstly use cross-correlation FFT to estimate the velocity roughly and then use the least burst error rule to estimate the accurate velocity. However, these methods may be struggled in scenario with multiple targets. Exploiting the fact that the Doppler signal of the moving targets is a chirp signal whose chirp rate is decided by the velocity of the target, the discrete Chirp-Fourier transform [35] is used to estimate the velocity of the targets by chirp-rate estimation of the Doppler frequency. However, it also needs a planar search in the frequency and chirp rate joint plane, resulting in a computational burden. The Hough Transform [36,37] and Lv's distribution [38,39] are also considered to solve the velocity estimation of the target for SFW radar. Both Hough transform and Lv’s distribution perform extremely high accuracy in estimation of the chirp rate. However, their computational burden is always mentioned in practical.

The fractional Fourier transform (FrFT) attracts a great deal of attention in recent years for its outstanding performance of chirp signal processing [40,41,42]. The fast discrete algorithm proposed by Ozaktas in [43] significantly reduces the computational complexity of the FrFT and makes it potential in signal processing. In this paper, we drew on the idea of estimating the velocity of the targets by chirp-rate estimation of the Doppler frequency based on FrFT and introduced a fast velocity estimation and compensation approach. The velocity estimation accuracy based on FrFT was quantitatively analyzed. An iterative algorithm based on golden section search (GSS), which we named iterative golden section search (IGSS), was proposed to alleviate the computational burden in the chirp rate domain. The proposed method showes a significant estimation accuracy and is suitable for the scenarios of multiple targets.

The main contributions are summarized as follows.

1)

We present a novel velocity estimation and compensation approach based on FrFT for SFW radar. It exploits the fact that the Doppler effect of the moving targets is a chirp signal whose chirp rate is decided by the velocity of the target. The The FrFT accurately determines this chirp rate by synthesizing the chirp signals in an optimal sequence. We have conducted a thorough quantitative analysis to validate the precision of our approach.

2)

We propose a rapid and precise iterative algorithm for the FrFT-based velocity estimation that minimizes computational demands without detracting from the estimation's accuracy. This innovative algorithm utilizes the golden section search (GSS) method, optimizing computational efficiency in pinpointing the FrFT's optimal order.

3)

To accommodate scenarios involving multiple targets, we develop an Improved Golden Section Search (IGSS). This enhancement introduces an additional looping mechanism to the GSS, enabling the IGSS to iteratively estimate velocities for individual scattering centers, thus refining the approach for more complex targeting situations.

The rest of the paper is organized as follows: Section 2 provides a detailed description of the beat signal of the SFW radar when detecting moving targets. Section 3 describes the framework of the proposed approach and derives the velocity accuracy that the FrFT can achieve. Section 4 introduces the velocity estimation algorithm proposed in this paper. Section 5 provides simulations to verify the proposed approach. Conclusions are formed in Section 6.

2. Principle of The SFCW Radar

Without loss of generality, we select continuous wave SFW radar as an example. Besides, the proposed approach is also suitable for stepped-frequency chirp waveform radar. According to the principle of the SFW radar, a frame of its transmitting signal can be represented as:

$$s_{tr}(t)=A_{tr}\exp [j2\pi (f_0+i\Delta f)(t-iT)],iT<t<(i+1)T,i=0,1\cdots I-1,$$

(1)

Where $j=\sqrt{-1}$ represents the imaginary number,$A_{tr}=\left|A_{tr}\right|\exp (j{\theta }_0)$ represents the complex amplitude of the transmitting signal, $f_0$ is the carrier frequency of the transmitting signal, $\Delta f$is the step bandwidth, $T$is the dwell time of each frequency bin, and $i$represents the index of the frequency bin. When there are K scatterers, with the distance$R_k$ and the velocity $v_k$, detected by the SFW radar, the target echoes received by radar can be denoted as:

$$s_{re}(t)={\displaystyle \sum _{k=1}^K{A}_{re}^k}\exp [\mathrm{j}2\pi (f_0+i\Delta f)(t-iT-\frac{2(R_k-v_{k}t)}{\mathrm{c}})],iT<t<(i+1)T,i=0,1\cdots I-1,$$

(2)

where $A_{\mathrm{re}}^k=\left|A_{tr}\right|\exp (j{\theta }_k)$ is the complex amplitude of the k-th scatterer and c represents the speed of electromagnetic wave. After mixing the receiving signal with the reference signal and sampling the output beat signal with a sampling rate of $f_s=1/T$, the obtained digital beat signal can be expressed as

$$\begin{array}{c}{s}_{beat}(i)={\displaystyle \sum _{k=1}^K{A}_{be}^k}\exp [j2\pi (f_0+i\Delta f)\frac{2(R_k-v_{k}i)}{\mathrm{c}}],iT<t<(i+1)T,i=0,1\cdots I-1\\ ={\displaystyle \sum _{k=1}^K{A}_{be}^k}\exp [j2\pi f_0\frac{2R_k}{\mathrm{c}}+j2\pi (\frac{2R_k\Delta f}{\mathrm{c}}-\frac{2v_k{f}_0}{\mathrm{c}})i-j2\pi \frac{2v_k\Delta f}{\mathrm{c}T}{i}^2]\end{array},$$

(3)

Where $A_{be}^k=(A_{tr}^{\ast }{A}_{re}^k)/2$ is the complex amplitude of the digital beat signal. According to equation (3), the beat signal of the SFW radar can be consider as chirp signal of $i$ with the initial frequency $f_{ini}=(2R_k\Delta f/\mathrm{c}-2v_k{f}_0/\mathrm{c})$ and the chirp rate ${\beta }_{beat}=4v_k\Delta f/(\mathrm{c}T)$. When the velocity of the target is low, the chirp rate ${\beta }_{beat}$ and the Doppler frequency $f_d=-2v_k{f}_0/\mathrm{c}$ can be ignored. Therefore, the range information of the target can be simply obtained through an inverse fast Fourier transform (iFFT). However, the initial frequency of the beat signal shifts seriously and the chirp rate increases significantly when the speed of the target is high. It results in serious cell shift and energy diffusion of the SFW radar. The iFFT of the beat signal when the target is located at $R_k=580m$ with different velocities is depicted in Figure 1. According to the simulation result, the peak value of the iFFT spectrum appears at 580 m when the velocity of the target is 10 m/s. However, the peak value of the iFFT spectrum appears at 520 m when the velocity of the target is 1000 m/s, shifting 60 m away from the true range cell. Besides, its value is almost half of that when the speed is slow. Therefore, motion compensation is a key problem in SFW radar when detecting high-speed targets.

3. Velocity Estimation Approach for SFW Radar by FRFT

From section 2, it can be concluded that the digital beat signal of the SFW radar is a chirp signal whose chirp rate is calculated as ${\beta }_{beat}=4v_k\Delta f/(\mathrm{c}T)$ when encountering a moving target. Therefore, the FrFT is considered to be used to estimate the velocity of the target. The block diagram of the velocity estimation approach proposed for SFW radar in this paper is shown in Figure 2. It takes advantage of the aggregation of the chirp signal in fractional Fourier domain with a specific transform order and estimates the velocity of the targets by searching for the optimal order. After velocity estimation, motion compensation and the iFFT is performed on the beat signal to ultimately obtain the actual range of the target. This approach translates the problem of velocity estimation into peak search of the FrFT in frequency chirp-rate jointly plane. However, the computation complexity of a planar search is high. Therefore, we adopted the GSS algorithm for the peak search in chirp rate domain to reduce the computational complexity of the proposed approach.

3.1. The FrFT Spectrum of the Beat Signal at Optimal Order

For the convenience of analysis, we first consider the case that there is only one ideal point target $k$ in the observation scene of the SFW radar. Then, the FrFT of the beat signal can be denoted as:

$$\begin{array}{cc}\hfill F_{\alpha }\{s_{beat}(i)\}(u)& =\frac{1}{I}{\displaystyle \sum _{i=1}^I{s}_{beat}(i)}{K}_{\alpha }(i,u)\hfill \\ & =\frac{A_{be}{A}_{\alpha }}{I}\exp (j2\pi f_0\frac{2R_k}{\mathrm{c}}+\mathrm{j}\pi u^2\cot \alpha ){\displaystyle \sum _{i=1}^I\exp [j2\pi (\frac{2R_k\Delta f}{\mathrm{c}}}-\hfill \\ & \frac{2v_k{f}_0}{\mathrm{c}}-u\csc \alpha )i-j2\pi (\frac{2v_k\Delta f}{\mathrm{c}T}-\frac{\cot \alpha }{2})i^2]\hfill \end{array},$$

(4)

In (4), $K_{\alpha }(i,u)$ is the FrFT Kernel, which can be expressed as

$$K_{\alpha }(i,u)=\{\begin{array}{c}{A}_{\alpha }\exp [j\pi (t^2\cot \alpha +u^2\cot \alpha -2tu\csc \alpha ),\alpha \ne n\pi \\ \Delta (t-u),\alpha =2n\pi \\ \Delta (t+u),\alpha =(2n+1)\pi \end{array},$$

(5)

$\alpha$ represents the transformation order of the FrFT, $A_{\alpha }=\sqrt{1-j\cot (\alpha )}$ represents the transformation gain of the FrFT and $\Delta (·)$ represents a unit impulse function. Shown in Figure 3, the FrFT can be considered as a rotation of the time frequency plane with the angle $\alpha$. When $\alpha =\pi /2$, the FrFT returns to traditional Fourier transform. When $\alpha ={\alpha }_{opt}$, the chirp signal would exhibit a delta function in the FrFT domain.

According to the property of FrFT, the optimal order of the FrFT in (4) is determined by ${\alpha }_{op}=\mathrm{arccot}(4v_k\Delta f/\mathrm{c}T)+\pi /2$ and the FrFT spectrum of the beat signal at ${\alpha }_{op}$ can be finally expressed as

$$\begin{array}{c}{F}_{{\alpha }_{op}}\{s_{beat}(i)\}(u)=\frac{A_{be}{A}_{{\alpha }_{op}}}{I}\exp (\mathrm{j}2\pi f_0\frac{2R_k}{\mathrm{c}}+\mathrm{j}\pi u^2\cot {\alpha }_{op}){\displaystyle \sum _{i=1}^I\exp [j2\pi (\frac{2R_k\Delta f}{\mathrm{c}}}-\\ \frac{2v_k{f}_0}{\mathrm{c}}-u\csc {\alpha }_{op})i]\end{array},$$

(6)

According to (6), the FrFT spectrum of the beat signal exhibits a delta function at the optimal order ${\alpha }_{op}=\mathrm{arccot}(4v_k\Delta f/(\mathrm{c}T))+\pi /2$. Therefore, the velocity estimation problem can be solved by the optimal order search of the FrFT. Due to the orthogonality of the FrFT matrice, the conclusion that${\Vert s_{beat}(i)\Vert }_2^2={\Vert F_{\alpha }\{s_{beat}(i)\}(u)\Vert }_2^2$ can be inferred. Therefore, the peak of the FrFT spectrum has the max power at the optimal order. Then, the optimal order search can be realized by max power search of the FrFT spectrum in chirp rate domain. Namely, the problem of optimal order search can be solved by the following optimization problem:

$$\widehat{{\alpha }_{op}}=\underset{\alpha }{\arg \max }\left\{\underset{u}{\max }\right|F_{\alpha }\{s_{beat}(i)\}(u){|}^2\},$$

(7)

After searching for the optimal order, the velocity estimation of the target can be calculated by

$$\widehat{v_k}=\frac{\mathrm{c}T\cot \widehat{{\alpha }_{op}}}{4\Delta f},$$

(8)

From equation (6), the range information of the target can be obtained with the frequency index of the FrFT spectrum at the optimal order. However, the discrete FrFT algorithm proposed by Ozaktas involves shifting and symmetry of the frequency, resulting in the estimation result of $u\csc {\alpha }_{op}$ corresponds to the center frequency of the beat signal rather than the initial frequency of the beat signal. Therefore, the proposed approach in Figure 2 estimates the range information of the target by motion compensation and iFFT, not the estimation of the initial frequency by FrFT.

3.2. Estimation Accuracy of the Fractional Fourier Transform

In this subsection, the estimation accuracy of the proposed approach is analyzed by the variance of the estimated velocity. According to equation (8), the estimation accuracy of the proposed approach depends on the estimation accuracy of the optimal order (or chirp rate). Given the noiseless case of a beat signal with a single ideal point target, $\left|F_{\alpha }\right\{s_{beat}(i)\}(u){|}^2$ attains its maximums at (α_op, u_op) in the frequency chirp rate jointly plane. Assuming $\left|F_{\alpha }\right\{s_{beat}(i)\}(u){|}^2$ is differentiable at its peak, then we can form the following equation.

$${\frac{\partial \left|F_{\alpha }\right\{s_{\mathrm{beat}}(i)\}(u){|}^2}{\partial \alpha }|}_{\begin{array}{c}{\alpha }_{\mathrm{op}}\\ u_{\mathrm{op}}\end{array}}=0,{\frac{\partial \left|F_{\alpha }\right\{s_{\mathrm{beat}}(i)\}(u){|}^2}{\partial u}|}_{\begin{array}{c}{\alpha }_{\mathrm{op}}\\ u_{\mathrm{op}}\end{array}}=0,$$

(9)

At this point, $\widehat{{\alpha }_{\mathrm{op}}}={\alpha }_{\mathrm{op}}$ is satisfied. Namely, the estimation result of the optimal order is equal to its true value. However, in practical applications, the beat signal output by the SFW radar is always accompanied by noise of varying intensities. When the noise is present, the FrFT spectrum is perturbed by the noise as

$$\begin{array}{l}\left|{\tilde{F}}_{\alpha }\right\{s_{\mathrm{beat}}(i)\}(u){|}^2=|F_{\alpha }\{s_{\mathrm{beat}}(i)\}(u)+F_{\alpha }\{w(i)\}(u){|}^2\\ =F_{\alpha }\{s_{\mathrm{beat}}(i)\}(u)F_{\alpha }^{\ast }\{s_{\mathrm{beat}}(i)\}(u)+F_{\alpha }\{s_{\mathrm{beat}}(i)\}(u)F_{\alpha }^{\ast }\{w(i)\}(u)+\\ F_{\alpha }\{w(i)\}(u)F_{\alpha }^{\ast }\{s_{\mathrm{beat}}(i)\}(u)+F_{\alpha }\{w(i)\}(u)F_{\alpha }^{\ast }\{w(i)\}(u)\end{array},$$

(10)

In equation (10), ${(•)}^{\ast }$ represents conjugate operation, $F_{\alpha }\{w(i)\}(u)$ is the FrFT spectrum of the noise, and $w(i)$ is the noise carried by beat signal. Usually, $F_{\alpha }\{w(i)\}(u)F_{\alpha }^{\ast }\{w(i)\}(u)$ can be ignored because the spectrum of the noise is expanded throughout the entire FrFT domain and is low-powered in comparison to the power of the signal. In this case, equation (9) can be approximately simplified as:

$$\left|{\tilde{F}}_{\alpha }\right\{s_{\mathrm{beat}}(i)\}(u){|}^2\approx |F_{\alpha }\{s_{\mathrm{beat}}(i)\}(u){|}^2+|\Delta {\tilde{F}}_{\alpha }(u){|}^2,$$

(11)

where $|\Delta {\tilde{F}}_{\alpha }(u){|}^2=F_{\alpha }\{s_{\mathrm{beat}}(i)\}(u)F_{\alpha }^{\ast }\{w(i)\}(u)+F_{\alpha }\{w(i)\}(u)F_{\alpha }^{\ast }\{s_{\mathrm{beat}}(i)\}(u)$ is the disturbance of the FrFT spectrum caused by noise. Assuming the disturbances caused by noise in frequency domain and chirp rate domain are $\Delta \alpha$ and $\Delta u$ respectively, the estimated results can be expressed as $\widehat{{\alpha }_{\mathrm{op}}}={\alpha }_{\mathrm{op}}+\Delta \alpha$ and $\widehat{u_{\mathrm{op}}}=u_{\mathrm{op}}+\Delta u$. Then, the first-order partial derivative of the FrFT spectral at its peak satisfies

$$\begin{array}{l}{\frac{\partial [|F_{\alpha }\{s_{\mathrm{beat}}(i)\}(u){|}^2+|\Delta {\tilde{F}}_{\alpha }(u){|}^2]}{\partial \alpha }|}_{\begin{array}{c}{\alpha }_{\mathrm{op}}+\Delta \alpha \\ u_{\mathrm{op}}+\Delta u\end{array}}=0\\ {\frac{\partial [|F_{\alpha }\{s_{\mathrm{beat}}(i)\}(u){|}^2+|\Delta {\tilde{F}}_{\alpha }(u){|}^2]}{\partial u}|}_{\begin{array}{c}{\alpha }_{\mathrm{op}}+\Delta \alpha \\ u_{\mathrm{op}}+\Delta u\end{array}}=0\end{array},$$

(12)

Ignoring the higher-order terms, the first-order Taylor series expansion of equation (11) can be expressed as

$$\begin{array}{l}{\frac{\partial \left|F_{\alpha }\right\{s_{\mathrm{beat}}(i)\}(u){|}^2}{\partial \alpha }|}_{\begin{array}{c}{\alpha }_{\mathrm{op}}\\ u_{\mathrm{op}}\end{array}}+{\frac{\partial |\Delta {\tilde{F}}_{\alpha }(u){|}^2}{\partial \alpha }|}_{\begin{array}{c}{\alpha }_{\mathrm{op}}\\ u_{\mathrm{op}}\end{array}}+\\ {\frac{{\partial }^2\left|F_{\alpha }\right\{s_{\mathrm{beat}}(i)\}(u){|}^2}{\partial {\alpha }^2}|}_{\begin{array}{c}{\alpha }_{\mathrm{op}}\\ u_{\mathrm{op}}\end{array}}\Delta \alpha +{\frac{{\partial }^2\left|F_{\alpha }\right\{s_{\mathrm{beat}}(i)\}(u){|}^2}{\partial \alpha \partial u}|}_{\begin{array}{c}{\alpha }_{\mathrm{op}}\\ u_{\mathrm{op}}\end{array}}\Delta u\approx 0\\ {\frac{\partial \left|F_{\alpha }\right\{s_{\mathrm{beat}}(i)\}(u){|}^2}{\partial u}|}_{\begin{array}{c}{\alpha }_{\mathrm{op}}\\ u_{\mathrm{op}}\end{array}}+{\frac{\partial |\Delta {\tilde{F}}_{\alpha }(u){|}^2}{\partial u}|}_{\begin{array}{c}{\alpha }_{\mathrm{op}}\\ u_{\mathrm{op}}\end{array}}+\\ {\frac{{\partial }^2\left|F_{\alpha }\right\{s_{\mathrm{beat}}(i)\}(u){|}^2}{\partial u^2}|}_{\begin{array}{c}{\alpha }_{\mathrm{op}}\\ u_{\mathrm{op}}\end{array}}\Delta u+{\frac{{\partial }^2\left|F_{\alpha }\right\{s_{\mathrm{beat}}(i)\}(u){|}^2}{\partial \alpha \partial u}|}_{\begin{array}{c}{\alpha }_{\mathrm{op}}\\ u_{\mathrm{op}}\end{array}}\Delta \alpha \approx 0\end{array},$$

(13)

Substitute equation (9) into equation (13), and make $s_{11}={\frac{{\partial }^2\left|F_{\alpha }\right\{s_{\mathrm{beat}}(i)\}(u){|}^2}{\partial {\alpha }^2}|}_{\begin{array}{c}{\alpha }_{\mathrm{op}}\\ u_{\mathrm{op}}\end{array}}$, $s_{12}=s_{21}={\frac{{\partial }^2\left|F_{\alpha }\right\{s_{\mathrm{beat}}(i)\}(u){|}^2}{\partial \alpha \partial u}|}_{\begin{array}{c}{\alpha }_{\mathrm{op}}\\ u_{\mathrm{op}}\end{array}}$, $s_{22}={\frac{{\partial }^2\left|F_{\alpha }\right\{s_{\mathrm{beat}}(i)\}(u){|}^2}{\partial u^2}|}_{\begin{array}{c}{\alpha }_{\mathrm{op}}\\ u_{\mathrm{op}}\end{array}}$, $s_1={\frac{\partial |\Delta {\tilde{F}}_{\alpha }(u){|}^2}{\partial \alpha }|}_{\begin{array}{c}{\alpha }_{\mathrm{op}}\\ u_{\mathrm{op}}\end{array}}$, $s_2={\frac{\partial |\Delta {\tilde{F}}_{\alpha }(u){|}^2}{\partial u}|}_{\begin{array}{c}{\alpha }_{\mathrm{op}}\\ u_{\mathrm{op}}\end{array}}$. Ultimately, equation (13) can be expressed as linear equations

$$\left[\begin{array}{c}\Delta \alpha \\ \Delta u\end{array}\right]=-{\left[\begin{array}{cc}{s}_{11}& s_{12}\\ s_{21}& s_{22}\end{array}\right]}^{-1}\left[\begin{array}{c}{s}_1\\ s_2\end{array}\right],$$

(14)

According to equation (14), the following equation can be obtained:

$$\Delta \alpha =\frac{s_1{s}_{22}-s_2{s}_{12}}{s_{11}{s}_{22}-s_{12}^2},$$

(15)

In equation (15), since $s_{11},s_{22},s_{12},s_{21}$ is only related to the beat signal itself and has nothing to do with the noise carried by the beat signal, they can be calculated directly

$$\begin{array}{l}{s}_{11}=-\frac{A_{be}^2{I}^6}{90}{\pi }^2{[{(\frac{4v_k\Delta f}{\mathrm{c}T})}^2+1]}^{\frac{5}{2}}\\ s_{22}=-\frac{A_{be}^2{I}^4}{90}{\pi }^2{[{(\frac{4v_k\Delta f}{\mathrm{c}T})}^2+1]}^{\frac{3}{2}}\\ s_{12}=-\frac{2A_{be}^2{I}^4}{3}{\pi }^2\frac{4v_k\Delta f}{\mathrm{c}T}(\frac{2R_k\Delta f}{\mathrm{c}}-\frac{2v_k{f}_0}{\mathrm{c}})[{(\frac{4v_k\Delta f}{\mathrm{c}T})}^2+1]\end{array},$$

(16)

The variance of the optimal order can be calculated by

$$\begin{array}{l}E\left\{{[\widehat{{\alpha }_{op}}-E(\widehat{{\alpha }_{op}})]}^2\right\}=E(\Delta {\alpha }^2)\\ =\frac{E\left\{{\left|s_1\right|}^2\right\}{\left|s_{22}\right|}^2+E\left\{{\left|s_2\right|}^2\right\}{\left|s_{12}\right|}^2}{{\left|s_{11}{s}_{22}-s_{12}^2\right|}^2}-\frac{E\left\{s_1{s}_2^{\ast }\right\}{s}_{22}{s}_{12}^{\ast }+E\left\{s_2{s}_1^{\ast }\right\}{s}_{12}{s}_{22}^{\ast }}{{\left|s_{11}{s}_{22}-s_{12}^2\right|}^2}\end{array}$$

(17)

As $E\left\{{\left|s_1\right|}^2\right\},E\left\{{\left|s_2\right|}^2\right\},E\left\{s_1^{\ast }{s}_2\right\},E\left\{s_2^{\ast }{s}_1\right\}$ can be calculated as $\begin{array}{c}E\left\{{\left|s_1\right|}^2\right\}=-\frac{A_{be}^2{I}^7{\sigma }^2}{90}{\pi }^2{[{(\frac{4v_k\Delta f}{\mathrm{c}T})}^2+1]}^3\\ E\left\{{\left|s_2\right|}^2\right\}=-\frac{2A_{be}^2{I}^5{\sigma }^2}{90}{\pi }^2{[{(\frac{4v_k\Delta f}{\mathrm{c}T})}^2+1]}^2\\ E\left\{s_1^{\ast }{s}_2\right\}=E\left\{s_2^{\ast }{s}_1\right\}=-\frac{2A_{be}^2{I}^5{\sigma }^2}{3}{\pi }^2(\frac{2R_k\Delta f}{\mathrm{c}}-\frac{2v_k{f}_0}{\mathrm{c}})\frac{4v_k\Delta f}{\mathrm{c}T}{[{(\frac{4v_k\Delta f}{\mathrm{c}T})}^2+1]}^{3/2}\end{array}$ (18)

Where ${\sigma }^2$ is the power of the noise. The variance of the optimal order can be finally expressed as

$$E(\Delta {\alpha }^2)=\frac{90{\sigma }^2}{A^2{[{(\frac{4v_k\Delta f}{\mathrm{c}T})}^2+1]}^2{\pi }^2{T}^4{I}^5}$$

(19)

Then, the variance of the chirp rate can be calculated as

$$\mathrm{var}\{\cot (\widehat{{\alpha }_{\mathrm{op}}})\}=\frac{90{\sigma }^2}{A^2{T}^4{I}^5},$$

(20)

According to equation (8) and equation (20), the variance of the velocity is calculated as:

$$\mathrm{var}(\widehat{v})=\frac{{\mathrm{c}}^{2}90{\sigma }^2}{16{(\Delta f)}^2{A}^2{T}^2{I}^5},$$

(21)

The detailed derivation of the $E\left\{{\left|s_1\right|}^2\right\},E\left\{{\left|s_2\right|}^2\right\},E\left\{s_1^{\ast }{s}_2\right\}$ and $E\left\{s_2^{\ast }{s}_1\right\}$ is shown in Appendix A.

4. Fast Estimation of The Chirp-Rate by Golden Section Search (GSS)

According to the principle of the proposed velocity estimation and compensation approach in section 3, its core is the estimation of the chirp rate of the beat signal. However, directly searching for the optimal order of the FrFT need traverse all values from $-\pi$ to $\pi$, which is large computational burden. Therefore, the GSS is introduced in searching for the optimal order to reduce computational complexity and improve the real-time performance of the SFW radar. Moreover, drawing lessons from the match pursuit algorithm, the IGSS algorithm is proposed in the paper to solve the problem of multi-scatterers.

4.1. Searching for optimal order by GSS

Some study has verified that the GSS shows well performance for the chirp rate estimation by FrFT [44]. The search process of the GSS is shown in Figure 4. The GSS starts by pre-determined search scope $[x_a,x_b]$ which is always initialized based on the prior knowledge. The gold ratio $r=(\sqrt{5}-1)/2$ is used to calculate the order of the first iteration point by

$$\begin{array}{l}{x}_1=x_b-r(x_b-x_a)\\ x_2=x_a+r(x_b-x_a)\end{array},$$

(22)

Assuming the objective function that requires searching for the maximum is denoted as $P(x)$, $P(x_1)$ and $P(x_2)$ are calculated to determine the next search scope. If $P(x_1)<P(x_2)$, the search scope is updated with $(x_1,x_b]$, and $x_3=x_1+r(x_b-x_1)$ is used to replace $x_1$ for the next judgment process. If $P(x_1)>P(x_2)$, the search scope is updated with $(x_a,x_2]$, and $x_3=x_a+r(x_2-x_a)$ is used to replace $x_2$. By continuously updating the search scope until the condition $\left|x_k-x_{k+1}\right|<\epsilon$ is satisfied, the index of the maximum can be determined as $x_{op}=\frac{x_k+x_{k+1}}{2}$.

Algorithm 1: iterative algorithm for estimating the velocity of the target based on GSS.

In this paper, the objective function can be considered as the peak of the FrFT spectrum$P(v)=\underset{u}{\max }{\left|F_{\alpha (v)}\{s_{be}(i)\}(u)\right|}^2$ with the argument $v$. $\alpha (v)=\frac{4v\Delta f}{\mathrm{c}T}$ is a function of $v$. The initial search scope $[v_a,v_b]$ of the GSS can be just preset with the possible speed of the targets. The pseudocode for the search algorithm is provided in Algorithm 1.

4.2. The Iterative GSS for Multi-scatterers

In the scenarios that there are multi-scatterers with different motion states, the beat signal output by the SFW radar is the superposition of several chirp signals with different initial frequencies and chirp rates. It can be expressed as

$$s_{be}(i)={\displaystyle \sum _{k=1}^K{A}_{be}^k}\exp [j2\pi f_0\frac{2R_k}{\mathrm{c}}+j2\pi (\frac{2R_k\Delta f}{\mathrm{c}}-\frac{2v_k{f}_0}{\mathrm{c}})i-\mathrm{j}2\pi \frac{2v_k\Delta f}{\mathrm{c}T}{i}^2],i=0,1\cdots I-1,$$

(23)

If the GSS mentioned in sub-section 4.1 is used to search for the optimal order, only the beat signal of the brightest scatterer can be detected. Therefore, an iterative GSS, which draws lessons from the match pursuit, is proposed in this sub-section to solve the problem of multi-scatterers. The iterative GSS first estimates the velocity $\widehat{v}$ and distance $\widehat{R}$ of the brightest scatterer $k$ by the ordinary GSS to estimate the velocity of the brightest scatterer. Then, the iFFT is used to calculate the range information of the scatterer $k$ after motion compensation.

Finally, the following equation

$$\widetilde{A_{be}^k}=\frac{1}{I}{\displaystyle \sum _{i=1}^I{s}_{beat}(i)}\exp [-j2\pi (\frac{2\widehat{R_k}\Delta f}{\mathrm{c}}-\frac{2\widehat{v_k}{f}_0}{\mathrm{c}})i+j2\pi \frac{2\widehat{v_k}\Delta f}{\mathrm{c}T}{i}^2]$$

(24)

is used to calculated the complex amplitude $\widetilde{A_{be}^k}$ of the beat signal corresponding to the k-th scatterer. After parameter estimation, the estimated parameters are used to reconstruct the beat signal corresponding to the k-th scatterer by$s_{beat}^k(i)=\widetilde{A_{be}^k}\exp [j2\pi (\frac{2\widehat{R_k}\Delta f}{\mathrm{c}}-\frac{2\widehat{v_k}{f}_0}{\mathrm{c}})i-j2\pi \frac{2\widehat{v_k}\Delta f}{\mathrm{c}T}{i}^2]$. The reconstructed signal is subtracted from the original beat signal $s_{be}(i)$, and the remained signal $s_{beat}^{remain}(i)=s_{beat}(i)-s_{beat}^k(i)$ is used as a new beat signal to estimate the parameters of the other scatterers. This process is repeated until the stopping criteria ${\Vert s_{\mathrm{beat}}(i)\Vert }_2^2<\mu$ or $m\le M$ is satisfied. Here, $\mu$ is a preset constant decided by the power of the noise, $m$ is the index of iterations, and $M$ is the total number of iterations. The pseudocode for the IGSS algorithm is provided in Algorithm 2.

Algorithm 2: Proposed Iterative GSS for Multi-scatterers.

4.3. Computational Cost of the IGSS

The computational cost of the proposed IGSS is decided by the number of iterations and the computational cost of the digital FrFT jointly. The fast digital FrFT proposed by Ozaktas et al includes interpolation, modulation, convolution. The core of this algorithm is the convolution of the chirp signal with the fractional kernel function, which is accelerated by the FFT. The total computational cost of the fast digital FrFT is

$$\frac{3(6I-5)}{2}{\log }_2(6I-5)$$

(25)

In GSS, the loop stops when $\left|v_a-v_b\right|<\epsilon$. According to [45], the number of iterations required by GSS is

$$Q=\lceil \frac{\log [\epsilon /(v_b-v_a)]}{\log r}\rceil ,$$

(26)

Where $\lceil •\rceil$ is the rounding function. According to (26), the number of iterations is mainly decided by the selected $\epsilon$. Therefore, an appropriate value of $\epsilon$ should be selected to decrease the total iterations based on the application. In (21), the theoretical variance of the velocity estimated by FrFT is given. To make the GSS converge and achieve this accuracy, the selected $\epsilon$ should satisfies that ${\epsilon }^2$ reaches an order of $O(I^{-5})$. To make the calculation simply, a convenient value of $\epsilon$ can be selected as

$$\epsilon =I^{-3},$$

(27)

With the $\epsilon$ selected by (27), the total number of iterations can be calculated as

$$Q=O(\lceil \frac{\log I}{\log r}\rceil )$$

(28)

Therefore, the computational complexity of the entire algorithm in the scenario with single scatter is approximately

$$\frac{(6I-5)}{2}{\log }_2(6I-5)O(\lceil \frac{\log I}{\log r}\rceil ),$$

(29)

and the computational complexity of the entire algorithm in the scenario with K scatters is approximately

$$K\frac{(6I-5)}{2}{\log }_2(6I-5)O(\lceil \frac{\log I}{\log r}\rceil ),$$

(29)

The computational costs of various motion compensation approaches are shown in Table 1. When compared to other approaches, the IGSS method undeniably offers an attractive computational cost.

5. Simulation Results

Extensive computer simulations have been conducted to verify the performance of the proposed approach. For the estimation accuracy of the velocity, 2000 times Monte Carlo simulations were conducted for each value of SNR. The standard deviation of the velocity has been calculated by

$$E=\sqrt{\frac{1}{2000}{\displaystyle {\sum }_{2000}{(v-\tilde{v})}^2}},$$

(25)

The SFW radar used in the simulation was a continuous wave radar, and the simulation parameters are shown in Table 2.

5.1. Mono-scatterer Case

The process of the GSS with a single scatterer at a speed of 1500m/s is shown in Figure 5. It can be concluded that the GSS gradually converges after 15 iterations under various SNR conditions. When SNR>0 dB, the estimation error of velocity after convergence is approximately 5.48 m/s; When SNR=-10 dB, the estimation error of velocity after convergence increases to 16.45 m/s. Namely, the estimation error of the GSS increases with the deterioration of SNR.

To give a detailed insight of the GSS, the standard deviation of the velocity calculated by equations (17) and (24) with different step points and stepped time-bandwidth products under SNR values ranging from -10 dB to 10 dB are shown in Figure 6. It can be concluded that the accuracy of the velocity is inversely proportional to both the number of step points I and the stepped time-bandwidth product. Therefore, we can improve the estimation accuracy of the GSS by increasing the number of steps, the dwell time of the stepped-frequency and the step bandwidth of the transmitting waveform. Besides, the figure makes it clear that the accuracy of the proposed GSS approach closely approximates the that of the theoretical value obtained from (17).

The comparison of the proposed approach with several different methods is shown in Figure 6. The methods, we selected for comparison include traditional velocimetry by Doppler frequency (without stepped frequency technology), the LVD-based methods [36], the Hough transform-based method [35], the ML-based method [28], and the cross-correlation FFT-based method [30]. The parameter of the radar is according to table Ⅱ (the step bandwidth of the traditional velocimetry is set to zero). From Figure 7, the proposed GSS approach has comparable performance to the ML-based method, and is superior to other methods.

5.2. Multi-scatterers Case

When there are two point-targets with velocities of 1500 m/s and 1000 m/s in the observation scenario, the convergence of the iterative GSS under SNR= 0dB condition is shown in Figure 8. It can be seen that the iterative GSS converges after 15 iterations when there are two point-targets with different velocities. The estimation accuracy of the iterative GSS under different SNR is shown in Figure 8(b). The estimation accuracy is almost equivalent to that of a single point-target. However, the estimation accuracy of target 2 is a slightly lower than that of target 1. This is because the target 1 is firstly estimated and compensated, and its residual influenced the estimation of target 2 by deteriorating the SNR. The HRRP image of a four-scatterers target with the velocity of 1000m/s obtained by iterative GSS is shown in Figure 9. The real range of the scatterers was obtained by iFFT with a velocity of 1000m/s, while the estimated range of the scatterers was obtained by the iterative GSS. From the results of the proposed iterative GSS, it can obtain the HRRP image of multi-scatterer targets.

The simulation model used in this paper is a SFW radar, whose detection range is limited by sampling rate (dwell time of stepped-frequency) of its transmitting waveform. In practice, we can combine the SFW with a chirp waveform (namely stepped-frequency chirp waveform radar) to improve the detection range and the SNR of the Doppler signal. The proposed iterative GSS is also suitable for the stepped-frequency chirp waveform radar with additional pulse compression before the velocity estimation. Here, we don’t provide the simulation of the stepped-frequency chirp waveform radar for its high similarity with SFW radar.

6. Conclusions

In this paper, we have put forward an accelerated method for estimating the velocity of moving targets using Stepped Frequency Waveform (SFW) radar, which harnesses the fractional Fourier transform (FrFT) complemented by an iterative golden section search (GSS) for the optimal FrFT order determination. Through a series of simulations, we have substantiated the robustness of our proposed methodology. The simulation outcomes indicate a marked improvement in velocity estimation precision over conventional FFT-based techniques. Moreover, our method approaches the estimative acuity of Maximum Likelihood (ML) methods, renowned for their supreme accuracy within Gaussian settings. Crucially, our approach achieves this while incurring a significantly lower computational cost, thereby presenting a more viable option for integration into radar systems. The iterative GSS algorithm we developed proves to be effective for scenarios involving multiple scatterers and succeeds in delivering a High Range Resolution Profile (HRRP) image for targets with up to four scatterers.