2D-Unitary ESPRIT Based Multi-Target Joint Range and Velocity Estimation Algorithm for FMCW Radar

: Millimeter-wave FMCW radar has been widely used in joint range-velocity estimation of multiple targets. However, most existing algorithms are unable to estimate the range-velocity information with high accuracy simultaneously and fail to discriminate the targets with either closely spaced ranges or closely spaced velocities in the 2D range-Doppler spectrum. In order to deal with these problems, this paper proposes a 2D-Unitary ESPRIT-based joint range and velocity estimation algorithm of multiple targets for FMCW radar. Firstly, The 1D-IF signal is constructed into a 2D virtual array signal, the virtual array signals are preprocessed by a 2D-spatial smoothing technique to generate a new matrix signal. Then, according to the 2D-Unitary ESPRIT algorithm, the 2D real-valued information of the target parameters is obtained from this matrix signal, and then a new complex-value matrix is constructed. Finally, the eigenvalue decomposition of this new complex-value matrix is performed, and the range-velocity estimates of multiple targets are, respectively, calculated from the real and imaginary parts of the eigenvalues, and paired automatically. The simulation results illustrate that the proposed algorithm not only provides highly accurate range-velocity estimates but also has high-resolution performance and achieves automatic pairing of the range-velocity estimates in multi-target scenarios, thus effectively improving the multi-target joint range and velocity estimation performance of FMCW radar.


Introduction
The millimeter-wave radar, known for its high resolution, finds widespread applications in various domains including autonomous driving [1,2] and small target detection [3].The radar system exhibits exceptional detection capabilities and a high probability of target recognition, enabling accurate determination of the target's nature.The millimeter-wave radar is mainly divided into pulsed radar and continuous wave radar according to different ways of radiating electromagnetic waves [4].The pulsed radar is commonly used for range estimation.The constant frequency continuous wave radar is generally adopted for velocity estimation.The frequency shift keying (FSK) [5] or the phase shift keying (PSK) continuous wave radar achieves range-velocity estimation of multi-targets simultaneously, but the estimation is not accurate when the targets are relatively stationary.Frequency Modulated Continuous Wave (FMCW) radar [6,7], which utilizes frequency modulation such as triangle wave and sawtooth wave, can provide range and velocity estimates of multiple targets simultaneously.
The FMCW radar transmits a chirp signal and receives the signals reflected by multiple moving targets.The information about the time delays and Doppler shifts due to targets' range and radial velocity are embedded in the difference frequency between the transmitted and received signals, and can be estimated by utilizing the signal processing techniques.Therefore, highly accurate frequency estimation of intermediate frequency (IF) signal is a critical step to obtain the range-velocity information of the moving targets.The classical 2D Fast Fourier Transform (2D-FFT) algorithm [8,9] can perform joint range-Doppler estimation, but the estimation accuracy is inefficient due to its fence effect and spectrum leakage.There are existing algorithms to improve the accuracy of one-dimension (1D) frequency estimation, such as the Rife algorithm [10], the modified Rife algorithm [11], and the Chirp Z-Transform (CZT) algorithm [12], which refine the frequency band to improve the frequency resolution.The Zoom-FFT algorithm [13] achieves a higher frequency estimation accuracy by shifting the spectral band and decreasing the sampling rate.However, the frequency estimation accuracy of the methods in [8][9][10][11][12][13] are limited by the Rayleigh limit [14,15].Super-resolution frequency estimation is achieved by utilizing the subspace techniques of the conventional Multiple Signal Classification (MUSIC) algorithm [16] and Estimating Signal Parameter via Rotational Invariance Techniques (ESPRIT) algorithm [17], which overcome the Rayleigh limit.The orthogonality of the noise subspace and signal subspace in the signal autocorrelation matrix is used by the MUSIC algorithm to construct a spatial-spectral function for spectral peak searching to obtain high-resolution spectral estimation.The rotational invariance of the signal subspace is employed by the ESPRIT algorithm to obtain higher accuracy without searching for the spectral peaks.However, the subspace-based methods [16,17] suffer from large computational effort and poor estimation performance in noisy environments.Meanwhile, there exists a mismatch pairing of the range-velocity estimates when there are multiple targets.Therefore, researchers have proposed algorithms to improve the accuracy of two-dimension (2D) frequency estimation [18][19][20][21][22][23][24][25][26].The 2D frequencies are estimated by the 2D-CZT algorithm [18] and the pairing of range-velocity estimates is solved by using the correlation of signals [19].The 2D-MUSIC algorithm [20] achieves 2D frequency estimation and automatic pairing by constructing a linear epidemic matrix of 2D frequencies.The 2D-RootMUSIC algorithm [21] performs two 1D root estimates of frequencies and then performs range-velocity pairing by using the noise subspace orthogonal to the range-velocity steering vector.The 2D-ESPRIT algorithm [22,23] derives the 2D frequency parameters by utilizing the rotation invariance and translation invariance of the 2D signal matrix.The Modified Matrix Enhancement and Matrix Pencil (MMEMP) algorithm [24,25] is an improvement of the Matrix Enhancement and Matrix Pencil (MEMP) algorithm [26], which allows joint estimation of 2D frequencies.However, the above methods [18][19][20][21][22][23][24][25][26] suffer from high computational complexity, limited accuracy of pairing, and low resolution in range-velocity estimation.
In order to deal with the problems described above, this paper proposes a 2D-Unitary ESPRIT-based joint range and velocity estimation algorithm for multiple targets.First, the IF signal is constructed into a 2D virtual array signal matrix and preprocessed by the 2D-spatial smoothing technique [27,28].Then, the preprocessed matrix signal is converted into a real signal matrix, and two real-valued matrixes associated with the range and velocity are extracted using the rotation invariance of this real signal subspace.Next, a new complex matrix is constructed by using the above two real-valued matrixes as its real and imaginary parts.Finally, the 2D-Unitary ESPRIT algorithm [29] is utilized on this complex matrix to obtain automatic pairing of the range-velocity estimates of multiple targets.
The main contributions of this article are as follows.The 2D-Unitary ESPRIT algorithm [29,30] is commonly utilized for 2D direction of arrival (DOA) estimation in array signals.In this paper, we employ this algorithm for range and velocity estimation in FMCW radar.The 1D IF frequency signal is constructed into a 2D virtual array signal, this virtual array signal is reconstructed as a new complex matrix using a 2D-spatial smoothing technique to apply the 2D-Unitary ESPRIT algorithm for 2D frequency estimation.It is worth noting that this part of the method is also applicable to the 2D-RootMUSIC algorithm.
The rest of this paper is organized as follows.In Section 2, the signal model of the FMCW radar is described and the range and velocity estimation principle is derived.In Section 3, the 2D-Unitary ESPRIT-based joint range and velocity estimation algorithm of multiple targets is proposed.In Section 4, simulation results are conducted to evaluate the performance of the proposed algorithm.Finally, conclusions are presented in Section 5.

Signal Model
The transmitted signal of the FMCW radar [31] can be expressed as: where a t , f 0 , µ, ϕ 0 and T c denote the complex amplitude, carrier frequency, frequency slope, initial phase and ramp duration of the transmitted signal, respectively.The received signal is reflected by D targets and can be expressed as follows: where a r(i) , τ i and ϕ 1i denote the strength, time delay and phase of the ith reflected signal, respectively, and n z (t) is the noise component.By mixing s TX (t) with s RX (t) and passing it through a low-pass filter (LPF), the analog IF signal of D targets is expressed as:

Principle of the Range and Velocity Estimation
Consider there are D targets moving at radial velocity v 0i (i = 1, 2, • • • , D) relative to the radar.For FMCW radar signals, the radial velocity of the target relative to the radar is positive when the target is far away from the radar.This is because a positive Doppler shift is observed when the target is moving away from the radar.Ignoring the target movement within a chirp period, the time delay corresponding to the ith target can be expressed as: where R 0i is the range between the ith target and the radar, c is the velocity of the electromagnetic wave in the air, T l is the chirp duration with T l = T c + T z , 1 m M, and M is the number of chirps.Substituting (4) into (3), ignoring the 1/c 2 term, the phase of the ith target can be expressed as: where τ 0i = 2R 0i /c , λ = c/ f 0 is the wavelength of the transmitted signal.
Employing the central linear expansion formula, i.e., t = T l /2 , (m − 1)T l = (M − 1)T l /2, the third term in (5) can be expressed as follows: Substituting ( 6) into (5), where the T 2 l term can be neglected, the phase in ( 5) is expressed as: where are frequencies related to the range and velocity of the ith target.The information of τ 0i and v 0i is simultaneously contained in f ri , while the information of v 0i is only contained in f vi .Therefore, after estimating f ri and f vi , τ 0i in f ri can be calculated from v 0i in f vi and in turn R 0i is obtained.
When the targets move slowly, i.e., c v 0i , f ri and f vi can be, respectively, simplified as follows: Therefore, the discrete signal can be expressed as: where t s = 1/ f s , f s is the sampling frequency, N = f s T c is the number of samples.When a moving target approaches the radar, the schematic diagram of the time delay and Doppler shift generated by the moving targets is shown in Figure 1.Therefore, the range and velocity of the ith target can be estimated as:

2D-Unitary ESPRIT Based Joint Range and Velocity Estimation Algorithm of Multiple Targets
In this section, a 2D virtual array is constructed by using a 2D-spatial smoothing preprocessing technique.Then, the 2D-Unitary ESPRIT-based multi-target joint range and velocity estimation algorithm is presented.

2D Virtual Array Construction Using 2D-Spatial Smoothing Preprtocessing
The signal in (10) is rewritten into a 2D matrix of M × N as follows: The angular frequency in the fast time domain and slow time domain are, respectively, defined as: Therefore, x(m, n) in ( 10) is rewritten as x m,n (w, u) for simplifying the analysis as follows: A smoothing submatrix of M s × N s is selected from the matrix X(m, n), as shown in Figure 2a.M s and N s are the window lengths along the slow and fast time domains, respectively, and the number of slides is m t = M − M s + 1 and n t = N − N s + 1, respectively.The submatrix of each smoothing is converted into a vector of M s N s × 1 shown in Figure 2b, which can be expressed as: The variable is expressed in vector form as: where A is the M s N s × D steering vector matrix as A = [a 1 , a 2 , . . .a D ] and a i is the M s N s × 1 steering vector as Finally, a matrix of M s N s × m t n t after smoothing is obtained as follows: x 0,N s . . .
which can be expressed as: where S is D × D signal diagonal matrix.C is m t n t × D steering matrix expressed as . , e ju i (n t −1) ] T .Therefore, the matrix Z has shift invariance, and the 2D-Unitary ESPRIT algorithm [32] can be employed to obtain the range-velocity estimates of multiple targets.

2D-Unitary ESPRIT Algorithm for Joint Range and Velocity Estimation
The 2D-Unitary ESPRIT algorithm is an extension of the 1D-Unitary ESPRIT algorithm [33] and is commonly used for 2D DOA estimation.
The steering vector of the angular frequency u of the matrix Z is written in conjugate centrosymmetric form when N s is even: Similarly, the steering vectors a M s (w) of the angular frequency w is obtained.Thus, the 2D-frequency vector of the signal matrix Z is written as follows: The real-valued realizations are obtained by using a bijective mapping between centro-Hermitian matrices and real matrices, where a complex matrix . The unitary matrix: where I K is the unit matrix and ∏ N is the exchange matrix with ones on its antidiagonal and zeros elsewhere.Therefore, the complex matrix A(w, u) is converted into a real matrix D(w, u) as follows: where are real-valued frequency vector.d N s (u) satisfies rotation invariance as follows: where J 1 and J 2 are selection matrices of (N s − 1) × N s written as J 1 = Π N s −1 0 and * .Thus, K 1 and K 2 denote the real and imaginary parts of Q H N s −1 J 2 Q N s as follows: substituting ( 26) and ( 27) into ( 25), ( 25) can be rewritten as: Thus, the relation is satisfied by d N s (u) as follows: Similarly, the above relations are satisfied by d M s (w) through replacing u with w and N s with M s .The vec operator is used to the real matrix D(w, u), i.e., d(w , u) = vec(D(w , u)), then the steering vector of M s N s × 1 satisfies the following relation: where K u1 and K u2 are M s (N s − 1) × M s N s matrices, written as follows: The relation in (30) is also satisfied by angular frequency w, so that the matrices K w1 and K w2 of (M s − 1)N s × M s N s are written as follows: where K 3 , K 4 are similar to those defined in ( 26) and ( 27), and N s is replaced by M s .Thus, the M s N s × D real steering matrix is denoted as in the same way as in (30), D satisfies as follows: where, Similarly, Ω w is obtained by simply replacing u with w.
The complex matrix in ( 20) is pre-multiplied by(Q H N s ⊗ Q H M s ) and then obtains as: Thus, the covariance matrix of YY = [Re{Y}, Im{Y}] = DS + N is written as: where S and N are signal and noise parts of the matrix YY, respectively.R s = E[SS H ] is the covariance matrix of S and σ 2 is the variance of N. The R YY is decomposed by the Eigenvalue Decomposition (EVD), and the eigenvectors corresponding to the D largest eigenval-ues are formed by the signal subspace E s , which is a matrix of M s N s × D. When the number of snapshots tends to infinity, we obtain E s = DT, where T is an unknown real matrix of D × D. Substituting D = E S T −1 into (35), the signal eigenvector relationship of angular frequency w i and u i are obtained as follows: Utilizing ( 39) and ( 40), the complex square matrix Ψ u + jΨ w is constructed.Then, the w i and u i are estimated from the real and imaginary parts of the eigenvalues of the matrix Ψ u + jΨ w , respectively.Then, f ri and f vi can be estimated according to ( 14) and ( 15), respectively.Finally, the range Ri and velocity vi can be estimated according to (11) and ( 12), respectively.Obviously, the pairing of the range and velocity is achieved automatically.The proposed algorithm is summarized in Algorithm 1. 1 Using 2D-spatial smoothing preprocessing, the signal matrix X ∈ R M×N in ( 13) is reconstructed into a matrix Z ∈ R M s N s ×m t n t in (20); 2 Compute the signal matrix Y in (37), calculate the covariance matrix R YY in (38), and perform EVD of R YY to obtain the signal subspace 39) and (40); 5 Calculate the largest D eigenvalues λ i of the complex matrix Ψ w + jΨ u , where i = 1, 2, • • • , D; 6 Calculate the range-Doppler angular frequency estimates: The paired f ri and f vi are calculated from w i and u i according to ( 14) and ( 15); 8 Ri and vi are calculated using (11) and (12).

Performance Metrics and Parameter Settings
Simulations show the effectiveness of the proposed algorithm and compare it with the existing algorithms including the 2D-FFT algorithm [8], the 2D-CZT algorithm [18], the 2D-RootMUSIC algorithm [21] and the MMEMP algorithm [25].The parameter settings are indicated in Table 1.In simulations, the additive random Gaussian white noise is considered, and the Signal to Noise Ratio (SNR) varies in the range of [−8 dB, 16 dB].All results are based on NU M = 500 independent Monte Carlo experiments.Root Mean Square Error (RMSE) is employed as a performance metric to measure the deviation between the estimated and true values.RMSE R of range estimates is expressed as follows: Similarly, RMSE v of velocity estimates is obtained.In simulations, two scenarios are considered: single-target scenario and multi-target scenario.The single-target scenario is utilized to compare the performance of the proposed algorithm with the existing algorithms in terms of range and velocity estimation accuracy, and the multi-target scenario is considered to compare the performance of the proposed algorithm with the existing algorithms in terms of resolution ability and pairing correctness, respectively.

Comparison of Estimation Accuracy in the Single-Target Scenario
In this simulation, we compare the estimation accuracy of various algorithms for targets with different ranges and velocities.Figure 3a,b show the RMSE of the range and velocity estimates of various algorithms versus SNR when R 0 = 11 m and v 0 = 5.8 m/s in a single-target scenario, respectively.As can be seen, the RMSE of the proposed algorithm and the MMEMP algorithm almost coincide and are smaller than the other algorithms.The RMSE of the 2D-FFT algorithm tends to be stable due to the fence effect.The RMSE of the 2D-RootMUSIC algorithm increases rapidly when the SNR decreases as it is sensitive to the noise.The RMSE of the 2D-CZT algorithm is lower than the 2D-FFT algorithm and the 2D-RootMUSIC algorithm, but still larger than the proposed algorithm and the MMEMP algorithm.Tables 2 and 3 show the RMSE values of the range and velocity estimates of various algorithms for targets with different ranges and velocities when SNR = 10 dB and SNR = −10 dB, respectively.As can be seen, the RMSE for different ranges and velocities (including positive and negative velocities) estimates of the proposed algorithm and the MMEMP algorithm are smaller than the other algorithms when SNR = 10 dB and SNR = −10 dB.Therefore, the proposed algorithm and the MMEMP algorithm have similar performance and outperform the other algorithms in terms of range and velocity estimation accuracy in the single target scenario.In this simulation, firstly, two targets with the same initial velocity and close range are, respectively, set as R 0 = [11.9,12] m and v 0 = [5.8,5.8] m/s to compare the range resolution of various algorithms.Secondly, two targets with the same initial range and close velocity are, respectively, set as R 0 = [11,11] m and v 0 = [5.6,5.8] m/s to compare the velocity resolution of various algorithms.

SNR/dB
Figures 4 and 5, respectively, show the scatter plots of the range resolution and velocity resolution of various algorithms when SNR = 10 dB for 10 independent Monte Carlo experiments.As can be seen, the proposed algorithm can distinguish two targets, while the other algorithms cannot.Therefore, the resolution ability of the proposed algorithm is higher than the other algorithms in this case.
Figure 6 shows the RMSE values of various algorithms versus SNR for range estimation of two targets with closer range.As can be seen from this figure, the RMSE of the proposed algorithm decreases as the SNR increases, which is less than the MMEMP algorithm.The 2D-FFT algorithm and the 2D-CZT algorithm fail to distinguish two targets with range difference 0.1 m due to the fence effect, and thus RMSE values of these algorithms are determined by the proximity of the estimated frequencies to the frequency points.Figure 6b shows that the 2D-RootMUSIC algorithm has the largest performance failure due to the estimated false targets.
Figure 7 shows the RMSE values of various algorithms versus SNR for velocity estimation of two targets with similar velocities.As can be seen from this figure, the RMSE values of the proposed algorithm reach 10 Figure 8 shows the scatter plot of the range-velocity pairing feasibility of various algorithms when SNR = 10 dB.As can be seen, the correct pairing of range-velocity estimates can be achieved by all algorithms except for the 2D-CZT algorithm.The reason is that the range and velocity spectral functions of the same target have the maximum correlation, which is utilized by the 2D-CZT algorithm for range-velocity pairing, but it has randomness in sawtooth wave frequency modulated signals.
Figure 9 shows the scatter plot of the range-velocity pairing results of various algorithms when SNR = 10 dB, and Figure 10 shows the correct rate of range-velocity pairing of various algorithms when SNR = 10 dB.As can be seen from Figures 9 and 10, the correct rate of range-velocity pairing can be achieved by the proposed algorithm, the MMEMP algorithm, and the 2D-FFT algorithm without wrong pairing and false detection.An arbitrary false signal is estimated by the 2D-RootMUSIC algorithm when there are targets

Conclusions
This paper proposes a 2D-Unitary ESPRIT-based multi-target joint range-velocity estimation algorithm for the FMCW radar.Firstly, the FMCW radar signal model and range-velocity estimation principle are analyzed.Then, the 2D-spatial smoothing technique is applied to preprocess the signal matrix to construct a 2D virtual array matrix.Finally, the 2D-Unitary ESPRIT algorithm is utilized in this virtual array matrix to obtain the paired range-velocity estimates of multiple targets.The simulation results show that the proposed algorithm has high accuracy in range and velocity estimation in the single-target scenario.In the multi-target scenario, the proposed algorithm has better pairing robustness than the existing algorithms.Moreover, the proposed algorithm has significantly better resolution than other algorithms when there are multiple targets with similar ranges or velocities, which fully reflects the superior performance of the proposed algorithm.The proposed algorithm enhances the range and velocity estimation performance of millimeter-wave FMCW radar, as well as improving its target discrimination capabilities, and thus can be applied to the fields of autonomous driving, target recognition, and related domains.

Figure 1 .
Figure 1.Schematic diagram of the time delay and Doppler shift generated by moving targets.

Algorithm 1 :
2D-Unitary ESPRIT based joint range and velocity estimation algorithm.Input : Discrete IF signal x(m, n); M s and N s .Output : Range-velocity estimates of D targets.

Figure 3 .
Figure 3. RMSE of range and velocity estimates of various algorithms versus SNR in a single-target scenario.(a) RMSE of range estimation versus SNR.(b) RMSE of velocity estimation versus SNR.

Figure 10 .
Figure 10.Correct rate of range-velocity pairing of various algorithms when SNR = 10 dB.

Table 2 .
RMSE of the range and velocity estimates when the SNR is 10 dB.

Table 3 .
RMSE of the range and velocity estimates when the SNR is −10 dB.