A Novel Unitary ESPRIT Algorithm for Monostatic FDA-MIMO Radar

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.


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.

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 Sensors 2020, 20, 827 3 of 17 where f 0 is the referenced frequency, and ∆f is the deviation of adjacent transmitting frequencies, where ∆f << f 0 . The baseband signals emitted by the different elements are mutually orthogonal. Assume that the emitted signal of m-th element satisfies where E is the total transmitted energy, T is the radar pulse duration, and φ m (t) denotes the baseband waveform. Suppose that all baseband signals are orthogonal to each other and normalized 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 , 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. 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(θ k ) = [1, e jπ sin θ k , · · · , e j(N−1)π sin θ k ] T ∈ C N×1 Then, the received data at the receiving array can be described as [11] β k e j2π f pk t b(θ k )a T (r k , θ k ) 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] β 2 e j2π f p2 t , . . . , β k e j2π f pk t 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 σ 2 I MN , where σ 2 and I MN denote the noise variance and M × N identity matrix, respectively.

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] π θ θ θ are selection matrices, and 0w denotes the w × w null matrix. The functions of J1 and J2 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 Under the assumption of K targets, the rotation invariance of each target can be written into a matrix form

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] where are selection matrices, and 0 w denotes the w × w null matrix. The functions of J 1 and J 2 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 Under the assumption of K targets, the rotation invariance of each target can be written into a matrix form Sensors 2020, 20, 827 where Ξ 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 where E s is composed of the K eigenvectors corresponding to the largest K eigenvalues of the covariance matrix of Y, and Θ is a non-singular matrix. Substituting Equation (16) into Equation (14), we can obtain the relationship between signal space of the subarrays where It can be noticed that the diagonal elements of Ξ R are the eigenvalues of Ψ R .

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] 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] where (.) H represents the conjugate transpose operator, and the sparse unitary matrix Qw is defined as 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 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 where 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 where Ξ T contains the ranges of all the targets. According to Equation (16), Equation (19) can be rewritten as where Ψ T = Θ −1 Ξ T Θ. 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.

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] 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] where (.) H represents the conjugate transpose operator, and the sparse unitary matrix Q w is defined as 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 where R Y and R Z are the covariance calculated by Y and Z, respectively. The signal subspaceÊ S corresponds to K eigenvectors of large eigenvalues of R Γ . The remaining MN − K eigenvectors of small eigenvalues can obtain the noise subspaceÊ N . Hence,Ê S andÊ 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 where d k = Q H MN 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 where D = [d 1 , d 2 , . . . , 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 where Θ R is the left eigenvector matrix of Σ R . By using the total least squares (TLS) method to solve Equation (31), DOA can be estimated as followŝ Similarly, the rotation invariance between Subarray 1 and Subarray 2 can be transformed into For K targets, Equation (33) can be integrated into matrix form 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 where Σ T = Θ −1 T Φ T Θ T . Θ T is the inverse of the left eigenvector matrix of Σ T . Substituting Equation (32) into Equation (37), the range estimation is solved with the TLS method.
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.

The Pairing of DOAs and Ranges
In this section, we analyze the speciality of Θ T and Θ R , and introduce the means to achieve pairing. Θ T and Θ R are eigenvectors of Σ T and Σ R , respectively. Since Σ T and Σ R are calculated byÊ S , there must be a random row of Θ R identical to a specific row of Θ T . Supposing that all of the K targets are independent, we notice that any two rows of Θ R are orthogonal because any two eigenvalues of Θ T are different. In this paper, considering algorithm complexity, we obtain the automatic pairing of Φ T and Φ R by decomposing the Σ T + jΣ R , which can be expressed as where Θ TR is the left eigenvector matrix. Hence, Φ T and Φ R can be automatically paired by the eigenvector matrix Θ TR . Considering the periodic phase ambiguity problem, we take a step to distinguish the phase ambiguity before calculating ranges.

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 where R H is the covariance of H in Equation (9), and L is the number of snapshots, and σ 2 denotes the noise power. Under the assumption of K targets, the unknown parameter to be estimated is Then, the Fisher information matrix (FIM), with respect to η, is [41] The expression of every block of F can be written as where C θ and C r are partial derivations of C with respect to θ and r, and The C θ and C r can be written as where We derive part of Equations (54) and (55) as and Then, every block of F is determined by Equations (48)-(51). Then, the CRB matrix can be obtained by where W = [C θ C r ], P = P 1 P 2 P 3 P 4 and P 1 = P 2 = P 3 = P 4 = R H C H R −1 Y CR H .

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) 2 } 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) 3 } flops. The complexity required to solve for Σ R is O{M(N − 1)(2K) 2 }, where K is the number of targets. Similarly, solving Σ T needs O{N(M − 1)(2K) 2 } flops. The eigenvalue decomposition in Equation (40) and the pairing of angle and range need O{4(2K) 3 + 2K 3 }. Here, we ignore the complexity of solving periodic ambiguity steps because they are too small. Thus, the complexity of the proposed algorithm is In [42], the complexity of the ESPRIT algorithm for estimation of angle and range is 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.

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 0 , namely the minimum frequency, is 3 GHz, and the frequency increment ∆f is 10 3 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.

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 Σ R and Σ 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

RMSE Versus SNR
In this section, the target is set to (45 • 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.

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 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. (a)

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. (a)

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.

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.
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.

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 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.

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.