# Three-Dimensional Coprime Array for Massive MIMO: Array Configuration Design and 2D DOA Estimation

Authors: Pan Gong; Tanveer Ahmed; Jianfeng Li
Journal: Wireless Communications and Mobile Computing
 (2020)

---

## Abstract

In massive multiple-input multiple-output (MIMO) systems, it is critical to obtain the accurate direction of arrival (DOA) estimation. Conventional three-dimensional array mainly focuses on the uniform array. Due to the dense arrangement of the sensors, the array aperture is limited and severe mutual coupling effects arise. In this paper, a coprime cubic array (CCA) configuration design is presented, which is composed of two uniform cubic subarrays and can extend the interelement spacing with a selection of three pairs of coprime integers. Compared with uniform cubic array (UCA), CCA achieves the larger array aperture and less MC effects. And the analytical expression of Cramer–Rao Bound (CRB) for CCA is derived which verifies that the proposed CCA geometry outperforms the conventional UCA in two-dimensional (2D) DOA estimation performance in massive MIMO systems. Meanwhile, we propose a computationally efficient 2D DOA estimation algorithm with high accuracy for CCA. Specifically, we utilize array mapping to extract two uniform arrays from the nonuniform array by exploiting the relation derived from the signal subspace and the two directional matrices. Then, we operate a reduced dimension process on the uniform arrays and convert the 2D spectrum peak searching (SPS) problem into one-dimensional (1D) one, which significantly reduces the computational complexity. In addition, we employ the polynomial root finding technique with a lower complexity instead of 1D SPS to further relieve the computational complexity. Simultaneously, with coprime property, the phase ambiguity problem is solved, which results from the large interelement spacing. Numerical simulation results demonstrate that the proposed algorithm is very computationally efficient without degradation of DOA estimation performance.


---

## Body

## 1. Introduction

Nowadays, massive multiple-input multiple-output (MIMO), known as large-scale MIMO, is considered as one of the promising technologies in the development of future wireless communication systems. Massive MIMO attracts considerable attention due to the high throughput, enhanced link reliability, and improved spectral efficiency [1, 2]. Also, massive MIMO can effectively reduce the latency and robustness to interference [3, 4], which are important factors in wireless communication systems.Direction of arrival (DOA) estimation plays an important role in massive MIMO, since precise DOA estimation is vital for the base station (BS) to conduct downlink precoding or beamforming [5]. DOA estimation is also widely utilized in engineering applications, like radar, sonar, wireless communication, and navigation [6–9]. Many classic DOA estimation algorithms [10–13] and plenty of derived algorithms [14–16] have been proposed to solve the DOA estimation problem. Most of these existing algorithms mainly utilize the conventional uniform array configurations [17, 18] with the interelement spacing no larger than typical half-wavelength to tackle the problem of phase ambiguity [19]. However, in the massive MIMO systems, with a large number of antennas, severe mutual coupling (MC) effects arise, which significantly affect the DOA estimation accuracy.In recent years, a new configuration called coprime array [20] has been a hot research direction and attracts much attention since it can significantly enhance the degree of freedom (DOF) [21], improve the resolution [22], and relieve MC effects [23]. This motivates many studies of DOA estimation using coprime array over the decade. In [24], a total angular search based algorithm by employing the multiple signals classification (MUSIC) algorithm was proposed. By exploiting the coprime property, the phase ambiguity problem is tackled. In [25], an ambiguity-free DOA estimation algorithm was proposed, which exploits the total information including self-information and mutual information to eliminate the ambiguity problem with high accuracy of DOA estimates. However, due to the total angular search, this algorithm results in large computational complexity. To deal with this obstacle, a partial spectral search (PSS) algorithm was introduced which only performs spectrum searching in a restricted sector and hence lowers the computational cost but has no effect on DOA estimation performance [26]. In the case of multiple sources, this algorithm will cause the target matching error when eliminating the phase ambiguity. To tackle the problem, an improved DOA estimation algorithm based on root-MUSIC was proposed in [27] for coprime linear array (CLA), which employs the relation between steering matrices and signal subspaces of two subarrays to achieve DOA estimation.The algorithms mentioned above [24–27] were presented for one-dimensional (1D) DOA estimation, whereas, practically, two-dimensional (2D) DOA estimation possesses more importance, and various studies have been introduced [28–30]. In [28], the PSS algorithm for coprime planar array (CPA) was presented, which can reduce the computational complexity. A reduced dimension notion to achieve 2D DOA estimation was proposed in [29]. A generalized CPA geometry with more flexible array layouts was constructed for 2D DOA estimation [30], which can attain more DOFs compared with CPA configuration. However, spectrum peak searching (SPS) is involved in these mentioned works [28–30] and hence causes expensive computational cost.In this paper, on the basis of the 1D CLA and 2D CPA, we construct a three-dimensional (3D) array geometry called coprime cubic array (CCA). In the context of massive MIMO systems [31, 32], MC effects become important which have a vital influence on DOA estimation accuracy. The geometry of CCA possesses the extended array aperture and can effectively reduce MC effects. Especially, in the cases with small elevations, the CCA obtains superior DOA estimation performance. Simultaneously, the Cramer–Rao Bound (CRB) is derived which demonstrates that the proposed CCA geometry outperforms the conventional uniform cubic array (UCA) configuration [33–35]. In addition, we propose an array mapping and reduced dimension based on computationally efficient MUSIC (AMRD-MUSIC) algorithm to estimate 2D DOAs. Different from the conventional multiple-dimensional MUSIC (MD-MUSIC) algorithm, where 2D SPS involves a tremendous computation burden, in the proposed algorithm, we utilize array mapping to extract two uniform arrays from the nonuniform array by exploiting the relation derived from the signal subspace and the two directional matrices. Then, we operate a reduced dimension process on the uniform arrays and convert the 2D spectrum peak searching (SPS) problem into 1D one, which significantly reduces the computational complexity. In addition, we employ the polynomial root finding technique with a lower complexity instead of 1D SPS to further reduce the computational complexity. Meanwhile, according to the coprime property, the phase ambiguity problem is solved which results from large interelement spacing and the proposed algorithm obtains paired angles automatically. Numerical simulation results demonstrate that the proposed AMRD-MUSIC algorithm can significantly reduce the computational complexity cost with no degradation of DOA estimation performance.Specifically, we summarize the main contributions of this paper.(1)
In massive MIMO systems, conventional 3D array mainly focuses on the uniform array, which limits the array aperture and suffers from severe MC effects due to the dense location of sensors. We construct a CCA configuration which can achieve the larger array aperture and less MC effects compared with UCA. In addition, the analytical expression of CRB with CCA is derived and simulations demonstrate that the proposed CCA performs better DOA estimation performance than the conventional UCA configuration.(2)
We propose an AMRD-MUSIC algorithm for 2D DOA estimation that can achieve almost the same DOA estimation performance as classic MD-MUSIC algorithm but with lower computational complexity. Through array mapping, we extract two uniform arrays from the nonuniform array and operate a reduced dimension on these two arrays to reduce 2D SPS into 1D one, which can effectively lower the computational complexity. To further reduce complexity, we utilize the polynomial root finding technique with a lower complexity instead of SPS. Therefore, computational complexity obtains a significant decrease.(3)
The proposed algorithm can achieve the full DOFs and obtain DOA estimates with pairing automatically.Section2 introduces the array configuration and signal model. In Section 3 and Section 4, we present the proposed algorithm and analyse the performance of the proposed algorithm, respectively. Sections 5 and 6 provide numerical simulations and conclusions, respectively.Notations 1.
We use⋅T as the transpose, and ⋅H is utilized as the conjugate transpose. ⊙ signifies Khatri–Rao product and ⊗ denotes Kronecker product. E⋅ is expectation. Dm⋅ is a diagonal matrix that is formed of the m-th row of the matrix. angle⋅ is a phase operator.


## 2. Array Configurations and Signal Model

### 2.1. Uniform Cubic Array

UCA geometry is designed withM×N×J sensors, where M denotes the number of sensors in the x-axis, N and J denote the number of sensors in the y-axis and z-axis, respectively. The interelement spacing of UCA is dx=dy=dz=λ/2, which denotes the interelement spacing of x-axis, y-axis, and z-axis, respectively. λ is wavelength. Figure 1(a) shows an example of UCA with M=N=J=3.Figure 1
Array configurations. (a) The structure of UCAM=N=J=3 and (b) the structure of CCA M1=3,N1=2,J1=3andM2=2,N2=3,J2=2.
(a)
(b)


### 2.2. Coprime Cubic Array

For UCA configuration, the array aperture is limited by the interelement spacing no larger than conventional half-wavelength. Simultaneously, MC effects exist due to the dense location of sensors. Based on CLA and CPA, we propose a CCA configuration, which can effectively increase the array aperture and decrease the MC effects.Definition 1 (coprime cubic array (CCA)).
The CCA consists of two subarrays. The first subarray is composed ofM1×N1×J1 sensors with the interelement spacing of dx1=M2λ/2 in x-axis direction, dy1=N2λ/2 in y-axis direction, and dz1=J2λ/2 in z-axis direction. For the second subarray, it contains M2×N2×J2 elements and dx2=M1λ/2, dy2=N1λ/2, dz2=J1λ/2, where M1,M2,N1,N2 and J1,J2 are coprime integer pairs and stand for the number of sensors in x-axis, y-axis, and z-axis, respectively. Therefore, there are TCCA=M1N1J1+M2N2J2−1 sensors in total. Figure 1(b) is an example of CCA geometry with M1=3,N1=2, J1=3 and M2=2, N2=3, J2=2.
In the following part, we compare the array performance in array aperture and MC effects of UCA and CCA, respectively.


### 2.3. Array Aperture Comparison

In this subsection, we compare the array aperture of UCA and CCA in the same cases. For a fair comparison, we assume that the number of sensors in UCA and CCA is basically the same. According to [36], we can compute the array apertures of these arrays, which are denoted as ΓUCA and ΓCCA:(1)ΓUCA=M−1N−1J−1×λ38,ΓCCA=maxM1−1M2N1−1N2J1−1J2λ38,M2−1M1N2−1N1J2−1J1λ38.For example, we set the UCA withM=4,N=4,J=6, so the total number of sensors in UCA is TUCA=M×N×J=4×4×6=96. One subarray of CCA is of M1=5,N1=2,J1=3 and the other subarray is of M2=4,N2=3,J2=5. The number of sensors in CCA is TCCA=M1×N1×J1+M2×N2×J2−1=89 in total. We can note that the CCA has fewer number of sensors than UCA. As a result, for the two given arrays, the array aperture can be obtained as ΓUCA=5.625λ3, ΓCCA=90λ3. So we have(2)ΓUCA<ΓCCA.It can be seen that the proposed CCA configuration can effectively enlarge the array aperture with fewer sensors.


### 2.4. Mutual Coupling Effects Comparison

Array signal models usually neglect the MC between the sensors. In practical application, the MC effects between the closely located sensors should be considered. Without loss of generality, we employ a coupling matrixC∈ℂT×T, where T denotes the total number of sensors in the array. In general, the expression of C is very complicated. According to [37], for given UCA and CCA configurations, C can be denoted approximately by a B-banded Toeplitz matrix as follows:(3)Ci,j=cdi−dj,di−dj≤B,0,di−dj>B,where the value of Ci,j is inversely proportional to sensor distance and satisfies 1=c0>c1>c2>⋯>cB>cB+1=0. B is a constant and di denotes the positions of sensors.According to [37], the coupling leakage can be exploited to evaluate the total MC(4)L=C−diagCFCF,where C−diagCF is the amount of MC. It is indicated that the smaller L is, the less MC is.We compare the coupling leakages of the given UCA ofM=4,N=4,J=6 and CCA of M1=5,N1=2,J1=3 and M2=4,N2=3,J2=5. The total number of sensors of UCA and CCA is TUCA=M×N×J=4×4×6=96 and TCCA=M1×N1×J1+M2×N2×J2−1=89, respectively. Hence, we have CUCA∈ℂ96×96 and CCCA∈ℂ89×89. So the coupling coefficients can be computed by the following:(5)cd=1,d=0,0.3ejπ/3⋅e−jd−1π/8d,1≤d≤B,0,d>B,where B=3 and d denotes the distance between the random two sensors. Then, we compute the coupling leakages of UCA and CCA as LTUCA=0.7588 and LTCCA=0.2829, respectively. It can see that the proposed CCA can effectively decrease the MC effects.(6)LCCA<LUCA.Table1 shows the comparison of array aperture and coupling leakage for UCA and CCA in the same conditions. From Table 1, we can see that the proposed CCA configuration has the least sensors but outperforms the UCA in array aperture and MC effects. In addition, the CCA achieves better DOA estimation performance that is shown in Section 5. Therefore, in the following part, we employ the CCA into consideration.Table 1
Performance of the array configurations.
Configuration
Total number of sensors
Array aperture
Coupling leakage
UCA
96
5.625
λ
3
0.7588
CCA
89
90
λ
3
0.2829


### 2.5. Signal Model

Assume that there areK far-field uncorrelated narrowband incident signals impinging on the CCA from θk,ϕkθk∈0°,90°,ϕk∈0°,180°,k=1,2,…,K, where θk and ϕk denote the DOA (the elevation angle and azimuth angle) of the k-th target [19]. And we define K<minM1N1J1,M2N2J2, uk=sinθkcosϕk∈−1,1, vk=sinθksinϕk∈0,1, and wk=cosθk∈0,1. In the following, we take the subarray that has Mi×Ni×Ji i=1,2 sensors to derive the proposed algorithm in the following part.We denote the received signal as follows [28]:(7)xit=Aist+nit,i=1,2,where st=s1t,s2t,…,sKt∈ℂK×L represents the signal matrix and L denotes the number of snapshots. ni∈ℂMiNiJi×L is the additive white Gaussian noise matrix, whose variance and mean are σn2 and zero, respectively. Ai∈ℂMiNiJi×K denotes directional matrix of the i-th subarray and is presented by the following:(8)Ai=ax,iu1⊗ay,iv1⊗az,iw1,…,ax,iuK⊗ay,ivK⊗az,iwK,where ax,iuk=1,ej2πdxiuk/λ,…,ej2πMi−1dxiuk/λT,ay,ivk=1,ej2πdyivk/λ,…,ej2πNi−1dyivk/λT and az,iwk=1,ej2πdziwk/λ,…,ej2πJi−1dziwk/λTk=1,2,...,K correspond to steering vectors of x-axis, y-axis, and z-axis of the subarray i.Denote the directional matrix in (7) as follows:(9)Ai=Ax,i⊙Ay,i⊙Az,i=Az,iD1Ay,iD1Ax,iAz,iD1Ay,iD2Ax,i⋮Az,iDNiAy,iDMiAx,i,where Ax,i denotes directional matrix of x-axis of the subarray i, and Ay,i,Az,i stand for the matrices of y-axis and z-axis, respectively, DniAy,in=1,2,…,N&i=1,2 denotes a diagonal matrix which is formed of the n-th row of the matrix Ay,i for the subarray i, and DmiAx,im=1,2,…,M&i=1,2 denotes a diagonal matrix which is formed of the m-th row of the matrix Ax,i for the subarray i.


## 2.1. Uniform Cubic Array

UCA geometry is designed withM×N×J sensors, where M denotes the number of sensors in the x-axis, N and J denote the number of sensors in the y-axis and z-axis, respectively. The interelement spacing of UCA is dx=dy=dz=λ/2, which denotes the interelement spacing of x-axis, y-axis, and z-axis, respectively. λ is wavelength. Figure 1(a) shows an example of UCA with M=N=J=3.Figure 1
Array configurations. (a) The structure of UCAM=N=J=3 and (b) the structure of CCA M1=3,N1=2,J1=3andM2=2,N2=3,J2=2.
(a)
(b)


## 2.2. Coprime Cubic Array

For UCA configuration, the array aperture is limited by the interelement spacing no larger than conventional half-wavelength. Simultaneously, MC effects exist due to the dense location of sensors. Based on CLA and CPA, we propose a CCA configuration, which can effectively increase the array aperture and decrease the MC effects.Definition 1 (coprime cubic array (CCA)).
The CCA consists of two subarrays. The first subarray is composed ofM1×N1×J1 sensors with the interelement spacing of dx1=M2λ/2 in x-axis direction, dy1=N2λ/2 in y-axis direction, and dz1=J2λ/2 in z-axis direction. For the second subarray, it contains M2×N2×J2 elements and dx2=M1λ/2, dy2=N1λ/2, dz2=J1λ/2, where M1,M2,N1,N2 and J1,J2 are coprime integer pairs and stand for the number of sensors in x-axis, y-axis, and z-axis, respectively. Therefore, there are TCCA=M1N1J1+M2N2J2−1 sensors in total. Figure 1(b) is an example of CCA geometry with M1=3,N1=2, J1=3 and M2=2, N2=3, J2=2.
In the following part, we compare the array performance in array aperture and MC effects of UCA and CCA, respectively.


## 2.3. Array Aperture Comparison

In this subsection, we compare the array aperture of UCA and CCA in the same cases. For a fair comparison, we assume that the number of sensors in UCA and CCA is basically the same. According to [36], we can compute the array apertures of these arrays, which are denoted as ΓUCA and ΓCCA:(1)ΓUCA=M−1N−1J−1×λ38,ΓCCA=maxM1−1M2N1−1N2J1−1J2λ38,M2−1M1N2−1N1J2−1J1λ38.For example, we set the UCA withM=4,N=4,J=6, so the total number of sensors in UCA is TUCA=M×N×J=4×4×6=96. One subarray of CCA is of M1=5,N1=2,J1=3 and the other subarray is of M2=4,N2=3,J2=5. The number of sensors in CCA is TCCA=M1×N1×J1+M2×N2×J2−1=89 in total. We can note that the CCA has fewer number of sensors than UCA. As a result, for the two given arrays, the array aperture can be obtained as ΓUCA=5.625λ3, ΓCCA=90λ3. So we have(2)ΓUCA<ΓCCA.It can be seen that the proposed CCA configuration can effectively enlarge the array aperture with fewer sensors.


## 2.4. Mutual Coupling Effects Comparison

Array signal models usually neglect the MC between the sensors. In practical application, the MC effects between the closely located sensors should be considered. Without loss of generality, we employ a coupling matrixC∈ℂT×T, where T denotes the total number of sensors in the array. In general, the expression of C is very complicated. According to [37], for given UCA and CCA configurations, C can be denoted approximately by a B-banded Toeplitz matrix as follows:(3)Ci,j=cdi−dj,di−dj≤B,0,di−dj>B,where the value of Ci,j is inversely proportional to sensor distance and satisfies 1=c0>c1>c2>⋯>cB>cB+1=0. B is a constant and di denotes the positions of sensors.According to [37], the coupling leakage can be exploited to evaluate the total MC(4)L=C−diagCFCF,where C−diagCF is the amount of MC. It is indicated that the smaller L is, the less MC is.We compare the coupling leakages of the given UCA ofM=4,N=4,J=6 and CCA of M1=5,N1=2,J1=3 and M2=4,N2=3,J2=5. The total number of sensors of UCA and CCA is TUCA=M×N×J=4×4×6=96 and TCCA=M1×N1×J1+M2×N2×J2−1=89, respectively. Hence, we have CUCA∈ℂ96×96 and CCCA∈ℂ89×89. So the coupling coefficients can be computed by the following:(5)cd=1,d=0,0.3ejπ/3⋅e−jd−1π/8d,1≤d≤B,0,d>B,where B=3 and d denotes the distance between the random two sensors. Then, we compute the coupling leakages of UCA and CCA as LTUCA=0.7588 and LTCCA=0.2829, respectively. It can see that the proposed CCA can effectively decrease the MC effects.(6)LCCA<LUCA.Table1 shows the comparison of array aperture and coupling leakage for UCA and CCA in the same conditions. From Table 1, we can see that the proposed CCA configuration has the least sensors but outperforms the UCA in array aperture and MC effects. In addition, the CCA achieves better DOA estimation performance that is shown in Section 5. Therefore, in the following part, we employ the CCA into consideration.Table 1
Performance of the array configurations.
Configuration
Total number of sensors
Array aperture
Coupling leakage
UCA
96
5.625
λ
3
0.7588
CCA
89
90
λ
3
0.2829


## 2.5. Signal Model

Assume that there areK far-field uncorrelated narrowband incident signals impinging on the CCA from θk,ϕkθk∈0°,90°,ϕk∈0°,180°,k=1,2,…,K, where θk and ϕk denote the DOA (the elevation angle and azimuth angle) of the k-th target [19]. And we define K<minM1N1J1,M2N2J2, uk=sinθkcosϕk∈−1,1, vk=sinθksinϕk∈0,1, and wk=cosθk∈0,1. In the following, we take the subarray that has Mi×Ni×Ji i=1,2 sensors to derive the proposed algorithm in the following part.We denote the received signal as follows [28]:(7)xit=Aist+nit,i=1,2,where st=s1t,s2t,…,sKt∈ℂK×L represents the signal matrix and L denotes the number of snapshots. ni∈ℂMiNiJi×L is the additive white Gaussian noise matrix, whose variance and mean are σn2 and zero, respectively. Ai∈ℂMiNiJi×K denotes directional matrix of the i-th subarray and is presented by the following:(8)Ai=ax,iu1⊗ay,iv1⊗az,iw1,…,ax,iuK⊗ay,ivK⊗az,iwK,where ax,iuk=1,ej2πdxiuk/λ,…,ej2πMi−1dxiuk/λT,ay,ivk=1,ej2πdyivk/λ,…,ej2πNi−1dyivk/λT and az,iwk=1,ej2πdziwk/λ,…,ej2πJi−1dziwk/λTk=1,2,...,K correspond to steering vectors of x-axis, y-axis, and z-axis of the subarray i.Denote the directional matrix in (7) as follows:(9)Ai=Ax,i⊙Ay,i⊙Az,i=Az,iD1Ay,iD1Ax,iAz,iD1Ay,iD2Ax,i⋮Az,iDNiAy,iDMiAx,i,where Ax,i denotes directional matrix of x-axis of the subarray i, and Ay,i,Az,i stand for the matrices of y-axis and z-axis, respectively, DniAy,in=1,2,…,N&i=1,2 denotes a diagonal matrix which is formed of the n-th row of the matrix Ay,i for the subarray i, and DmiAx,im=1,2,…,M&i=1,2 denotes a diagonal matrix which is formed of the m-th row of the matrix Ax,i for the subarray i.


## 3. 2D DOA Estimation Algorithm

Conventional MD-MUSIC algorithm can obtain DOA estimates with high accuracy by selecting a refined search step, which causes expensive computational complexity. To tackle the problems, an AMRD-MUSIC algorithm is proposed.


### 3.1. MD-MUSIC Spectrum

We concatenate the array measurementsx1t and x2t as follows:(10)xt=x1tx2t=A1A2st+n1tn2t=Ast+nt,where A=A1T,A2TT and nt=n1Tt,n2TtT. A1∈ℂM1N1J1×K and A2∈ℂM2N2J2×K Then, the covariance matrix of the total array can be computed by the following: (11)R^=1L∑t=1LxtxHt=UsΣsUsH+UnΣnUnH,where Un∈ℂM1N1J1+M2N2J2×M1N1J1+M2N2J2−K refers to noise subspace and Us∈ℂM1N1J1+M2N2J2×K signifies signal subspace. Σs∈ℂK×K is a diagonal matrix with the K largest eigenvalues of R^ and Σn contains the remaining eigenvalues. Here, the symbol of ( ^ ) means the estimation of the covariance matrix R.Then, according to 1D MUSIC spectrum [12], the MD-MUSIC spectrum can be solved as follows: (12)fu,v,w=1aHu,v,wUnUnHau,v,w=1a1Hu,v,wa2Hu,v,wUnUnHa1u,v,wa2u,v,w.


### 3.2. Array Mapping

In this subsection, we utilize array mapping to extract two uniform arrays from the nonuniform array by exploiting the relation derived from the signal subspace and the two directional submatrices.As the signal subspace spans the same space with the steering matrix doing, we have(13)Us=A1A2T,where T∈ℂK×K represents a nonsingular matrix. Then, the signal subspace can be decomposed into two parts:(14)Us=Us1Us2=A1A2T,where Us1∈ℂM1N1J1×K and Us2∈ℂM2N2J2×K.Then, the transformation matrices can be defined as follows [38]:(15)H1=Us2Us1+=A2TT−1A1+=A2A1+,H2=Us1Us2+=A1TT−1A2+=A1A2+,where ⋅+ denotes the pseudoinverse operation. From (15), we have the following:(16)A2=H1A1,A1=H2A2,and the steering vectors of a1u,v,w and a2u,v,w also satisfy this relation in (16), so we have the following:(17)a2u,v,w=H1a1u,v,w,a1u,v,w=H2a2u,v,w.By utilizing (17), we transform (13) into the following: (18)f1u,v,w=1a1Hu,v,wa1Hu,v,wH1HUnUnHa1u,v,wH1a1u,v,w=1a1Hu,v,w∏1∏1Ha1u,v,w,(19)f2u,v,w=1a2Hu,v,wH2Ha2Hu,v,wUnUnHH2a2u,v,wa2u,v,w=1a2Hu,v,w∏2∏2Ha2u,v,w,where ∏1=IM1N1J1H1HUn and ∏2=H2HIM2N2J2Un.


### 3.3. Reduced Dimension Processing

In this part, we employ the double reduced dimension processing on (16) and (17), respectively, to transform the 2D SPS into 1D one. Therefore, the computational complexity can be significantly reduced.Asaiu,v,w=ax,iu⊗ay,iv⊗az,iw=ax,iu⊗ay,iv⊗az,iw=ax,iu⊗INiJiay,iv⊗az,iw, then we have(20)a1u,v,w=ax,1u⊗IN1J1ay,1v⊗az,1w,(21)a2u,v,w=ax,2u⊗IN2J2ay,2v⊗az,2w,where IN1J1 and IN2J2 are both identity matrices.Then, we define the functionQiu,v,wi=1,2 as follows:(22)Qiu,v,w=aiHu,v,w∏i∏iHaiu,v,w,which transforms the 2D SPS to find the minimum value of (22). According to (20) and (21), we have the following:(23)Qiu,v,w=ax,iu⊗ay,iv⊗az,iwH∏i∏iHax,iu⊗ay,iv⊗az,iw=ax,iu⊗INiJiay,iv⊗az,iwH∏i∏iHax,iu⊗INiJiay,iv⊗az,iwH=ay,iv⊗az,iwHax,iu⊗INiJiH∏i∏iHax,iu⊗INiJiay,iv⊗az,iw=ay,iv⊗az,iwHGiuay,iv⊗az,iw=fiHv,wGiufiv,w,where Giu=ax,iu⊗INiJiH∏i∏iHax,iu⊗INiJi, fiv,w=ay,iv⊗az,iw, and fiv,w∈ℂNiJi×1. Note that (23) is an optimization problem. Based on eHfiHv,w=1, where e=1,0,…,0T∈ℂNiJi, we can eliminate the trivial solution fiHv,w=0. This optimization problem can be reconstructed as follows:(24)minu,v,wfiHv,wGiufiv,ws.t.eHfiHv,w=1.Then, we can construct the cost function as follows:(25)Lu,v,w=fiHv,wGiufiv,w−λeHfiHv,w−1,where λ is a constant. Then, we have the following:(26)∂∂fvi,wiLu,v,w=2Giufiv,w+λe.According to (26), we have fiv,w=ξGi−1ue, where ξ is a constant. Combining (24), we have the following:(27)ξ=1eHGi−1ue.Then,fiv,w can be denoted as follows:(28)fiv,w=Gi−1ueeHGi−1ue.Based on (24) and (28), we can obtain the following:(29)u^k=argminu1eHGi−1ue=argmaxeHGi−1ue.According to (29), the 2D SPS is transformed into 1D one by utilizing a reduced dimension process. To further reduce complexity, we employ the polynomial rooting technique with the lower computational complexity to obtain the final DOA estimates instead of the 1D SPS process.


### 3.4. Polynomial Root Finding Process

Definezx,i=e−j2πdiu/λ,i=1,2 and denote ax,iu as follows:(30)ax,iu=azx,i=1,zx,i,zx,i2,…,zx,iMi−1T.Then, we rewriteGiu as follows:(31)Gzx,i=aHzx,i⊗INiJiH∏i∏iHazx,i⊗INiJi.Then, we employ polynomial root finding technique instead of 1D SPS.(32)aHzx,i⊗INiJiH∏zx,i∏Hzx,iazx,i⊗INiJi=0.Ifzx,i cannot match with one source signal and if(33)Rank∏i≥NiJi,the matrix ∏i is invertible [38].In order to achieve the polynomial (32), zx,i satisfies(34)Dzx,i=det∏zx,i=0.Then, we can obtain the finalu^ik by finding K roots which are closest to the unit circle of the polynomial Dzx,i,zx,ikk=1,2,…,K,(35)u^ik=anglezx,ikλ2πdi.In order to obtain the estimation ofv^ik or w^ik, the similar reduced dimensional polynomial root finding technique can be employed again. We can reform the whole steering matrix as Ai′=ay,iv1⊗ax,iu1⊗az,iw1,…,ay,ivK⊗ax,iuK⊗az,iwK, which contains the same elements as the matrix Ai and only differs in elements positions. Then, a reconstructed received signal can be obtained. Then, we operate the proposed array mapping, reduced dimension, and polynomial rooting technique again to attain the estimate of v^ik(36)v^ik=anglezy,ikλ2πdi.Similarly, we can achievew^ik. Phase ambiguity problem, resulting from the larger element spacing than half-wavelength, will be solved by coprime property in the following section.


### 3.5. Ambiguity Elimination

Firstly, the phase ambiguity problem is illustrated. Then, the coprime property is exploited to tackle the problem.Assume that only one signal located withθτ,ϕτ impinges on the CCA subarray and the subarray consists of Mi×Ni×Ji sensors. θa,ϕa is assumed as one of the ambiguous DOAs of the real angle θτ,ϕτ.According to [38],(37)mod2πλdxiuτ+2kuπ=2πλdxiuτ,mod2πλdyivτ+2kvπ=2πλdyivτ,mod2πλdziwτ+2kwπ=2πλdziwτ,where the mod operation depends on the period of the exponential function 2π. So it can be noted that the phase difference between the real and ambiguous DOAs is 2kπk∈ℤ [39] and ku,kv,kw∈ℤ are integers. And dxi=Njλ/2, dyi=Mjλ/2, dzi=Jjλ/2i,j∈1,2,i≠j denote the element spacing among three axes, respectively. uτ=sinθτcosϕτ,vτ=sinθτsinϕτ, wτ=cosθτ.Then, we employ the obtained estimates to recover all estimates [28],(38)2πdxiuτ−uaλ=2kuπ,(39)2πdyivτ−vaλ=2kvπ,(40)2πdziwτ−waλ=2kwπ,where ua=sinθacosϕa, va=sinθasinϕa, and wa=cosθa are defined and they denote the ambiguous values.Asθ∈0,π/2 and ϕ∈0,π, we have ua∈−1,1, va∈0,1, wa∈0,1. Simultaneously, ua2+vb2∈0,1 and ua2+va2+wa2=1 demand to satisfy each estimate.As the interelement spacing of conventional UCA is limited to half-wavelength,ku, kv, kw can only be the value of zero, which means no ambiguity problem exists. However, in CCA, due to the interelement spacing larger than half-wavelength, there are other pairs to make (37)–(39) hold besides the infeasible ku, kv, and kw pairs, whereas only one pair is the candidate related to the actual DOA. Coprime property is utilized to solve the phase ambiguity problem and achieve true DOA.Assumeθa,ϕa as one of the ambiguous DOAs that appears in one of the subarrays. Simultaneously, define θt,ϕt as the corresponding actual DOA. There exists the following:(41)uτ−ua=2kuiMj,vτ−va=2kviNj,wτ−wa=2kwiJj.It exists thatkui∈−Mj,Mj, kvi∈−Nj/2,Nj/2, kwi∈−Jj/2,Jj/2, i,j∈1,2,i≠j [38], so we can obtain(42)ku1M2=ku2M1,kv1N2=kv2N1,kw1J2=kw2J1.SinceN1, N2, M1,M2 and J1,J2 are coprime integers pairs, respectively, (41) can hold only in the case of ku1=ku2=0, kv1=kv2=0, kw1=kw2=0. It is indicated that the ambiguity problem is solved and actual DOA estimates can be achieved.In practice, we only need two variables ofu^k,v^k,γ^k to achieve DOAs θk,ϕk estimates, where u^k,v^k, γ^k are the estimates of uk,vk,γkk=1,2,…,K.The final estimate ofθk,ϕk can be calculated by the following:(43)θ^k=arcsinu^k2+v^k2+arccosw^k2,ϕ^k=arctanu^k/v^k+arccosu^k/1−w^k22.


### 3.6. Detailed Steps of the Proposed Algorithm

The steps are provided in detail as follows:Step 1: compute the covariance of the received signal data to obtain the MD-MUSIC spectrum function (12).Step 2: according to the signal subspace, constructH^1 and H^2 by (15), with which we can obtain the reconstructed spectrum function according to (19).Step 3: by array mapping, we extract the two uniform arrays (19) from the nonuniform array (12) according to the relation derived from two directional matrices.Step 4: according to (20) and (21), transform the 2D SPS process into 1D SPS.Step 5: according to (32), convert 1D PSS to polynomial root finding process and obtain the K roots of u^ik. Similarly, we can achieve v^ik and w^ik.Step 6: by (41), all the ambiguous DOA estimates are recovered, and by searching for the nearest ones, the actual estimations of uk, vk and wk can be achieved.Step 7: according to (43), the final DOAs estimates can be obtained.


## 3.1. MD-MUSIC Spectrum

We concatenate the array measurementsx1t and x2t as follows:(10)xt=x1tx2t=A1A2st+n1tn2t=Ast+nt,where A=A1T,A2TT and nt=n1Tt,n2TtT. A1∈ℂM1N1J1×K and A2∈ℂM2N2J2×K Then, the covariance matrix of the total array can be computed by the following: (11)R^=1L∑t=1LxtxHt=UsΣsUsH+UnΣnUnH,where Un∈ℂM1N1J1+M2N2J2×M1N1J1+M2N2J2−K refers to noise subspace and Us∈ℂM1N1J1+M2N2J2×K signifies signal subspace. Σs∈ℂK×K is a diagonal matrix with the K largest eigenvalues of R^ and Σn contains the remaining eigenvalues. Here, the symbol of ( ^ ) means the estimation of the covariance matrix R.Then, according to 1D MUSIC spectrum [12], the MD-MUSIC spectrum can be solved as follows: (12)fu,v,w=1aHu,v,wUnUnHau,v,w=1a1Hu,v,wa2Hu,v,wUnUnHa1u,v,wa2u,v,w.


## 3.2. Array Mapping

In this subsection, we utilize array mapping to extract two uniform arrays from the nonuniform array by exploiting the relation derived from the signal subspace and the two directional submatrices.As the signal subspace spans the same space with the steering matrix doing, we have(13)Us=A1A2T,where T∈ℂK×K represents a nonsingular matrix. Then, the signal subspace can be decomposed into two parts:(14)Us=Us1Us2=A1A2T,where Us1∈ℂM1N1J1×K and Us2∈ℂM2N2J2×K.Then, the transformation matrices can be defined as follows [38]:(15)H1=Us2Us1+=A2TT−1A1+=A2A1+,H2=Us1Us2+=A1TT−1A2+=A1A2+,where ⋅+ denotes the pseudoinverse operation. From (15), we have the following:(16)A2=H1A1,A1=H2A2,and the steering vectors of a1u,v,w and a2u,v,w also satisfy this relation in (16), so we have the following:(17)a2u,v,w=H1a1u,v,w,a1u,v,w=H2a2u,v,w.By utilizing (17), we transform (13) into the following: (18)f1u,v,w=1a1Hu,v,wa1Hu,v,wH1HUnUnHa1u,v,wH1a1u,v,w=1a1Hu,v,w∏1∏1Ha1u,v,w,(19)f2u,v,w=1a2Hu,v,wH2Ha2Hu,v,wUnUnHH2a2u,v,wa2u,v,w=1a2Hu,v,w∏2∏2Ha2u,v,w,where ∏1=IM1N1J1H1HUn and ∏2=H2HIM2N2J2Un.


## 3.3. Reduced Dimension Processing

In this part, we employ the double reduced dimension processing on (16) and (17), respectively, to transform the 2D SPS into 1D one. Therefore, the computational complexity can be significantly reduced.Asaiu,v,w=ax,iu⊗ay,iv⊗az,iw=ax,iu⊗ay,iv⊗az,iw=ax,iu⊗INiJiay,iv⊗az,iw, then we have(20)a1u,v,w=ax,1u⊗IN1J1ay,1v⊗az,1w,(21)a2u,v,w=ax,2u⊗IN2J2ay,2v⊗az,2w,where IN1J1 and IN2J2 are both identity matrices.Then, we define the functionQiu,v,wi=1,2 as follows:(22)Qiu,v,w=aiHu,v,w∏i∏iHaiu,v,w,which transforms the 2D SPS to find the minimum value of (22). According to (20) and (21), we have the following:(23)Qiu,v,w=ax,iu⊗ay,iv⊗az,iwH∏i∏iHax,iu⊗ay,iv⊗az,iw=ax,iu⊗INiJiay,iv⊗az,iwH∏i∏iHax,iu⊗INiJiay,iv⊗az,iwH=ay,iv⊗az,iwHax,iu⊗INiJiH∏i∏iHax,iu⊗INiJiay,iv⊗az,iw=ay,iv⊗az,iwHGiuay,iv⊗az,iw=fiHv,wGiufiv,w,where Giu=ax,iu⊗INiJiH∏i∏iHax,iu⊗INiJi, fiv,w=ay,iv⊗az,iw, and fiv,w∈ℂNiJi×1. Note that (23) is an optimization problem. Based on eHfiHv,w=1, where e=1,0,…,0T∈ℂNiJi, we can eliminate the trivial solution fiHv,w=0. This optimization problem can be reconstructed as follows:(24)minu,v,wfiHv,wGiufiv,ws.t.eHfiHv,w=1.Then, we can construct the cost function as follows:(25)Lu,v,w=fiHv,wGiufiv,w−λeHfiHv,w−1,where λ is a constant. Then, we have the following:(26)∂∂fvi,wiLu,v,w=2Giufiv,w+λe.According to (26), we have fiv,w=ξGi−1ue, where ξ is a constant. Combining (24), we have the following:(27)ξ=1eHGi−1ue.Then,fiv,w can be denoted as follows:(28)fiv,w=Gi−1ueeHGi−1ue.Based on (24) and (28), we can obtain the following:(29)u^k=argminu1eHGi−1ue=argmaxeHGi−1ue.According to (29), the 2D SPS is transformed into 1D one by utilizing a reduced dimension process. To further reduce complexity, we employ the polynomial rooting technique with the lower computational complexity to obtain the final DOA estimates instead of the 1D SPS process.


## 3.4. Polynomial Root Finding Process

Definezx,i=e−j2πdiu/λ,i=1,2 and denote ax,iu as follows:(30)ax,iu=azx,i=1,zx,i,zx,i2,…,zx,iMi−1T.Then, we rewriteGiu as follows:(31)Gzx,i=aHzx,i⊗INiJiH∏i∏iHazx,i⊗INiJi.Then, we employ polynomial root finding technique instead of 1D SPS.(32)aHzx,i⊗INiJiH∏zx,i∏Hzx,iazx,i⊗INiJi=0.Ifzx,i cannot match with one source signal and if(33)Rank∏i≥NiJi,the matrix ∏i is invertible [38].In order to achieve the polynomial (32), zx,i satisfies(34)Dzx,i=det∏zx,i=0.Then, we can obtain the finalu^ik by finding K roots which are closest to the unit circle of the polynomial Dzx,i,zx,ikk=1,2,…,K,(35)u^ik=anglezx,ikλ2πdi.In order to obtain the estimation ofv^ik or w^ik, the similar reduced dimensional polynomial root finding technique can be employed again. We can reform the whole steering matrix as Ai′=ay,iv1⊗ax,iu1⊗az,iw1,…,ay,ivK⊗ax,iuK⊗az,iwK, which contains the same elements as the matrix Ai and only differs in elements positions. Then, a reconstructed received signal can be obtained. Then, we operate the proposed array mapping, reduced dimension, and polynomial rooting technique again to attain the estimate of v^ik(36)v^ik=anglezy,ikλ2πdi.Similarly, we can achievew^ik. Phase ambiguity problem, resulting from the larger element spacing than half-wavelength, will be solved by coprime property in the following section.


## 3.5. Ambiguity Elimination

Firstly, the phase ambiguity problem is illustrated. Then, the coprime property is exploited to tackle the problem.Assume that only one signal located withθτ,ϕτ impinges on the CCA subarray and the subarray consists of Mi×Ni×Ji sensors. θa,ϕa is assumed as one of the ambiguous DOAs of the real angle θτ,ϕτ.According to [38],(37)mod2πλdxiuτ+2kuπ=2πλdxiuτ,mod2πλdyivτ+2kvπ=2πλdyivτ,mod2πλdziwτ+2kwπ=2πλdziwτ,where the mod operation depends on the period of the exponential function 2π. So it can be noted that the phase difference between the real and ambiguous DOAs is 2kπk∈ℤ [39] and ku,kv,kw∈ℤ are integers. And dxi=Njλ/2, dyi=Mjλ/2, dzi=Jjλ/2i,j∈1,2,i≠j denote the element spacing among three axes, respectively. uτ=sinθτcosϕτ,vτ=sinθτsinϕτ, wτ=cosθτ.Then, we employ the obtained estimates to recover all estimates [28],(38)2πdxiuτ−uaλ=2kuπ,(39)2πdyivτ−vaλ=2kvπ,(40)2πdziwτ−waλ=2kwπ,where ua=sinθacosϕa, va=sinθasinϕa, and wa=cosθa are defined and they denote the ambiguous values.Asθ∈0,π/2 and ϕ∈0,π, we have ua∈−1,1, va∈0,1, wa∈0,1. Simultaneously, ua2+vb2∈0,1 and ua2+va2+wa2=1 demand to satisfy each estimate.As the interelement spacing of conventional UCA is limited to half-wavelength,ku, kv, kw can only be the value of zero, which means no ambiguity problem exists. However, in CCA, due to the interelement spacing larger than half-wavelength, there are other pairs to make (37)–(39) hold besides the infeasible ku, kv, and kw pairs, whereas only one pair is the candidate related to the actual DOA. Coprime property is utilized to solve the phase ambiguity problem and achieve true DOA.Assumeθa,ϕa as one of the ambiguous DOAs that appears in one of the subarrays. Simultaneously, define θt,ϕt as the corresponding actual DOA. There exists the following:(41)uτ−ua=2kuiMj,vτ−va=2kviNj,wτ−wa=2kwiJj.It exists thatkui∈−Mj,Mj, kvi∈−Nj/2,Nj/2, kwi∈−Jj/2,Jj/2, i,j∈1,2,i≠j [38], so we can obtain(42)ku1M2=ku2M1,kv1N2=kv2N1,kw1J2=kw2J1.SinceN1, N2, M1,M2 and J1,J2 are coprime integers pairs, respectively, (41) can hold only in the case of ku1=ku2=0, kv1=kv2=0, kw1=kw2=0. It is indicated that the ambiguity problem is solved and actual DOA estimates can be achieved.In practice, we only need two variables ofu^k,v^k,γ^k to achieve DOAs θk,ϕk estimates, where u^k,v^k, γ^k are the estimates of uk,vk,γkk=1,2,…,K.The final estimate ofθk,ϕk can be calculated by the following:(43)θ^k=arcsinu^k2+v^k2+arccosw^k2,ϕ^k=arctanu^k/v^k+arccosu^k/1−w^k22.


## 3.6. Detailed Steps of the Proposed Algorithm

The steps are provided in detail as follows:Step 1: compute the covariance of the received signal data to obtain the MD-MUSIC spectrum function (12).Step 2: according to the signal subspace, constructH^1 and H^2 by (15), with which we can obtain the reconstructed spectrum function according to (19).Step 3: by array mapping, we extract the two uniform arrays (19) from the nonuniform array (12) according to the relation derived from two directional matrices.Step 4: according to (20) and (21), transform the 2D SPS process into 1D SPS.Step 5: according to (32), convert 1D PSS to polynomial root finding process and obtain the K roots of u^ik. Similarly, we can achieve v^ik and w^ik.Step 6: by (41), all the ambiguous DOA estimates are recovered, and by searching for the nearest ones, the actual estimations of uk, vk and wk can be achieved.Step 7: according to (43), the final DOAs estimates can be obtained.


## 4. Performance Analysis of the Algorithm

### 4.1. DOF Comparison

The proposed AMRD-MUSIC algorithm can obtainM1N1J1+M2N2J2−1 DOFs, which indicates that the proposed algorithm can achieve M1N1J1+M2N2J2−2 signals by employing M1N1J1+M2N2J2−1 antennas.(44)DOFCCAAMRD−MUSIC=M1N1J1+M2N2J2−1.The algorithms in [24, 26] process the two subarrays separately so that they have a great loss in DOFs since they only exploit the autoinformation. The number of achievable DOFs by using the method in [24, 26] can be denoted as follows:(45)DOFCCAalgorithms in24,26=minM1N1J1,M2N2J2,where M1=N1=J1>M2=N2=J2.Thus, the proposed algorithms can detect many more signals than conventional algorithms with the same sensors. From (44) and (45), we can note the relationship of DOFs of different algorithms as follows:(46)DOFCCAalgorithms in24,26<DOFCCAARMD−MUSIC.Figure2 depicts the maximum number of achievable DOFs to illustrate the signals explicitly without ambiguous DOA estimations, where the number of snapshots is set to be L=500 with noise-free, the total number of signals is K=28, and the total number of sensors is T=3×2×3+2×3×2−1=29. The elevation angles and azimuth angles are set to be θ=20∘,20∘,20∘,20∘,20∘,20∘,20∘,40∘,40∘,40∘,40∘,40∘,40∘,40∘,60∘,60∘,60∘,60∘,60∘,60∘,60∘,80∘,80∘,80∘,80∘,80∘,80∘,80∘ and ϕ=15∘,25∘,35∘,45∘,55∘,65∘,75∘,15∘,25∘,35∘,45∘,55∘,65∘,75∘,15∘,25∘,35∘,45∘,55∘,65∘,75∘,15∘,25∘,35∘,45∘,55∘,65∘,75∘ respectively. From Figure 2, we can obviously see that the 28 signals can be detected uniquely due to the employment of total correlation matrix information including autocorrelation matrix and mutual correlation matrix, where no ambiguous angles arise; however, the methods in [24, 26] fail to detect 28 signals and they can detect 11 signals at most in this scenario.Figure 2
Maximum achievable DOFs without ambiguous angles by using the AMRD-MUSIC algorithm.


### 4.2. Complexity Analysis

In this subsection, we compute the complexity of the classical MD-MUSIC algorithm, multiple-dimensional ESPRIT (MD-ESPRIT) algorithm and the proposed AMRD-MUSIC algorithm for the same array configuration of CCA with the same sensors.For the MD-MUSIC algorithm, the total complexity is composed of three parts, including computing covariance, eigenvalues decomposition, and spectrum peak searching. The complexity of these three parts is presented asT2L, T3 and 180/τ90/τTT−K, respectively, where L represents the number of snapshots T=T1+T2,T1=M1N1J1,T2=M2N2J2 and τ denotes the searching step. So the total computational complexity can be denoted as T2L+T3+180/τ90/τTT−K.And we compute the complexity of the MD-ESPRIT algorithm. For the two subarrays, complexity consists of three parts, computing the covariance asT2L, eigenvalue decomposition as T13+T23, computing pseudoinverse and eigenvalue decomposition after reconstruction as T1−N1J13+T2−N2J23+T1−M1J13+T2−M2J23+T1−N1M13+T2−N2M23. So the total complexity of these operations is denoted as T2L+T13+T23+T1−N1J13+T2−N2J23+T1−M1J13+T2−M2J23M2J23+T1−N1M13+T2−N2M23.In the proposed AMRD-MUSIC algorithm, the complexity can be denoted asT12L+2KT1T2+2T1T22+2N1J1T1T+N1J12T1+T−N1J3+2N2J2T2T+N2J22T2+T23, and the complexity of the other subarray is approximately similar to it.We compare the complexity of the conventional MD-MUSIC algorithm, the MD-ESPRIT algorithm, and the proposed AMRD-MUSIC algorithm with the same CCA configuration. For clarity, we list the computational complexity of these algorithms with CCA in Table2 with L=200,K=2,M1=3,M2=2,N1=2,N2=3,J1=3,J2=2,J1=3,J2=2, and the searching step τ is set to be 0.02.Table 2
Computational complexity of different algorithms.
Algorithms
Complex multiplication
MD-MUSIC
5.4432
×
10
9
MD-ESPRIT
1.9212
×
10
5
AMRD-MUSIC
1.20096
×
10
5We compare the computational complexity versus the number of snapshots and sensors respectively of three algorithms including MD-MUSIC, MD-ESPRIT, and ARMD-MUSIC. Figures3 and 4 show the complexity comparison with the increasing number of snapshots and sensors, respectively. It can be seen that the proposed AMRD-MUSIC algorithm achieves the lowest computational complexity than other algorithms. Figure 5 shows the complexity comparison versus the searching step. We can see that the complexity of MD-MUSIC reduces with the searching step increasing. However, the proposed AMRD-MUSIC algorithm outperforms the other algorithms.Figure 3
Comparison of total complexity versus snapshots.Figure 4
Comparison of total complexity versus the number of elements.Figure 5
Comparison of total complexity versus searching step.


### 4.3. Advantages of the Proposed Algorithm

In summary, the advantages of the proposed algorithm are listed as follows:(1)
The proposed algorithm achieves significant complexity decrease by reducing 2D SPS into 1D one and utilizing polynomial root finding technique instead of 1D SPS.(2)
The proposed algorithm possesses approximately the same DOA estimation performance as the MD-MUSIC algorithm and behaves better in its estimation than the classic MD-ESPRIT algorithm.(3)
The proposed algorithm can obtain full DOFs due to the utilization of the total matrix information.(4)
The proposed algorithm can achieve DOA angles with pairing automatically.


### 4.4. Cramer–Rao Bound

We provide the Cramer–Rao Bound (CRB) [40–42] analytical derivation of DOA estimation for CCA. The signal covariance can be obtained by the following:(47)Rs=1L∑t=1LstsHt,and we define the matrix parameter vector η=θ1,θ2,…,θK,ϕ1,ϕ2,…,ϕKT.Then, we can denote theξ0,ζ0-th element of the Fisher information matrix (FIM) as follows [42]:(48)FIMξ0,ζ0=ς0tr∂p¯∂ηξ0P−1∂p¯∂ηζ0P−1=ς0∂p¯∂ηξ0HPT⊗P−1∂p¯∂ηζ0,where(49)ς0=M1N1J1+M2N2J2,(50)P=APsAH+Iς0,(51)p¯=vecP=A∗⊙AΞ+vecIς0,and we denote ∂p¯/∂η=∂p¯/∂θ1,∂p¯/∂θ2,…,∂p¯/∂θK,∂p¯/∂ϕ1,∂p¯/∂ϕ2,…,∂p¯/∂ϕK; then we have the following:(52)FIM=ς0∂p¯∂ηHPT⊗P−1∂p¯∂η.By utilizing (52), we can calculate the following:(53)∂p¯∂η=Ad1∏,Ad2∏,where(54)Ad1=D1∗⊙A+A∗⊙D1,Ad2=D2∗⊙A+A∗⊙D2,D1=∂au1,v1∂θ1,∂au2,v2∂θ2,…,∂auK,vK∂θK,D2=∂au1,v1∂ϕ1,∂au2,v2∂ϕ2,…,∂auK,vK∂ϕK,auk,vk=ax1uk⊗ay1vk⊗az1wkax2uk⊗ay2vk⊗az2wk,k=1,2,…,K,∏=diagΞ.The matrixPT⊗P is positive definite since the P in (50) is positive definite.By letting(55)F=PT⊗P−1/2AdR0,the CRB can be obtained as follows:(56)CRB=1ς0FHF−1,where Ad=Ad1,Ad2,R0=diagRs,Rs.


## 4.1. DOF Comparison

The proposed AMRD-MUSIC algorithm can obtainM1N1J1+M2N2J2−1 DOFs, which indicates that the proposed algorithm can achieve M1N1J1+M2N2J2−2 signals by employing M1N1J1+M2N2J2−1 antennas.(44)DOFCCAAMRD−MUSIC=M1N1J1+M2N2J2−1.The algorithms in [24, 26] process the two subarrays separately so that they have a great loss in DOFs since they only exploit the autoinformation. The number of achievable DOFs by using the method in [24, 26] can be denoted as follows:(45)DOFCCAalgorithms in24,26=minM1N1J1,M2N2J2,where M1=N1=J1>M2=N2=J2.Thus, the proposed algorithms can detect many more signals than conventional algorithms with the same sensors. From (44) and (45), we can note the relationship of DOFs of different algorithms as follows:(46)DOFCCAalgorithms in24,26<DOFCCAARMD−MUSIC.Figure2 depicts the maximum number of achievable DOFs to illustrate the signals explicitly without ambiguous DOA estimations, where the number of snapshots is set to be L=500 with noise-free, the total number of signals is K=28, and the total number of sensors is T=3×2×3+2×3×2−1=29. The elevation angles and azimuth angles are set to be θ=20∘,20∘,20∘,20∘,20∘,20∘,20∘,40∘,40∘,40∘,40∘,40∘,40∘,40∘,60∘,60∘,60∘,60∘,60∘,60∘,60∘,80∘,80∘,80∘,80∘,80∘,80∘,80∘ and ϕ=15∘,25∘,35∘,45∘,55∘,65∘,75∘,15∘,25∘,35∘,45∘,55∘,65∘,75∘,15∘,25∘,35∘,45∘,55∘,65∘,75∘,15∘,25∘,35∘,45∘,55∘,65∘,75∘ respectively. From Figure 2, we can obviously see that the 28 signals can be detected uniquely due to the employment of total correlation matrix information including autocorrelation matrix and mutual correlation matrix, where no ambiguous angles arise; however, the methods in [24, 26] fail to detect 28 signals and they can detect 11 signals at most in this scenario.Figure 2
Maximum achievable DOFs without ambiguous angles by using the AMRD-MUSIC algorithm.


## 4.2. Complexity Analysis

In this subsection, we compute the complexity of the classical MD-MUSIC algorithm, multiple-dimensional ESPRIT (MD-ESPRIT) algorithm and the proposed AMRD-MUSIC algorithm for the same array configuration of CCA with the same sensors.For the MD-MUSIC algorithm, the total complexity is composed of three parts, including computing covariance, eigenvalues decomposition, and spectrum peak searching. The complexity of these three parts is presented asT2L, T3 and 180/τ90/τTT−K, respectively, where L represents the number of snapshots T=T1+T2,T1=M1N1J1,T2=M2N2J2 and τ denotes the searching step. So the total computational complexity can be denoted as T2L+T3+180/τ90/τTT−K.And we compute the complexity of the MD-ESPRIT algorithm. For the two subarrays, complexity consists of three parts, computing the covariance asT2L, eigenvalue decomposition as T13+T23, computing pseudoinverse and eigenvalue decomposition after reconstruction as T1−N1J13+T2−N2J23+T1−M1J13+T2−M2J23+T1−N1M13+T2−N2M23. So the total complexity of these operations is denoted as T2L+T13+T23+T1−N1J13+T2−N2J23+T1−M1J13+T2−M2J23M2J23+T1−N1M13+T2−N2M23.In the proposed AMRD-MUSIC algorithm, the complexity can be denoted asT12L+2KT1T2+2T1T22+2N1J1T1T+N1J12T1+T−N1J3+2N2J2T2T+N2J22T2+T23, and the complexity of the other subarray is approximately similar to it.We compare the complexity of the conventional MD-MUSIC algorithm, the MD-ESPRIT algorithm, and the proposed AMRD-MUSIC algorithm with the same CCA configuration. For clarity, we list the computational complexity of these algorithms with CCA in Table2 with L=200,K=2,M1=3,M2=2,N1=2,N2=3,J1=3,J2=2,J1=3,J2=2, and the searching step τ is set to be 0.02.Table 2
Computational complexity of different algorithms.
Algorithms
Complex multiplication
MD-MUSIC
5.4432
×
10
9
MD-ESPRIT
1.9212
×
10
5
AMRD-MUSIC
1.20096
×
10
5We compare the computational complexity versus the number of snapshots and sensors respectively of three algorithms including MD-MUSIC, MD-ESPRIT, and ARMD-MUSIC. Figures3 and 4 show the complexity comparison with the increasing number of snapshots and sensors, respectively. It can be seen that the proposed AMRD-MUSIC algorithm achieves the lowest computational complexity than other algorithms. Figure 5 shows the complexity comparison versus the searching step. We can see that the complexity of MD-MUSIC reduces with the searching step increasing. However, the proposed AMRD-MUSIC algorithm outperforms the other algorithms.Figure 3
Comparison of total complexity versus snapshots.Figure 4
Comparison of total complexity versus the number of elements.Figure 5
Comparison of total complexity versus searching step.


## 4.3. Advantages of the Proposed Algorithm

In summary, the advantages of the proposed algorithm are listed as follows:(1)
The proposed algorithm achieves significant complexity decrease by reducing 2D SPS into 1D one and utilizing polynomial root finding technique instead of 1D SPS.(2)
The proposed algorithm possesses approximately the same DOA estimation performance as the MD-MUSIC algorithm and behaves better in its estimation than the classic MD-ESPRIT algorithm.(3)
The proposed algorithm can obtain full DOFs due to the utilization of the total matrix information.(4)
The proposed algorithm can achieve DOA angles with pairing automatically.


## 4.4. Cramer–Rao Bound

We provide the Cramer–Rao Bound (CRB) [40–42] analytical derivation of DOA estimation for CCA. The signal covariance can be obtained by the following:(47)Rs=1L∑t=1LstsHt,and we define the matrix parameter vector η=θ1,θ2,…,θK,ϕ1,ϕ2,…,ϕKT.Then, we can denote theξ0,ζ0-th element of the Fisher information matrix (FIM) as follows [42]:(48)FIMξ0,ζ0=ς0tr∂p¯∂ηξ0P−1∂p¯∂ηζ0P−1=ς0∂p¯∂ηξ0HPT⊗P−1∂p¯∂ηζ0,where(49)ς0=M1N1J1+M2N2J2,(50)P=APsAH+Iς0,(51)p¯=vecP=A∗⊙AΞ+vecIς0,and we denote ∂p¯/∂η=∂p¯/∂θ1,∂p¯/∂θ2,…,∂p¯/∂θK,∂p¯/∂ϕ1,∂p¯/∂ϕ2,…,∂p¯/∂ϕK; then we have the following:(52)FIM=ς0∂p¯∂ηHPT⊗P−1∂p¯∂η.By utilizing (52), we can calculate the following:(53)∂p¯∂η=Ad1∏,Ad2∏,where(54)Ad1=D1∗⊙A+A∗⊙D1,Ad2=D2∗⊙A+A∗⊙D2,D1=∂au1,v1∂θ1,∂au2,v2∂θ2,…,∂auK,vK∂θK,D2=∂au1,v1∂ϕ1,∂au2,v2∂ϕ2,…,∂auK,vK∂ϕK,auk,vk=ax1uk⊗ay1vk⊗az1wkax2uk⊗ay2vk⊗az2wk,k=1,2,…,K,∏=diagΞ.The matrixPT⊗P is positive definite since the P in (50) is positive definite.By letting(55)F=PT⊗P−1/2AdR0,the CRB can be obtained as follows:(56)CRB=1ς0FHF−1,where Ad=Ad1,Ad2,R0=diagRs,Rs.


## 5. Simulation Results

Numerical simulations are provided to evaluate the performance of the proposed array configuration of CCA and algorithm in this part. For the CCA geometry, it consists of two uniform cubic subarrays of3×2×3 and 2×3×2 sensors, where dx1=2λ/2, dy1=3λ/2, dz1=2λ/2, dx2=3λ/2, dy2=2λ/2, and dz2=3λ/2. The total number of these two subarrays is computed as T=3×2×3+2×3×2−1=29. The UCA consists of 3×2×5=30 sensors. Assume K=2 signals impinging on the arrays from θ1,ϕ1=20∘,30∘ and θ2,ϕ2=50∘,70∘. From Figure 6, we can clearly see that the AMRD-MUSIC algorithm is able to detect the DOA targets successfully.Figure 6
Scattering figure with AMRD-MUSIC for CCA.The root mean square error (RMSE), as the performance metric of the DOA estimation, is defined as follows:(57)RMSE=1PK∑p=1P∑k=1Kθk−θ^kp2+ϕk−ϕ^kp2,where K is the number of sources. θ^kp,ϕ^kpk=1,2,…,K are estimates of θk,ϕk of the p-th Monte Carlo trial, respectively. In this paper, we set P=500.


### 5.1. DOA Estimation Performance Comparison of Different Array Configurations with the Same Algorithm

We employ simulations to illustrate the DOA estimation performance of UCA and CCA configurations by using the same algorithm, where.K=2,θ1,ϕ1=15∘,30∘,θ2,ϕ2=20∘,40∘.Figures7 and 8 exhibit the RMSE performance comparison between CCA and UCA by utilizing the same algorithm. It captures explicitly that the proposed AMRD-MUSIC algorithm with CCA possesses better DOA estimation performance because of the extended array aperture and less MC effects. And CRB comparison shows that CCA can obtain the lower bound of theoretical DOA estimation performance compared with UCA geometry.Figure 7
RMSE of different arrays versus snapshots.Figure 8
RMSE of different arrays versus SNR.


### 5.2. DOA Estimation Performance Comparison of Different Algorithms with the Same Array

In this subsection, we compare the RMSE performance of different algorithms with the same array of CCA, including the proposed AMRD-MUSIC algorithm, conventional MD-MUSIC algorithm, and MD-ESPRIT algorithm.Figures9 and 10 present the 2D angle estimation comparisons versus SNR and L, respectively. The two figures clearly show that the proposed AMRD-MUSIC algorithm has approximately the same DOA estimation performance as the MD-MUSIC algorithm and outperforms the classic MD-ESPRIT algorithm owing to the utilization of total information including autocorrelation matrix and mutual correlation matrix, while the proposed algorithm has a lower computational complexity as shown in Figure 2.Figure 9
RMSE comparison of different algorithms versus SNR.Figure 10
RMSE comparison of different algorithms versus snapshots.


### 5.3. Comparison of DOA Estimation Performance versus the Number of Sensors

In this subsection, we employ the proposed approach to a CCA geometry which consists of two subarrays including3×N×3 and 2×3×2 sensors, where N=2,4,5. The RMSE performance versus SNR and snapshots are displayed in Figures 11 and 12, respectively. It is obviously shown that the RMSE estimation performance of the proposed AMRD-MUSIC algorithm is improving with the increase of SNR or snapshots due to the diversity gain.Figure 11
RMSE performance of different number of sensors versus SNR.Figure 12
RMSE performance of different number of sensors versus snapshots.


## 5.1. DOA Estimation Performance Comparison of Different Array Configurations with the Same Algorithm

We employ simulations to illustrate the DOA estimation performance of UCA and CCA configurations by using the same algorithm, where.K=2,θ1,ϕ1=15∘,30∘,θ2,ϕ2=20∘,40∘.Figures7 and 8 exhibit the RMSE performance comparison between CCA and UCA by utilizing the same algorithm. It captures explicitly that the proposed AMRD-MUSIC algorithm with CCA possesses better DOA estimation performance because of the extended array aperture and less MC effects. And CRB comparison shows that CCA can obtain the lower bound of theoretical DOA estimation performance compared with UCA geometry.Figure 7
RMSE of different arrays versus snapshots.Figure 8
RMSE of different arrays versus SNR.


## 5.2. DOA Estimation Performance Comparison of Different Algorithms with the Same Array

In this subsection, we compare the RMSE performance of different algorithms with the same array of CCA, including the proposed AMRD-MUSIC algorithm, conventional MD-MUSIC algorithm, and MD-ESPRIT algorithm.Figures9 and 10 present the 2D angle estimation comparisons versus SNR and L, respectively. The two figures clearly show that the proposed AMRD-MUSIC algorithm has approximately the same DOA estimation performance as the MD-MUSIC algorithm and outperforms the classic MD-ESPRIT algorithm owing to the utilization of total information including autocorrelation matrix and mutual correlation matrix, while the proposed algorithm has a lower computational complexity as shown in Figure 2.Figure 9
RMSE comparison of different algorithms versus SNR.Figure 10
RMSE comparison of different algorithms versus snapshots.


## 5.3. Comparison of DOA Estimation Performance versus the Number of Sensors

In this subsection, we employ the proposed approach to a CCA geometry which consists of two subarrays including3×N×3 and 2×3×2 sensors, where N=2,4,5. The RMSE performance versus SNR and snapshots are displayed in Figures 11 and 12, respectively. It is obviously shown that the RMSE estimation performance of the proposed AMRD-MUSIC algorithm is improving with the increase of SNR or snapshots due to the diversity gain.Figure 11
RMSE performance of different number of sensors versus SNR.Figure 12
RMSE performance of different number of sensors versus snapshots.


## 6. Conclusions

In this paper, a CCA geometry in massive MIMO systems is constructed, which consists of two uniform cubic subarrays and obtains the extension of interelement spacing by selecting three pairs of coprime integers. Meanwhile, CRB for CCA is derived and the proposed CCA geometry is verified to outperform the conventional UCA geometry in 2D DOA estimation performance in massive MIMO systems. In addition, we propose a computationally efficient 2D DOA estimation algorithm via array mapping and reduced dimension for CCA. In the proposed algorithm, we utilize the array mapping and map the nonuniform array into two uniform arrays and operate the reduced dimension process on the two uniform arrays to transform the 2D SPS problem into 1D one. In addition, the computational complexity gets further reduced by converting the 1D SPS problem into a polynomial root finding problem. Numerical simulation results verify the effectiveness of the proposed algorithm that can reduce the computational complexity cost with no degradation of DOA estimation performance.



---

*Source: 2686257-2020-01-13.xml*