# Parallel Factor-Based Model for Two-Dimensional Direction Estimation

Authors: Nizar Tayem; Khaqan Majeed; Ahmed A. Hussain
Journal: International Journal of Antennas and Propagation
 (2017)

---

## Abstract

Two-dimensional (2D) Direction-of-Arrivals (DOA) estimation for elevation and azimuth angles assuming noncoherent, mixture of coherent and noncoherent, and coherent sources using extended three parallel uniform linear arrays (ULAs) is proposed. Most of the existing schemes have drawbacks in estimating 2D DOA for multiple narrowband incident sources as follows: use of large number of snapshots, estimation failure problem for elevation and azimuth angles in the range of typical mobile communication, and estimation of coherent sources. Moreover, the DOA estimation for multiple sources requires complex pair-matching methods. The algorithm proposed in this paper is based on first-order data matrix to overcome these problems. The main contributions of the proposed method are as follows: (1) it avoids estimation failure problem using a new antenna configuration and estimates elevation and azimuth angles for coherent sources; (2) it reduces the estimation complexity by constructing Toeplitz data matrices, which are based on a single or few snapshots; (3) it derives parallel factor (PARAFAC) model to avoid pair-matching problems between multiple sources. Simulation results demonstrate the effectiveness of the proposed algorithm.


---

## Body

## 1. Introduction

Antenna arrays have long been in use in the field of wireless communications in both civilian and military applications. A fundamental problem in array signal processing is the Direction of Arrival (DOA) estimation in two dimensions (i.e., elevation and azimuth) [1–4] for source localization of multiple sources impinging on this array of antennas. Several methods were proposed in the literature for DOA estimation. These methods vary based on estimation accuracy, the geometry of the antenna array, the complexity of the underlying algorithm, computational cost, and so forth.In the conventional methods of DOA estimation, the observation space is decomposed into a signal subspace and a noise subspace. Two classic subspace techniques, multiple signal classification (MUSIC [5]) and estimation of signal parameters via rotational invariance technique (ESPRIT [6]), are widely used. These and other similar techniques are computationally complex as they employ either eigenvalue decomposition (EVD) or the singular-value decomposition (SVD) of the data matrix. These techniques also suffer from estimation accuracy and fail in the presence of highly correlated and coherent signals. A spatial smoothing (SS) technique was introduced in [7, 8] as a preprocessing step to improve the system performance and estimation accuracy. However, this results in the increase in computational complexity of the algorithm.Improved techniques which are simpler and less complex were reported in the literature [9–13]. These techniques do not rely on either EVD or SVD. However, some of these methods [10–12] require several identical subarrays configuration (and consequently suffer heavy losses of the array aperture) and encounter estimation failure problem when elevation angles are between 70° and 90° (typical mobile communications elevation angle range). Some cumulant-based methods [14, 15] avoid aperture loss problem. But these methods are computationally intensive, require parameter pairing, and also do not solve the estimation failure problem. A trilinear decomposition-based [16] blind 2D DOA estimation algorithm for an L-shaped array is described in [17]. One drawback of this method is that it requires large number of snapshots. The method in [18] works only for noncoherent sources and also requires large number of snapshots.The algorithm in [19] uses a fourth-order cumulants-based method and performs EVD of the Toeplitz matrices constructed from cumulant elements of received signals to obtain 2D angle parameters without the need for parameter pairing. This method, however, suffers from estimation accuracy at low SNR and smaller number of snapshots. A recent work reported in [20] proposes an oblique projection based approach [21, 22], which is an extended orthogonal projection of the measurement onto a low-rank subspace along a nonorthogonal subspace, by applying some cross-correlation between the received array data. The oblique projection is used to estimate the elevation and azimuth angles by isolating the coherent signals from the noncoherent ones, thereby alleviating the effect of additive noise, while avoiding the computationally intensive EVD and parameter pairing problem. Another recent work reported in [23] investigates the problem of tracking the 2D DOA of multiple moving targets by associating the estimated azimuth and elevation angles of different targets at two successive time instants. This method does not require pair-matching of the estimated azimuth and elevation angles and avoids the computationally expensive EVD. However, it works only for noncoherent signals. Another 2-destination algorithm of interest was proposed in [24]. In this algorithm, the elevation angle is first estimated based on the polynomial root methods (such as fast Root-MUSIC or ESPRIT) using a 1D uniform linear subarray along the z-axis. Next, it obtains the azimuth angle estimate using the elevation angle estimated earlier and a 2D uniform linear subarray along the x-axis based on a subspace method (such as 2D MUSIC). While this algorithm does not require pair-matching, it has relatively high estimation errors at low SNR and requires a large number of snapshots.Array geometries also have a significant impact on 2D DOA estimation accuracy. Several array geometries exist in the literature, such as uniform linear array (ULA), uniform circular array (UCA), rectangular array, parallel array [25], and L-shaped array [26, 27]. The L-shaped array [26, 27] is reported to have higher accuracy compared with other geometries. However, since the L-shaped array consists of two orthogonal linear subarrays, two electric angles are separately estimated from each subarray. Failure to perform pair-matching of these angles will result in incorrect 2D DOA estimation of elevation and azimuth angles, and hence severe performance degradation occurs. The work reported in [28] solves the problem of pair-matching in 2D DOA estimation using L-shaped arrays. Another type of array geometry called cylindrical conformal array is proposed in [29, 30]. This type of array, however, suffers from polarization diversity of element patterns which results in DOA estimation difficulty. This is due to the coupling between the angle information and the polarization parameter which are incorporated in the snapshot data model.The problem of parameter pairing (or pair-matching) for the L-shaped array in [18] and for the conformal array in [31] is successfully alleviated using a parallel factor (PARAFAC) model [32]. PARAFAC, first introduced in psychometrics, is a method of multiple data analysis and is used in many fields such as statistics, arithmetic complexity, and chemometrics. The PARAFAC model [32] is widely used for low-rank decomposition of three-way (TWA) and higher order arrays.The method in [33] employs cross-correlation information of the received signals for constructing a data matrix and uses linear operations which reduce the computational complexity significantly. However, it has several drawbacks compared with the 2D DOA estimation method proposed in this paper. The work in [33] does not solve the pair-matching problem nor does it work for coherent sources. It also requires a large number of snapshots for estimation accuracy.In this paper, we propose a novel 2D DOA estimation algorithm that employs a new three parallel uniform antenna arrays. The proposed method uses PARAFAC model to avoid pair-matching problem. Compared with existing methods [15, 17], the proposed method works for both coherent and noncoherent sources, has lower computational complexity, requires very few snapshots (as low as 1), and has no estimation failure problem especially in typical mobile communication range. The main contributions of the proposed method are as follows: (1) to avoid estimation failure problem using a new antenna configuration and estimate elevation and azimuth angles for coherent sources; (2) to reduce the estimation complexity by constructing Toeplitz data matrices, which are based on a single or few snapshots; and (3) to derive parallel factor (PARAFAC) model to avoid pair-matching problems between multiple sources. Simulation results demonstrate the effectiveness of the proposed algorithm. These advantages make the proposed algorithm a realistic candidate for real-time hardware implementation for high speed wireless communications applications. Simulation results are shown against the Cramer-Rao Bound as the reference.The remaining portion of the paper is organized as follows: The system model and proposed algorithm are presented in Section2. Simulation results and complexity analysis of the proposed method are described in Section 3. Finally, conclusions are drawn in Section 4.


## 2. System Model and Proposed Algorithm

The proposed methodology considers three extended parallel-shaped arrays, where Toeplitz matrices are constructed from the received signals on these arrays to estimate the elevation and azimuth angles. The proposed algorithm assumes noncoherent, coherent, or a mixture of noncoherent and coherent sources to estimate the angles. The following notations are used to represent variables throughout the paper. The operations·H, ·T, ·†, ·F, R·, ⊙, ⊗, and ⊚ represent conjugate transpose, transpose, pseudoinverse, Frobenius norm, real part of a complex number, Khatri-Rao product, Kronecker product, and Hadamard product, respectively. The scalar is denoted by υ, constant byV, vector byv, matrix byV, and ijth member of a matrixV by υi,j. A three-way array (TWA) is denoted byV and its ijkth member is given by υi,j,k. The operations diag(v) and diag-1(V) denote conversion of the vectorvinto a diagonal matrix and the diagonal matrixV into a vector, respectively. The following subsections describe the proposed extended parallel-shaped array model and formulation of the problem to estimate the elevation and azimuth angles in the presence of noncoherent and/or coherent sources.


### 2.1. Proposed Extended Array Geometry and Signal Model

The proposed extended parallel-shaped array geometry is shown in Figure1. The figure shows three parallel arrays placed on y-axis, y-z plane, and x-y plane. The middle array contains one extra element as compared to the other arrays. There are 2N antenna elements in the middle array. The arrays above and below the middle array contain 2N-1 elements each. These arrays are partitioned into four subarrays such that each subarray contains 2N-1 elements. These are denoted as X, Y, Z, and W as shown in the figure. The distance between adjacent elements on a single array is d meters. The array on the x-axis is separated by d meters from the arrays in the x-z and x-y planes.Figure 1
Geometry of the proposed extended parallel-shaped array.ConsiderK narrowband noncoherent and/or coherent sources in far field of the extended parallel-shaped array. Let θk and ϕk denote the elevation and azimuth angles, respectively, for the kth source. The 2N-1×1 signal vectors received at the X, Y, Z, and W subarrays at the tth snapshot are given, respectively, as(1)xt=x-N-1t⋯x0t⋯xN-1tT,yt=y-N-2t⋯y0t⋯yNtT,zt=z-N-1t⋯z0t⋯zN-1tT,wt=w-N-1t⋯w0t⋯wN-1tT,where t=1,2,…,L and superscript T denotes the transpose. The 2N-1×1-dimensional signal vectors received at the antenna arrays can be represented as(2)xt=Aθ,ϕst+nxt,(3)yt=Aθ,ϕΦyθ,ϕst+nyt,(4)zt=Aθ,ϕΦzθst+nzt,(5)wt=Aθ,ϕΦwθ,ϕst+nwt,where the matrix Aθ,ϕ is defined as (6)Aθ,ϕ=u1-N-1u2-N-1⋯uK-N-1⋮⋮⋱⋮u1-1u2-1⋯uK-111⋯1u11u21⋯uK1⋮⋮⋱⋮u1N-1u2N-1⋯uKN-1.The dimensions of Aθ,ϕ are 2N-1×K. The entry uk in (6) is defined as (7)uk=exp-j2πdsin⁡θksin⁡ϕkλ,where θk and ϕk are the elevation and azimuth angles, respectively, defined earlier for the kth source. The signal vector st in (2) to (5) is defined as(8)st=s1ts2⋯sKtT.The matrix Φy in (3) is diagonal and contains information about the elevation and azimuth angles. Its dimensions are K×K and is defined as(9)Φy=diagq1,q2,…,qK,where qK is given as(10)qK=exp-j2πdsin⁡θksin⁡ϕkλ.Similarly, Φz and Φw in (4) and (5), respectively, are defined as(11)Φz=diagr1,r2,…,rK,where rK is given as(12)rK=exp-j2πdcos⁡θkλ,Φw=diagv1,v2,…,vK,where vK is given as(13)vK=exp-j2πdsin⁡θkcos⁡ϕkλ.The2N-1×1-dimensional vectors nxt, nyt, nzt, and nwt appearing in (2), (3), (4), and (5), respectively, are noise vectors. These vectors are assumed independent of each other and each element of the noise vector (nxt,nyt,nzt,ornwt) is also independent of other elements in that vector. The noise is assumed to be Additive White Gaussian (AWG) with mean 0 and variance σ2. The noise vectors are also independent of the source signals. The goal now is to estimate the diagonal matrices Φy, Φz, and Φw, which contain information about the elevation and azimuth angles. The PARAFAC model is employed to estimate elevation and azimuth angles for multiple sources and solve the pair-matching problem, which occurs in the presence of two or more sources. The following subsection introduces the basics of the PARAFAC model [32].


### 2.2. PARAFAC Model

Consider a TWAD_ with loading matricesA,B, andC. The array D_ is defined as(14)D_=A⊗BCT.The P-component trilinear decomposition of the array D_ is defined as(15)di,j,k=∑p=1Pai,pbj,pck,p,where i=1,2,…,Id,j=1,2,…,Jd,k=1,2,…,Kd, and di,j,k is the i,j,kth element of the array D_. The element ai,P is the i,Pth element of the Id×P-dimensional matrixA. Similarly, bi,P and ck,P are j,Pth and k,Pth elements of the Jd×P-dimensional and Kd×P-dimensional matricesB andC, respectively. The matrices A=a1a2⋯aP, B=b1b2⋯bP, and C=c1c2⋯cP are called the mode or loading matrices for a given PARAFAC model. The PARAFAC analysis of D_ in (14) is represented by (15).In the proposed method, PARAFAC model is used to estimate the elevation and azimuth angles correctly for each signal source without the problem of pair-matching.


### 2.3. Problem Formulation

In the proposed method, the Toeplitz matrices from multiple signals received by the extended parallel-shaped antenna array are first constructed. The advantage of constructing Toeplitz matrix is its capability to estimate the angles regardless of the signals being noncoherent and/or coherent. The Toeplitz matrices are constructed using the signals received at the subarraysX,Y,Z, and W from a single snapshot. As inferred from the array geometry in Figure 1, the array elements are indexed from -N-1 to N for the array in the middle and -N-1 to N-1 for the other two arrays. Using the index 0 element as a reference, we define the Toeplitz matrix Rxt for the subarray X as(16)Rxt=x0tx-1t⋯x-N-1tx1tx0t⋯x-N-2t⋮⋮⋱⋮xN-1txN-2t⋯x0t.The reference element for subarray Y is at index 1 since it is defined as subarray X shifted right by one element. Similarly, for subarray Y, the Toeplitz matrix is defined as(17)Ryt=y1ty0t⋯y-N-2ty2ty1t⋯y-N-3t⋮⋮⋱⋮yNtyN-1t⋯y1t.The reference antenna elements for Z and W subarrays both are at index 0 and the Toeplitz matrices for subarrays Z and W are defined, respectively, as(18)Rzt=z0tz-1t⋯z-N-1tz1tz0t⋯z-N-2t⋮⋮⋱⋮zN-1tzN-2t⋯z0t,(19)Rwt=w0tw-1t⋯w-N-1tw1tw0t⋯w-N-2t⋮⋮⋱⋮wN-1twN-2t⋯w0t.The dimensions of the matrices Rxt, Ryt, Rzt, and Rwt are N×N each. These matrices are obtained at the tth snapshot of the received signals.The Toeplitz matrixRxt in (16) can be decomposed as follows:(20)Rxt=Eθ,ϕStEHθ,ϕ+Nxt,where Nxt is a noise matrix with dimensions N×N and Eθ,ϕ and St are defined as (21)Eθ,ϕ=11⋯1u1u2⋯uKu12u22⋯uK2⋮⋮⋱⋮u1N-1u2N-1⋯uKN-1,(22)St=s1t0⋯00s2t⋯0⋮⋮⋱⋮00⋯sKt,where uk is defined in (7) and St is a diagonal matrix containing the signals from K sources. The dimensions of Eθ,ϕ and St are N×K and K×K, respectively. Since the number of antenna elements N is usually greater than the number of sources K so the rank of Eθ,ϕ is K, which is full column rank. The diagonal matrix St is of full rank whether the K sources are coherent or not. Consider the case of noncoherent sources with diagonal signal matrix as shown in (22). The dot product of any two columns is zero; that is, s1t0⋯0T,0s2t⋯0T=0. Now, consider two coherent sources such that s2t=αs1t. The dot product corresponding to these two columns is again zero; that is, s1t0⋯0T,0αs1t⋯0T=0. This indicates that the matrix St is of full rank, which is equal to K, the number of sources. So the rank is actually the minimum of N and K, which indicates that the rank is equal to the number of sources K. As inferred from (20), the rank of Toeplitz matrix Rxt is K also. Consider the matrices U and V with dimensions N×K and K×K, respectively, where the matrix U is of full column rank. The theorem states that rUVUH=rV, where rV denotes the rank of the matrix V, where Nyt, Nzt, and Nwt are the noise matrices for the subarrays Y, Z, and W, respectively. Therefore, the proposed method can estimate elevation and azimuth angles for noncoherent and/or coherent sources.Similarly, the Toeplitz matricesRyt, Rzt, and Rwt can be decomposed, respectively, as(23)Ryt=Eθ,ϕΦyStEHθ,ϕ+Nyt,Rzt=Eθ,ϕΦzStEHθ,ϕ+Nzt,Rwt=Eθ,ϕΦwStEHθ,ϕ+Nwt.The Toeplitz matrices defined above are rearranged to form anN×N×4 TWA F_t. This is represented as(24)F_t=F:,:,1;tF:,:,2;tF:,:,3;tF:,:,4;t=RxtRytRztRwt+NxtNytNztNwt=Eθ,ϕStEHθ,ϕEθ,ϕΦyStEHθ,ϕEθ,ϕΦzStEHθ,ϕEθ,ϕΦwStEHθ,ϕ+NxtNytNztNwt,where this TWA is in the same form as defined in (14) and the noise matrices Nxt, Nyt, Nzt, and Nwt are as defined earlier.The TWAF_ is rearranged into three slices of 2D matrices F1, F2, and F3 as follows:(25)F1=A⊙BCT+N1,F2=A⊙BCT+N2,F3=A⊙BCT+N3,where the matrices N1, N2, and N3 are the rearranged noise matrices, A=Eθ,ϕ, and B=EHθ,ϕ. The matrix C contains information about the elevation and azimuth angles. The structure of C is shown in the following:(26)C=c1Tc2Tc3Tc4T=diag-1Sdiag-1ΦySdiag-1ΦzSdiag-1ΦwS.The cost functions for the matricesA,B, andC to be minimized are defined as follows:(27)∫1A,B,C;F1=F1=A⊙BCTF2,∫2A,B,C;F2=F2=B⊙CATF2,∫3A,B,C;F3=F3=C⊙ABTF2.The alternating least squares method is applied to the cost function∫1A,B,C;F1 to find the matrix C, which contains the information about elevation and azimuth angles. This solution is given as follows: (28)C^=argminCF1-A^⊙B^C^TF2=A^⊙B^†F1T,where A^ and B^ are obtained, similarly, as follows:(29)A^=B^⊙C^†F2T,B^=C^⊙A^†F3T.The matricesA^, B^, and C^ contain the result of one iteration. Note that only a few iterations are required for the algorithm to converge.The COMFAC [34] algorithm is used to fit PARAFAC model defined in (15) to the TWA F_t. The COMFAC MATLAB function has the following form: A,B,C,∙,i=comfacF,P,∙,∙,∙,∙, where F_ (index t is removed for simplification) and P are the decomposing TWA and corresponding factor number, respectively. The factor number P represents the number of sources, which are K. The outputsA,B, andC represent the identification results (matrices), while i represents the iteration number required for the low-rank decomposition. The “∙” represents other options.After obtaining the identification matricesA^, B^, and C^ by performing PARAFAC analysis on F_, the diagonal matrix Φ^y is estimated from the identification matrix C^ as shown below(30)φ^yk=c^2kc^1k,where φ^yk is the estimated kth entry of the diagonal matrix Φ^y and c^1k and c^2k represent the kth entries of the row vectors c^1T and c^2T, respectively. Similarly, the kth diagonal entries of Φ^z and Φ^w are estimated as shown in the following:(31)φ^zk=c^3kc^1k,(32)φ^wk=c^4kc^1k,where c^3k and c^4k denote the kth entries of the row vectors c^3T and c^4T, respectively.The azimuth angleϕ^k is then estimated for the kth source using the following:(33)ϕ^k=tan-1∠φ^yk∠φ^wkand elevation angle θ^k for the kth source using the following:(34)θ^k=tan-1∠φ^wk∠φ^zk×cos⁡ϕ^k.The estimated azimuth and elevation angles in (33) and (34), respectively, are for the kth source. So, for K sources, the following pairs are obtained: θ1,ϕ1,θ2,ϕ2,…,θK,ϕK. The estimated identification matrices C^, A^, and B^ obtained earlier have the same column permutation matrix [31]. The pair-matching problem is automatically avoided since the kth column of A^ corresponds to the kth column of B^. Moreover, the elevation and azimuth angles estimated in (34) and (33), respectively, do not result in failure in practical mobile elevation angle range (70° to 90°).The estimation performance is improved by using more than one snapshot of the received signals. For this purpose, consider a TWAG_, which contains L snapshots of the received signals. The TWA G_ is constructed such that the L snapshots of the Toeplitz matrices defined earlier are concatenated together as shown in the following:(35)G_=F_1F_2⋯F_L=Rx1Rx2⋯RxLRy1Ry2⋯RyLRz1Rz2⋯RzLRw1Rw2⋯RwL.The PARAFAC analysis is then performed on TWA G followed by the estimation of azimuth and elevation angles in (33) and (34), respectively. The performance comparison using one and multiple snapshots (up to 25) is investigated in the simulation results section.


### 2.4. Summary of the Proposed Algorithm

As described earlier, the proposed algorithm considers the extended parallel-shaped array shown in Figure1. The problem of noncoherent and/or coherent sources is resolved by constructing the Toeplitz matrices and the problem of pair-matching between two or more sources is avoided using the PARAFAC model. The following steps summarize the proposed algorithm to estimate DOA for noncoherent and/or coherent sources without the problem of pair-matching.Step 1.
Construct the Toeplitz matricesRxt, Ryt, Rzt, and Rwt shown in (16), (17), (18), and (19), respectively, from single snapshot of the received signals.Step 2.
The Toeplitz matrices used to form theN×N×4 TWA F_t given in (24) for one snapshot or G_ in (35) for L snapshots are rearranged according to the PARAFAC model shown in (14).Step 3.
Apply alternating least squares (ALS) given in (28) and (29) to the TWA F_t to obtain the output matrix C, which has the same structure as shown in (26).Step 4.
Estimate the diagonal matricesΦ^y, Φ^z, and Φ^w using (30), (31), and (32), respectively. These matrices contain the elevation and azimuth angle information of the K sources.Step 5.
Estimate the azimuth angleϕ^k and elevation angle θ^k using (33) and (34), respectively.


### 2.5. Complexity and Cramer-Rao Bound Analysis of the Proposed Method

The computational complexity of the proposed method is compared with those of 2D DOA using cumulant-based method [15], novel 2D DOA estimation algorithm with L-shaped array [17], universal 2D DOA estimation [24], and joint elevation azimuth and elevation angle estimation [28]. For N antenna elements and K sources and considering only the major processing operations, the computational complexity of the proposed method is in the order of O4N+1Lp+n3K3+12N+12K, where Lp denotes the number of snapshots of the received signals used in the proposed method. The complexity of the novel L-shaped method in [17] is O4N-12Ln+n3K3+12N-12K, where Ln denotes the number of snapshots of the received signals. The cumulant-based method in [15] has complexity in the order of O212N+12Lc+n3K3+122N+12K, where Lc denotes the number of snapshots of the received signals. The complexity of the universal 2D DOA estimation in [24] is OLuN2+2N3+N2N-K+N+8N2Lu  +  OLuN2+2N3+V/ΔN-KNLu3, where Lu denotes the number of snapshots of the received signals, V is 180 grid number, and Δ is a number between 0 and 1. The first part represents the complexity to estimate elevation angle using Root-MUSIC algorithm and second part to estimate azimuth angle using MUSIC algorithm. The joint L-shaped method in [28] has complexity in the order of O2N+12Lj2+4/32N+13+20001×10N3+20001N-1N+1 (extracted from [28]), where Lj denotes the number of snapshots of the received signals. This implies that the methods in [15, 17, 24, 28] require a lot more computations compared to the proposed method. Moreover, the received signals’ snapshots used in the proposed method are far less than those used in the above methods; that is, Lp≪Ln,Lc,Lu, and Lj. This is described in Simulation Results.Cramer-Rao Bound (CRB) for 2-dimensional DOA estimation is derived in a similar manner as shown in [19, 35–37]. The received signals in (2), (3), (4), and (5) can be concatenated as shown in the following:(36)xtytztwt=Aθ,ϕAθ,ϕΦyθ,ϕAθ,ϕΦzθAθ,ϕΦwθ,ϕst+nxtnytnztnwt=Mst+Nt.Using the signal model above, CRB can be expressed as(37)CRB=σ22LRVHΠM⊥V⊙P^T-1,where V=[∂m1/∂θ1,∂m1/∂θ1,…, ∂mK/∂θK, ∂m1/∂θ1, ∂m2/∂θ2,…, ∂mK/∂θK], mi is the ith column of M, and(38)ΠM⊥=I8N-4-MMHM-1MH.The matrix P^ in (37) is given as(39)P^=P^sP^sP^sP^s,where P^s is as follows:(40)P^s=1L∑l=1LstsHt.The matrix P^s is not a diagonal matrix in the presence of coherent sources.


## 2.1. Proposed Extended Array Geometry and Signal Model

The proposed extended parallel-shaped array geometry is shown in Figure1. The figure shows three parallel arrays placed on y-axis, y-z plane, and x-y plane. The middle array contains one extra element as compared to the other arrays. There are 2N antenna elements in the middle array. The arrays above and below the middle array contain 2N-1 elements each. These arrays are partitioned into four subarrays such that each subarray contains 2N-1 elements. These are denoted as X, Y, Z, and W as shown in the figure. The distance between adjacent elements on a single array is d meters. The array on the x-axis is separated by d meters from the arrays in the x-z and x-y planes.Figure 1
Geometry of the proposed extended parallel-shaped array.ConsiderK narrowband noncoherent and/or coherent sources in far field of the extended parallel-shaped array. Let θk and ϕk denote the elevation and azimuth angles, respectively, for the kth source. The 2N-1×1 signal vectors received at the X, Y, Z, and W subarrays at the tth snapshot are given, respectively, as(1)xt=x-N-1t⋯x0t⋯xN-1tT,yt=y-N-2t⋯y0t⋯yNtT,zt=z-N-1t⋯z0t⋯zN-1tT,wt=w-N-1t⋯w0t⋯wN-1tT,where t=1,2,…,L and superscript T denotes the transpose. The 2N-1×1-dimensional signal vectors received at the antenna arrays can be represented as(2)xt=Aθ,ϕst+nxt,(3)yt=Aθ,ϕΦyθ,ϕst+nyt,(4)zt=Aθ,ϕΦzθst+nzt,(5)wt=Aθ,ϕΦwθ,ϕst+nwt,where the matrix Aθ,ϕ is defined as (6)Aθ,ϕ=u1-N-1u2-N-1⋯uK-N-1⋮⋮⋱⋮u1-1u2-1⋯uK-111⋯1u11u21⋯uK1⋮⋮⋱⋮u1N-1u2N-1⋯uKN-1.The dimensions of Aθ,ϕ are 2N-1×K. The entry uk in (6) is defined as (7)uk=exp-j2πdsin⁡θksin⁡ϕkλ,where θk and ϕk are the elevation and azimuth angles, respectively, defined earlier for the kth source. The signal vector st in (2) to (5) is defined as(8)st=s1ts2⋯sKtT.The matrix Φy in (3) is diagonal and contains information about the elevation and azimuth angles. Its dimensions are K×K and is defined as(9)Φy=diagq1,q2,…,qK,where qK is given as(10)qK=exp-j2πdsin⁡θksin⁡ϕkλ.Similarly, Φz and Φw in (4) and (5), respectively, are defined as(11)Φz=diagr1,r2,…,rK,where rK is given as(12)rK=exp-j2πdcos⁡θkλ,Φw=diagv1,v2,…,vK,where vK is given as(13)vK=exp-j2πdsin⁡θkcos⁡ϕkλ.The2N-1×1-dimensional vectors nxt, nyt, nzt, and nwt appearing in (2), (3), (4), and (5), respectively, are noise vectors. These vectors are assumed independent of each other and each element of the noise vector (nxt,nyt,nzt,ornwt) is also independent of other elements in that vector. The noise is assumed to be Additive White Gaussian (AWG) with mean 0 and variance σ2. The noise vectors are also independent of the source signals. The goal now is to estimate the diagonal matrices Φy, Φz, and Φw, which contain information about the elevation and azimuth angles. The PARAFAC model is employed to estimate elevation and azimuth angles for multiple sources and solve the pair-matching problem, which occurs in the presence of two or more sources. The following subsection introduces the basics of the PARAFAC model [32].


## 2.2. PARAFAC Model

Consider a TWAD_ with loading matricesA,B, andC. The array D_ is defined as(14)D_=A⊗BCT.The P-component trilinear decomposition of the array D_ is defined as(15)di,j,k=∑p=1Pai,pbj,pck,p,where i=1,2,…,Id,j=1,2,…,Jd,k=1,2,…,Kd, and di,j,k is the i,j,kth element of the array D_. The element ai,P is the i,Pth element of the Id×P-dimensional matrixA. Similarly, bi,P and ck,P are j,Pth and k,Pth elements of the Jd×P-dimensional and Kd×P-dimensional matricesB andC, respectively. The matrices A=a1a2⋯aP, B=b1b2⋯bP, and C=c1c2⋯cP are called the mode or loading matrices for a given PARAFAC model. The PARAFAC analysis of D_ in (14) is represented by (15).In the proposed method, PARAFAC model is used to estimate the elevation and azimuth angles correctly for each signal source without the problem of pair-matching.


## 2.3. Problem Formulation

In the proposed method, the Toeplitz matrices from multiple signals received by the extended parallel-shaped antenna array are first constructed. The advantage of constructing Toeplitz matrix is its capability to estimate the angles regardless of the signals being noncoherent and/or coherent. The Toeplitz matrices are constructed using the signals received at the subarraysX,Y,Z, and W from a single snapshot. As inferred from the array geometry in Figure 1, the array elements are indexed from -N-1 to N for the array in the middle and -N-1 to N-1 for the other two arrays. Using the index 0 element as a reference, we define the Toeplitz matrix Rxt for the subarray X as(16)Rxt=x0tx-1t⋯x-N-1tx1tx0t⋯x-N-2t⋮⋮⋱⋮xN-1txN-2t⋯x0t.The reference element for subarray Y is at index 1 since it is defined as subarray X shifted right by one element. Similarly, for subarray Y, the Toeplitz matrix is defined as(17)Ryt=y1ty0t⋯y-N-2ty2ty1t⋯y-N-3t⋮⋮⋱⋮yNtyN-1t⋯y1t.The reference antenna elements for Z and W subarrays both are at index 0 and the Toeplitz matrices for subarrays Z and W are defined, respectively, as(18)Rzt=z0tz-1t⋯z-N-1tz1tz0t⋯z-N-2t⋮⋮⋱⋮zN-1tzN-2t⋯z0t,(19)Rwt=w0tw-1t⋯w-N-1tw1tw0t⋯w-N-2t⋮⋮⋱⋮wN-1twN-2t⋯w0t.The dimensions of the matrices Rxt, Ryt, Rzt, and Rwt are N×N each. These matrices are obtained at the tth snapshot of the received signals.The Toeplitz matrixRxt in (16) can be decomposed as follows:(20)Rxt=Eθ,ϕStEHθ,ϕ+Nxt,where Nxt is a noise matrix with dimensions N×N and Eθ,ϕ and St are defined as (21)Eθ,ϕ=11⋯1u1u2⋯uKu12u22⋯uK2⋮⋮⋱⋮u1N-1u2N-1⋯uKN-1,(22)St=s1t0⋯00s2t⋯0⋮⋮⋱⋮00⋯sKt,where uk is defined in (7) and St is a diagonal matrix containing the signals from K sources. The dimensions of Eθ,ϕ and St are N×K and K×K, respectively. Since the number of antenna elements N is usually greater than the number of sources K so the rank of Eθ,ϕ is K, which is full column rank. The diagonal matrix St is of full rank whether the K sources are coherent or not. Consider the case of noncoherent sources with diagonal signal matrix as shown in (22). The dot product of any two columns is zero; that is, s1t0⋯0T,0s2t⋯0T=0. Now, consider two coherent sources such that s2t=αs1t. The dot product corresponding to these two columns is again zero; that is, s1t0⋯0T,0αs1t⋯0T=0. This indicates that the matrix St is of full rank, which is equal to K, the number of sources. So the rank is actually the minimum of N and K, which indicates that the rank is equal to the number of sources K. As inferred from (20), the rank of Toeplitz matrix Rxt is K also. Consider the matrices U and V with dimensions N×K and K×K, respectively, where the matrix U is of full column rank. The theorem states that rUVUH=rV, where rV denotes the rank of the matrix V, where Nyt, Nzt, and Nwt are the noise matrices for the subarrays Y, Z, and W, respectively. Therefore, the proposed method can estimate elevation and azimuth angles for noncoherent and/or coherent sources.Similarly, the Toeplitz matricesRyt, Rzt, and Rwt can be decomposed, respectively, as(23)Ryt=Eθ,ϕΦyStEHθ,ϕ+Nyt,Rzt=Eθ,ϕΦzStEHθ,ϕ+Nzt,Rwt=Eθ,ϕΦwStEHθ,ϕ+Nwt.The Toeplitz matrices defined above are rearranged to form anN×N×4 TWA F_t. This is represented as(24)F_t=F:,:,1;tF:,:,2;tF:,:,3;tF:,:,4;t=RxtRytRztRwt+NxtNytNztNwt=Eθ,ϕStEHθ,ϕEθ,ϕΦyStEHθ,ϕEθ,ϕΦzStEHθ,ϕEθ,ϕΦwStEHθ,ϕ+NxtNytNztNwt,where this TWA is in the same form as defined in (14) and the noise matrices Nxt, Nyt, Nzt, and Nwt are as defined earlier.The TWAF_ is rearranged into three slices of 2D matrices F1, F2, and F3 as follows:(25)F1=A⊙BCT+N1,F2=A⊙BCT+N2,F3=A⊙BCT+N3,where the matrices N1, N2, and N3 are the rearranged noise matrices, A=Eθ,ϕ, and B=EHθ,ϕ. The matrix C contains information about the elevation and azimuth angles. The structure of C is shown in the following:(26)C=c1Tc2Tc3Tc4T=diag-1Sdiag-1ΦySdiag-1ΦzSdiag-1ΦwS.The cost functions for the matricesA,B, andC to be minimized are defined as follows:(27)∫1A,B,C;F1=F1=A⊙BCTF2,∫2A,B,C;F2=F2=B⊙CATF2,∫3A,B,C;F3=F3=C⊙ABTF2.The alternating least squares method is applied to the cost function∫1A,B,C;F1 to find the matrix C, which contains the information about elevation and azimuth angles. This solution is given as follows: (28)C^=argminCF1-A^⊙B^C^TF2=A^⊙B^†F1T,where A^ and B^ are obtained, similarly, as follows:(29)A^=B^⊙C^†F2T,B^=C^⊙A^†F3T.The matricesA^, B^, and C^ contain the result of one iteration. Note that only a few iterations are required for the algorithm to converge.The COMFAC [34] algorithm is used to fit PARAFAC model defined in (15) to the TWA F_t. The COMFAC MATLAB function has the following form: A,B,C,∙,i=comfacF,P,∙,∙,∙,∙, where F_ (index t is removed for simplification) and P are the decomposing TWA and corresponding factor number, respectively. The factor number P represents the number of sources, which are K. The outputsA,B, andC represent the identification results (matrices), while i represents the iteration number required for the low-rank decomposition. The “∙” represents other options.After obtaining the identification matricesA^, B^, and C^ by performing PARAFAC analysis on F_, the diagonal matrix Φ^y is estimated from the identification matrix C^ as shown below(30)φ^yk=c^2kc^1k,where φ^yk is the estimated kth entry of the diagonal matrix Φ^y and c^1k and c^2k represent the kth entries of the row vectors c^1T and c^2T, respectively. Similarly, the kth diagonal entries of Φ^z and Φ^w are estimated as shown in the following:(31)φ^zk=c^3kc^1k,(32)φ^wk=c^4kc^1k,where c^3k and c^4k denote the kth entries of the row vectors c^3T and c^4T, respectively.The azimuth angleϕ^k is then estimated for the kth source using the following:(33)ϕ^k=tan-1∠φ^yk∠φ^wkand elevation angle θ^k for the kth source using the following:(34)θ^k=tan-1∠φ^wk∠φ^zk×cos⁡ϕ^k.The estimated azimuth and elevation angles in (33) and (34), respectively, are for the kth source. So, for K sources, the following pairs are obtained: θ1,ϕ1,θ2,ϕ2,…,θK,ϕK. The estimated identification matrices C^, A^, and B^ obtained earlier have the same column permutation matrix [31]. The pair-matching problem is automatically avoided since the kth column of A^ corresponds to the kth column of B^. Moreover, the elevation and azimuth angles estimated in (34) and (33), respectively, do not result in failure in practical mobile elevation angle range (70° to 90°).The estimation performance is improved by using more than one snapshot of the received signals. For this purpose, consider a TWAG_, which contains L snapshots of the received signals. The TWA G_ is constructed such that the L snapshots of the Toeplitz matrices defined earlier are concatenated together as shown in the following:(35)G_=F_1F_2⋯F_L=Rx1Rx2⋯RxLRy1Ry2⋯RyLRz1Rz2⋯RzLRw1Rw2⋯RwL.The PARAFAC analysis is then performed on TWA G followed by the estimation of azimuth and elevation angles in (33) and (34), respectively. The performance comparison using one and multiple snapshots (up to 25) is investigated in the simulation results section.


## 2.4. Summary of the Proposed Algorithm

As described earlier, the proposed algorithm considers the extended parallel-shaped array shown in Figure1. The problem of noncoherent and/or coherent sources is resolved by constructing the Toeplitz matrices and the problem of pair-matching between two or more sources is avoided using the PARAFAC model. The following steps summarize the proposed algorithm to estimate DOA for noncoherent and/or coherent sources without the problem of pair-matching.Step 1.
Construct the Toeplitz matricesRxt, Ryt, Rzt, and Rwt shown in (16), (17), (18), and (19), respectively, from single snapshot of the received signals.Step 2.
The Toeplitz matrices used to form theN×N×4 TWA F_t given in (24) for one snapshot or G_ in (35) for L snapshots are rearranged according to the PARAFAC model shown in (14).Step 3.
Apply alternating least squares (ALS) given in (28) and (29) to the TWA F_t to obtain the output matrix C, which has the same structure as shown in (26).Step 4.
Estimate the diagonal matricesΦ^y, Φ^z, and Φ^w using (30), (31), and (32), respectively. These matrices contain the elevation and azimuth angle information of the K sources.Step 5.
Estimate the azimuth angleϕ^k and elevation angle θ^k using (33) and (34), respectively.


## 2.5. Complexity and Cramer-Rao Bound Analysis of the Proposed Method

The computational complexity of the proposed method is compared with those of 2D DOA using cumulant-based method [15], novel 2D DOA estimation algorithm with L-shaped array [17], universal 2D DOA estimation [24], and joint elevation azimuth and elevation angle estimation [28]. For N antenna elements and K sources and considering only the major processing operations, the computational complexity of the proposed method is in the order of O4N+1Lp+n3K3+12N+12K, where Lp denotes the number of snapshots of the received signals used in the proposed method. The complexity of the novel L-shaped method in [17] is O4N-12Ln+n3K3+12N-12K, where Ln denotes the number of snapshots of the received signals. The cumulant-based method in [15] has complexity in the order of O212N+12Lc+n3K3+122N+12K, where Lc denotes the number of snapshots of the received signals. The complexity of the universal 2D DOA estimation in [24] is OLuN2+2N3+N2N-K+N+8N2Lu  +  OLuN2+2N3+V/ΔN-KNLu3, where Lu denotes the number of snapshots of the received signals, V is 180 grid number, and Δ is a number between 0 and 1. The first part represents the complexity to estimate elevation angle using Root-MUSIC algorithm and second part to estimate azimuth angle using MUSIC algorithm. The joint L-shaped method in [28] has complexity in the order of O2N+12Lj2+4/32N+13+20001×10N3+20001N-1N+1 (extracted from [28]), where Lj denotes the number of snapshots of the received signals. This implies that the methods in [15, 17, 24, 28] require a lot more computations compared to the proposed method. Moreover, the received signals’ snapshots used in the proposed method are far less than those used in the above methods; that is, Lp≪Ln,Lc,Lu, and Lj. This is described in Simulation Results.Cramer-Rao Bound (CRB) for 2-dimensional DOA estimation is derived in a similar manner as shown in [19, 35–37]. The received signals in (2), (3), (4), and (5) can be concatenated as shown in the following:(36)xtytztwt=Aθ,ϕAθ,ϕΦyθ,ϕAθ,ϕΦzθAθ,ϕΦwθ,ϕst+nxtnytnztnwt=Mst+Nt.Using the signal model above, CRB can be expressed as(37)CRB=σ22LRVHΠM⊥V⊙P^T-1,where V=[∂m1/∂θ1,∂m1/∂θ1,…, ∂mK/∂θK, ∂m1/∂θ1, ∂m2/∂θ2,…, ∂mK/∂θK], mi is the ith column of M, and(38)ΠM⊥=I8N-4-MMHM-1MH.The matrix P^ in (37) is given as(39)P^=P^sP^sP^sP^s,where P^s is as follows:(40)P^s=1L∑l=1LstsHt.The matrix P^s is not a diagonal matrix in the presence of coherent sources.


## 3. Simulation Results

The performance of the proposed method is evaluated in this section using simulations. The performance is compared with the novel L-shaped method in [17], cumulant-based method in [15], universal 2D DOA estimation in [24], and joint elevation and azimuth angles estimation using L-shaped array in [28]. The performance of the proposed method is also evaluated by considering noncoherent and/or coherent sources. Furthermore, the improvement in performance is also investigated by using more than one snapshot of the received signals. The number of snapshots considered in the proposed method is much less than the ones considered in [15, 17, 24, 28] as mentioned earlier. The extended parallel-shaped array model proposed in Figure 1 assumes the spacing between adjacent antenna elements as d=λ/2. The curves shown in the following subsections are obtained by averaging over several runs of the proposed algorithm.


### 3.1. Effect on Standard Deviation (STD) by Increasing SNR

Two noncoherent sourcesK=2 are considered here, which are located at θ1=70°,ϕ1=50° and θ2=80°,ϕ2=60°. The number of antenna elements on the positive y-axis is N=20. For the purpose of simulation, the number of snapshots considered for the proposed algorithm is L=25, and for the method in [15, 17, 24, 28] they are L=100, L=400, L=200, and L=200, respectively. Figure 2 shows the effect of standard deviation (STD) of the estimated elevation and azimuth angles against increasing SNR for sources 1 and 2. The comparison of the curves shows that the proposed method outperforms the methods [15, 17, 24, 28], at low SNR. The proposed method shows approximately 57% and 16% performance improvement at 5 dB as compared to the methods in [17] and [15], respectively. The performance is also better than the methods in [24, 28]. Moreover, the proposed method performs well even at much lower number of snapshots. The performance of the proposed method and the method in [28] is comparable with CRB. However, the method in [28] requires high number of snapshots.Figure 2
Standard deviation versus SNR for source 1 located atθ1=70°,ϕ1=50° and source 2 at θ2=80°,ϕ2=60°.


### 3.2. Effect on Standard Deviation (STD) by Increasing the Number of Snapshots

The effect of increasing the number of snapshots for the proposed method is investigated in the presence of two noncoherent sources. The same effect is also investigated by considering two coherent sources. The sources are located atθ1=70°,ϕ1=50° and θ2=80°,ϕ2=60°. The number of antenna elements on the positive y-axis is N=20. The curves are plotted for the following number of snapshots: L=1, L=2, L=4, L=8, and L=16. Figures 3 and 4 show the effect of increasing snapshots for noncoherent sources. Similarly, Figures 5 and 6 show the effect of increasing snapshots for coherent sources. The analysis of Figures 3, 4, 5, and 6 shows that the performance is improved by increasing the number of snapshots for both the noncoherent and coherent sources. However, the number of snapshots used here is far less than those used in the aforementioned methods.Figure 3
Standard deviation versus SNR for source 1 (noncoherent) located atθ1=70°,ϕ1=50° for the proposed method.Figure 4
Standard deviation versus SNR for source 2 (noncoherent) located atθ2=80°,ϕ2=60° for the proposed method.Figure 5
Standard deviation versus SNR for source 1 (coherent) located atθ1=70°,ϕ1=50° for the proposed method.Figure 6
Standard deviation versus SNR for source 2 (coherent) located atθ2=80°,ϕ2=60° for the proposed method.


### 3.3. Effect on Standard Deviation (STD) in the Presence of Noncoherent and Coherent Sources Together

The performance of the proposed method is evaluated here by considering a mixture of one noncoherent and two coherent sources. The number of antenna elements on the positivey-axis considered here is N=20. The number of snapshots is L=20. The two coherent sources are located at θ1=55°,ϕ1=40° and θ2=70°,ϕ2=50° and a noncoherent source is located at θ3=80°,ϕ3=60°. Figure 7 shows that the estimation performances are almost same for the mixture of coherent and noncoherent sources.Figure 7
Standard deviation versus SNR for the proposed method in the presence of two coherent sources located atθ1=55°,ϕ1=40° and θ2=70°,ϕ2=50° and a noncoherent source located at θ3=80°,ϕ3=60°.


### 3.4. Performance Comparison Using Scatter Plot

The scatter plot is obtained to make comparison between the proposed technique and the method in [17]. The number of antenna elements on the positive y-axis is N=18. The number of snapshots used is L=5. Figures 8 and 9 show the scatter plot for three noncoherent sources located at θ1=50°,ϕ1=40°, θ2=60°,ϕ2=55°, and θ3=70°,ϕ3=65° with SNR = 15 dB and 25 dB, respectively. The scatter plot is also obtained for the proposed method in the presence of noncoherent and coherent sources separately. The number of antenna elements on the positive y-axis is N=18. The number of snapshots used is L=5. Figure 10 shows the scatter plot for three noncoherent sources located at θ1=50°,ϕ1=40°, θ2=60°,ϕ2=55°, and θ3=70°,ϕ3=65° with SNR = 20 dB. Similarly, Figure 11 shows the scatter plot for three coherent sources located at θ1=50°,ϕ1=40°, θ2=60°,ϕ2=55°, and θ3=70°,ϕ3=65° with SNR = 20 dB. Again, Figures 10 and 11 show that the performance is still better than the method in [17] regardless of the noncoherent or coherent sources.Figure 8
Scatter plot for three noncoherent sources located atθ1=55°,ϕ1=40°, θ2=60°,ϕ2=55°, and θ3=70°,ϕ3=65° with SNR = 15 dB.Figure 9
Scatter plot for three noncoherent sources located atθ1=50°,ϕ1=40°, θ2=60°,ϕ2=55°, and θ3=70°,ϕ3=65° with SNR = 25 dB.Figure 10
Scatter plot for the proposed method in the presence of three noncoherent sources located atθ1=50°,ϕ1=40°, θ2=60°,ϕ2=55°, and θ3=70°,ϕ3=65° with SNR = 20 dB.Figure 11
Scatter plot for the proposed method in the presence of three coherent sources located atθ1=50°,ϕ1=40°, θ2=60°,ϕ2=55°, and θ3=70°,ϕ3=65° with SNR = 20 dB.In Figures12 and 13, the performance of the proposed method at low SNR is presented. The number of antenna elements on the positive y-axis is N=18. The numbers of snapshots used are L=1andL=5. Two coherent sources are located at θ1=40°,ϕ1=50° and θ2=60°,ϕ2=70°, and SNR = 0 dB. From the scatter plot in Figures 12 and 13, we observe that the performance of the proposed method degrades in estimating the azimuth and elevation angle at low SNR and at fewer snapshots. This concludes that in practical applications at low SNR several snapshots should be considered to get a reasonably accurate estimation.Figure 12
Scatter plot for the proposed method in the presence of two coherent sources located atθ1=40°,ϕ1=50° and θ2=60°,ϕ2=70° with SNR = 0 dB and single snapshot.Figure 13
Scatter plot for the proposed method in the presence of two coherent sources located atθ1=40°,ϕ1=50° and θ2=60°,ϕ2=70° with SNR = 0 dB and five snapshots.


### 3.5. Effect of Varying Elevation Angleθ

The comparison is made between the method in [15] and the proposed method by varying the elevation angle θ from 70° to 88° in steps of 3°. The azimuth angle ϕ is fixed at 20°. The number of antenna elements on the positive y-axis is fixed at N=20. The number of signal snapshots is L=25. The STD is obtained for both the methods at three different values of SNR, which are 5 dB, 15 dB, and 25 dB. The algorithm is run for 500 trials. Table 1 shows the comparison between two methods at different SNR values by varying the elevation angle θ. The comparison of STD shows that the proposed method outperforms the method in [15] by approximately 50%.Table 1
Comparison between the cumulant-based method in [15] and the proposed method.
Azimuth and elevation angles
Cumulant-based method in [15] standard deviation (degrees)
Proposed method standard deviation (degree)
θ
,
ϕ
SNR = 5 dB
SNR = 15 dB
SNR = 25 dB
SNR = 5 dB
SNR = 15 dB
SNR = 25 dB
7
0
∘
,
2
0
∘
0.405
0.121
0.035
0.338
0.117
0.023
7
3
∘
,
2
0
∘
0.732
0.245
0.073
0.420
0.132
0.045
7
6
∘
,
2
0
∘
0.794
0.237
0.078
0.447
0.135
0.044
7
9
∘
,
2
0
∘
0.705
0.226
0.074
0.413
0.138
0.045
8
2
∘
,
2
0
∘
0.743
0.248
0.077
0.486
0.142
0.047
8
5
∘
,
2
0
∘
0.734
0.227
0.068
0.470
0.143
0.046
8
8
∘
,
2
0
∘
0.729
0.232
0.076
0.459
0.145
0.047


## 3.1. Effect on Standard Deviation (STD) by Increasing SNR

Two noncoherent sourcesK=2 are considered here, which are located at θ1=70°,ϕ1=50° and θ2=80°,ϕ2=60°. The number of antenna elements on the positive y-axis is N=20. For the purpose of simulation, the number of snapshots considered for the proposed algorithm is L=25, and for the method in [15, 17, 24, 28] they are L=100, L=400, L=200, and L=200, respectively. Figure 2 shows the effect of standard deviation (STD) of the estimated elevation and azimuth angles against increasing SNR for sources 1 and 2. The comparison of the curves shows that the proposed method outperforms the methods [15, 17, 24, 28], at low SNR. The proposed method shows approximately 57% and 16% performance improvement at 5 dB as compared to the methods in [17] and [15], respectively. The performance is also better than the methods in [24, 28]. Moreover, the proposed method performs well even at much lower number of snapshots. The performance of the proposed method and the method in [28] is comparable with CRB. However, the method in [28] requires high number of snapshots.Figure 2
Standard deviation versus SNR for source 1 located atθ1=70°,ϕ1=50° and source 2 at θ2=80°,ϕ2=60°.


## 3.2. Effect on Standard Deviation (STD) by Increasing the Number of Snapshots

The effect of increasing the number of snapshots for the proposed method is investigated in the presence of two noncoherent sources. The same effect is also investigated by considering two coherent sources. The sources are located atθ1=70°,ϕ1=50° and θ2=80°,ϕ2=60°. The number of antenna elements on the positive y-axis is N=20. The curves are plotted for the following number of snapshots: L=1, L=2, L=4, L=8, and L=16. Figures 3 and 4 show the effect of increasing snapshots for noncoherent sources. Similarly, Figures 5 and 6 show the effect of increasing snapshots for coherent sources. The analysis of Figures 3, 4, 5, and 6 shows that the performance is improved by increasing the number of snapshots for both the noncoherent and coherent sources. However, the number of snapshots used here is far less than those used in the aforementioned methods.Figure 3
Standard deviation versus SNR for source 1 (noncoherent) located atθ1=70°,ϕ1=50° for the proposed method.Figure 4
Standard deviation versus SNR for source 2 (noncoherent) located atθ2=80°,ϕ2=60° for the proposed method.Figure 5
Standard deviation versus SNR for source 1 (coherent) located atθ1=70°,ϕ1=50° for the proposed method.Figure 6
Standard deviation versus SNR for source 2 (coherent) located atθ2=80°,ϕ2=60° for the proposed method.


## 3.3. Effect on Standard Deviation (STD) in the Presence of Noncoherent and Coherent Sources Together

The performance of the proposed method is evaluated here by considering a mixture of one noncoherent and two coherent sources. The number of antenna elements on the positivey-axis considered here is N=20. The number of snapshots is L=20. The two coherent sources are located at θ1=55°,ϕ1=40° and θ2=70°,ϕ2=50° and a noncoherent source is located at θ3=80°,ϕ3=60°. Figure 7 shows that the estimation performances are almost same for the mixture of coherent and noncoherent sources.Figure 7
Standard deviation versus SNR for the proposed method in the presence of two coherent sources located atθ1=55°,ϕ1=40° and θ2=70°,ϕ2=50° and a noncoherent source located at θ3=80°,ϕ3=60°.


## 3.4. Performance Comparison Using Scatter Plot

The scatter plot is obtained to make comparison between the proposed technique and the method in [17]. The number of antenna elements on the positive y-axis is N=18. The number of snapshots used is L=5. Figures 8 and 9 show the scatter plot for three noncoherent sources located at θ1=50°,ϕ1=40°, θ2=60°,ϕ2=55°, and θ3=70°,ϕ3=65° with SNR = 15 dB and 25 dB, respectively. The scatter plot is also obtained for the proposed method in the presence of noncoherent and coherent sources separately. The number of antenna elements on the positive y-axis is N=18. The number of snapshots used is L=5. Figure 10 shows the scatter plot for three noncoherent sources located at θ1=50°,ϕ1=40°, θ2=60°,ϕ2=55°, and θ3=70°,ϕ3=65° with SNR = 20 dB. Similarly, Figure 11 shows the scatter plot for three coherent sources located at θ1=50°,ϕ1=40°, θ2=60°,ϕ2=55°, and θ3=70°,ϕ3=65° with SNR = 20 dB. Again, Figures 10 and 11 show that the performance is still better than the method in [17] regardless of the noncoherent or coherent sources.Figure 8
Scatter plot for three noncoherent sources located atθ1=55°,ϕ1=40°, θ2=60°,ϕ2=55°, and θ3=70°,ϕ3=65° with SNR = 15 dB.Figure 9
Scatter plot for three noncoherent sources located atθ1=50°,ϕ1=40°, θ2=60°,ϕ2=55°, and θ3=70°,ϕ3=65° with SNR = 25 dB.Figure 10
Scatter plot for the proposed method in the presence of three noncoherent sources located atθ1=50°,ϕ1=40°, θ2=60°,ϕ2=55°, and θ3=70°,ϕ3=65° with SNR = 20 dB.Figure 11
Scatter plot for the proposed method in the presence of three coherent sources located atθ1=50°,ϕ1=40°, θ2=60°,ϕ2=55°, and θ3=70°,ϕ3=65° with SNR = 20 dB.In Figures12 and 13, the performance of the proposed method at low SNR is presented. The number of antenna elements on the positive y-axis is N=18. The numbers of snapshots used are L=1andL=5. Two coherent sources are located at θ1=40°,ϕ1=50° and θ2=60°,ϕ2=70°, and SNR = 0 dB. From the scatter plot in Figures 12 and 13, we observe that the performance of the proposed method degrades in estimating the azimuth and elevation angle at low SNR and at fewer snapshots. This concludes that in practical applications at low SNR several snapshots should be considered to get a reasonably accurate estimation.Figure 12
Scatter plot for the proposed method in the presence of two coherent sources located atθ1=40°,ϕ1=50° and θ2=60°,ϕ2=70° with SNR = 0 dB and single snapshot.Figure 13
Scatter plot for the proposed method in the presence of two coherent sources located atθ1=40°,ϕ1=50° and θ2=60°,ϕ2=70° with SNR = 0 dB and five snapshots.


## 3.5. Effect of Varying Elevation Angleθ

The comparison is made between the method in [15] and the proposed method by varying the elevation angle θ from 70° to 88° in steps of 3°. The azimuth angle ϕ is fixed at 20°. The number of antenna elements on the positive y-axis is fixed at N=20. The number of signal snapshots is L=25. The STD is obtained for both the methods at three different values of SNR, which are 5 dB, 15 dB, and 25 dB. The algorithm is run for 500 trials. Table 1 shows the comparison between two methods at different SNR values by varying the elevation angle θ. The comparison of STD shows that the proposed method outperforms the method in [15] by approximately 50%.Table 1
Comparison between the cumulant-based method in [15] and the proposed method.
Azimuth and elevation angles
Cumulant-based method in [15] standard deviation (degrees)
Proposed method standard deviation (degree)
θ
,
ϕ
SNR = 5 dB
SNR = 15 dB
SNR = 25 dB
SNR = 5 dB
SNR = 15 dB
SNR = 25 dB
7
0
∘
,
2
0
∘
0.405
0.121
0.035
0.338
0.117
0.023
7
3
∘
,
2
0
∘
0.732
0.245
0.073
0.420
0.132
0.045
7
6
∘
,
2
0
∘
0.794
0.237
0.078
0.447
0.135
0.044
7
9
∘
,
2
0
∘
0.705
0.226
0.074
0.413
0.138
0.045
8
2
∘
,
2
0
∘
0.743
0.248
0.077
0.486
0.142
0.047
8
5
∘
,
2
0
∘
0.734
0.227
0.068
0.470
0.143
0.046
8
8
∘
,
2
0
∘
0.729
0.232
0.076
0.459
0.145
0.047


## 4. Conclusion

In this paper, 2D-DOA estimation method for estimating the elevation and azimuth angles is proposed. The proposed method estimates the 2D DOA for noncoherent and/or coherent sources employing a new antenna array configuration. It has lower computational complexity compared to the existing schemes since it requires only single or few snapshots and is based on first-order data matrices. In addition, it has no failure estimation especially in typical mobile communication. Furthermore, a parallel factor (PARAFAC) model has been derived to avoid pair-matching problems. These advantages make the proposed algorithm a realistic candidate for high speed wireless communications applications.



---

*Source: 1813497-2017-07-03.xml*