# Computationally Efficient Ambiguity-Free Two-Dimensional DOA Estimation Method for Coprime Planar Array: RD-Root-MUSIC Algorithm

Authors: Luo Chen; Changbo Ye; Baobao Li
Journal: Mathematical Problems in Engineering
 (2020)

---

## Abstract

While the two-dimensional (2D) spectral peak search suffers from expensive computational burden in direction of arrival (DOA) estimation, we propose a reduced-dimensional root-MUSIC (RD-Root-MUSIC) algorithm for 2D DOA estimation with coprime planar array (CPA), which is computationally efficient and ambiguity-free. Different from the conventional 2D DOA estimation algorithms based on subarray decomposition, we exploit the received data of the two subarrays jointly by mapping CPA to the full array of the CPA (FCPA), which contributes to the enhanced degrees of freedom (DOFs) and improved estimation performance. In addition, due to the ambiguity-free characteristic of the FCPA, the extra ambiguity elimination operation can be avoided. Furthermore, we convert the 2D spectral search process into 1D polynomial rooting via reduced-dimension transformation, which substantially reduces the computational complexity while preserving the estimation accuracy. Finally, numerical simulations demonstrate the superiority of the proposed algorithm.


---

## Body

## 1. Introduction

Two-dimensional direction of arrival (2D DOA) estimation has been extensively utilized in radar, sonar, wireless communication, and other fields [1–3]. Numerous DOA estimation algorithms, e.g., multiple signal classification (MUSIC) [4], estimation of signal parameters via rotational invariance techniques (ESPRIT) [5–7], propagator method (PM) [8], and PARAllel FACtor (PARAFAC) technique, have been applied to various planar arrays, such as uniform planar arrays (UPAs) [9–11], L-shaped arrays [12, 13], uniform circular arrays [5], and two parallel linear arrays [14]. However, the distance between adjacent elements in these traditional arrays is limited to no larger than half-wavelength to avoid spatial aliasing, which bring undesired serious mutual coupling effect. Besides, since the DOA estimation accuracy has positive correlation with the array aperture, the limited interelement spacing brings a negative impact on the estimation performance [15, 16].In recent years, sparse arrays such as coprime arrays [15, 17, 18], nested arrays [19], and minimum redundancy arrays [20] have been proposed to tackle this issue. As a typical sparse array, the coprime arrays have inherent superiorities over the conventional compact arrays, including enlarged array aperture, increased DOFs, and reduced mutual coupling [21], employed the traditional 2D-MUSIC algorithm to CPA by exploring the transformation relation between true and ambiguous estimates, and further proposed a 2D partial spectrum search (2D-PSS) method [22], which considerably relieves the computational burden of 2D total spectrum search (TSS). By combining the reduced-dimensional MUSIC (RD-MUSIC) [23] method with the PSS method, the reduced-dimension transformation is performed to further reduce complexity [24]. A generalized CPA array structure was designed in [25] based on the mechanism of ambiguity elimination method, which provides a more flexible array layout and significant increase in DOFs. The aforementioned methods [22, 24, 25] can be categorized as decomposition algorithms, which process the received data of each subarray separately and then the estimates are combined to determine the final DOAs, whereas the mutual information between the two subarrays is unfortunately neglected. The root-MUSIC method to CPA with low complexity is applied [26], while the estimation performance of this cascade approach depends heavily on the initial estimates, especially at low SNRs. An ambiguity-free MUSIC (AF-MUSIC) method was proposed in [27], where the output of two subarrays were stacked and processed jointly to avoid the ambiguous problem. Although the entire information of CPA is fully exploited, the 2D spectrum search leads to heavy computational burden.In the above research studies for 2D DOA estimation methods with CPA, they either treat the two subarrays as individual arrays, which suffers performance degradation due to the loss of mutual information, or 2D spectrum search is required leading to expensive computational cost, or extra ambiguity elimination process is involved. To address these issues, we propose a computationally efficient ambiguity-free algorithm via reduced-dimensional polynomial rooting technique. Specifically, we first map the CPA into full array of the CPA (FCPA) using an extraction matrix, based on the characteristics of array configuration, which enables the sufficient utilization of the entire received data of the CPA. Meanwhile, the number of achievable DOFs is enhanced benefiting from the utilization of full information, as compared to the conventional decomposition algorithms. Furthermore, we transform the 2D spectrum search into 2D polynomial root-finding process and further perform reduced-dimension transformation to convert the 2D root-finding operation into two 1D one, which substantially reduces complementation complexity as well as the computational burden. In addition, extra ambiguity elimination can be avoided owing to the inherent ambiguity-free characteristic of the FCPA.We summarize the major contributions of our work below:(1)
We construct the FCPA corresponding to CPA, which processes ambiguity-free characteristic and thereby extra ambiguity elimination operation can be avoided(2)
We exploit the received data of the two subarrays jointly, where improved DOA estimation performance as well as enhanced achievable DOFs can be achieved(3)
We propose a reduced-dimensional polynomial root-finding algorithm with CPA for 2D DOA estimation, which transforms the 3D spectrum search into 1D polynomial rooting and hence reduces the complexity significantly while preserving the estimation accuracyWe outline this paper as follows. Section2 introduces the data model of CPA and its corresponding FCPA. The proposed algorithm is elaborated in Section 3, and we analyze the complexity and DOFs in Section 4. Section 5 provides simulation results to corroborate the effectiveness of the proposed algorithm, and Section 6 concludes this paper.


### 1.1. Notations

Bold uppercase (lowercase) characters represent matrices (vectors).⋅T, ⋅H, ⋅−1, and ⋅∗ denote the transpose, conjugate transpose, inverse, and conjugate operation, respectively. ⊗ and ⊙ are Kronecker product and Khatri–Rao product, respectively. Rank⋅ means the rank of the matrix. angle⋅ represents the phase operator. det⋅ denotes the determinant of the matrix.


## 1.1. Notations

Bold uppercase (lowercase) characters represent matrices (vectors).⋅T, ⋅H, ⋅−1, and ⋅∗ denote the transpose, conjugate transpose, inverse, and conjugate operation, respectively. ⊗ and ⊙ are Kronecker product and Khatri–Rao product, respectively. Rank⋅ means the rank of the matrix. angle⋅ represents the phase operator. det⋅ denotes the determinant of the matrix.


## 2. Preliminaries

### 2.1. Data Model with CPA

Assume thatK far-field narrowband uncorrelated signals impinge on the CPA with DOAs θk,ϕk,k=1,2,…,K, where θk and ϕk are the elevation and azimuth angles of the kth signal, respectively. The CPA consists of two uniform planar arrays with M1×M1 and M2×M2 sensors. The spacing between adjacent elements of subarray 1 with M1×M1 sensors is d1=M2λ/2, and subarray 2 with M2×M2 sensors has the interelement spacing d2=M1λ/2, where M1 and M2 are coprime integers and λ is the wavelength. The total number of elements is T=M12+M22−1 since the two subarrays share the same element at the origin. Define a transformation as uk=sinθksinϕk and vk=sinθkcosϕk for simplification. A CPA configuration is displayed in Figure 1 as an example, where M1=2,M2=3, and T=12.Figure 1
CPA configuration withM1=2 and M2=3.For theith (i = 1, 2) subarray, the received signal can be expressed by [22](1)Xi=AiS+Ni,where S represents the source matrix and S=s1,s2,…,sKT∈ℂK×L, sk denotes source vector and sk=sk1,…,skLT∈ℂL×1, L is the number of snapshots, ℂ represents a complex set, Ai is the steering matrix of the ith (i = 1, 2) subarray, and Ai=ayiu1⊗axiv1,…,ayiuK⊗axivK∈ℂMi2×K, ayiuk and axivk represent the steering vectors along the y-axis and x-axis, respectively, the specific forms can be expressed as ayiuk=1,ej2πdiuk/λ,…,ej2πMi−1diuk/λT∈ℂMi×1 and axivk=1,ej2πdivk/λ,…,ej2πMi−1divk/λT∈ℂMi×1, Ni is white Gaussian noise with mean value zero and variance σ2 of the ith (i = 1, 2) subarray.The output of the whole CPA can be stacked as [27](2)X=X1X2=A1A2S+N1N2=AS+N,where A represents the direction matrix of the whole CPA, A=A1T,A2TT∈ℂM12+M22×K, and N denotes the white Gaussian noise of the whole arrays and N=N1T,N2TT,N1∈ℂM12×L,N2∈ℂM22×L.In practice, the covariance matrix ofX can be calculated using L snapshots by(3)R^=1LXXH.Perform eigenvalue decomposition (EVD) andR^ can be decomposed by(4)R^=EsΛsEsH+EnΛnEnH,where Λs is a diagonal matrix whose diagonal elements are the K largest eigenvalues, Λn is a diagonal matrix composed of the rest eigenvalues, Es represents the signal subspace spanned of the eigenvectors corresponding to the K largest eigenvalues, and En is the noise subspace composed of the remaining eigenvectors.


### 2.2. Full Array of CPA and Extraction

The CPA can be extracted from a large nonuniform planar array, which can be denoted by the full array of CPA (FCPA). The sensor location of the FCPA can be expressed as(5)ΩFCPA=ΩxFCPA,ΩyFCPA,where ΩxFCPA and ΩyFCPA represent the location sets of x-axis and y-axis, respectively. ΩxFCPA=m1d1∪m2d2, ΩyFCPA=m1d1∪m2d2, 0≤m1≤M1−1, 0≤m2≤M2−1, and m1,m2∈ℤ. Accordingly, the number of elements in FCPA is TF=M1+M2−12.Figure2 illustrates the FCPA corresponding to the CPA shown in Figure 1, where ΩxFCPA=0,2,3,4d and ΩyFCPA=0,2,3,4d, d=λ/2. It can be observed that the FCPA contains all elements of the CPA and has four additional elements with sensor number 7, 10, 12, and 13, respectively, which demonstrates that the CPA can be regarded as an extraction from FCPA.Figure 2
The FCPA configuration corresponding to CPA.According to the correspondence of the CPA and FCPA, we introduce an extraction matrixG∈0,1 to characterize the mapping relation as(6)A=GAF,where G∈ℤM12+M22×M1+M2−12 represents the extraction matrix, ℤ denotes a set of integers, AF is the steering matrix of the FCPA, AF=aFyu1⊗aFxv1,…,aFyuK⊗aFxvK∈ℂM1+M2−12×K, aFyuk and aFxvk represent the steering vector of the FCPA along the y-axis and x-axis, the specific forms are aFyuk=1,ej2πdFy2uk/λ,…,ej2πdFyM1+M2−1uk/λTaFxvk=1,ej2πdFx2vk/λ,…,ej2πdFxM1+M2−1vk/λT respectively, and dFyi∈ΩyFCPA,dFxi∈ΩxFCPA1≤i≤M1+M2−1 denote the location sets of elements on the y-axis and x-axis, respectively.To demonstrate the extraction more specifically, each element in the CPA and FCPA is labeled according to their order in the steering vectors, i.e.,1∼M12 for the subarray 1 and M12+1∼M12+M22 for the subarray 2 in CPA, 1∼M1+M2−12 for the FCPA. If the ith sensor in the CPA and the jth sensor in the FCPA overlap; then, gij=1, otherwise gij=0, where gij denotes the (i, j)th element of G. For the FCPA given in Figure 2, G is a 13×16 matrix with 3 columns of all zeros.Definition 1.
(extraction efficiency). The extraction efficiency is the proportion of nonzero elements in the extraction matrix.
The sensor location of CPA can be given by(7)ΩCPA=ΩxCPA,0∪0,ΩyCPA,where ΩxCPA=m1d∪m2d2, ΩyCPA=m1d1∪m2d2, 0≤m1≤M1−1, 0≤m2≤M2−1, and m1,m2∈ℤ. Then, we construct a uniform planar array that has the same array aperture as the CPA with the location set:(8)ΩUPA=ΩxUPA,0∪0,ΩyUPA,where ΩxUPA=md, ΩyUPA=md, d=λ/2, 0≤m≤M1M2−1, and m∈ℤ.
The locations of the sensors in the FCPA corresponding to the CPA can be expressed as(9)ΩFCPA=ΩxFCPA,ΩyFCPA,where ΩxCPA⊆ΩxFCPA⊆ΩxUPA and ΩyCPA⊆ΩyFCPA⊆ΩyUPA represent the location set of FCPA on the x- and y-axes, respectively.
For the CPA study in this paper, the corresponding FCPA can be constructed by the following four forms:(10)ΩxFCPA=0,2,3,4d,ΩyFCPA=0,2,3,4d,ΩxFCPA_1=0,2,3,4d,ΩyFCPA_1=0,1,2,3,4d,ΩxFCPA_2=0,1,2,3,4d,ΩyFCPA_2=0,2,3,4d,ΩxFCPA_3=0,1,2,3,4d,ΩyFCPA_3=0,1,2,3,4d.
According to Definition1, the extraction efficiency of the above schemes is 0.0625, 0.05, 0.05, and 0.04, respectively. It is clearly seen that the FCPA we designed in (5) has the highest extraction efficiency.


## 2.1. Data Model with CPA

Assume thatK far-field narrowband uncorrelated signals impinge on the CPA with DOAs θk,ϕk,k=1,2,…,K, where θk and ϕk are the elevation and azimuth angles of the kth signal, respectively. The CPA consists of two uniform planar arrays with M1×M1 and M2×M2 sensors. The spacing between adjacent elements of subarray 1 with M1×M1 sensors is d1=M2λ/2, and subarray 2 with M2×M2 sensors has the interelement spacing d2=M1λ/2, where M1 and M2 are coprime integers and λ is the wavelength. The total number of elements is T=M12+M22−1 since the two subarrays share the same element at the origin. Define a transformation as uk=sinθksinϕk and vk=sinθkcosϕk for simplification. A CPA configuration is displayed in Figure 1 as an example, where M1=2,M2=3, and T=12.Figure 1
CPA configuration withM1=2 and M2=3.For theith (i = 1, 2) subarray, the received signal can be expressed by [22](1)Xi=AiS+Ni,where S represents the source matrix and S=s1,s2,…,sKT∈ℂK×L, sk denotes source vector and sk=sk1,…,skLT∈ℂL×1, L is the number of snapshots, ℂ represents a complex set, Ai is the steering matrix of the ith (i = 1, 2) subarray, and Ai=ayiu1⊗axiv1,…,ayiuK⊗axivK∈ℂMi2×K, ayiuk and axivk represent the steering vectors along the y-axis and x-axis, respectively, the specific forms can be expressed as ayiuk=1,ej2πdiuk/λ,…,ej2πMi−1diuk/λT∈ℂMi×1 and axivk=1,ej2πdivk/λ,…,ej2πMi−1divk/λT∈ℂMi×1, Ni is white Gaussian noise with mean value zero and variance σ2 of the ith (i = 1, 2) subarray.The output of the whole CPA can be stacked as [27](2)X=X1X2=A1A2S+N1N2=AS+N,where A represents the direction matrix of the whole CPA, A=A1T,A2TT∈ℂM12+M22×K, and N denotes the white Gaussian noise of the whole arrays and N=N1T,N2TT,N1∈ℂM12×L,N2∈ℂM22×L.In practice, the covariance matrix ofX can be calculated using L snapshots by(3)R^=1LXXH.Perform eigenvalue decomposition (EVD) andR^ can be decomposed by(4)R^=EsΛsEsH+EnΛnEnH,where Λs is a diagonal matrix whose diagonal elements are the K largest eigenvalues, Λn is a diagonal matrix composed of the rest eigenvalues, Es represents the signal subspace spanned of the eigenvectors corresponding to the K largest eigenvalues, and En is the noise subspace composed of the remaining eigenvectors.


## 2.2. Full Array of CPA and Extraction

The CPA can be extracted from a large nonuniform planar array, which can be denoted by the full array of CPA (FCPA). The sensor location of the FCPA can be expressed as(5)ΩFCPA=ΩxFCPA,ΩyFCPA,where ΩxFCPA and ΩyFCPA represent the location sets of x-axis and y-axis, respectively. ΩxFCPA=m1d1∪m2d2, ΩyFCPA=m1d1∪m2d2, 0≤m1≤M1−1, 0≤m2≤M2−1, and m1,m2∈ℤ. Accordingly, the number of elements in FCPA is TF=M1+M2−12.Figure2 illustrates the FCPA corresponding to the CPA shown in Figure 1, where ΩxFCPA=0,2,3,4d and ΩyFCPA=0,2,3,4d, d=λ/2. It can be observed that the FCPA contains all elements of the CPA and has four additional elements with sensor number 7, 10, 12, and 13, respectively, which demonstrates that the CPA can be regarded as an extraction from FCPA.Figure 2
The FCPA configuration corresponding to CPA.According to the correspondence of the CPA and FCPA, we introduce an extraction matrixG∈0,1 to characterize the mapping relation as(6)A=GAF,where G∈ℤM12+M22×M1+M2−12 represents the extraction matrix, ℤ denotes a set of integers, AF is the steering matrix of the FCPA, AF=aFyu1⊗aFxv1,…,aFyuK⊗aFxvK∈ℂM1+M2−12×K, aFyuk and aFxvk represent the steering vector of the FCPA along the y-axis and x-axis, the specific forms are aFyuk=1,ej2πdFy2uk/λ,…,ej2πdFyM1+M2−1uk/λTaFxvk=1,ej2πdFx2vk/λ,…,ej2πdFxM1+M2−1vk/λT respectively, and dFyi∈ΩyFCPA,dFxi∈ΩxFCPA1≤i≤M1+M2−1 denote the location sets of elements on the y-axis and x-axis, respectively.To demonstrate the extraction more specifically, each element in the CPA and FCPA is labeled according to their order in the steering vectors, i.e.,1∼M12 for the subarray 1 and M12+1∼M12+M22 for the subarray 2 in CPA, 1∼M1+M2−12 for the FCPA. If the ith sensor in the CPA and the jth sensor in the FCPA overlap; then, gij=1, otherwise gij=0, where gij denotes the (i, j)th element of G. For the FCPA given in Figure 2, G is a 13×16 matrix with 3 columns of all zeros.Definition 1.
(extraction efficiency). The extraction efficiency is the proportion of nonzero elements in the extraction matrix.
The sensor location of CPA can be given by(7)ΩCPA=ΩxCPA,0∪0,ΩyCPA,where ΩxCPA=m1d∪m2d2, ΩyCPA=m1d1∪m2d2, 0≤m1≤M1−1, 0≤m2≤M2−1, and m1,m2∈ℤ. Then, we construct a uniform planar array that has the same array aperture as the CPA with the location set:(8)ΩUPA=ΩxUPA,0∪0,ΩyUPA,where ΩxUPA=md, ΩyUPA=md, d=λ/2, 0≤m≤M1M2−1, and m∈ℤ.
The locations of the sensors in the FCPA corresponding to the CPA can be expressed as(9)ΩFCPA=ΩxFCPA,ΩyFCPA,where ΩxCPA⊆ΩxFCPA⊆ΩxUPA and ΩyCPA⊆ΩyFCPA⊆ΩyUPA represent the location set of FCPA on the x- and y-axes, respectively.
For the CPA study in this paper, the corresponding FCPA can be constructed by the following four forms:(10)ΩxFCPA=0,2,3,4d,ΩyFCPA=0,2,3,4d,ΩxFCPA_1=0,2,3,4d,ΩyFCPA_1=0,1,2,3,4d,ΩxFCPA_2=0,1,2,3,4d,ΩyFCPA_2=0,2,3,4d,ΩxFCPA_3=0,1,2,3,4d,ΩyFCPA_3=0,1,2,3,4d.
According to Definition1, the extraction efficiency of the above schemes is 0.0625, 0.05, 0.05, and 0.04, respectively. It is clearly seen that the FCPA we designed in (5) has the highest extraction efficiency.


## 3. The Proposed Algorithm

### 3.1. 2D-MUSIC Algorithm

Derive from the orthogonal relationship between the noise subspace and the steering vector, and the spectral function of CPA can be represented by [26, 27](11)Pu,v=1aHu,vEnEnHau,v,where au,v denotes the steering vector of CPA, au,v=a1Tu,v,a2Tu,vT∈ℂM12+M22×1, aiu,v=ayiu⊗axiv, ayiu=1,ej2πdiu/λ,…,ej2πMi−1diu/λT, axiv=1,ej2πdiv/λ,…,ej2πMi−1div/λT (i = 1, 2), and En is the noise subspace of CPA.According to (6), we have au,v=GaFu,v. Then, (11) can be rewritten as [26, 27](12)Pu,v=1aFHu,vGHEnEnHGaFu,v=1aFHu,vEFnEFnHaFu,v,where aFu,v is the steering vector of the FCPA, aFu,v=aFyu⊗aFxv∈ℂM1+M2−12×1, aFyu=1,ej2πdFy2u/λ,…,ej2πdFyM1+M2−1u/λT∈ℂM1+M2−1×1, aFxv=1,ej2πdFx2v/λ,…,ej2πdFxM1+M2−1v/λT∈ℂM1+M2−1×1. EFn denotes the noise space of the FCPA, and EFn=GHEn.Although the autopaired 2D DOA estimates can be obtained via performing spectral search on (12), it suffers from tremendously expensive computational cost. To tackle this issue, we first performed reduced-dimension transformation and then exploited 1D polynomial root-finding technique to estimate u and v.


### 3.2. Reduced-Dimensional Polynomial Root-Finding Process

Construct the polynomial based on (12) as(13)Vu,v=aFyu⊗aFxvHEFnEFnHEFyu⊗aFxv=aFxHvaFyu⊗IM1+M2−1HEFnEFnHaFyu⊗IM1+M2−1aFxv=aFxHvQuaFxv,(14)Vu,v=aFyu⊗aFxvHEFnEFnHaFyu⊗aFxv=aFyHuIM1+M2−1⊗aFxvHEFnEFnHIM1+M2−1⊗aFxvaFyu=aFyHuQvaFyu,where Qu=aFyu⊗IM1+M2−1HEFnEFnHaFyu⊗IM1+M2−1 and Qv=IM1+M2−1⊗aFxvHEFnEFnHIM1+M2−1⊗aFxv.According to the relation of rank of the matrix product, the following constraint(15)0<RankEFnEFnH≤RankEFnhas to be satisfied, and then we have(16)detEFnEFnH≠0.We can conclude thatdetQu is nonzero polynomial from (16); thus, Qu is a factor of Vu,v. As Qu depends only on the variable u, the roots of detQu=0 can make the following equation hold:(17)Vu,v=aFxHvQuaFxv=0.It is noteworthy that only 1D polynomial is involved to achieve the estimates ofu. Similarly, we can get the estimates of v. Consequently, the problem of obtaining paired estimates of u and v from the 2D polynomial is transformed into two 1D root-finding process. Then, we reconstruct (13) and (14) as(18)detQu=detaFyu⊗IM1+M2−1HEFnEFnHaFyu⊗IM1+M2−1=0,(19)detQv=detIM1+M2−1⊗aFxvHEFnEFnHIM1+M2−1⊗aFxv=0.Define(20)z1=ej2πdu/λ,z2=ej2πdv/λ,where d=λ/2.Define the steering vector of UPA that has the same array aperture as FCPA along thex-axis as aExv=1,ej2πdv/λ,…,ej2πminM1,M2maxM1,M2−1dv/λT∈ℂminM1,M2maxM1,M2−1+1×1. We assume that M1<M2 for simplification and then aExv can be expressed as aExv=1,ej2πdv/λ,…,ej2πM1M2−1dv/λT∈ℂM1M2−1+1×1. Based on the correspondence between FCPA and UPA with same array aperture, we have(21)aFxv=G1aExv,where G1∈ℤM1+M2−1×M1M2−1+1. To be specific, g1ij=1 holds when the ith sensor in the aFxv and the jth sensor in the aExv overlap, otherwise gij=0, where gij is the (i, j)th element of G1. For the FCPA displayed in Figure 2, G1 is a 4×5 matrix with one columns of all zeros. Similarly, we can obtain the relation aFyu=G1aEyu.Correspondingly, the steering vectors can be rewritten as(22)aFyu=G1aEyu=G11,ej2πdu/λ,…,ej2πM1M2−1du/λT=G11,z1,…,z1M1M2−1T=∞Fyz1,(23)aFxv=G1aExv=G11,ej2πdv/λ,…,ej2πM1M2−1dv/λT=∞11,z2,…,z2M1M2−1T=aFxz2.Without loss of generality, substitutingz1M1M2−1aFyTz1−1⊗IM1+M2−1H for aFyu⊗IM1+M2−1H and substituting z2M1M2−1IM1+M2−1⊗aFxTz2−1H for IM1+M2−1⊗aFxvH, i.e.,(24)detQz1=detz1M1M2−1aFyTz1−1⊗IM1+M2−1HEnEnHaFyz1−1⊗IM1+M2−1=0,(25)detQz2=detz2M1M2−1IM1+M2−1⊗aFxTz2−1HEnEnHIM1+M2−1⊗aFxz2−1=0.Considering thatdetQz1 and detQz2 are polynomials of even degree, u^k and v^i can be obtained from the K roots distributed closest to the unit circle corresponding to (24) and (25), and the roots are denoted by z^11,…,z^1k,…,z^1K and z^21,…,z^2i,…,z^2K, respectively, i.e.,(26)u^k=anglez^1kλ2πd,k=1,…,K,(27)v^i=anglez^2iλ2πd,i=1,…,K.For the conventional DOA estimation methods with CPA, the ambiguity elimination operation is required since the interelement spacing in the two subarrays is larger than half-wavelength. The FCPA is an unambiguous array which has at least one sensor pair with separation no larger than half-wavelength according to (5), and we can obtain the true DOA estimates directly after the parameter pairing without extra ambiguity elimination process.


### 3.3. Parameter Pairing and DOA Estimation

In this part, we determine the pairing ofu^k and v^i since the two root-finding procedures are conducted separately. Construct the cost function for pairing as(28)Vk,i=argi=1,…,KminaHu^k,v^iEnEnHau^k,v^ik=1,…,K,where au^k,v^i represents the steering vector reconstructed u^k and v^i, which can be obtained according to (1).For eachu^k, we can obtain the value of i and k that minimize Vk,i1≤i≤K, and we define the paired index as i′ and k. Finally, the 2D DOAs can be calculated by(29)θ^k=arcsinu^k2+v^i′2,1≤k≤K,(30)ϕ^k=arctanu^kv^i′,1≤k≤K,where v^i′ is reconstructed by i′1≤i′≤K.


### 3.4. The Procedure of the Proposed Algorithm

We summarize the major steps of the proposed algorithm as follows:Step1: calculateR^ of the received data X and perform EVD to obtain the noise space EnStep 2: reconstruct the spectral functionPu,v according to (6)Step 3: construct the polynomialVu,v and conduct reduced-dimension transformation according to (12)Step 4: calculateu^k and v^i according to (26) and (27)Step 5: perform parameter matching to obtainθ^k and ϕ^k according to (28)–(30)


## 3.1. 2D-MUSIC Algorithm

Derive from the orthogonal relationship between the noise subspace and the steering vector, and the spectral function of CPA can be represented by [26, 27](11)Pu,v=1aHu,vEnEnHau,v,where au,v denotes the steering vector of CPA, au,v=a1Tu,v,a2Tu,vT∈ℂM12+M22×1, aiu,v=ayiu⊗axiv, ayiu=1,ej2πdiu/λ,…,ej2πMi−1diu/λT, axiv=1,ej2πdiv/λ,…,ej2πMi−1div/λT (i = 1, 2), and En is the noise subspace of CPA.According to (6), we have au,v=GaFu,v. Then, (11) can be rewritten as [26, 27](12)Pu,v=1aFHu,vGHEnEnHGaFu,v=1aFHu,vEFnEFnHaFu,v,where aFu,v is the steering vector of the FCPA, aFu,v=aFyu⊗aFxv∈ℂM1+M2−12×1, aFyu=1,ej2πdFy2u/λ,…,ej2πdFyM1+M2−1u/λT∈ℂM1+M2−1×1, aFxv=1,ej2πdFx2v/λ,…,ej2πdFxM1+M2−1v/λT∈ℂM1+M2−1×1. EFn denotes the noise space of the FCPA, and EFn=GHEn.Although the autopaired 2D DOA estimates can be obtained via performing spectral search on (12), it suffers from tremendously expensive computational cost. To tackle this issue, we first performed reduced-dimension transformation and then exploited 1D polynomial root-finding technique to estimate u and v.


## 3.2. Reduced-Dimensional Polynomial Root-Finding Process

Construct the polynomial based on (12) as(13)Vu,v=aFyu⊗aFxvHEFnEFnHEFyu⊗aFxv=aFxHvaFyu⊗IM1+M2−1HEFnEFnHaFyu⊗IM1+M2−1aFxv=aFxHvQuaFxv,(14)Vu,v=aFyu⊗aFxvHEFnEFnHaFyu⊗aFxv=aFyHuIM1+M2−1⊗aFxvHEFnEFnHIM1+M2−1⊗aFxvaFyu=aFyHuQvaFyu,where Qu=aFyu⊗IM1+M2−1HEFnEFnHaFyu⊗IM1+M2−1 and Qv=IM1+M2−1⊗aFxvHEFnEFnHIM1+M2−1⊗aFxv.According to the relation of rank of the matrix product, the following constraint(15)0<RankEFnEFnH≤RankEFnhas to be satisfied, and then we have(16)detEFnEFnH≠0.We can conclude thatdetQu is nonzero polynomial from (16); thus, Qu is a factor of Vu,v. As Qu depends only on the variable u, the roots of detQu=0 can make the following equation hold:(17)Vu,v=aFxHvQuaFxv=0.It is noteworthy that only 1D polynomial is involved to achieve the estimates ofu. Similarly, we can get the estimates of v. Consequently, the problem of obtaining paired estimates of u and v from the 2D polynomial is transformed into two 1D root-finding process. Then, we reconstruct (13) and (14) as(18)detQu=detaFyu⊗IM1+M2−1HEFnEFnHaFyu⊗IM1+M2−1=0,(19)detQv=detIM1+M2−1⊗aFxvHEFnEFnHIM1+M2−1⊗aFxv=0.Define(20)z1=ej2πdu/λ,z2=ej2πdv/λ,where d=λ/2.Define the steering vector of UPA that has the same array aperture as FCPA along thex-axis as aExv=1,ej2πdv/λ,…,ej2πminM1,M2maxM1,M2−1dv/λT∈ℂminM1,M2maxM1,M2−1+1×1. We assume that M1<M2 for simplification and then aExv can be expressed as aExv=1,ej2πdv/λ,…,ej2πM1M2−1dv/λT∈ℂM1M2−1+1×1. Based on the correspondence between FCPA and UPA with same array aperture, we have(21)aFxv=G1aExv,where G1∈ℤM1+M2−1×M1M2−1+1. To be specific, g1ij=1 holds when the ith sensor in the aFxv and the jth sensor in the aExv overlap, otherwise gij=0, where gij is the (i, j)th element of G1. For the FCPA displayed in Figure 2, G1 is a 4×5 matrix with one columns of all zeros. Similarly, we can obtain the relation aFyu=G1aEyu.Correspondingly, the steering vectors can be rewritten as(22)aFyu=G1aEyu=G11,ej2πdu/λ,…,ej2πM1M2−1du/λT=G11,z1,…,z1M1M2−1T=∞Fyz1,(23)aFxv=G1aExv=G11,ej2πdv/λ,…,ej2πM1M2−1dv/λT=∞11,z2,…,z2M1M2−1T=aFxz2.Without loss of generality, substitutingz1M1M2−1aFyTz1−1⊗IM1+M2−1H for aFyu⊗IM1+M2−1H and substituting z2M1M2−1IM1+M2−1⊗aFxTz2−1H for IM1+M2−1⊗aFxvH, i.e.,(24)detQz1=detz1M1M2−1aFyTz1−1⊗IM1+M2−1HEnEnHaFyz1−1⊗IM1+M2−1=0,(25)detQz2=detz2M1M2−1IM1+M2−1⊗aFxTz2−1HEnEnHIM1+M2−1⊗aFxz2−1=0.Considering thatdetQz1 and detQz2 are polynomials of even degree, u^k and v^i can be obtained from the K roots distributed closest to the unit circle corresponding to (24) and (25), and the roots are denoted by z^11,…,z^1k,…,z^1K and z^21,…,z^2i,…,z^2K, respectively, i.e.,(26)u^k=anglez^1kλ2πd,k=1,…,K,(27)v^i=anglez^2iλ2πd,i=1,…,K.For the conventional DOA estimation methods with CPA, the ambiguity elimination operation is required since the interelement spacing in the two subarrays is larger than half-wavelength. The FCPA is an unambiguous array which has at least one sensor pair with separation no larger than half-wavelength according to (5), and we can obtain the true DOA estimates directly after the parameter pairing without extra ambiguity elimination process.


## 3.3. Parameter Pairing and DOA Estimation

In this part, we determine the pairing ofu^k and v^i since the two root-finding procedures are conducted separately. Construct the cost function for pairing as(28)Vk,i=argi=1,…,KminaHu^k,v^iEnEnHau^k,v^ik=1,…,K,where au^k,v^i represents the steering vector reconstructed u^k and v^i, which can be obtained according to (1).For eachu^k, we can obtain the value of i and k that minimize Vk,i1≤i≤K, and we define the paired index as i′ and k. Finally, the 2D DOAs can be calculated by(29)θ^k=arcsinu^k2+v^i′2,1≤k≤K,(30)ϕ^k=arctanu^kv^i′,1≤k≤K,where v^i′ is reconstructed by i′1≤i′≤K.


## 3.4. The Procedure of the Proposed Algorithm

We summarize the major steps of the proposed algorithm as follows:Step1: calculateR^ of the received data X and perform EVD to obtain the noise space EnStep 2: reconstruct the spectral functionPu,v according to (6)Step 3: construct the polynomialVu,v and conduct reduced-dimension transformation according to (12)Step 4: calculateu^k and v^i according to (26) and (27)Step 5: perform parameter matching to obtainθ^k and ϕ^k according to (28)–(30)


## 4. Performance Analysis

### 4.1. Complexity Analysis

Herein, we compare the computational complexity of the proposed algorithm, 2D-PSS [22], RD-MUSIC [23], AF-MUSIC [27], and 2D-ROOT [26] methods in this section. For the proposed algorithm, calculating the covariance matrix requires OT2L and the complexity of eigenvalue decomposition is OT3. The root-finding operation costs O22M1M2−13 and parameter matching process requires OK2T−KT+1. Consequently, the total complexity is OT2L+T3+22M1M2−13+K2T−KT+1.Table1 lists the total complexity of above algorithms, where Δ=0.0001 is the spectral search interval. In addition, Figure 3 displays the complexity comparison versus number of sensors, where K=2andL=500, while the complexity comparison with different number of snapshots is shown in Figure 4, where K=2,M1=2,andM2=3. It is clearly seen that the proposed algorithm owns the approximate low complexity to the 2D-ROOT method, which is significantly lower than that of the 2D-PSS, RD-MUSIC, and AF-MUSIC methods, as the spectral search process with heavy computational burden is transformed into computationally efficient polynomial root-finding.Table 1
Complexity of different algorithms.
AlgorithmComplexityProposedOT2L+T3+22M1M2−123+K2T−KT+12D-PSSOM14+M24L+M16+M26+4M12M12−K/M22Δ2+4M22M22−K/M12Δ2RD-MUSICOM14+M24L+M16+M26+22M13M13−K+M14+M13/Δ+22M23M23−K+M24+M23/Δ2D-ROOTOT2L+T3+8M16+KM13+M26+KM23AF-MUSICOT2L+T3+8TT−K/Δ2Figure 3
Complexity comparison versus number of elements.Figure 4
Complexity comparison versus snapshots.


### 4.2. Cramer–Rao Bound

In this part, the derivation of Cramer–Rao Bound (CRB) of the 2D DOA estimation with CPA is given as the performance comparison metric.DefineA=A1A2, where A1=ay1θ1,ϕ1⊗ax1θ1,ϕ1,…,ay1θK,ϕK⊗ax1θK,ϕK and A2=ay2θ1,ϕ1⊗ax2θ1,ϕ1,…,ay2θK,ϕK⊗ax2θK,ϕK. According to [28], the CRB of CPA can be given by(31)CRB=σ22LReDHΠA⊥D⊕P^T−1,where D=∂a1/∂θ1,…,∂aK/∂θK,∂a1/∂ϕ1,…,∂aK/∂ϕK, ΠA⊥=IM12+M22−1−AAHA−1AH, P^=P^sP^sP^sP^s, P^s=SSH/L, ak denotes the kth column of A, and σ2 is the variance of the received noise.


### 4.3. Achievable DOFs

As for the conventional decomposition-based 2D DOA estimation methods, such as RD-MUSIC, 2D-PSS, and other algorithms, the number of achievable DOFs isminM12−1,M22−1. According [26], M12+M22−2 signals at most can be resolved by utilizing AF-MUSIC algorithm which processes the received data of the two subarrays jointly, whereas extremely high computational complexity is involved due to the 2D spectral peak search. RD root-finding technique is employed in the proposed algorithm to deal with the complexity and the achievable DOFs can be obtained from (15). If u or v does not match any of the incident signals, we have the constraint(32)maxM1,M2≤RankEFnEFnH≤RankEFn,where RankEFn=minM1+M2−12,M12+M22−1−K. Consequently, the maximum number of signals which can be identified by the proposed algorithm is(33)K≤M12+M22−1−maxM1,M2.It is clear that the proposed algorithm can greatly improve the DOFs compared with the traditional 2D DOA estimation methods.


### 4.4. Advantages

Based on the above discussion, the advantages of the proposed algorithm can be listed as follows:(1)
It can achieve superior estimation performance and higher achievable DOFs than the decomposition methods, owing to the utilization of entire received data. Moreover, additional ambiguity elimination is no longer required due to the ambiguity-free characteristic of the FCPA.(2)
It outperforms the AF-MUSIC algorithm, RD-MUSIC algorithm, 2D-PSS algorithm, and 2D-ROOT algorithm in estimation performance.(3)
The proposed algorithm owns much lower computational complexity than the RD-MUSIC, 2D-PSS, and AF-MUSIC methods, which is approximately as low as 2D-ROOT algorithm.


## 4.1. Complexity Analysis

Herein, we compare the computational complexity of the proposed algorithm, 2D-PSS [22], RD-MUSIC [23], AF-MUSIC [27], and 2D-ROOT [26] methods in this section. For the proposed algorithm, calculating the covariance matrix requires OT2L and the complexity of eigenvalue decomposition is OT3. The root-finding operation costs O22M1M2−13 and parameter matching process requires OK2T−KT+1. Consequently, the total complexity is OT2L+T3+22M1M2−13+K2T−KT+1.Table1 lists the total complexity of above algorithms, where Δ=0.0001 is the spectral search interval. In addition, Figure 3 displays the complexity comparison versus number of sensors, where K=2andL=500, while the complexity comparison with different number of snapshots is shown in Figure 4, where K=2,M1=2,andM2=3. It is clearly seen that the proposed algorithm owns the approximate low complexity to the 2D-ROOT method, which is significantly lower than that of the 2D-PSS, RD-MUSIC, and AF-MUSIC methods, as the spectral search process with heavy computational burden is transformed into computationally efficient polynomial root-finding.Table 1
Complexity of different algorithms.
AlgorithmComplexityProposedOT2L+T3+22M1M2−123+K2T−KT+12D-PSSOM14+M24L+M16+M26+4M12M12−K/M22Δ2+4M22M22−K/M12Δ2RD-MUSICOM14+M24L+M16+M26+22M13M13−K+M14+M13/Δ+22M23M23−K+M24+M23/Δ2D-ROOTOT2L+T3+8M16+KM13+M26+KM23AF-MUSICOT2L+T3+8TT−K/Δ2Figure 3
Complexity comparison versus number of elements.Figure 4
Complexity comparison versus snapshots.


## 4.2. Cramer–Rao Bound

In this part, the derivation of Cramer–Rao Bound (CRB) of the 2D DOA estimation with CPA is given as the performance comparison metric.DefineA=A1A2, where A1=ay1θ1,ϕ1⊗ax1θ1,ϕ1,…,ay1θK,ϕK⊗ax1θK,ϕK and A2=ay2θ1,ϕ1⊗ax2θ1,ϕ1,…,ay2θK,ϕK⊗ax2θK,ϕK. According to [28], the CRB of CPA can be given by(31)CRB=σ22LReDHΠA⊥D⊕P^T−1,where D=∂a1/∂θ1,…,∂aK/∂θK,∂a1/∂ϕ1,…,∂aK/∂ϕK, ΠA⊥=IM12+M22−1−AAHA−1AH, P^=P^sP^sP^sP^s, P^s=SSH/L, ak denotes the kth column of A, and σ2 is the variance of the received noise.


## 4.3. Achievable DOFs

As for the conventional decomposition-based 2D DOA estimation methods, such as RD-MUSIC, 2D-PSS, and other algorithms, the number of achievable DOFs isminM12−1,M22−1. According [26], M12+M22−2 signals at most can be resolved by utilizing AF-MUSIC algorithm which processes the received data of the two subarrays jointly, whereas extremely high computational complexity is involved due to the 2D spectral peak search. RD root-finding technique is employed in the proposed algorithm to deal with the complexity and the achievable DOFs can be obtained from (15). If u or v does not match any of the incident signals, we have the constraint(32)maxM1,M2≤RankEFnEFnH≤RankEFn,where RankEFn=minM1+M2−12,M12+M22−1−K. Consequently, the maximum number of signals which can be identified by the proposed algorithm is(33)K≤M12+M22−1−maxM1,M2.It is clear that the proposed algorithm can greatly improve the DOFs compared with the traditional 2D DOA estimation methods.


## 4.4. Advantages

Based on the above discussion, the advantages of the proposed algorithm can be listed as follows:(1)
It can achieve superior estimation performance and higher achievable DOFs than the decomposition methods, owing to the utilization of entire received data. Moreover, additional ambiguity elimination is no longer required due to the ambiguity-free characteristic of the FCPA.(2)
It outperforms the AF-MUSIC algorithm, RD-MUSIC algorithm, 2D-PSS algorithm, and 2D-ROOT algorithm in estimation performance.(3)
The proposed algorithm owns much lower computational complexity than the RD-MUSIC, 2D-PSS, and AF-MUSIC methods, which is approximately as low as 2D-ROOT algorithm.


## 5. Simulations

In this section, we perform 500 Monte Carlo simulations to validate the performance of the proposed algorithm. Assume thatK uncorrelated far-field narrowband signals impinge on the CPA. Define root mean square error (RMSE) by(34)RMSE=1K∑k=1K1500∑i=1500ϕ^k,i−ϕk2+θ^k,i−θk2,where ϕ^k,i and θ^k,i are the estimates of the kth signal in the ith trial corresponding to the true azimuth ϕk and elevation θk, respectively.


### 5.1. Scatter Figures

The scatter figures of the proposed algorithm underK sources are as follows, where K=6,SNR=30dB,L=1000, θ1,θ2,θ3,θ4,θ5,θ6=5∘,15∘,25∘,35∘,45∘,55∘, and ϕ1,ϕ2,ϕ3,ϕ4,ϕ5,ϕ6=5∘,15∘,25∘,35∘,45∘,55∘. The CPA is composed of two UPAs with M1×M1=2×2 and M2×M2=3×3. As illustrated in Figure 5, the proposed algorithm can effectively distinguish all sources incident on the CPA.Figure 5
(a) Elevation estimates of the proposed algorithm. (b) Azimuth estimates of the proposed algorithm.
(a)(b)


### 5.2. RMSE Results Versus Snapshots

In this part, we present the DOA estimation performance of the proposed algorithm with different number of snapshots in Figure6, where θ1,ϕ1=20∘,30∘, θ2,ϕ2=40∘,50∘, M1×M1=2×2, and M2×M2=3×3. The result demonstrates that the estimation accuracy improves as the number of snapshots increases, owing to the more accurate covariance.Figure 6
RMSE performance of the proposed algorithm versus snapshots.


### 5.3. RMSE Results Versus Number of Sensors

Herein, we provide the RMSE results of the proposed algorithm versus number of sensors in Figure7, where θ1,ϕ1=20∘,30∘, θ2,ϕ2=40∘,50∘, and L=200. As the number of array elements increases, the diversity gain of the receiving antenna increases. It is illustrated clearly that the increased number of sensors leads to improved DOA estimation performance.Figure 7
RMSE performance of the proposed algorithm versus number of sensors.


### 5.4. RMSE Comparison of Different Algorithms

Figure8 exhibits the RMSE comparison of the proposed algorithm, 2D-ROOT algorithm, 2D-PSS algorithm, RD-MUSIC, and AF-MUSIC algorithm, where M1×M1=2×2 and M2×M2=3×3. It is indicated explicitly in Figure 8 that the proposed algorithm and AF-MUSIC algorithm, benefiting from the utilization of the received data of the entire CPA, outperform the decomposition-based RD-MUSIC, 2D-PSS, and 2D-ROOT algorithms. Furthermore, the proposed algorithm yields superior estimation performance to the AF-MUSIC algorithm, as the proposed algorithm directly performs two root-finding operations based on the spectral function (11), while the cascading process in the AF-MUSIC algorithm may result in performance degradation.Figure 8
(a) RMSE comparison of different algorithms versus snapshots. (b) RMSE comparison of different algorithms versus SNR.
(a)(b)


## 5.1. Scatter Figures

The scatter figures of the proposed algorithm underK sources are as follows, where K=6,SNR=30dB,L=1000, θ1,θ2,θ3,θ4,θ5,θ6=5∘,15∘,25∘,35∘,45∘,55∘, and ϕ1,ϕ2,ϕ3,ϕ4,ϕ5,ϕ6=5∘,15∘,25∘,35∘,45∘,55∘. The CPA is composed of two UPAs with M1×M1=2×2 and M2×M2=3×3. As illustrated in Figure 5, the proposed algorithm can effectively distinguish all sources incident on the CPA.Figure 5
(a) Elevation estimates of the proposed algorithm. (b) Azimuth estimates of the proposed algorithm.
(a)(b)


## 5.2. RMSE Results Versus Snapshots

In this part, we present the DOA estimation performance of the proposed algorithm with different number of snapshots in Figure6, where θ1,ϕ1=20∘,30∘, θ2,ϕ2=40∘,50∘, M1×M1=2×2, and M2×M2=3×3. The result demonstrates that the estimation accuracy improves as the number of snapshots increases, owing to the more accurate covariance.Figure 6
RMSE performance of the proposed algorithm versus snapshots.


## 5.3. RMSE Results Versus Number of Sensors

Herein, we provide the RMSE results of the proposed algorithm versus number of sensors in Figure7, where θ1,ϕ1=20∘,30∘, θ2,ϕ2=40∘,50∘, and L=200. As the number of array elements increases, the diversity gain of the receiving antenna increases. It is illustrated clearly that the increased number of sensors leads to improved DOA estimation performance.Figure 7
RMSE performance of the proposed algorithm versus number of sensors.


## 5.4. RMSE Comparison of Different Algorithms

Figure8 exhibits the RMSE comparison of the proposed algorithm, 2D-ROOT algorithm, 2D-PSS algorithm, RD-MUSIC, and AF-MUSIC algorithm, where M1×M1=2×2 and M2×M2=3×3. It is indicated explicitly in Figure 8 that the proposed algorithm and AF-MUSIC algorithm, benefiting from the utilization of the received data of the entire CPA, outperform the decomposition-based RD-MUSIC, 2D-PSS, and 2D-ROOT algorithms. Furthermore, the proposed algorithm yields superior estimation performance to the AF-MUSIC algorithm, as the proposed algorithm directly performs two root-finding operations based on the spectral function (11), while the cascading process in the AF-MUSIC algorithm may result in performance degradation.Figure 8
(a) RMSE comparison of different algorithms versus snapshots. (b) RMSE comparison of different algorithms versus SNR.
(a)(b)


## 6. Conclusion

In this paper, we propose a computationally efficient 2D DOA estimation algorithm for CPA by exploiting the RD polynomial root-finding technique. The proposed algorithm first maps CPA into FCPA and exploits the received data of two subarrays jointly, where the mutual information loss is avoided and simultaneously the improved estimation performance as well as enhanced DOFs can be achieved. In particular, the FCPA we constructed is an ambiguity-free array with high extraction efficiency. Furthermore, we convert the 2D total spectral search into one 1D polynomial root-finding process via reduced-dimension transformation, which significantly reduces the computational cost and simultaneously preserves the estimation accuracy. Simulations demonstrates the superiority of the proposed approach in regard to complexity, achievable DOFs and DOA estimation performance.



---

*Source: 2794387-2020-07-10.xml*