Parallel Factor-Based Model for Two-Dimensional Direction Estimation

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.

[ProQuest: [...] denotes non US-ASCII text; see PDF]

Nizar Tayem 1 and Khaqan Majeed 1 and Ahmed A. Hussain 1

Academic Editor:Angelo Liseno

Department of Electrical Engineering, Prince Mohammad Bin Fahd University (PMU), Al-Khobar, Saudi Arabia

Received 8 February 2017; Revised 22 May 2017; Accepted 1 June 2017; 3 July 2017

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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 Section 2. 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 , ^(·)[dagger] , _(·)F , R(·), [ecedil]9;, [ecedil]7;, and [ecedil]a; 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 [upsilon], constant by V , vector by v , matrix by V , and ijth member of a matrix V by _[upsilon]i,j . A three-way array (TWA) is denoted by V and its ijkth member is given by _{[upsilon]i,j,k} . The operations diag(v ) and ^diag-1 (V ) denote conversion of the vector v into a diagonal matrix and the diagonal matrix V 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 Figure 1. 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.

[figure omitted; refer to PDF]

Consider K narrowband noncoherent and/or coherent sources in far field of the extended parallel-shaped array. Let _θk and _[varphi]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 [figure omitted; refer to PDF] 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 [figure omitted; refer to PDF] where the matrix A(θ,[varphi]) is defined as [figure omitted; refer to PDF] The dimensions of A(θ,[varphi]) are (2N-1)×K. The entry _uk in (6) is defined as [figure omitted; refer to PDF] where _θk and _[varphi]k are the elevation and azimuth angles, respectively, defined earlier for the kth source. The signal vector s(t) in (2) to (5) is defined as [figure omitted; refer to PDF] 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 [figure omitted; refer to PDF] where _qK is given as [figure omitted; refer to PDF] Similarly, _Φz and _Φw in (4) and (5), respectively, are defined as [figure omitted; refer to PDF] where _rK is given as [figure omitted; refer to PDF] where _vK is given as [figure omitted; refer to PDF]

The (2N-1)×1-dimensional vectors _nx (t), _ny (t), _nz (t), and _nw (t) 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 (_nx (t),_ny (t),_nz (t),or _nw (t)) 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 TWA D_ with loading matrices A , B , and C . The array D_ is defined as [figure omitted; refer to PDF] The P-component trilinear decomposition of the array D_ is defined as [figure omitted; refer to PDF] where i=1,2,...,_Id , j=1,2,...,_Jd , k=1,2,...,_Kd , and _di,j,k is the (i,j,k)th element of the array D_. The element _ai,P is the (i,P)th element of the (_Id ×P)-dimensional matrix A . Similarly, _bi,P and _ck,P are (j,P)th and (k,P)th elements of the (_Jd ×P)-dimensional and (_Kd ×P)-dimensional matrices B and C , 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 subarrays X,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 _Rx (t) for the subarray X as [figure omitted; refer to PDF] 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 [figure omitted; refer to PDF] 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 [figure omitted; refer to PDF] The dimensions of the matrices _Rx (t), _Ry (t), _Rz (t), and _Rw (t) are N×N each. These matrices are obtained at the tth snapshot of the received signals.

The Toeplitz matrix _Rx (t) in (16) can be decomposed as follows: [figure omitted; refer to PDF] where _Nx (t) is a noise matrix with dimensions N×N and E(θ,[varphi]) and S(t) are defined as [figure omitted; refer to PDF] where _uk is defined in (7) and S(t) is a diagonal matrix containing the signals from K sources. The dimensions of E(θ,[varphi]) and S(t) 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(θ,[varphi]) is K, which is full column rank. The diagonal matrix S(t) 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, (^{(_s1 (t)0[...]0)T} ,^{(0_s2 (t)[...]0)T} )=0. Now, consider two coherent sources such that _s2 (t)=α_s1 (t). The dot product corresponding to these two columns is again zero; that is, (^{(_s1 (t)0[...]0)T} ,^{(0α_s1 (t)[...]0)T} )=0. This indicates that the matrix S(t) 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 _Rx (t) 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 r(UV^UH )=r(V), where r(V) denotes the rank of the matrix V, where _Ny (t), _Nz (t), and _Nw (t) 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 matrices _Ry (t), _Rz (t), and _Rw (t) can be decomposed, respectively, as [figure omitted; refer to PDF]

The Toeplitz matrices defined above are rearranged to form an N×N×4 TWA F_(t). This is represented as [figure omitted; refer to PDF] where this TWA is in the same form as defined in (14) and the noise matrices _Nx (t), _Ny (t), _Nz (t), and _Nw (t) are as defined earlier.

The TWA F_ is rearranged into three slices of 2D matrices _F1 , _F2 , and _F3 as follows: [figure omitted; refer to PDF] where the matrices _N1 , _N2 , and _N3 are the rearranged noise matrices, A=E(θ,[varphi]), and B=^EH (θ,[varphi]). The matrix C contains information about the elevation and azimuth angles. The structure of C is shown in the following: [figure omitted; refer to PDF]

The cost functions for the matrices A , B , and C to be minimized are defined as follows: [figure omitted; refer to PDF]

The alternating least squares method is applied to the cost function _∫1 (A,B,C;_F1 ) to find the matrix C, which contains the information about elevation and azimuth angles. This solution is given as follows: [figure omitted; refer to PDF] where A^ and B^ are obtained, similarly, as follows: [figure omitted; refer to PDF]

The matrices A^, 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)=comfac(F,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 outputs A , B , and C 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 matrices A^, B^, and C^ by performing PARAFAC analysis on F_, the diagonal matrix _Φ^y is estimated from the identification matrix C^ as shown below [figure omitted; refer to PDF] where _{[straight phi]^y} (k) is the estimated kth entry of the diagonal matrix _Φ^y and _c^1 (k) and _c^2 (k) 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: [figure omitted; refer to PDF] where _c^3 (k) and _c^4 (k) denote the kth entries of the row vectors ^{c^3T} and ^{c^4T} , respectively.

The azimuth angle _[varphi]^k is then estimated for the kth source using the following: [figure omitted; refer to PDF] and elevation angle _θ^k for the kth source using the following: [figure omitted; refer to PDF]

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 ,_[varphi]1 ),(_θ2 ,_[varphi]2 ),...,(_θK ,_[varphi]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 TWA G_, 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: [figure omitted; refer to PDF] 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 Figure 1. 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 matrices _Rx (t), _Ry (t), _Rz (t), and _Rw (t) shown in (16), (17), (18), and (19), respectively, from single snapshot of the received signals.

Step 2.

The Toeplitz matrices used to form the N×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 _[varphi]^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 O(4(N+1)_Lp +n(3^K3 +12^(N+1)2 K)), 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 O(4^(N-1)2_Ln +n(3^K3 +12^(N-1)2 K)), where _Ln denotes the number of snapshots of the received signals. The cumulant-based method in [15] has complexity in the order of O(21^(2N+1)2_Lc +n(3^K3 +12^(2N+1)2 K)), where _Lc denotes the number of snapshots of the received signals. The complexity of the universal 2D DOA estimation in [24] is O(_Lu^N2 +2^N3 +^N2 (N-K)+N+8^N2_Lu ) + O(_Lu^N2 +2^N3 +(V/Δ)(N-K)N^Lu3 ), 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 O(^(2N+1)2Lj2 +(4/3)^(2N+1)3 +20001 × 10^N3 +20001(N-1)(N+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: [figure omitted; refer to PDF] Using the signal model above, CRB can be expressed as [figure omitted; refer to PDF] where V=[∂_m1 /∂_θ1 ,∂_m1 /∂_θ1 ,..., ∂_mK /∂_θK , ∂_m1 /∂_θ1 , ∂_m2 /∂_θ2 ,..., ∂_mK /∂_θK ], _mi is the ith column of M, and [figure omitted; refer to PDF] The matrix P^ in (37) is given as [figure omitted; refer to PDF] where _P^s is as follows: [figure omitted; refer to PDF] 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 sources (K=2) are considered here, which are located at (_θ1 =^70° ,_[varphi]1 =^50° ) and (_θ2 =^80° ,_[varphi]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° ,_[varphi]1 =^50° ) and source 2 at (_θ2 =^80° ,_[varphi]2 =^60° ).

[figure omitted; refer to PDF]

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° ,_[varphi]1 =^50° ) and (_θ2 =^80° ,_[varphi]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° ,_[varphi]1 =^50° ) for the proposed method.

[figure omitted; refer to PDF]

Figure 4: Standard deviation versus SNR for source 2 (noncoherent) located at (_θ2 =^80° ,_[varphi]2 =^60° ) for the proposed method.

[figure omitted; refer to PDF]

Figure 5: Standard deviation versus SNR for source 1 (coherent) located at (_θ1 =^70° ,_[varphi]1 =^50° ) for the proposed method.

[figure omitted; refer to PDF]

Figure 6: Standard deviation versus SNR for source 2 (coherent) located at (_θ2 =^80° ,_[varphi]2 =^60° ) for the proposed method.

[figure omitted; refer to PDF]

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 positive y-axis considered here is N=20. The number of snapshots is L=20. The two coherent sources are located at (_θ1 =^55° ,_[varphi]1 =^40° ) and (_θ2 =^70° ,_[varphi]2 =^50° ) and a noncoherent source is located at (_θ3 =^80° ,_[varphi]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° ,_[varphi]1 =^40° ) and (_θ2 =^70° ,_[varphi]2 =^50° ) and a noncoherent source located at (_θ3 =^80° ,_[varphi]3 =^60° ).

[figure omitted; refer to PDF]

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° ,_[varphi]1 =^40° ), (_θ2 =^60° ,_[varphi]2 =^55° ), and (_θ3 =^70° ,_[varphi]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° ,_[varphi]1 =^40° ), (_θ2 =^60° ,_[varphi]2 =^55° ), and (_θ3 =^70° ,_[varphi]3 =^65° ) with SNR = 20 dB. Similarly, Figure 11 shows the scatter plot for three coherent sources located at (_θ1 =^50° ,_[varphi]1 =^40° ), (_θ2 =^60° ,_[varphi]2 =^55° ), and (_θ3 =^70° ,_[varphi]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° ,_[varphi]1 =^40° ), (_θ2 =^60° ,_[varphi]2 =^55° ), and (_θ3 =^70° ,_[varphi]3 =^65° ) with SNR = 15 dB.

[figure omitted; refer to PDF]

Figure 9: Scatter plot for three noncoherent sources located at (_θ1 =^50° ,_[varphi]1 =^40° ), (_θ2 =^60° ,_[varphi]2 =^55° ), and (_θ3 =^70° ,_[varphi]3 =^65° ) with SNR = 25 dB.

[figure omitted; refer to PDF]

Figure 10: Scatter plot for the proposed method in the presence of three noncoherent sources located at (_θ1 =^50° ,_[varphi]1 =^40° ), (_θ2 =^60° ,_[varphi]2 =^55° ), and (_θ3 =^70° ,_[varphi]3 =^65° ) with SNR = 20 dB.

[figure omitted; refer to PDF]

Figure 11: Scatter plot for the proposed method in the presence of three coherent sources located at (_θ1 =^50° ,_[varphi]1 =^40° ), (_θ2 =^60° ,_[varphi]2 =^55° ), and (_θ3 =^70° ,_[varphi]3 =^65° ) with SNR = 20 dB.

[figure omitted; refer to PDF]

In Figures 12 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=1 and L=5. Two coherent sources are located at (_θ1 =^40° ,_[varphi]1 =^50° ) and (_θ2 =^60° ,_[varphi]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° ,_[varphi]1 =^50° ) and (_θ2 =^60° ,_[varphi]2 =^70° ) with SNR = 0 dB and single snapshot.

[figure omitted; refer to PDF]

Figure 13: Scatter plot for the proposed method in the presence of two coherent sources located at (_θ1 =^40° ,_[varphi]1 =^50° ) and (_θ2 =^60° ,_[varphi]2 =^70° ) with SNR = 0 dB and five snapshots.

[figure omitted; refer to PDF]

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 [varphi] 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)
( θ , [varphi] )	SNR = 5 dB	SNR = 15 dB	SNR = 25 dB	SNR = 5 dB	SNR = 15 dB	SNR = 25 dB
( 7 0 [composite function] , 2 0 [composite function] )	0.405	0.121	0.035	0.338	0.117	0.023
( 7 3 [composite function] , 2 0 [composite function] )	0.732	0.245	0.073	0.420	0.132	0.045
( 7 6 [composite function] , 2 0 [composite function] )	0.794	0.237	0.078	0.447	0.135	0.044
( 7 9 [composite function] , 2 0 [composite function] )	0.705	0.226	0.074	0.413	0.138	0.045
( 8 2 [composite function] , 2 0 [composite function] )	0.743	0.248	0.077	0.486	0.142	0.047
( 8 5 [composite function] , 2 0 [composite function] )	0.734	0.227	0.068	0.470	0.143	0.046
( 8 8 [composite function] , 2 0 [composite function] )	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.

References

[1] T. Lv, F. Tan, H. Gao, S. Yang, "A beamspace approach for 2-D localization of incoherently distributed sources in massive MIMO systems,", Signal Processing , vol. 121, pp. 30-45, 2016.

[2] K.-J. Xu, W.-K. Nie, D.-Z. Feng, X.-J. Chen, D.-Y. Fang, "A multi-direction virtual array transformation algorithm for 2D DOA estimation,", Signal Processing , vol. 125, pp. 122-133, 2016.

[3] M. F. Khan, M. Tufail, "Computationally efficient 2D beamspace matrix pencil method for direction of arrival estimation,", Digital Signal Processing: A Review Journal , vol. 20, no. 6, pp. 1526-1534, 2010.

[4] A. B. Gershman, M. Rübsamen, M. Pesavento, "One- and two-dimensional direction-of-arrival estimation: an overview of search-free techniques,", Signal Processing , vol. 90, no. 5, pp. 1338-1349, 2010.

[5] R. O. Schmidt, "Multiple emitter location and signal parameter estimation,", IEEE Transactions on Antennas and Propagation , vol. 34, no. 3, pp. 276-280, 1986.

[6] R. Roy, T. Kailath, "ESPRIT-estimation of signal parameters via rotational invariance techniques,", IEEE Transactions on Acoustics, Speech and Signal Processing , vol. 37, no. 7, pp. 984-995, 1989.

[7] S. U. Pillai, B. H. Kwon, "Forward/backward spatial smoothing techniques for coherent signal identification,", IEEE Transactions on Acoustics, Speech, and Signal Processing , vol. 37, no. 1, pp. 8-15, 1989.

[8] W. Du, R. L. Kirlin, "Improved spatial smoothing techniques for DOA estimation of coherent signals,", IEEE Transactions on Signal Processing , vol. 39, no. 5, pp. 1208-1210, 1991.

[9] S. Marcos, A. Marsal, M. Benidir, "The propagator method for source bearing estimation,", Signal Processing , vol. 42, no. 2, pp. 121-138, 1995.

[10] Y. Wu, G. Liao, H. C. So, "A fast algorithm for 2-D direction-of-arrival estimation,", Signal Processing , vol. 83, no. 8, pp. 1827-1831, 2003.

[11] F.-J. Chen, S. Kwong, C.-W. Kok, "ESPRIT-like two-dimensional DOA estimation for coherent signals,", IEEE Transactions on Aerospace and Electronic Systems , vol. 46, no. 3, pp. 1477-1484, 2010.

[12] Y. Zhang, Z. Ye, X. Xu, J. Cui, "Estimation of two-dimensional direction-of-arrival for uncorrelated and coherent signals with low complexity,", IET Radar, Sonar and Navigation , vol. 4, no. 4, pp. 507-519, 2010.

[13] N. Tayem, M. Omer, A. Abul Hussain, "DOA estimation method using R matrix of the QR factorized data and its prototype implementation on NI-PXI platform," in Proceedings of the 33rd Annual IEEE Military Communications Conference (MILCOM '14), pp. 333-337, Baltimore, Md, USA, October 2014.

[14] P. Chevalier, L. Albera, A. Ferreol, P. Comon, "On the virtual array concept for higher order array processing,", IEEE Transactions on Signal Processing , vol. 53, no. 4, pp. 1254-1271, 2005.

[15] J. L. Liang, "Joint azimuth and elevation direction finding using cumulant,", IEEE Sensors Journal , vol. 9, no. 4, pp. 390-398, 2009.

[16] R. Bro, N. D. Sidiropoulos, G. B. Giannakis, "Optimal joint azimuth-elevation and signal-array response estimation using parallel factor analysis," in Proceedings of the 32nd Asilomar Conference on Signals, Systems & Computers, vol. 2, pp. 1594-1598, November 1998.

[17] Z. Xiaofei, L. Jianfeng, X. Lingyun, "Novel two-dimensional DOA estimation with L-shaped array,", EURASIP Journal on Advances in Signal Processing , vol. 1, pp. 1-7, 2011.

[18] D. Liu, J. Liang, "L-shaped array-based 2-D DOA estimation using parallel factor analysis," in Proceedings of the 2010 8th World Congress on Intelligent Control and Automation, WCICA 2010, pp. 6949-6952, chn, July 2010.

[19] H. Chen, C.-P. Hou, Q. Wang, L. Huang, W.-Q. Yan, "Cumulants-based toeplitz matrices reconstruction method for 2-D coherent DOA estimation,", IEEE Sensors Journal , vol. 14, no. 8, pp. 2824-2832, 2014.

[20] H. Tao, J. Xin, J. Wang, N. Zheng, A. Sano, "Two-dimensional direction estimation for a mixture of noncoherent and coherent signals,", IEEE Transactions on Signal Processing , vol. 63, no. 2, pp. 318-333, 2015.

[21] R. T. Behrens, L. L. Scharf, "Signal processing applications of oblique projection operators,", IEEE Transactions on Signal Processing , vol. 42, no. 6, pp. 1413-1424, 1994.

[22] R. Boyer, G. Bouleux, "Oblique projections for direction-of-arrival estimation with prior knowledge,", IEEE Transactions on Signal Processing , vol. 56, no. 4, pp. 1374-1387, 2008.

[23] G. Wang, J. Xin, J. Wang, N. Zheng, A. Sano, "Subspace-based two-dimensional direction estimation and tracking of multiple targets,", IEEE Transactions on Aerospace and Electronic Systems , vol. 51, no. 2, pp. 1386-1402, 2015.

[24] J. Jia-Jia, D. Fa-Jie, W. Xian-Quan, "A universal two-dimensional direction of arrival estimation method without parameter match,", IEEE Sensors Journal , vol. 16, no. 9, pp. 3141-3146, 2016.

[25] N. Tayem, H. M. Kwon, "Azimuth and elevation angle estimation with no failure and no eigen decomposition,", Signal Processing , vol. 86, no. 1, pp. 8-16, 2006.

[26] N. Tayem, H. M. Kwon, "L-shape 2-dimensional arrival angle estimation with propagator method," in Proceedings of the 2005 IEEE 61st Vehicular Technology Conference, vol. 1, pp. 6-10, 2005.

[27] T. Shu, X. Liu, J. Lu, "Comments on "L-Shape 2-dimensional arrival angle estimation with propagator method",", IEEE Transactions on Antennas and Propagation , vol. 56, no. 5, pp. 1502-1503, 2008.

[28] J. Liang, D. Liu, "Joint elevation and azimuth direction finding using L-shaped array,", Institute of Electrical and Electronics Engineers. Transactions on Antennas and Propagation , vol. 58, no. 6, pp. 2136-2141, 2010.

[29] L. Zou, J. Lasenby, Z. He, "Direction and polarisation estimation using polarised cylindrical conformal arrays,", IET Signal Processing , vol. 6, no. 5, pp. 395-403, 2012.

[30] W. Si, L. Wan, L. Liu, Z. Tian, "Fast estimation of frequency and 2-D doas for cylindrical conformal array antenna using state-space and propagator method,", Progress in Electromagnetics Research , vol. 137, pp. 51-71, 2013.

[31] L. Wan, W. Si, L. Liu, Z. Tian, N. Feng, "High accuracy 2D-DOA estimation for conformal array using PARAFAC,", International Journal of Antennas and Propagation , vol. 2014, 2014.

[32] R. A. Harshman Foundations of the PARAFAC procedure: Model and conditions for an "explanatory" multi-mode factor analysis

[33] N. Tayem, K. Majeed, A. A. Hussain, "Two-Dimensional DOA Estimation Using Cross-Correlation Matrix With L-Shaped Array,", IEEE Antennas and Wireless Propagation Letters , vol. 15, pp. 1077-1080, 2016.

[34] N. D. Sidiropoulos, "COMFAC: Matlab Code for LS Fitting of the Complex PARAFAC Model in 3-D," http://www.telecom.tuc.gr/~nikos

[35] T. Xia, Y. Zheng, Q. Wan, X. Wang, "Decoupled estimation of 2-D angles of arrival using two parallel uniform linear arrays,", IEEE Transactions on Antennas and Propagation , vol. 55, no. 9, pp. 2627-2632, 2007.

[36] P. Stoica, A. Nehorai, "MUSIC, maximum likelihood, and cramér-rao bound: further results and comparisons,", IEEE Transactions on Acoustics, Speech, and Signal Processing , vol. 38, no. 12, pp. 2140-2150, 1990.

[37] P. Stoica, A. Nehorai, "Performance study of conditional and unconditional direction-of-arrival estimation,", IEEE Transactions on Acoustics, Speech, and Signal Processing , vol. 38, no. 10, pp. 1783-1795, 1990.