DOA Estimation for Coprime Linear Array Based on MI-ESPRIT and Lookup Table

In order to improve the angle measurement performance of a coprime linear array, this paper proposes a novel direction-of-arrival (DOA) estimation algorithm for a coprime linear array based on the multiple invariance estimation of signal parameters via rotational invariance techniques (MI-ESPRIT) and a lookup table method. The proposed algorithm does not require a spatial spectrum search and uses a lookup table to solve ambiguity, which reduces the computational complexity. To fully use the subarray elements, the DOA estimation precision is higher compared with existing algorithms. Moreover, the algorithm avoids the matching error when multiple signals exist by using the relationship between the signal subspace of two subarrays. Simulation results verify the effectiveness of the proposed algorithm.


Introduction
DOA estimation is an important problem in array signal processing, and is widely used in radar, communication, sonar, and other detection equipment [1][2][3][4][5][6]. Traditional subspace-based methods, which include the multiple signal classification (MUSIC) algorithm [7] and the estimation of signal parameters via rotational invariance techniques (ESPRIT) [8,9], have been verified as efficient estimation techniques by using the eigenvalue decomposition of the received covariance matrix. Previous studies focused on the uniform array, such as the uniform linear array, uniform circular array, etc. Many DOA estimation algorithms have been proposed [10][11][12][13]. However, as the number of array elements and array aperture are restricted in practice, the uniform array is not the optimal array structure. The sparse array has attracted considerable attention because it obtains a larger array aperture without increasing the number of array sensors, thus producing better DOA estimation performance. The earliest sparse array is the minimum redundancy array (MRA). However, MRA cannot obtain the closed form expression, which makes it impossible to apply in practice. Subsequently, Vaidyanathan proposed the nested array and coprime array [14,15]. The application of the nested array is restricted by the mutual coupling of array elements. The coprime array is a non-uniform array system whose inter-element spacing is larger than half-wavelength. The coprime array processes spatial signals in a sparser array structure, which results in good angle measurement performance. Compared with the uniform array, it has a larger array aperture when the number of array elements is the same, and fewer array sensors are required when the array aperture is the same. Additionally, the mutual coupling of array elements is reduced, which weakens the influence on the DOA estimation performance.

DOA Estimation Based on the MI-ESPRIT Algorithm
The ESPRIT algorithm achieves DOA estimation by taking advantage of only a single displacement invariance in the sensor array. However, there are many situations where the subarray possesses multiple invariance structures. In order to make full use of the displacement invariance of the subarray and improve angle measurement precision, we used the MI-ESPRIT algorithm to estimate DOAs for the CLA. According to the reference [29], in order to obtain the multiple invariance structures, the th i subarray is divided into p arrays, and each array owns h elements. There are 1 h − overlapping elements between adjacent arrays. Therefore, the number of arrays and the number of array elements satisfy the following relationship: Since the distance between adjacent elements is unequal, there is no general steering vector expression. So, we constructed the output vector model from the perspective of two subarrays. For the ith subarray with M i elements, the output vector can be modeled as: where S(t) is the signal vector and N(t) is the additive white Gaussian noise vector. The steering vector A i can be expressed as: The array manifold of the kth signal is: where d i is the inter-element spacing and i = 1, 2. Then, we calculated the covariance matrix R X i X i of the output vector X i (t): The eigen-decomposition of R X i X i can be expressed as: where U S i is the signal subspace of the ith subarray.

DOA Estimation Based on the MI-ESPRIT Algorithm
The ESPRIT algorithm achieves DOA estimation by taking advantage of only a single displacement invariance in the sensor array. However, there are many situations where the subarray possesses multiple invariance structures. In order to make full use of the displacement invariance of the subarray and improve angle measurement precision, we used the MI-ESPRIT algorithm to estimate DOAs for the CLA. According to the reference [29], in order to obtain the multiple invariance structures, the ith subarray is divided into p arrays, and each array owns h elements. There are h − 1 overlapping elements between adjacent arrays. Therefore, the number of arrays and the number of array elements satisfy the following relationship: where M i denotes the number of the ith subarray's elements, as shown in Figure 2.
For the CLA, a unique non-singular matrix T exists: According to the MI structure of the subarray, we constructed a singular value decomposition (SVD)-like matrix i E From Equation (10), we see that 1 i Q and 1 i Q exist in the following relationship:  (12) can be modified as:  The signal subspace U S i and the steering vector of the ith subarray span the same space, i.e., For the CLA, a unique non-singular matrix T exists: According to the MI structure of the subarray, we constructed a singular value decomposition where U S ij = U S i (j : j + h − 1, :), (j = 1, 2, · · · , p) means extracting the jth row to the (j + h − 1)th row of U S i as a new matrix.
Define matrix E i 1 and E i 2 from E i : From Equation (10), we see that Q i 1 and Q i 1 exist in the following relationship: where ) is a diagonal matrix and ψ i contains the direction information of incoming signals. Thus, we obtain: Define (12) can be modified as: where E i 1 + denotes the Moore-Penrose pseudo-inverse of E i 1 .
Since T is non-singular matrix, Ω i and ψ i have the same eigenvalues. After completing eigen-decomposition on matrix Ω i , we obtain K eigenvalues λ i 1 , λ i 2 , · · · , λ i K . According to λ i k = e −j2πd i sin(θ k )/λ , the DOAs of K signals can be estimated: Sensors 2018, 18, 3043 5 of 11

Solving Ambiguity Based on the LUT
Since the inter-element spacing of each subarray was larger than a half-wavelength, the DOA estimation results were ambiguous. In Section 3.1, for every incoming signal, only an estimated value was obtained based on the MI-ESPRIT algorithm. First, we needed to calculate all estimated values according to the coprime property, which includes real DOA and ambiguous DOA. Then, by using the solving ambiguity method, real DOA can be obtained. The earliest solving ambiguity method obtains the real DOA by searching the nearest value from all DOA estimated values of two subarrays, which is computationally complex. Additionally, matching errors may occur when the incoming wave contains multiple signals. The reference [27] proposed a search-free solving ambiguity method, but solving the modular equations requires an iterative process, which still requires considerable calculation. Therefore, a solving ambiguity method with ultra-high speed is urgently required. In engineering applications, we usually take advantage of the LUT method to reduce the computational complexity.

Construct the LUT
For any signal, real DOAθ r and ambiguous DOAθ a satisfy: where M is the number of subarray elements and l m is a non-zero integer number between −(M − 1) and (M − 1). For subarray 1, M = M 2 . For subarray 2, M = M 1 .
In order to simplify the procedure of constructing the LUT, we performed the transformation as u = sin(θ), −1 ≤ u ≤ 1. Thus, Equation (15) can be written as: In the transformation domain, it can be seen from Equation (16) that the estimated values are uniformly distributed.
By calculating Equations (14)- (17), we obtain all DOA estimation values of the two subarrays. According to the reference [23], real DOAs uniquely exist, which satisfies the estimated values of the two subarrays. We considered one signal impinging on the CLA. Denote θ m 1 as the estimated values set of subarray 1 and denote θ m 2 as the estimated values set of subarray 2, where m 1 = 1, 2, · · · M 2 , m 2 = 1, 2, · · · , M 1 . By solving the minimum of the following formula, we can obtain the real DOA:θ For the given CLA, the number of subarray elements is determined. In the transformation domain, the estimated values of subarray 1 are uniformly distributed in For any incoming signal θ k ∈ [−π/2, π/2], every subinterval always has an estimated value, and all DOA estimated values can be obtained according to any subinterval estimated value. Select any subinterval of two subarrays as the reference interval, as shown in the solid line of Figure 3. So, we chose the reference interval of the two subarrays in the transformation domain to construct the LUT. In the reference interval, we set d s as the step. The smaller the d s , the higher the angle measurement precision. However, too small d s increases the size of the table, so we generally chose d s = 0.01 in practice. As the two subarrays traverse the entire reference interval at step d s , we can obtain the corresponding incident angle: where I 1 = 1, · · · , (R(2/M 2 d s ) + 1) and I 2 = 1, · · · , (R(2/M 1 d s ) + 1), R(·) represents the rounding symbol.
For any incoming signal [ 2,2] k θ π π ∈ − , every subinterval always has an estimated value, and all DOA estimated values can be obtained according to any subinterval estimated value. Select any subinterval of two subarrays as the reference interval, as shown in the solid line of Figure 3. So, we chose the reference interval of the two subarrays in the transformation domain to construct the LUT.
In the reference interval, we set s d as the step. The smaller the s d , the higher the angle measurement precision. However, too small s d increases the size of the table, so we generally chose 0.01 s d = in practice. As the two subarrays traverse the entire reference interval at step s d , we can obtain the corresponding incident angle:  However, there is no consistent one-to-one matching relationship between the estimated values of the two subarrays. In order to perform DOA estimation based on the above LUT, it was necessary to find the pairing relationship between the estimated values of the two subarrays.
This problem has been mentioned by the reference [25]. We defined 1 H and 2 H according to the signal subspace of the two subarrays.

H U U A TT A
which means Firstly, the estimated values set to {θ 1 (I 1 )} and {θ 2 (I 2 )}, which correspond to θ 1 (I 1 ) and θ 2 (I 2 ), are calculated based on Equations (16) and (17), respectively. Then, substituting {θ 1 (I 1 )} and {θ 2 (I 2 )} into Equations (18) and (19), respectively, we obtain the real DOAs corresponding to {θ 1 (I 1 )} and {θ 2 (I 2 )}. Finally, the real DOA of each pair I 1 and I 2 is stored in the table. We then constructed the LUT by traversing I 1 and I 2 . Since the LUT is only constructed in the reference interval, the table is smaller.

Solving Real DOA Based on the LUT
The established LUT in the above section assumes there only one signal exists. When K(K ≥ 2) signals impinge on the CLA, each subarray can have K estimated values based on MI-ESPRIT. However, there is no consistent one-to-one matching relationship between the estimated values of the two subarrays. In order to perform DOA estimation based on the above LUT, it was necessary to find the pairing relationship between the estimated values of the two subarrays.
This problem has been mentioned by the reference [25]. We defined H 1 and H 2 according to the signal subspace of the two subarrays.
For any set of estimated valuesθ 1k andθ 2k of the two subarrays, we first calculated their corresponding u 1k and u 2k in the transformation domain, and then obtained their index values in the table.
Substituting I 1k and I 2k into the LUT, the real DOAθ k can be solved. The complete calculation procedure of the algorithm is shown in Figure 4. The proposed method realizes DOA estimation of multiple signals, by only using one signal to construct the LUT, which is relatively simple and convenient for practical application. In practical engineering, the proposed algorithm in this paper could be implemented on the hardware platform composed of the digital signal processor (DSP) and field-programmable gate array (FPGA), where the LUT is built in the microwave darkroom in advance and stored in the RAM module of the FPGA chip.
Substituting 1k I and 2k I into the LUT, the real DOA ˆk θ can be solved. The complete calculation procedure of the algorithm is shown in Figure 4. The proposed method realizes DOA estimation of multiple signals, by only using one signal to construct the LUT, which is relatively simple and convenient for practical application. In practical engineering, the proposed algorithm in this paper could be implemented on the hardware platform composed of the digital signal processor (DSP) and field-programmable gate array (FPGA), where the LUT is built in the microwave darkroom in advance and stored in the RAM module of the FPGA chip.

Simulation Analysis
Consider a CLA consisting of two uniform linear subarrays with M 1 = 13 and M 2 = 11 elements, and inter-element spacing of d 1 = M 2 λ/2 and d 2 = M 1 λ/2, respectively. In MI-ESPRIT, we take h = 4, and then the subarray 1 can be divided into p M 1 = 10 arrays and subarray 2 can be divided into p M 2 = 8 arrays. Two signals imping on the CLA from −20 • , 20 • . In order to verify the performance of the proposed algorithm, we performed comparison simulations between the previously proposed methods [23,25,28] and our method. In the reference [23], the rough step was d S1 = 0.1 • and the fine step was d S2 = 0.2 • . Moreover, in order to illustrate the advantage of CLA, a comparison with the uniform linear array (ULA) was added to the simulation, in which the number of array elements was consistent with the CLA, i.e., M 1 + M 2 − 1. We performed 1000 independent Monte Carlo simulations. We recorded 512 snapshots and Figure 5a shows that the DOA estimation root mean square error (RMSE) of the six algorithms varied with the signal-to-noise ratio (SNR). When the SNR = 10 dB, the DOA estimation RMSE of the six algorithms versus snapshots is shown in Figure 5b.
comparison with the uniform linear array (ULA) was added to the simulation, in which the number of array elements was consistent with the CLA, i.e., 1 2 1 M M + − . We performed 1000 independent Monte Carlo simulations. We recorded 512 snapshots and Figure 5a shows that the DOA estimation root mean square error (RMSE) of the six algorithms varied with the signal-to-noise ratio (SNR). When the SNR 10 = dB, the DOA estimation RMSE of the six algorithms versus snapshots is shown in Figure 5b. With increasing SNR or snapshots, the DOA estimation RMSE of the algorithms decreased rapidly. Compared with the ULA, the DOA estimation precision of CLA was significantly better. Compared to both CLA and ULA, the DOA estimation precision of MI-ESPRIT was better than ESPRIT. This is because the MI-ESPRIT algorithm makes full use of the subarray elements, so the angle measurement precision of the proposed algorithm is better than that of the reference [28], which was based on the ESPRIT algorithm. The DOA estimation precision of the reference [25] was between that of the reference [28] and the proposed algorithm. It can be seen from Figure 5 that the DOA estimation precision of the algorithm introduced by the reference [23] is basically consistent with our proposed method when SNR is low. With increasing SNR, the DOA estimation precision of the algorithm of the reference [23] was poor compared to our proposed method. Because the angle measurement accuracy in the reference [23] is closely related to the precision of the fine spectrum search step dS2, if we continued to reduce dS2, the DOA estimation precision of the reference's algorithm [23] may be better at high SNR, but the high computational complexity caused by the fine spectrum search would be unfeasible. In summary, the proposed algorithm possesses the best DOA estimation accuracy.
In order to verify the angular resolution of the proposed method, other simulation conditions remained unchanged and the DOAs of two signals were reduced to −10°, 10°; −5°, 5° and −0.5°, 0.5°. The SNR 5 = dB, SNR 10 = dB, and SNR 15 = dB DOA estimation results are shown in Table 1. It can be seen from Table 1 that the proposed method could still distinguish two signals, although the DOA spacing was 1°. This is because the MI-ESPRIT algorithm has higher DOA estimation precision. With increasing SNR or snapshots, the DOA estimation RMSE of the algorithms decreased rapidly. Compared with the ULA, the DOA estimation precision of CLA was significantly better. Compared to both CLA and ULA, the DOA estimation precision of MI-ESPRIT was better than ESPRIT. This is because the MI-ESPRIT algorithm makes full use of the subarray elements, so the angle measurement precision of the proposed algorithm is better than that of the reference [28], which was based on the ESPRIT algorithm. The DOA estimation precision of the reference [25] was between that of the reference [28] and the proposed algorithm. It can be seen from Figure 5 that the DOA estimation precision of the algorithm introduced by the reference [23] is basically consistent with our proposed method when SNR is low. With increasing SNR, the DOA estimation precision of the algorithm of the reference [23] was poor compared to our proposed method. Because the angle measurement accuracy in the reference [23] is closely related to the precision of the fine spectrum search step d S2 , if we continued to reduce d S2 , the DOA estimation precision of the reference's algorithm [23] may be better at high SNR, but the high computational complexity caused by the fine spectrum search would be unfeasible. In summary, the proposed algorithm possesses the best DOA estimation accuracy.
In order to verify the angular resolution of the proposed method, other simulation conditions remained unchanged and the DOAs of two signals were reduced to −10 • , 10 • ; −5 • , 5 • and −0.5 • , 0.5 • . The SNR = 5 dB, SNR = 10 dB, and SNR = 15 dB DOA estimation results are shown in Table 1. It can be seen from Table 1 that the proposed method could still distinguish two signals, although the DOA spacing was 1 • . This is because the MI-ESPRIT algorithm has higher DOA estimation precision. Additionally, in the process of constructing the lookup table, we set the step d s = 0.01 in the reference interval, which guarantees high estimation accuracy. Higher DOA estimation precision means higher angular resolution. Therefore, the angular resolution of the proposed method can achieve 1 • . The computational complexity of the various algorithms are analyzed in Table 2. The CLA consists of two uniform linear subarrays with M 1 and M 2 sensors. Consider that K signals impinge on the CLA and T snapshots are used, and the number of searches in the reference [23] is set as S. It can be seen from Table 2 that the computational complexity of each algorithm is mainly created by two parts. The first part obtains the estimated values, which include covariance matrix estimation, eigenvalue decomposition, and solving estimated values, so the computational burden of the four algorithms is different. The second part is the solving ambiguity; the computational complexities in the previously formulated algorithms [23,25,28] are the same, i.e., K 2 M 1 M 2 , which is larger than the complexity in our method. In the proposed algorithm, by simply indexing in the LUT, the real DOA can be obtained. In general, the computational burden of the first part is larger than the second part. Table 2. Comparison of computation complexity.

Algorithm
Computation Complexity In the reference [23], the computational burden is caused by the spectrum peak search process and S is usually much larger than other variables, so CLA-Decom-MUSIC [23] has the highest computational complexity. CLA-Root-MUSIC [25] estimates the covariance matrix and performs the eigenvalue decomposition by combining the two subarrays, which increases the computational complexity, and the polynomial root finding is very time-consuming in practice. Therefore, CLA-Root-MUSIC [25] also has higher computational complexity, which second to CLA-Decom-MUSIC [23]. Both CLA-ESPRIT [28] and our proposed algorithm are essentially based on the ESPRIT algorithm, whose computational complexity is lower than CLA-Decom-MUSIC [23] and CLA-Root-MUSIC [25]. For the ESPRIT algorithm, because the computational burden of solving estimated values can be ignored, its computational burden of the first part mainly includes covariance matrix estimation and eigenvalue decomposition of both the covariance matrix R X i X i and matrix Ω i , and the corresponding complexities are (M 1 2 + M 2 2 )T, (M 1 3 + M 2 3 ) and (3(M 1 + M 2 ) + 4K)K 2 , respectively. The MI-ESPRIT algorithm uses the multiple invariance structure of ESPRIT, which only adds the linear transformation of the matrix compared to ESPRIT. The added computational burden can also be neglected. Therefore, the proposed algorithm has approximately similar computational complexity in terms of the first part compared with CLA-ESPRIT [28]. However, we simplified the solving ambiguity by using a LUT, which substantially reduces the computational burden, i.e., K K 2 M 1 M 2 . Therefore, the proposed algorithm is more efficient. According to the above simulation conditions, M 1 = 13 remains unchanged and M 2 is 8, 10, 12, 14, and 16. We compared the processing time of the four different algorithms in Figure 6. The processing time was determined by a PC (Lenovo manufactory, Beijing, China) with AMD Phenom™ IIX6 1055T Processor 2.8 GHz CPU and 8 GB RAM by running the MATLAB codes in the same environment. It can be seen from Figure 6 that the proposed algorithm has the highest computational efficiency.
Sensors 2018, 18, x FOR PEER REVIEW 10 of 12 processing time was determined by a PC (Lenovo manufactory, Beijing, China) with AMD Phenom™ IIX6 1055T Processor 2.8 GHz CPU and 8 GB RAM by running the MATLAB codes in the same environment. It can be seen from Figure 6 that the proposed algorithm has the highest computational efficiency.

Conclusions
CLA has been widely studied due to its superior DOA estimation performance. This paper proposed a novel DOA estimation algorithm for CLA based on MI-ESPRIT and LUT. MI-ESPRIT

Conclusions
CLA has been widely studied due to its superior DOA estimation performance. This paper proposed a novel DOA estimation algorithm for CLA based on MI-ESPRIT and LUT. MI-ESPRIT fully uses the subarray's elements, which improves the angle measurement accuracy. Then, according to the property of the CLA, the phase ambiguity was solved using the LUT, which reduced the computational complexity. At the same time, using the relationship between the signal subspace of two subarrays, matching errors were avoided when in the presence of multiple signals. Compared with the existing algorithms, the proposed method not only has higher DOA estimation accuracy, but also has lower computational complexity. Additionally, our DOA estimation method, which is based on the LUT, has broad application prospects in practice. However, in our study, the coprime array achieved DOA estimation by decomposing CLA into two uniform linear subarrays, which sacrifices the degrees of freedom, i.e., a reduction in the number of sources that the CLA can resolve. By increasing the number of CLA sensors, we could obtain increased DOFs. In addition, the virtualization array sensor method, using the Khatri-Rao transformation, could also be applied to the CLA to yields the virtual array structure. Based on the extended virtual array aperture, the DOFs can also be increased. Therefore, future research efforts will aim to improve the DOF of CLA.
Author Contributions: W.Z. proposed the original writing idea of the full text; K.C. designed the simulation experiment and implemented the algorithm in MATLAB with W.Z.; X.C. and T.X. wrote the manuscript under the guidance of N.Y.; W.Z. and N.Y. revised the manuscript.
Funding: This research received no external funding.