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.
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 -axis. Next, it obtains the azimuth angle estimate using the elevation angle estimated earlier and a 2D uniform linear subarray along the -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 , , , , , , , 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 by V, vector by v, matrix by V, and th member of a matrix V by . A three-way array (TWA) is denoted by V and its th member is given by . The operations diag(v) and (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 -axis, plane, and plane. The middle array contains one extra element as compared to the other arrays. There are antenna elements in the middle array. The arrays above and below the middle array contain elements each. These arrays are partitioned into four subarrays such that each subarray contains elements. These are denoted as , , , and as shown in the figure. The distance between adjacent elements on a single array is meters. The array on the -axis is separated by meters from the arrays in the and planes.
Consider narrowband noncoherent and/or coherent sources in far field of the extended parallel-shaped array. Let and denote the elevation and azimuth angles, respectively, for the th source. The signal vectors received at the , , , and subarrays at the th snapshot are given, respectively, aswhere and superscript denotes the transpose. The -dimensional signal vectors received at the antenna arrays can be represented aswhere the matrix is defined as The dimensions of are . The entry in (6) is defined as where and are the elevation and azimuth angles, respectively, defined earlier for the th source. The signal vector in (2) to (5) is defined asThe matrix in (3) is diagonal and contains information about the elevation and azimuth angles. Its dimensions are and is defined aswhere is given asSimilarly, and in (4) and (5), respectively, are defined aswhere is given aswhere is given as
The -dimensional vectors , , , and 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 () is also independent of other elements in that vector. The noise is assumed to be Additive White Gaussian (AWG) with mean 0 and variance . The noise vectors are also independent of the source signals. The goal now is to estimate the diagonal matrices , , and , 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 with loading matrices A, B, and C. The array is defined asThe -component trilinear decomposition of the array is defined aswhere , and is the element of the array . The element is the element of the -dimensional matrix A. Similarly, and are and elements of the -dimensional and -dimensional matrices B and C, respectively. The matrices , , and are called the mode or loading matrices for a given PARAFAC model. The PARAFAC analysis of 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 , and from a single snapshot. As inferred from the array geometry in Figure 1, the array elements are indexed from to for the array in the middle and to for the other two arrays. Using the index 0 element as a reference, we define the Toeplitz matrix for the subarray asThe reference element for subarray is at index 1 since it is defined as subarray shifted right by one element. Similarly, for subarray , the Toeplitz matrix is defined asThe reference antenna elements for and subarrays both are at index 0 and the Toeplitz matrices for subarrays and are defined, respectively, asThe dimensions of the matrices , , , and are each. These matrices are obtained at the th snapshot of the received signals.
The Toeplitz matrix in (16) can be decomposed as follows:where is a noise matrix with dimensions and and are defined as where is defined in (7) and is a diagonal matrix containing the signals from sources. The dimensions of and are and , respectively. Since the number of antenna elements is usually greater than the number of sources so the rank of is , which is full column rank. The diagonal matrix is of full rank whether the 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, . Now, consider two coherent sources such that . The dot product corresponding to these two columns is again zero; that is, . This indicates that the matrix is of full rank, which is equal to , the number of sources. So the rank is actually the minimum of and , which indicates that the rank is equal to the number of sources . As inferred from (20), the rank of Toeplitz matrix is also. Consider the matrices and with dimensions and , respectively, where the matrix is of full column rank. The theorem states that , where denotes the rank of the matrix , where , , and are the noise matrices for the subarrays , , and , respectively. Therefore, the proposed method can estimate elevation and azimuth angles for noncoherent and/or coherent sources.
Similarly, the Toeplitz matrices , , and can be decomposed, respectively, as
The Toeplitz matrices defined above are rearranged to form an TWA . This is represented aswhere this TWA is in the same form as defined in (14) and the noise matrices , , , and are as defined earlier.
The TWA is rearranged into three slices of 2D matrices , , and as follows:where the matrices , , and are the rearranged noise matrices, , and . The matrix contains information about the elevation and azimuth angles. The structure of is shown in the following:
The cost functions for the matrices A, B, and C to be minimized are defined as follows:
The alternating least squares method is applied to the cost function to find the matrix , which contains the information about elevation and azimuth angles. This solution is given as follows: where and are obtained, similarly, as follows:
The matrices , , and 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 . The COMFAC MATLAB function has the following form: , where (index is removed for simplification) and are the decomposing TWA and corresponding factor number, respectively. The factor number represents the number of sources, which are . The outputs A, B, and C represent the identification results (matrices), while represents the iteration number required for the low-rank decomposition. The “” represents other options.
After obtaining the identification matrices , , and by performing PARAFAC analysis on , the diagonal matrix is estimated from the identification matrix as shown belowwhere is the estimated th entry of the diagonal matrix and and represent the th entries of the row vectors and , respectively. Similarly, the th diagonal entries of and are estimated as shown in the following:where and denote the th entries of the row vectors and , respectively.
The azimuth angle is then estimated for the th source using the following:and elevation angle for the th source using the following:
The estimated azimuth and elevation angles in (33) and (34), respectively, are for the th source. So, for sources, the following pairs are obtained: . The estimated identification matrices , , and obtained earlier have the same column permutation matrix [31]. The pair-matching problem is automatically avoided since the th column of corresponds to the th column of . 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 , which contains snapshots of the received signals. The TWA is constructed such that the snapshots of the Toeplitz matrices defined earlier are concatenated together as shown in the following:The PARAFAC analysis is then performed on TWA 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 , , , and shown in (16), (17), (18), and (19), respectively, from single snapshot of the received signals.
Step 2. The Toeplitz matrices used to form the TWA given in (24) for one snapshot or in (35) for 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 to obtain the output matrix , which has the same structure as shown in (26).
Step 4. Estimate the diagonal matrices , , and using (30), (31), and (32), respectively. These matrices contain the elevation and azimuth angle information of the sources.
Step 5. Estimate the azimuth angle and elevation angle 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 antenna elements and sources and considering only the major processing operations, the computational complexity of the proposed method is in the order of , where 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 , where denotes the number of snapshots of the received signals. The cumulant-based method in [15] has complexity in the order of , where denotes the number of snapshots of the received signals. The complexity of the universal 2D DOA estimation in [24] is + , where denotes the number of snapshots of the received signals, 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 (extracted from [28]), where 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, , and . 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:Using the signal model above, CRB can be expressed aswhere , , , , , is the th column of , andThe matrix in (37) is given aswhere is as follows:The matrix 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 . 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 are considered here, which are located at and . The number of antenna elements on the positive -axis is . For the purpose of simulation, the number of snapshots considered for the proposed algorithm is , and for the method in [15, 17, 24, 28] they are , , , and , 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.
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 and . The number of antenna elements on the positive -axis is . The curves are plotted for the following number of snapshots: , , , , and . 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.
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 -axis considered here is . The number of snapshots is . The two coherent sources are located at and and a noncoherent source is located at . Figure 7 shows that the estimation performances are almost same for the mixture of coherent and noncoherent sources.
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 -axis is . The number of snapshots used is . Figures 8 and 9 show the scatter plot for three noncoherent sources located at , , and 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 -axis is . The number of snapshots used is . Figure 10 shows the scatter plot for three noncoherent sources located at , , and with SNR = 20 dB. Similarly, Figure 11 shows the scatter plot for three coherent sources located at , , and 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.
In Figures 12 and 13, the performance of the proposed method at low SNR is presented. The number of antenna elements on the positive -axis is . The numbers of snapshots used are . Two coherent sources are located at and , 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.
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 to in steps of . The azimuth angle is fixed at . The number of antenna elements on the positive -axis is fixed at . The number of signal snapshots is . 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%.
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.
Conflicts of Interest
The authors declare that they have no conflicts of interest.