A Computationally Efficient and Virtualization-Free Two-Dimensional DOA Estimation Method for Nested Planar Array: RD-Root-MUSIC Algorithm

To address the problem of expensive computation in traditional two-dimensional (2D) direction of arrival (DOA) estimation, in this paper, we propose a 2D DOA estimation method based on a reduced dimension and root-finding MUSIC algorithm for nested planar arrays (NPAs). Specifically, the algorithm proposed in this paper transforms the problem based on 2D spectral peak search into two one-dimensional estimation problems by reducing the dimension, and then transforms the one-dimensional estimation problem into a problem of polynomial root finding. Finally the parameters are paired to realize the 2D DOA estimation. The proposed algorithm not only performs two root finding operations directly according to the 2D spectral function transformation, avoiding the performance degradation caused by intermediate operations, but can also fully exploit the enlarged array aperture offered by NPAs with reduced computational complexity and no need for virtualization. The superiorities of the proposed algorithm in terms of estimation accuracy and complexity are verified by simulations.


Introduction
Two-dimensional direction of arrival (2D-DOA) estimation has attracted extensive attention over the years and has been applied to wireless communication, sonar, radar, and many other fields [1][2][3]. There have been increasing number of algorithms proposed for two-dimensional DOA estimation such as estimation of signal parameters via rotational invariance techniques (ESPRIT) [4][5][6], multiple signal classification (MUSIC) [7], parallel factor (PARAFAC), and the propagator method (PM) [8] technique based on various array geometries, for example, L-shaped arrays [9,10], two parallel linear arrays [11], uniform planar arrays (UPA) [12][13][14], and uniform circular arrays [4]. In addition, traditional 2D-DOA methods mainly include spectral peak search [7,[15][16][17][18][19] and rotation invariance algorithm [4][5][6][7][8]. Among them, although the 2D MUSIC algorithm based on spectral peak search [7] has high DOA estimation accuracy, 2D total spectrum search (TSS) has great computational cost. Rotation invariance algorithm such as 2D ESPRIT can directly calculate DOA estimation, avoiding the complex calculation of 2D spectral peak search, but cannot achieve the estimation performance of spectral peak search. Of course, there are many new DOA estimation methods based on information geometry [20][21][22], which will not be discussed in this paper. Ref. [23] proposes a MUSIC-based dimensionality reduction and root-finding algorithm for uniform planar arrays (RD-ROOT-UPA), which further reduces the amount of computation by reducing the dimension of the spectral function and then finding the root. However, in order to avoid spatial aliasing, these traditional arrays require adjacent elements to be no more than half a wavelength apart, which induces a heavy 1.
We use NPAs instead of UPAs and CPAs to further expand the array aperture and obtain better estimation results; 2.
We use the reduced dimension and root-finding MUSIC algorithm for 2D DOA estimation. The proposed algorithm changes the 2D spectral search into two 1D estimation problems and then transforms the 1D estimation problem into a root finding problem of a 1D polynomial, which reduces the complexity while maintaining the estimation accuracy; 3.
We regard the nested array as extracted from a uniform array with the same array aperture and use the polynomial root finding method to deal with the problem, avoiding the virtualization process of NPA.
We summarize this paper below. Section 2 recommends the data model of NPA. The proposed algorithm is expounded in Sections 3 and 4, which analyze the complexity and DOFs of the proposed algorithm. Section 5 gives simulation results to verify the efficiency of the proposed algorithm, and Section 6 is a conclusion of the whole paper.

Notations
Bold uppercase characters denote matrices, while the lowercase characters mean vectors. (·) T and (·) H represent the transpose and conjugate transpose; (·) * and (·) −1 Sensors 2022, 22, 5220 3 of 13 denote conjugate operation and inverse, respectively. and ⊗ are the Khatri-Rao product and Kronecker product, respectively. angle(·) represents the phase operator. Rank(·) is the rank of the matrix. det(·) means the determinant of the matrix.

Data Model
Consider an NPA, which is a 2D sparse planar array extended configuration of a twolevel nested linear array (NLA) with M 1 and M 2 . In other words, each row and column of the NPA is a two-level NLA, which is composed of two uniform linear arrays in series. The element spacing of the first level uniform linear array is d 1 = λ/2 with M 1 elements, while that of the second level is d 2 = (M 1 + 1)λ/2 with M 2 elements, where λ is the wavelength. Thus, the total number of array sensors is Specifically, the location set of sensors in the NPA can be represented by Figure 1 displays the structural layout of an NPA, taking M 1 = 2, M 2 = 3, T = 25 as an example.

Notations
Bold uppercase characters denote matrices, while the lowercase characters mean vectors. () T  and () H  represent the transpose and conjugate transpose; ()   and 1 () −  denote conjugate operation and inverse, respectively. and  are the Khatri-Rao product and Kronecker product, respectively.
() angle  represents the phase operator.
() Rank  is the rank of the matrix. det( )  means the determinant of the matrix.

Data Model
Consider an NPA, which is a 2D sparse planar array extended configuration of a twolevel nested linear array (NLA) with 1 M and 2 M . In other words, each row and column of the NPA is a two-level NLA, which is composed of two uniform linear arrays in series. The element spacing of the first level uniform linear array is M elements, while that of the second level is is the wavelength. Thus, the total number of array sensors is Specifically, the location set of sensors in the NPA can be represented by , mm . Figure 1 displays the structural layout of an NPA, taking For the NPA, the received signal can be modeled as [15]   Suppose K narrowband far-field uncorrelated signals impact the NPA. The K signals are (θ k , φ k ), k = 1, 2, · · · , K, where φ k and θ k are the azimuth and elevation angles of the kth signal, respectively. Define u k = sin θ k sin φ k and v k = sin θ k cos φ k for simplification.
For the NPA, the received signal can be modeled as [15] where S = [s 1 , s 2 , · · · , s K ] T ∈ C K×L is the signal matrix, S k represents the source vector and s k = [s k (1), s k (2), · · · , s k (L)] T ∈ C L×1 , L indicates the number of snapshots, A denotes the direction matrix of the array, and A = [a y ( and a y (u k ) denote the direction vectors along the x-axis and y-axis, respectively. They are indicated as a y (u k ) = [1, e j2πd 1 u k /λ , · · · , e j2πM 1 d 1 u k /λ , e j2π(M 1 d 1 +d 2 )u k /λ , · · · , and a x (v k ) = [1, e j2πd 1 v k /λ , · · · , e j2πM 1 d 1 u k /λ , . N denotes the received noise of the array, which is white Gaussian noise with zero mean and variance σ 2 . Actually, the covariance matrix of the received signal X with L snapshots can be approximated byR The Eigenvalue decomposition (EVD) of the covariance matrixR is expressed bŷ whereΛ s andΛ n are diagonal matrixes where the diagonal elements are denoted by the K largest Eigenvalues and (T − K) smallest Eigenvalues ofR, respectively. Signal subspacê E s contains the Eigenvectors of the K largest Eigenvalues and noise subspaceÊ n comprises the (T − K) smallest remaining Eigenvectors.

The Proposed Algorithm
In this section, different from other methods on NPAs that require virtualization with high computational complexity, we first transform the 2D spectral search into two 1D spectral search problems. Then, we exploit the 1D polynomial root finding method to estimate the two parameters. Finally, after completing the parameter matching, we can estimate the values of elevation and azimuth to complete DOA estimation.

D-MUSIC Spectrum Function
According to the orthogonal property between the signal steering vectors and the noise subspace, we form the spectrum function of the NPA as [29,32] represents the direction vector of NPA, is the noise subspace of NPA. The DOA estimation can be obtained by a 2D spatial spectral search according to (4). Although it can be paired automatically, it requires expensive computation. To overcome this problem, we transform the 2D spectral search into two 1D spectral search problems, and then estimate the two parameters by exploiting the 1D polynomial root finding method.

Root-Finding Process of Reduced-Dimensional Polynomial
Construct polynomials based on the spectrum function (4) as and where According to the rank relationship of the matrix product, the condition below must be met Sensors 2022, 22, 5220 5 of 13 and then we have Obviously, on the basis of (8), we can conclude that det{Q(u)} is a nonzero polynomial; hence, Q(u) is a factor of V(u, v). As Q(u) is only related to u, the roots of det{Q(u)} = 0 can establish the following formula: The estimates of u and v can be expressed by the problem of 1D polynomial root finding. Therefore, the pairwise estimates of u and v obtained from 2D spectral research are changed into 1D root finding processes. After that, reconstruct (5) and (6) as Define where d = λ/2. Define the direction vector along the x-axis for the UPA with the same array aperture size as the NPA: the same array aperture between the NPA and UPA, we have where G ∈ Z (M 1 +M 2 )×((M 1 +1)M 2 ) . Specifically, g ij = 1 holds when the ith sensor in the a x (v) overlaps with the jth sensor in the a Ex (v), or else g ij = 0, where g ij is the (i, j)th element of G. For the NPA shown in Figure 1, G is a 5 × 9 matrix, in which four columns are zeros. Furthermore, we can also obtain the relation a y (u) = Ga Ey (u). Accordingly, the direction vectors can be further expressed as Without losing generality, we can substitute z Namely, According to the above formula,û k andv i are acquired from the K-root distribution nearest to the unit circle based on (16) and (17), which are expressed asẑ 11 , · · · ,ẑ 1k , · · · ,ẑ 1K andẑ 21 , · · · ,ẑ 2k , · · · ,ẑ 2K , respectively, i.e., Among the traditional DOA estimation methods with sparse arrays, the array element separation of the two subarrays is larger than half wavelength; as a result, it is necessary to eliminate the ambiguity. According to the structure of the NPA, since it has at least one set of adjacent sensors spaced no more than half a wavelength apart, it is an unambiguous array and can avoid additional disambiguation operations.
Sinceû k andv i in the above operation are obtained by finding the roots of two 1D polynomials, they are not paired with each other, so we need to pair these unpaired parameters. After the parameters are paired, we can directly obtain the real DOA estimation.

Parament Pairing and DOA Estimation
Since the two root-finding processes ofû k andv i are carried out separately, in this part, we pair them up. Therefore, a cost function is given as where a(û k ,v i ) denotes the direction vector constructed byû k andv i , which can be obtained on the basis of (1). For anyû k , we can calculate the values of i and k when V k,i (1 ≤ i ≤ K) takes the minimum value, and then we define this i paired with k as i . Finally, the elevation and azimuth angles can be calculated bŷ

The Main Steps of the Proposed Algorithm
We provide a summary of the main steps of the proposed algorithm below: Step 1: Calculate the covariance matrixR of the received signal X; the noise subspace E n is obtained by EVD; Step 2: Reduce the dimension of constructed polynomial V(u, v) according to (4); Step 3: Estimateû k andv i by the 1D polynomial root finding method; Step 4: Match the parametersû k andv i based on (20); Step 5: Calculateθ k andφ k on the basis of (21) and (22).

Complexity Analysis
In this part, we compare the computational complexity of the proposed algorithm with other algorithms, including the ESPRIT algorithm for NPAs [34] and some other MUSIC based algorithms such as the 2D-PSS [15] and RD-ROOT-UPA [23] methods for UPAs. For the proposed algorithm, it takes O T 2 L to calculate the covariance matrix and the computational complexity of EVD is O T 3 . The polynomial root-finding process needs O 2(2((M 1 + 1)M 2 ) 2 ) 3 and the parameter matching process costs O K 2 (T − K)(T + 1) . Therefore, the total complexity According to the computational complexity of the above algorithms listed in Table 1 in which ∆ = 0.0001 is the spectral search interval, M = M 1 + M 2 represents the number of elements in each row and column of the uniform planar array (UPA). Figure 2 shows the computational complexity of different methods with the increase of sensors, where K = 2 and L = 200. It can be observed from Figure 2 that the complexity of the proposed algorithm is significantly lower than that of 2D-PSS, because the computational efficiency of polynomial root finding is much higher than that of spectral search. Meanwhile, it can be observed that the complexity of the proposed algorithm is similar to that of the ESPRIT algorithm, but higher accuracy can be obtained by finding the root of the spectral function. Finally, although the complexity of the proposed algorithm is higher than that of the RD-ROOT-UPA algorithm, the larger array aperture of the NPA compared to the UPA gives better estimation accuracy with an equivalent number of sensors.

Algorithm Computational Complexity
Proposed

Cramer-Rao Bound
The Cramer-Rao bound (CRB) is a lower bound of the error variance of a parameter estimator. In this part, we derive the CRB of 2D DOA estimation for NPAs.
According to [35], the CRB of an NPA can be given by where

Achievable DOFs
For the traditional 2D DOA estimation algorithms of a UPA, such as 2D-PSS, the DOF is represents the number of elements in each row and column of the UPA. According to [32], the RD-ROOT-UPA algorithm can distinguish ( 1) MM − signals at most, which reduces the dimension and directly finds the roots of the spectral function of the UPA. The proposed algorithm utilizes an RD root finding technique to reduce the computational complexity, and we can obtain the achievable DOFs on the basis of (7) as follows:

Cramer-Rao Bound
The Cramer-Rao bound (CRB) is a lower bound of the error variance of a parameter estimator. In this part, we derive the CRB of 2D DOA estimation for NPAs.
According to [35], the CRB of an NPA can be given by ;P S = SS H /L, a k represents the kth column of A, and the variance of the received noise is σ 2 .

Achievable DOFs
For the traditional 2D DOA estimation algorithms of a UPA, such as 2D-PSS, the DOF is M 2 − 1, where M = M 1 + M 2 represents the number of elements in each row and column of the UPA. According to [32], the RD-ROOT-UPA algorithm can distinguish M(M − 1) signals at most, which reduces the dimension and directly finds the roots of the spectral function of the UPA. The proposed algorithm utilizes an RD root finding technique to reduce the computational complexity, and we can obtain the achievable DOFs on the basis of (7) as follows: Therefore, the proposed algorithm can identify the maximum number of signals as We provide the achievable DOFs of different methods in Table 2. Compared with the traditional 2D DOA estimation methods, it is clear that the proposed algorithm is able to achieve relatively high DOFs. Table 2. Achievable DOFs of different methods.

Algorithm
Achievable DOFs

Advantages
On account of the above analysis, the main advantages of the proposed algorithm are summarized as below: 1.
The proposed algorithm transforms the problem of 2D spectral search into two 1D polynomial root finding problems, which greatly reduces the complexity compared with the traditional 2D-MUSIC algorithms such as the 2D-PSS algorithm; 2.
Compared with other 2D DOA estimation method for NPAs, the proposed algorithm avoids virtualization steps and spatial smoothing, so it has lower complexity than ESPRIT-NPA; 3.
The proposed algorithm fully utilizes the expanded aperture and DOF of the NPA, allowing it to achieve a higher DOF and array aperture compared to RD-ROOT-UPA and 2D-PSS algorithms based on uniform arrays.

Simulations
In this section, assuming that there are K uncorrelated far-field narrowband signals impinging on the NPA, we provide extensive simulations to verify the performance of the proposed algorithm in terms of the root mean square error (RMSE) under various conditions. According to [33], the RMSE is expressed as follows: where C = 500 is the number of trials;φ k,i andθ k,i are the estimated values of the kth signal in the ith test that correspond with the true φ k and θ k , respectively.

Scatter Figure
In this part, we draw the scatter figure of K sources estimated by the proposed algorithm in Figure 3

RMSE Results of the Proposed Algorithm Versus Snapshots
Here, we present the RMSE result of the proposed algorithm with different snapsh in Figure 4

RMSE Results of the Proposed Algorithm Versus Number of Sensors
Herein, Figure 5 shows the relationship between the RMSE results of the propo algorithm and the number of sensors, where

RMSE Results of the Proposed Algorithm Versus Snapshots
Here, we present the RMSE result of the proposed algorithm with different snaps in Figure

RMSE Results of the Proposed Algorithm Versus Number of Sensors
Herein, Figure 5

RMSE Results of the Proposed Algorithm Versus Number of Sensors
Herein, Figure 5 shows the relationship between the RMSE results of the proposed algorithm and the number of sensors, where (θ 1 , φ 1 ) = (20 • , 30 • ), (θ 2 , φ 2 ) = (40 • , 50 • ), and L = 200. As depicted in Figure 5, with the number of array elements increasing, the DOA estimation performance improves because of diversity gain of the receiving antenna increases.         Figure 6 shows that all the algorithms obtain better estimation results with the increase of snapshots; when the SNR is 10 dB, the performance of the 2D-PSS algorithm is better than that of RD-ROOT-UPA and ESPRIT-NPA because the performance of spectral peak search will be better than that of spectral function root finding method. It will also be better than the ESPRIT method based on rotation invariance. The proposed algorithm achieves the best estimation performance because the proposed method can exploit the extended array aperture and increased DOF offered by NPAs.

RMSE Comparison of Different Algorithms
The result of Figure 7 shows that the proposed algorithm and ESPRRIT-NPA algorithm outperform other algorithms for UPAs when the SNR is low, profiting from the larger array aperture of NPA, and the proposed algorithm performs two root finding processes directly according to the Formula (4), avoiding the performance degradation caused by intermediate operations. The 2D-PSS algorithm obtains better performance with spectral peak search than the RD-ROOT-UPA algorithm which directly finds the root of the same uniform array. As the SNR increases gradually, the performance of ESPRRIT-NPA becomes worse than that of other algorithms, but the proposed algorithm is still the best.

Conclusions
In this paper, we proposed a 2D DOA estimation method for NPAs with high computational efficiency and no virtualization by using the root finding method of reduced dimension polynomials. The proposed algorithm converts the 2D spectral peak search into two 1D polynomial root finding problems, and then after parameters pairing, 2D DOA estimation is realized, which is directly processed by the root-finding MUSIC method, avoiding the process of virtualization and spatial smoothing as well as the huge computational load brought by spectral peak search while maintaining estimation accuracy. Simulation results verified that the proposed method outperforms many exiting methods in terms of computational complexity and estimation accuracy. In future research, we will consider using this dimensionality reduction and root-finding method on improved arrays with larger array apertures, or consider combining this method with virtual array technology to obtain greater DOFs.  Figure 6 shows that all the algorithms obtain better estimation results with the increase of snapshots; when the SNR is 10 dB, the performance of the 2D-PSS algorithm is better than that of RD-ROOT-UPA and ESPRIT-NPA because the performance of spectral peak search will be better than that of spectral function root finding method. It will also be better than the ESPRIT method based on rotation invariance. The proposed algorithm achieves the best estimation performance because the proposed method can exploit the extended array aperture and increased DOF offered by NPAs.
The result of Figure 7 shows that the proposed algorithm and ESPRRIT-NPA algorithm outperform other algorithms for UPAs when the SNR is low, profiting from the larger array aperture of NPA, and the proposed algorithm performs two root finding processes directly according to the Formula (4), avoiding the performance degradation caused by intermediate operations. The 2D-PSS algorithm obtains better performance with spectral peak search than the RD-ROOT-UPA algorithm which directly finds the root of the same uniform array. As the SNR increases gradually, the performance of ESPRRIT-NPA becomes worse than that of other algorithms, but the proposed algorithm is still the best.

Conclusions
In this paper, we proposed a 2D DOA estimation method for NPAs with high computational efficiency and no virtualization by using the root finding method of reduced dimension polynomials. The proposed algorithm converts the 2D spectral peak search into two 1D polynomial root finding problems, and then after parameters pairing, 2D DOA estimation is realized, which is directly processed by the root-finding MUSIC method, avoiding the process of virtualization and spatial smoothing as well as the huge computational load brought by spectral peak search while maintaining estimation accuracy. Simulation results verified that the proposed method outperforms many exiting methods in terms of computational complexity and estimation accuracy. In future research, we will consider using this dimensionality reduction and root-finding method on improved arrays with larger array apertures, or consider combining this method with virtual array technology to obtain greater DOFs.