A Novel Unitary ESPRIT Algorithm for Monostatic FDA-MIMO Radar

Abstract

A novel unitary estimation of signal parameters via rotational invariance techniques (ESPRIT) algorithm, for the joint direction of arrival (DOA) and range estimation in a monostatic multiple-input multiple-output (MIMO) radar with a frequency diverse array (FDA), is proposed. Firstly, by utilizing the property of Centro-Hermitian of the received data, the extended real-valued data is constructed to improve estimation accuracy and reduce computational complexity via unitary transformation. Then, to avoid the coupling between the angle and range in the transmitting array steering vector, the DOA is estimated by using the rotation invariance of the receiving subarrays. Thereafter, an automatic pairing method is applied to estimate the range of the target. Since phase ambiguity is caused by the phase periodicity of the transmitting array steering vector, a removal method of phase ambiguity is proposed. Finally, the expression of Cramér–Rao Bound (CRB) is derived and the computational complexity of the proposed algorithm is compared with the ESPRIT algorithm. The effectiveness of the proposed algorithm is verified by simulation results.

1. Introduction

The multiple-input multiple-output (MIMO) radar [1,2], which utilizes multiple antennas to simultaneously transmit diverse waveforms and receive the reflected signals in similar ways, has many potential advantages [3]. Unlike the conventional phased-array (PA) radar, MIMO radar has many superiorities based on its spatial diversity and waveform diversity, such as improving the system performance with higher degrees-of-freedom (DOFs) [4,5]. Normally, MIMO radar can be classified into two types based on the spatial location of antenna elements. One is named distributed MIMO radar, where the transmitting or receiving array elements are placed in different positions, with a relatively large spacing of elements [6]. The other is the collocated MIMO radar, where both transmitting and receiving antennas are arranged close to each other [7]. According to whether the receiver and transmitter are located in the same place, the collocated MIMO radar is further categorized into the monostatic MIMO radar [8,9] and bistatic MIMO radar [10,11]. A monostatic MIMO radar is superior in its excellent maneuverability and synchronization between the transmitter and receiver. In the contemporary defense system, the monostatic radar system is the most mainstream and common sensor unit in the modern radar network system. A monostatic MIMO radar is addressed in this paper.

In radar systems, the target localization is one of the main issues. The angle and range estimation are the important parts for target localization [12,13,14,15,16,17,18,19]. However, in the beam scanning of both the PA radar and MIMO radar, the beam pointing is angle-dependent and range-independent. The frequency diverse array (FDA) radar, in which the direction of focus changes with range, is proposed in [20]. For the FDA radar, the frequencies of each transmitting array element are different, which can lead to a range-angle-dependent beampattern [21]. Then, joint angle and range can be estimated simultaneously in the FDA radar [22,23,24]. Nevertheless, the resolution of parameter estimation is influenced by the maximum frequency increment and the array aperture [25]. The maximum number of distinguishable targets is determined by the number of DOFs and the frequency increment. The range-dependent beampattern produced by a linear frequency modulation structure is analyzed [26]. In order to estimate angle and range, a double-pulse method was introduced to obtain the range and angle separately, where the antenna successively transmits two pulses with non-zero and zero frequency increments [27]. A transmitting subaperture scheme in FDA radar was investigated, where the uniform linear array (ULA) is divided into multiple overlapping subarrays [28]. To decouple the angle and range, a special FDA radar with different frequency increments was proposed [29]. A method based on nonuniform frequency increment was proposed to decouple the FDA beampattern [30]. There are some other interesting investigations reported in [31,32,33,34] about decoupling range and angle in the FDA radar.

In the last few years, FDA-MIMO radar was regarded as a novel radar system combining the advantages of FDA radar with MIMO radar [35,36,37]. The FDA-MIMO radar can decouple angle and range by exploiting the high DOFs of MIMO technique and the angle-range-dependent beampattern of FDA radar [28,37]. Furthermore, Xu et al. applied a maximum likelihood estimator to obtain an unambiguous angle and range estimation [38]. The sparse reconstruction-based algorithm was utilized for the target location with an FDA-MIMO radar [39]. Multiple signal classification (MUSIC) algorithms were also introduced to the FDA-MIMO radar [40,41]. However, these algorithms with peak-searching require a large computational complexity. Although a method based on estimation of signal parameters via rotational invariance techniques (ESPRIT) algorithm has been introduced to obtain the range and angle estimation and reduce the complexity [42], there are still some challenges, such as a low estimation performance, high complexity and phase ambiguity.

In this paper, an improved unitary ESPRIT method for angle and range estimation in monostatic FDA-MIMO radar is presented. Then, a method is proposed to achieve the automatic pairing of angle and range. The reason for phase ambiguity is analyzed, and a simple and efficient solution is proposed. The Cramér–Rao Bound (CRB) for target parameter of monostatic FDA-MIMO radar is derived. Computational complexity analysis is provided to be compared with that of [42].

The paper is summarized as follows. Section 2 introduces the signal model of the monostatic FDA radar. We propose an improved unitary ESPRIT algorithm with an automatic pairing, and investigate a phase judgment in Section 3. In Section 4, we derive the CRB of the monostatic FDA-MIMO radar and offer the specific analysis of computational complexity. Several simulation results are provided to indicate the effectiveness of the presented algorithm in Section 5. The conclusion of the paper is indicated in Section 6.

2. Signal Model of Monostatic FDA-MIMO Radar

The signal model of monostatic FDA-MIMO radar is shown in Figure 1, which is configured with M-antenna-transmitting ULA and N-antenna-receiving ULA. The transmitting and receiving arrays are placed together to guarantee that the DOA is the same as the direction of departure (DOD). The element interval of the transmitting array is equal to that of receiving array, and they are both marked with d. The element interval d is set to half of the maximal wavelength. The first element of transmitting array is assumed as the reference. The frequency of the m-th signal emitted by the transmitting array is calculated by

$$f_m=f_0+(m-1)\mathsf{\Delta }f,\hspace{1em}m=1,2,\cdots ,m$$

(1)

where f₀ is the referenced frequency, and Δf is the deviation of adjacent transmitting frequencies, where Δf << f₀. The baseband signals emitted by the different elements are mutually orthogonal. Assume that the emitted signal of m-th element satisfies

$$s_m(t)=\sqrt{\frac{E}{M}}{\varphi }_m(t)e^{j2\pi f_{m}t}, 0\le t\le t$$

(2)

where E is the total transmitted energy, T is the radar pulse duration, and ${\varphi }_m(t)$ denotes the baseband waveform. Suppose that all baseband signals are orthogonal to each other and normalized

$${\displaystyle \underset{T}{\int }{\varphi }_m(t){\varphi }_n^{*}(t-\tau )e^{j2\pi (m-n)\mathsf{\Delta }ft}}dt=\{\begin{array}{c}0,\\ 1,\end{array}\begin{array}{c}\hspace{1em}m\ne n,\\ \hspace{1em}m=n\end{array}\begin{array}{c}\forall \tau \\ \end{array}$$

(3)

where (.)* stands for the conjugate operator and τ is the time delay. There are K-independent targets distributed in the far-field, and the ranges of all targets are much larger than the FDA-MIMO radar aperture. The angle and range of k-th target are represented by (θ_k, r_k). Due to the limitation of the maximum unambiguity range, r_k should be less than c/(2Δf), and c denotes the speed of signal propagation, and c/(2Δf) represents the maximum unambiguity range [32]. The expression of steering vector of transmitting array is [33]

$$a(r_k,{\theta }_k)=r(r_k)\odot d({\theta }_k)\in {ℂ}^{M\times 1}$$

(4)

$$r(r_k)=[1,e^{-j4\pi \mathsf{\Delta }f{r}_k/c},\cdots ,e^{-j(m-1)4\pi \mathsf{\Delta }f{r}_k/c}{]}^T\in {ℂ}^{M\times 1}$$

(5)

$$d({\theta }_k)=[1,e^{j\pi \sin {\theta }_k},\cdots ,e^{j(m-1)\pi \sin {\theta }_k}{]}^T\in {ℂ}^{M\times 1}$$

(6)

where r(r_k) and d(θ_k) denote the range-dependent part and the angle-dependent part of transmitting the joint-steering vector of the k-th target, respectively. $\odot$ stands for Hadamard product operator. (.)^T denotes the transpose operator. The coupling relationship between angle and range can be shown in Equation (4). The receiving spatial steering vector can be given by [34]

$$b({\theta }_k)=[1,e^{j\pi \sin {\theta }_k},\cdots ,e^{j(N-1)\pi \sin {\theta }_k}{]}^T\in {ℂ}^{N\times 1}$$

(7)

Then, the received data at the receiving array can be described as [11]

$$X(t)={\displaystyle \sum _{k=1}^K{\beta }_k{e}^{j2\pi f_{pk}t}b({\theta }_k)a^T(r_k,{\theta }_k)}\left[\begin{array}{c}{s}_1(t)\\ ⋮\\ s_M(t)\end{array}\right]+n(t)$$

(8)

where β_k and f_pk represent the amplitude and Doppler frequency of the k-th target, respectively. n(t) represents an N × 1 complex Gaussian white noise vector with zero mean. According to Equation (3), the received data after the matched filtering is expressed as [37]

$$Y(t)=CH(t)+N(t)$$

(9)

$$C=[c_1(r_1,{\theta }_1),\cdots ,c_K(r_K,{\theta }_K)]=\left[a\right(r_1,{\theta }_1)\otimes b({\theta }_1),\cdots ,a(r_K,{\theta }_K)\otimes b({\theta }_K\left)\right]$$

(10)

$$H(t)=\sqrt{\frac{E}{M}}{\left[{\beta }_1{e}^{j2\pi f_{p1}t},{\beta }_2{e}^{j2\pi f_{p2}t},\dots ,{\beta }_k{e}^{j2\pi f_{pk}t}\right]}^T$$

(11)

where C is the joint transmit-receive steering matrix, and ⊗ denotes the Kronecker product. H(t) is the signal matrix after matched filters. N(t) represents the noise vector after matching filter with the transmitted signal, and the noise covariance matrix is σ²I_MN, where σ² and I_MN denote the noise variance and M × N identity matrix, respectively.

3. Unitary ESPRIT Algorithm for Ange and Range Estimation

3.1. Rotation Invariance of Subarrays

In this section, the complex-valued rotation invariance relation is introduced with respect to the transmitting array and receiving array, respectively. As shown in Figure 2, the transmitting array and receiving array of monostatic FDA-MIMO radar are divided into two overlapping subarrays,. Assuming that the two adjacent subarrays are identical, there exists a rotation invariance between Subarray 1 and Subarray 2, or Subarray 3 and Subarray 4. According to Equation (4), due to the existence of a coupling relationship between DOA and the range, it is necessary to obtain DOA information by the rotation invariance relationship of the receiving subarrays, and substitute the estimated DOA into the rotation invariance of transmitting arrays to get the range information. The complex-valued invariance relationship of the Subarray 3 and Subarray 4 can be expressed as [42]

$$J_{2}b({\theta }_k)=e^{j\pi \sin {\theta }_k}{J}_{1}b({\theta }_k)$$

(12)

where $J_1=[\begin{array}{cc}{I}_{N-1}& 0_{(N-1)\times 1}\end{array}]$ and $J_2=[\begin{array}{cc}{0}_{(N-1)\times 1}& I_{N-1}\end{array}]$ are selection matrices, and 0_w denotes the w × w null matrix. The functions of J₁ and J₂ are to select the first and last N − 1 rows of a matrix, respectively. Based on Equation (12), the invariance relationship in the joint steering vector can be expressed as

$$(I_M\otimes J_2)c(r_k,{\theta }_k)=e^{j\pi \sin {\theta }_k}(I_M\otimes J_1)c(r_k,{\theta }_k)$$

(13)

Under the assumption of K targets, the rotation invariance of each target can be written into a matrix form

$$(I_M\otimes J_2)C=(I_M\otimes J_1)C{\mathsf{\Xi }}_R$$

(14)

$${\mathsf{\Xi }}_R=\left[\begin{array}{ccc}{e}^{j\pi \sin {\theta }_1}& & \\ & \ddots & \\ & & e^{j\pi \sin {\theta }_k}\end{array}\right]$$

(15)

where ${\mathsf{\Xi }}_R$ contains the angle information of all targets. It can be shown that the columns in C span the same signal subspace as the column vectors in the signal subspace E_s [43], we can obtain the following relationship

$$E_S=C\mathsf{\Theta }$$

(16)

where E_s is composed of the K eigenvectors corresponding to the largest K eigenvalues of the covariance matrix of Y, and $\mathsf{\Theta }$ is a non-singular matrix. Substituting Equation (16) into Equation (14), we can obtain the relationship between signal space of the subarrays

$$(I_M\otimes J_2)E_S=(I_M\otimes J_1)E_S{\mathsf{\Psi }}_R$$

(17)

where ${\mathsf{\Psi }}_R={\mathsf{\Theta }}^{-1}{\mathsf{\Xi }}_R\mathsf{\Theta }$. It can be noticed that the diagonal elements of ${\mathsf{\Xi }}_R$ are the eigenvalues of Ψ_R.

Next, the complex-valued rotation invariance between Subarray 1 and Subarray 2 is considered to estimate range. The transmitting steering vectors of Subarray 1 and Subarray 2 satisfy the following equation

$$J_{4}a(r_k,{\theta }_k)={\mathrm{e}}^{j(\pi \sin {\theta }_k-4\pi \frac{\mathsf{\Delta }f}{c}{r}_k)}{J}_{3}a(r_k,{\theta }_k)$$

(18)

where $J_3=[\begin{array}{cc}{I}_{M-1}& 0_{(m-1)\times 1}\end{array}]$ and $J_4=[\begin{array}{cc}{0}_{(m-1)\times 1}& I_{M-1}\end{array}]$ stand for selection matrices to select the first and last M − 1 rows, respectively. For K targets, this relationship is extended to the joint steering vector which can be expressed as

$$(I_N\otimes J_4)C=(I_N\otimes J_3)C{\mathsf{\Xi }}_T$$

(19)

$${\mathsf{\Xi }}_T=\left[\begin{array}{ccc}{e}^{j(\pi \sin {\theta }_1-4\pi \frac{\mathsf{\Delta }f}{c}{r}_1)}& & \\ & \ddots & \\ & & e^{j(\pi \sin {\theta }_K-4\pi \frac{\mathsf{\Delta }f}{c}{r}_K)}\end{array}\right]$$

(20)

where ${\mathsf{\Xi }}_T$ contains the ranges of all the targets. According to Equation (16), Equation (19) can be rewritten as

$$(I_N\otimes J_4)E_S=(I_N\otimes J_3)E_S{\mathsf{\Psi }}_T$$

(21)

where ${\mathsf{\Psi }}_T={\mathsf{\Theta }}^{-1}{\mathsf{\Xi }}_T\mathsf{\Theta }$. Because the calculation of the covariance matrix of Y and the acquisition of E_S are based on complex-valued data, DOAs and ranges are estimated with relatively high complexity. Hence, based on the idea of a unitary ESPRIT algorithm, a novel unitary ESPRIT algorithm is proposed to reduce complexity and improve estimation accuracy.

3.2. Unitary ESPRIT in FDA-MIMO Radar

In this section, a complex-valued invariance is transformed into a real-valued invariance, and the DOAs and ranges are estimated by using the unitary ESPRIT algorithm. As there is no central Hermitian symmetric characteristic in Y, an extended receiving data matrix with the symmetric structure is defined as [44,45]

$$Z=\left[\begin{array}{cc}Y& {\mathsf{\Pi }}_{\mathrm{MN}}{Y}^{*}\end{array}{\mathsf{\Pi }}_L\right]$$

(22)

where Π_MN is an M × N exchange matrix with ones on its anti-diagonal and zeros elsewhere. Over the construction of Equation (22), Z is a generalized Centro-Hermitian matrix. The complex matrix Z is transformed into the real-valued matrix Γ by utilizing the unitary transformation. It can be expressed as [43]

$$\mathsf{\Gamma }=Q_{MN}^{H}Z{Q}_{2L}=Q_{MN}^H\left[\begin{array}{cc}Y& {\mathsf{\Pi }}_{\mathrm{MN}}{Y}^{*}\end{array}{\mathsf{\Pi }}_L\right]Q_{2L}$$

(23)

where (.)^H represents the conjugate transpose operator, and the sparse unitary matrix Q_w is defined as

$$\{\begin{array}{c}{Q}_w=\frac{1}{\sqrt{2}}\left[\begin{array}{cc}{I}_w& j{I}_w\\ {\mathsf{\Pi }}_w& -j{\mathsf{\Pi }}_w\end{array}\right]\\ Q_w=\frac{1}{\sqrt{2}}\left[\begin{array}{ccc}{I}_w& 0& j{I}_w\\ 0^T& \sqrt{2}& 0^T\\ {\mathsf{\Pi }}_w& 0& -j{\mathsf{\Pi }}_w\end{array}\right]\end{array}\begin{array}{ccc}& & \begin{array}{c}w \mathrm{is} \mathrm{even}\\ \\ w \mathrm{is} \mathrm{odd}\end{array}\end{array}$$

(24)

Compared with Equation (9), Equation (22) is competent in doubling the number of snapshots. Then, the real-valued covariance R_Γ of the extended received data can be acquired by using the maximum likelihood estimation

$$\begin{array}{l}{R}_{\mathsf{\Gamma }}= \frac{1}{2L}\mathsf{\Gamma }{\mathsf{\Gamma }}^H\\ \begin{array}{cc}& \end{array}=\frac{1}{2}\left(Q_{MN}^H{R}_Z{Q}_{MN}\right)\\ \begin{array}{cc}& \end{array}=\frac{1}{2}\left[Q_{MN}^H\left(R_Y+{\mathsf{\Pi }}_{MN}{R}_Y^{*}{\mathsf{\Pi }}_{MN}\right)Q_{MN}\right]\end{array}$$

(25)

where R_Y and R_Z are the covariance calculated by Y and Z, respectively. The signal subspace ${\hat{E}}_S$ corresponds to K eigenvectors of large eigenvalues of R_Γ. The remaining MN − K eigenvectors of small eigenvalues can obtain the noise subspace ${\hat{E}}_N$. Hence, ${\hat{E}}_S$ and ${\hat{E}}_N$ are both real-valued. Due to the unitary transformation in Equation (23), the complex-valued invariance relation in Equation (13) is transformed into real-valued invariance relation as follows

$$K_2{d}_k=\tan \left(\frac{\pi \sin {\theta }_k}{2}\right)K_1{d}_k$$

(26)

$$K_1=\mathrm{Re}\{Q_{(N-1)m}^H(I_M\otimes J_2)Q_{MN}\}$$

(27)

$$K_2=\mathrm{Im}\{Q_{(N-1)m}^H(I_M\otimes J_2)Q_{MN}\}$$

(28)

where $d_k=Q_{MN}^H{c}_k$ is the real-valued steering vector. Re{.} and Im{.} denote the real and imaginary parts of a complex number, respectively. Considering K-independent targets, Equation (26) is expressed as

$$K_{2}D=K_{1}D{\mathsf{\Phi }}_R$$

(29)

$${\mathsf{\Phi }}_R=\left[\begin{array}{ccc}\tan (\frac{\pi \sin {\theta }_1}{2})& & \\ & \ddots & \\ & & \tan (\frac{\pi \sin {\theta }_k}{2})\end{array}\right]$$

(30)

where D = [d₁, d₂, …, d_k], and Φ_R is a real diagonal matrix whose diagonal elements contain the desired angle information. According to Equation (16), Equation (29) can be rewritten as

$$K_2{\hat{E}}_S=K_1{\hat{E}}_S{\mathsf{\Sigma }}_R$$

(31)

where ${\mathsf{\Sigma }}_R={\mathsf{\Theta }}_R^{-1}{\mathsf{\Phi }}_R{\mathsf{\Theta }}_R$. ${\mathsf{\Theta }}_R$ is the left eigenvector matrix of ${\mathsf{\Sigma }}_R$. By using the total least squares (TLS) method to solve Equation (31), DOA can be estimated as follows

$${\widehat{\theta }}_k=\arcsin \left(\frac{2\arctan {[({\mathsf{\Phi }}_R)]}_k}{\pi }\right)$$

(32)

Similarly, the rotation invariance between Subarray 1 and Subarray 2 can be transformed into

$$K_4{d}_k=\tan \left(\frac{\pi (\sin {\theta }_k-4\mathsf{\Delta }f\frac{r_k}{c})}{2}\right)K_3{d}_k$$

(33)

$$K_3=\mathrm{Re}\{Q_{(m-1)N}^H(I_N\otimes J_4)Q_{MN}\}$$

(34)

$$K_4=\mathrm{Im}\{Q_{(m-1)N}^H(I_N\otimes J_4)Q_{MN}\}$$

(35)

For K targets, Equation (33) can be integrated into matrix form

$$K_{4}D=K_{3}D{\mathsf{\Phi }}_T$$

(36)

$${\mathsf{\Phi }}_T=\left[\begin{array}{ccc}\tan (\frac{\pi \sin {\theta }_1-4\pi \mathsf{\Delta }f\frac{r_1}{c}}{2})& & \\ & \ddots & \\ & & \tan (\frac{\pi \sin {\theta }_k-4\pi \mathsf{\Delta }f\frac{r_k}{c}}{2})\end{array}\right]$$

(37)

where Φ_T is a real diagonal matrix whose diagonal elements contain the desired range of information. In the same way as Equation (31), we can obtain the rotational invariance of the signal subspace

$$K_4{\widehat{E}}_S=K_3{\widehat{E}}_S{\mathsf{\Sigma }}_T$$

(38)

where ${\mathsf{\Sigma }}_T={\mathsf{\Theta }}_T^{-1}{\mathsf{\Phi }}_T{\mathsf{\Theta }}_T$. ${\mathsf{\Theta }}_T$ is the inverse of the left eigenvector matrix of ${\mathsf{\Sigma }}_T$. Substituting Equation (32) into Equation (37), the range estimation is solved with the TLS method.

$${\widehat{r}}_k=\frac{2\arctan {[({\mathsf{\Phi }}_R)]}_k-2\arctan {[({\mathsf{\Phi }}_T)]}_k}{4\pi \mathsf{\Delta }f}c$$

(39)

Due to the correlation of Φ_R and Φ_T in Equation (39), the ranges will be miscalculated without the pairing of Φ_R and Φ_T. Hence, we employ an automatic pairing method to implement correct range estimation.

3.3. The Pairing of DOAs and Ranges

In this section, we analyze the speciality of ${\mathsf{\Theta }}_T$ and ${\mathsf{\Theta }}_R$, and introduce the means to achieve pairing. ${\mathsf{\Theta }}_T$ and ${\mathsf{\Theta }}_R$ are eigenvectors of ${\mathsf{\Sigma }}_T$ and ${\mathsf{\Sigma }}_R$, respectively. Since ${\mathsf{\Sigma }}_T$ and ${\mathsf{\Sigma }}_R$ are calculated by ${\hat{E}}_S$, there must be a random row of ${\mathsf{\Theta }}_R$ identical to a specific row of ${\mathsf{\Theta }}_T$. Supposing that all of the K targets are independent, we notice that any two rows of ${\mathsf{\Theta }}_R$ are orthogonal because any two eigenvalues of ${\mathsf{\Theta }}_T$ are different. In this paper, considering algorithm complexity, we obtain the automatic pairing of Φ_T and Φ_R by decomposing the ${\mathsf{\Sigma }}_T+j{\mathsf{\Sigma }}_R$, which can be expressed as

$${\mathsf{\Sigma }}_T+j{\mathsf{\Sigma }}_R={\mathsf{\Theta }}_{TR}^{-1}\{{\mathsf{\Phi }}_T+j{\mathsf{\Phi }}_R\}{\mathsf{\Theta }}_{TR}$$

(40)

where ${\mathsf{\Theta }}_{TR}$ is the left eigenvector matrix. Hence, Φ_T and Φ_R can be automatically paired by the eigenvector matrix ${\mathsf{\Theta }}_{TR}$. Considering the periodic phase ambiguity problem, we take a step to distinguish the phase ambiguity before calculating ranges.

3.4. The Solution of Periodic Ambiguity of Transmitter

In this section, we analyze the periodic phase ambiguity and adopt an ambiguity judgment method to obtain the correct range estimation. Since the period of tan in Equation (26) is π, and sinθ_k ∈ (−1, 1), DOAs can be estimated by Equation (32) without periodic ambiguity. However, there is phase ambiguity in range estimation due to sinθ_k − 4Δfr_k/c ∈ (−3, 1) and r_k ∈ (0, c/2Δf). Therefore, r_k obtained by Equation (39) is misestimated, because the tangent of π(sinθ_k − 4Δfr_k/c)/2 is equal to the tangent of π(sinθ_k − 4Δfr_k/c)/2 + π when sinθ_k − 4Δfr_k/c ∈ (−3, −1). Note that π(sinθ_k − 4Δfr_k/c)/2 < πsinθ_k/2< π(sinθ_k − 4Δfr_k/c)/2 + π. Since Φ_T and Φ_R are relevant to π(sinθ_k − 4Δfr_k/c)/2 and πsinθ_k/2, respectively, we determine the range of π(sinθ_k − 4Δfr_k/c)/2 by comparing arctan[Φ_T]_k and arctan[Φ_R]_k. If arctan[Φ_T]_k > arctan[(Φ_R)]_k, due to arctan[Φ_T]_k ∈ (−π/2, π/2), there is a phase ambiguity, as π(sinθ_k − 4Δfr_k/c)/2 is considered to lie in (−π/2, π/2) when π(sinθ_k − 4Δfr_k/c)/2 ∈ (−3π/2, −π/2). We can use a phase shift π to solve the periodic phase ambiguity problem. Hence, the true phase value of Φ_T can be calculated as

$$\frac{\pi }{2}(\sin {\widehat{\theta }}_k-4\mathsf{\Delta }f{r}_k/c)=\arctan ({\mathsf{\Phi }}_T)-\pi$$

(41)

r_k can be given by

$${\widehat{r}}_k=\frac{\arctan {[({\mathsf{\Phi }}_R)]}_k-\arctan {[({\mathsf{\Phi }}_T)]}_k+\pi }{2\pi \mathsf{\Delta }f}c$$

(42)

Otherwise, there is no phase ambiguity and the true phase value of Φ_T can be given by

$$\frac{\pi }{2}(\sin {\widehat{\theta }}_k-4\mathsf{\Delta }f{r}_k/c)=\arctan ({\mathsf{\Phi }}_T)$$

(43)

r_k can be given by

$${\widehat{r}}_k=\frac{\arctan {[({\mathsf{\Phi }}_R)]}_k-\arctan {[({\mathsf{\Phi }}_T)]}_k}{2\pi \mathsf{\Delta }f}c$$

(44)

The main steps of the proposed algorithm are summarized in Algorithm 1.

Algorithm 1: A Novel Unitary ESPRIT for Monostatic FDA-MIMO Radar
1:	Construct the extended received data matrix Z via (22).
2:	Take the unitary transformation using (23) and obtain RΓ via (25).
3:	Perform eigenvalue decomposition of RΓ and return ${\hat{E}}_S$.
4:	Calculate ${\mathsf{\Sigma }}_T$ and ${\mathsf{\Sigma }}_R$ via (31) and (38), respectively.
5:	Obtain the automatically paired ΦT and ΦR via (40).
6:	Compute the DOAs by (32) and ranges by (42) or (44).

4. CRB and Complexity Analysis

4.1. CRB

In this section, we analyze the CRBs of angle and range for the monostatic FDA-MIMO radar. According to Equation (9), the concrete expression of R_Y is written as

$$R_Y=\frac{1}{L}Y{Y}^H=C{R}_H{C}^H+{\sigma }^{2}I$$

(45)

where $R_H$ is the covariance of H in Equation (9), and L is the number of snapshots, and σ² denotes the noise power. Under the assumption of K targets, the unknown parameter to be estimated is

$$\eta ={[{\theta }^T,r^T]}^T={[{\theta }_1,\dots ,{\theta }_k,r_1,\dots ,r_k]}^T$$

(46)

Then, the Fisher information matrix (FIM), with respect to η, is [41]

$$F=\left[\begin{array}{cc}{F}_{\theta \theta }& F_{\theta r}\\ F_{r\theta }& F_{\mathrm{rr}}\end{array}\right]$$

(47)

The expression of every block of F can be written as

$$F_{\theta \theta }=\frac{2L}{{\sigma }^2}\mathrm{Re}\left\{\left(C_{\theta }^H{\mathsf{\Pi }}_C^{\perp }{C}_{\theta }\right)\odot {\left(R_H{C}^H{R}_Y^{-1}C{R}_H\right)}^T\right\}$$

(48)

$$F_{\theta r}=\frac{2L}{{\sigma }^2}\mathrm{Re}\left\{\left(C_{\theta }^H{\mathsf{\Pi }}_C^{\perp }{C}_r\right)\odot {\left(R_H{C}^H{R}_Y^{-1}C{R}_H\right)}^T\right\}$$

(49)

$$F_{r\theta }=\frac{2L}{{\sigma }^2}\mathrm{Re}\left\{\left(C_r^H{\mathsf{\Pi }}_C^{\perp }{C}_{\theta }\right)\odot {\left(R_H{C}^H{R}_Y^{-1}C{R}_H\right)}^T\right\}$$

(50)

$$F_{rr}=\frac{2L}{{\sigma }^2}\mathrm{Re}\left\{\left(C_r^H{\mathsf{\Pi }}_C^{\perp }{C}_r\right)\odot {\left(R_H{C}^H{R}_Y^{-1}C{R}_H\right)}^T\right\}$$

(51)

where C_θ and C_r are partial derivations of C with respect to θ and r, and ${\mathsf{\Pi }}_C^{\perp }=I-C{\left(C^{H}C\right)}^{-1}{C}^H$. The C_θ and C_r can be written as

$$C_{\theta }=\left[c_{1\theta },c_{2\theta },\dots ,c_{K\theta }\right]$$

(52)

$$C_r=\left[c_{1r},c_{2r},\dots ,c_{\mathrm{Kr}}\right]$$

(53)

where

$$c_{k\theta }=\frac{\partial c({\theta }_k,r_k)}{\partial {\theta }_k}=\frac{\partial a({\theta }_k,r_k)}{\partial {\theta }_k}\otimes b({\theta }_k)+a({\theta }_k,r_k)\otimes \frac{\partial b({\theta }_k)}{\partial {\theta }_k}\begin{array}{cc}& \end{array}k=1,2,\dots K$$

(54)

$$c_{kr}=\frac{\partial c({\theta }_k,r_k)}{\partial r_k}=\frac{\partial a({\theta }_k,r_k)}{\partial r_k}\otimes b({\theta }_k),\begin{array}{cc}& \end{array}k=1,2,\dots K$$

(55)

We derive part of Equations (54) and (55) as

$$\frac{\partial a({\theta }_k,r_k)}{\partial {\theta }_k}=j\pi \cos ({\theta }_k)\left[\begin{array}{ccc}0& & \\ & \ddots & \\ & & M-1\end{array}\right]a({\theta }_k,r_k)$$

(56)

$$\frac{\partial b({\theta }_k)}{\partial {\theta }_k}=j\pi \cos ({\theta }_k)\left[\begin{array}{ccc}0& & \\ & \ddots & \\ & & N-1\end{array}\right]b({\theta }_k)$$

(57)

and

$$\frac{\partial a({\theta }_k,r_k)}{\partial r_k}=-j4\pi \frac{\mathsf{\Delta }f}{c}\left[\begin{array}{ccc}0& & \\ & \ddots & \\ & & M-1\end{array}\right]a({\theta }_k,r_k)$$

(58)

Then, every block of F is determined by Equations (48)–(51). Then, the CRB matrix can be obtained by

$$CRB=F^{-1}=\frac{{\sigma }^2}{2L}\mathrm{Re}{\left\{\left(W^H{\mathsf{\Pi }}_C^{\perp }W\right)\odot P^T\right\}}^{-1}$$

(59)

where $W=[C_{\theta } C_r]$, $P=\left[\begin{array}{cc}{P}_1& P_2\\ P_3& P_4\end{array}\right]$ and $P_1=P_2=P_3=P_4=R_H{C}^H{R}_Y^{-1}C{R}_H$.

4.2. Complexity

To compare the presented algorithm and the ESPRIT algorithm of [42], the specific analysis of the computational complexity is provided. In this paper, we transform received data from the complex domain to the real domain by a unitary transformation. Hence, the calculation of eigenvalue decomposition and generalized inverse depend on the real domain. The calculation of a complex product is equivalent to the calculation of four real products. The concentration of computational complexity in the presented algorithm is based on calculating the covariance matrix, utilizing the eigenvalue decomposition, obtaining signal subspace, solving the solution for angle and range, and achieving pairing for angle and range. The calculation of R_Γ needs O{2L(MN)²} flops, where M and N denote the number of transmitting and receiving array elements, respectively, and L is the number of snapshots. The eigenvalue decomposition of R_Γ, to obtain the signal subspace and the noise subspace, needs O{(MN)³} flops. The complexity required to solve for ${\mathsf{\Sigma }}_R$ is O{M(N − 1)(2K)²}, where K is the number of targets. Similarly, solving ${\mathsf{\Sigma }}_T$ needs O{N(M − 1)(2K)²} flops. The eigenvalue decomposition in Equation (40) and the pairing of angle and range need O{4(2K)³ + 2K³}. Here, we ignore the complexity of solving periodic ambiguity steps because they are too small. Thus, the complexity of the proposed algorithm is

$$O\left\{2L{\left(MN\right)}^2+{(MN)}^3+4K^2(2MN-m-N)+34K^3\right\}$$

(60)

In [42], the complexity of the ESPRIT algorithm for estimation of angle and range is

$$O\left\{2L{\left(MN\right)}^2+{(2MN)}^3+4K^2(5MN-2m-2N)+31K^3\right\}$$

(61)

By the comparison of Equations (60) and (61), the computational complexity of the proposed algorithm is much lower than [42]. Furthermore, later in the simulation, we give the comparison results regarding complexity.

5. Simulation Results

In this section, we provide several simulation results to evaluate the performance of the proposed algorithm for angle and range estimation in monostatic FDA-MIMO radar with ULA. The ESPRIT algorithm in the same model is chosen for comparison [42]. In all simulations, assume that the reference carrier frequency f₀, namely the minimum frequency, is 3 GHz, and the frequency increment Δf is 10³ Hz. According to the relationship of Equation (1), the maximum frequency and the number of bins depend on the number of transmitting arrays. The noise is assumed to be the uniform complex white Gaussian noise. The reflection coefficient of the target is set to 1. The number of Monte Carlo experiments is set to 500.

5.1. Estimated Results

In this section, the SNR is set to 10 dB, and the number of snapshots is 50, and the number of transmitting array elements M and receiving array elements N are both set to eight. Figure 3a,b shows the unpaired and paired estimation of range, respectively, obtained by the proposed algorithm, where the two-dimensional parameters of the target are set to (45°, 40 km) and (30°, 10 km). It is noted that an incorrect range estimation is shown in Figure 3a, which is caused by the mismatch between the eigenvalues of ${\mathsf{\Sigma }}_R$ and ${\mathsf{\Sigma }}_T$. It is seen in Figure 3b that the pairing method can obtain the correct range estimation. Figure 4a,b shows the estimation results of angle and range acquired by the proposed algorithm and the ESPRIT algorithm, where the targets in Figure 4a are the same as those of Figure 3, and the targets in Figure 4b are assumed to be (45°, 40 km) and (−30°, 70 km). As there is no period ambiguity for (45°, 40 km) and (30°, 10 km), the proposed algorithm and the ESPRIT algorithm can both obtain accurate estimation, as shown in Figure 4a. It is shown in Figure 4b that the proposed algorithm can effectively solve the period ambiguity, because (−30°, 70 km) satisfies sinθ_k − 4Δfr_k/c ∈ (−3, −1).

5.2. RMSE Versus SNR

In this section, the target is set to (45°, 40 km) and (30°, 10 km), and the number of snapshots is 50. The number of transmitting array elements M and receiving array elements N are set to M = N = 4 and M = N = 8, respectively. The SNR increases from 0 dB to 20 dB, with each step being 2 dB. Figure 5a,b shows the root mean square errors (RMSEs) of the proposed algorithm and the ESPRIT algorithm with the different SNR. Meanwhile, the CRBs of angle and range in the monostatic FDA-MIMO radar are chosen for the assessment of the performance of the proposed algorithm. The RMSEs of angle and range are respectively defined as

$$RMS{E}_{\theta }=\sqrt{\frac{1}{G}\frac{1}{K}{\displaystyle \sum _{g=1}^G{\displaystyle \sum _{k=1}^K{({\theta }_k-{\widehat{\theta }}_k)}^2}}}$$

(62)

$$RMS{E}_r=\sqrt{\frac{1}{G}\frac{1}{K}{\displaystyle \sum _{g=1}^G{\displaystyle \sum _{k=1}^K{(r_k-{\widehat{r}}_k)}^2}}}$$

(63)

where G represents the number of Monte Carlo experiments. We can observe that the RMSEs of the proposed algorithm are closer to the CRBs. This indicates that the performance of the proposed method is better than the ESPRIT algorithm with the identical SNR.

5.3. RMSE Versus Number of Snapshots

In the simulation, the target is set to (45°, 40 km) and (30°, 10 km), and the SNR is 0 dB. The number of transmitting array elements M and receiving array elements N are both set to eight and four, respectively. We set the initial number of snapshots to be 50, and observe the effect of the number of snapshots on the RMSEs by intervals of 100. Figure 6a, b shows the RMSEs of angle and range versus the number of snapshots, respectively. This indicates that the performance of the proposed algorithm is better than that of the ESPRIT algorithm with the same number of snapshots.

5.4. Computational Complexity

In this part, the runtime of the proposed algorithm is compared with that of the ESPRIT algorithm. The target is set to (45°, 40 km), and (30°, 10 km), the SNR is set to 0 dB, and the number of snapshots is 50. The number of transmitting array elements is equal to that of the receiving array elements, i.e., M = N, and the transmitting array number M is changed in this simulation. The required runtime of the two algorithms is shown in Figure 7. The runtime of the proposed algorithm is less than that of the ESPRIT algorithm.

We can summarize this with two situations, according to the existence of periodic ambiguity. In the case of periodic ambiguity, the ESPRIT algorithm cannot obtain the correct estimation of target parameters, but the proposed algorithm can accurately get the angles and ranges of the target. In the absence of periodic ambiguity, the proposed algorithm and ESPRIT algorithm can both obtain a correct estimation of target parameters. The estimation accuracy of the proposed algorithm is higher than that of the ESPRIT algorithm, and the running time is shorter than that of the ESPRIT algorithm. Due to the extended receiving data and the unitary transformation operation of the proposed algorithm, the number of snapshots is as twice as the original number, and the complex data is transformed into the real data, which greatly reduces the computational complexity.

6. Conclusions

In this paper, a novel unitary ESPRIT algorithm is proposed for the angle and range estimation in a monostatic FDA-MIMO radar. In the proposed method, the angle and range are estimated by using the rotation invariance between the specific subarrays. Then, we make a specific analysis of periodic ambiguity and propose a method to solve that. Additionally, the computational complexity of the proposed algorithm is compared with that of the ESPRIT algorithm. The theoretical performance of the proposed algorithm is verified by computer simulation. In future work, we will focus on how to estimate the parameters of targets when mutual coupling errors exist in the FDA-MIMO radar, how to use the proposed algorithm in more general array structures, and how to use the proposed algorithm to estimate parameters in a colored noise environment.