Azimuth / Elevation Directional Finding with Automatic Pair Matching

We addressed the problem of two-dimensional (2D) direction-of-arrival (DOA) elevation and azimuth angles estimation for multiple uncorrelated signals using L-shaped antenna array configuration.The key points of the proposedmethod are the following: (1) it obtains azimuth and elevation angles through construction of three cross-correlation matrices from the collected data of the received signals; this implies that the noise reduces significantly in the reconstructed data matrices; (2) it derives a parallel factor analysis (PARAFAC) model and applies trilinear least squares method to avoid pair matching problem between 2D DOA azimuth and elevation angles for multiple sources; (3) it does not require spectral peak searching; and (4) it has better 2D DOA estimation compared with signal parameters via rotational invariance technique and fourth-order signal parameters via rotational invariance technique. Simulation results demonstrate the estimation accuracy and the effectiveness of the proposed method.


Introduction
The problem of estimating the two-dimensional directionof-arrival (2D DOA) azimuth and elevation (, ) angles of multiple incident sources plays an important role in many practical applications in the fields of wireless communication and signal processing such as radar, sonar, and multipleinput-multiple-output (MIMO) systems.Several 2D DOA estimation methods have been proposed in the literature considering different geometries of the antenna arrays such as circular array [1], parallel uniform linear array [2], uniform rectangular array [3], and L-shaped array [4,5].In particular, the L-shaped array which is composed of two orthogonal uniform linear arrays with one placed on the -axis and the other on the -axis has received considerable attention due to its geometric configuration, ease of implementation, use of less number of antennas, and higher estimation accuracy compared with other geometrical arrays.
Conventional methods of DOA estimation rely on the decomposition of the observation space into a signal subspace and a noise subspace.Two widely used subspace techniques, MUSIC [6] and ESPRIT [7], are computationally complex as they require either eigenvalue decomposition (EVD) of the sample covariance matrix or the singular value decomposition (SVD) of the received data matrix.Improved techniques which are simpler and less complex have been reported in the literature [6][7][8][9][10][11] which do not rely on either EVD or SVD.However, some of these methods [6,7] suffer from heavy losses of the array aperture and encounter an estimation failure problem.To avoid aperture loss problem many cumulant-based methods [12][13][14][15][16][17] have been proposed.But these methods are computationally intensive and require parameter pairing.
The 2D DOA methods using L-shaped array [3,4] fail to resolve the pair matching problem for multiple incident sources which result in incorrect 2D estimation of azimuth and elevation (, ) angles and hence severe performance degradation.Many schemes have been proposed to solve the problem of pair matching [11,[18][19][20][21].These pair matching techniques have high computational cost and complexity and they do not always provide accurate pairing results.
A parallel factor analysis (PARAFAC) [22] model is a method that transfers low rank matrix decomposition to three-way arrays (TWAs); it was first introduced in psychometrics and flow injection analysis.It has been widely used in subspace domain in array signal processing and wireless communication areas [23,24].Several methods have applied PARAFAC model to alleviate the problem of parameter pairing (or pair matching) for 2D DOA estimation for azimuth and elevation angles [25][26][27][28][29].A trilinear decomposition-based blind 2D DOA estimation algorithm employing PARAFAC data model for parallel shaped array has been proposed in [17] to achieve automatic pair matching.One drawback of this method is that it requires high number of snapshots and has high computational complexity.
In this paper, we propose a novel 2D DOA estimation method that employs an L-shaped uniform antenna array based on a new computationally efficient cross-correlation with automatic pair matching based on PARAFAC data model.The constructed data consists of three cross-correlation matrices which contain information about azimuth and elevation angles of multiple uncorrelated signals that are spatially independent; this implies that the noise significantly reduces in the constructed data matrices.Further, we derive a PARAFAC model and apply trilinear least squares method to avoid pair matching problem.Compared with existing methods, the proposed method works with less than a hundred snapshots and has very accurate estimation and has lower computational complexity.
The rest of the paper is structured as follows.The signal model and proposed method are presented and developed in Section 2. Analysis of computational complexity is presented in Section 3. Simulation results are presented in Section 4. Conclusions are drawn in Section 5.

Signal Model and Proposed Method
The geometry of the proposed L-shaped array is shown in Figure 1.The distance between adjacent antenna elements is  where  = /2 with  being the wavelength of the incident waveform.The arrays are divided into four subarrays: x, y, z, and w.Each linear subarray consists of ( − 1) antenna elements.The element located at point (0, 0, 0) is considered as a reference element.Hence, the total number of antenna elements in both -axis and -axis is (2 − 1).Consider  narrowband noncoherent sources in far-field of the antenna array.
( − 1) × 1 collected signal vectors received at x, y, z, and w subarrays are defined as follows: where  represents the snapshot index and superscript  represents the transpose operation.The received ( − 1) × 1 signal vectors in (1) can be represented as follows: where the array manifold vector for th source a  (

𝑇
,   = 2 sin   cos   /.() is the complex envelope vector of  incident sources and at snapshot , n  (), n  (), n  (), and n  () are the Gaussian white noise vectors of zero mean and variance  2 .A  () and A  (, ) are array manifold matrices.
The matrices Φ  () in (3) and Φ  (, ) in ( 5) are  ×  diagonal matrices containing information about the azimuth angle   and elevation angle   which can be presented as follows: In the proposed method, we first construct three crosscorrelation matrices (z subarray, x subarray), (w subarray, x subarray), and (z subarray, y subarray) as follows: where where the superscript  represents the conjugate and transpose operations.Note that {n  (), n  ()} and {n  (), n  ()} are spatially independent of each other.Therefore, where 0 matrix has a dimension of ( − 1) × ( − 1) with all entries zero.Additive noises are not correlated with incident signals.

For uncorrelated sources s(𝑡
, {z(), x()}, {w(), x()} and {z(), y()} are wide-sense stationary sequences.As a consequence, the correlation matrix of the signal sources R  = [s()s  ()] is a diagonal matrix where its entries represent the power of signal sources.The matrices The cross-correlation matrices in (10) are concatenated to form a new data matrix R as follows: Parallel factor (PARAFAC) model is applied on the data in (11) along with trilinear least squares method to jointly estimate the correct pair of azimuth   and elevation   angles for each signal source.
Using the definition of PARAFAC model, the outer product (a ∘ b ∘ c) of the three vectors a ∈ C ×1 , b ∈ C ×1 and c ∈ C ×1 can be expressed in a third-order tensor form as Q ∈ C ×× with typical elements defined as   =       .Q can be expressed as a sum of tensor product: where a  , b  , and c  are th columns of the following load matrices A, B, and C.These matrices for a given PARAFAC model can be defined as follows: The PARAFAC decomposition in (12) can also be represented in another matrix form of 3D tensor matrix Q ∈ C ×× and can be represented using three slice matrices Q (1) , Q (2) , Q (3) as follows: where is the Khatri-Rao product based on column wise Kronecker products and ⊗ is Kronecker product.
The PARAFAC decomposition in ( 14) is considered to be essentially unique to arbitrary permutation and scaling under the condition ( [31,32]): where   ,   ,   denote the maximum number of arbitrary linearly independent columns of matrices A, B, and C, respectively.The arbitrary permutation and scaling implies International Journal of Antennas and Propagation that there exists a triplet matrix ( Ã, B, C) related to (A, B, C) as follows: where Π is a permutation matrix and {Δ 1 , Δ 2 , Δ 3 } are arbitrary diagonal matrices satisfying Δ 1 Δ 2 Δ 3 = I.
On the basis of PARAFAC theorem, three-way array (TWA) ( − 1) × ( − 1) × 3 can be constructed using (11) Let C = A   (); according to the definition of Khatri-Rao product the matrix in (17) can be transformed as follows: where Λ −1 (R  ) represents the row vector data built from the diagonal elements of the diagonal matrix R  .The uniqueness of ( 18), (19), and (20) will be guaranteed if the following inequality holds: For different DOAs and independent sources C and A  have Vandermonde structure with minimum rank equal to the number of sources, and D also is a nonsingular matrix whose rank equals the number of incident sources.This implies that the minimum rank of rank(D) + rank(A) + rank(C) = 3; for multiple incident sources  ≥ 2 and, therefore, (22) will be always guaranteed.
One of the methods of solving PARAFAC model in ( 18), (19), and ( 20) is trilinear alternative least squares (TALS) approach [24,[31][32][33].TALS method can be applied to solve the matrices D, C, and A  and then estimate the azimuth and elevation angles.There are three basic steps behind TALS method: (a) update one of the matrices D, C, and A  each time using least squares (LS) method, (b) continue updating the remaining matrices based on the LS results from the previous estimation step, and (c) repeat previous steps (a) and (b) until convergence of the LS cost function.The detailed procedure of estimating D, C, and A  using TALS is as follows.
Define the cost functions for finding the matrices D, C, and A  as where ‖ ⋅ ‖  stands for Frobenius norm.Given the estimation of matrices Ĉ and Â , the matrix D can be found from (23) as follows: where (J) † represents the pseudoinverse of matrix J.A  can be also obtained by minimizing the cost function in (24) and keeping Ĉ and D fixed.
Similarly, an estimation of matrix Ĉ can be obtained as (1) .
The process in ( 27), (28), and (29) will continue until matrices D, Ĉ, and Â converge.PARAFAC along with TALS method can be initialized to speed up the convergence by exploiting an ESPRIT method idea on the concatenated data formed in (11).
The estimated matrix D in (27) will be sufficient for 2D DOA estimation of azimuth and elevation angles.The TALS method guarantees the convergence but it is slow.For fast implementation of alternative least squares method to solve PARAFAC model in (23), (24), and (25), the COMFAC algorithm is employed which speeds up the least square fitting by utilizing a compressed version of the three-way data into smaller matrix dimensions.COMFAC MATLAB function will be used to estimate the matrices D, Ĉ, and Â (as described in [24]) as follows: where R is the input data,  is the number of sources,  represents the number of iterations, and where φ () is the estimated th entry of the diagonal matrix Φ ().Similarly, th diagonal entry of Φ (, ) can be obtained as follows: The azimuth angle φ and elevation angle θ for th source can be estimated as follows: The estimated azimuth and elevation angles in (33) are for th source.In case of multiple sources, the following pairs are obtained: ( θ1 , φ1 ), ( θ2 , φ2 ), . . ., ( θ , φ ).The estimated matrices D, Ĉ, and Â have the same column permutation matrix.This implies automatic pair matching since th column of the steering matrix A  corresponds to th column of matrix D.
The procedure of the 2D DOA proposed method is summarized as follows.
Step 2. Estimate the cross-correlation matrices R , R , and R from multiple snapshots in (10) as follows: Step 3. Concatenate the estimated cross-correlation matrices { R , R , R } according to (34).
Step 7. Obtain the estimated diagonal matrices φ () and φ () from the identification matrix D.
Step 8. Estimate the 2D DOA azimuth and elevation angles according to (33).

Analysis of the Computational Complexity of the Proposed Method
The computational complexity of the proposed method is compared with that of 2D DOA fourth-order cumulant method [17] and novel 2D DOA with L-shaped array [30].For  total snapshots,  number of antenna elements,  number of iterations, and  number of sources, considering major processing operations like forming the sample covariance or cross-correlation matrices and applying the alternative least squares method, the total computational complexity of the proposed method is in the order of (3( − 1) + (3 3 ) + 9( − 1) 2 ).The complexity of the novel L-shaped array method is in the order of (4(−1)+(3 3 )+12(−1) 2 ) and the complexity of the fourth-order cumulant method is in the order of (21(2 + 1) 2  + (3 3 ) + 12(2 + 1) 2 ).
Upon comparison, it can be seen that the proposed method requires slightly less computational complexity compared to the method in [17] and significantly less computational complexity compared with the method in [30].

Simulation Results
The performance of the proposed method is presented in the section and compared with the novel L-shaped method in [30] and cumulant-based method in [17].The performance is measured in terms of root mean square error (RMSE) for the azimuth and elevation angles estimation.We consider 21 antenna elements in total for the first three experiments.The distance between the adjacent elements is taken to be half the wave length of the incoming signal, and the number of uncorrelated sources is taken as ( = 1,  = 2 and  = 4).Several simulation experiments have been conducted to evaluate the performance of the proposed method.The RMSE for the joint 2D DOA estimation azimuth and elevation angles is defined as follows: where  represents the source index, [  ] represents the expectation value of a random variable   , and ( θ , φ ) are the pair of the estimated elevation and azimuth angles.
In the first experiment, we consider two uncorrelated sources with direction-of-arrival azimuth and elevation angles (, ) = (75 ∘ , 60 ∘ ) and (80 ∘ , 70 ∘ ), SNR range is set from −5 to 30 dB, and the number of snapshots is 200.Monte-Carlo trials of 500 are used.The RMSE values for source 1 and source 2 are shown in Figures 2 and 3 versus SNR for both novel L-shaped and cumulant methods and compared with the proposed method.We observe that the proposed International Journal of Antennas and Propagation Cumulant method for source 1 Novel L-shaped array for source 1 Proposed method for source 1  method has better performance which is indicated through lower RMSE especially at low SNR.For a given RMSE value of 0.15 degrees for source 1 and 0.25 degrees for source 2, it is clear that the proposed method is 5 dB better compared with novel L-shaped for source 1 and around 2.5 dB for source 2. It is also 7.5 and 8 dB better when compared with cumulant method for source 1 and source 2, respectively.We observe also that proposed method has better performance compared Azimuth angle (degrees) Figure 5: Scatter plot for azimuth and elevation for four uncorrelated sources at (65 ∘ , 60 ∘ ), (80 ∘ , 75 ∘ ), (90 ∘ , 80 ∘ ), and (100 ∘ , 90 ∘ ) at SNR = 10 dB using novel L-shaped method [30].
with [17,30] even at low SNR.It can be deduced from Figures 2 and 3 that the performance of the proposed method will be affected due to noise.For example, for the proposed method the RMSE at SNR of −5 dB for source 1 is about 1.3 degrees and for source 2 is about 1.5 degrees.
In the second experiment, four uncorrelated sources with DOAs at (, ) = (65 ∘ , 60 ∘ ), (80 ∘ , 75 ∘ ), (90 ∘ , 80 ∘ ), and (100 ∘ , 90 ∘ ) are considered and the number of snapshots is set to 200.Monte-Carlo trials of 200 are conducted.SNR is set to 10 dB.Figures 4, 5, and 6 illustrate joint azimuth and elevation scatter diagrams of 2D DOA estimation for the proposed method, novel L-shaped, and cumulant method, respectively.It is shown that the four incoming sources can be clearly observed by all methods.The proposed method in Figure 4 gives better and precise estimation compared to the other two methods in [17,30].It is also observed that
the cumulant method in Figure 6 has the worst azimuth and elevation estimation for the four sources.
In the third experiment, the effect of the number of snapshots on the performance of the proposed method is evaluated.We consider two uncorrelated sources with DOAs at (, ) = (45 ∘ , 50 ∘ ) and (65 ∘ , 80 ∘ ), SNR is set to 10 dB, and the number of snapshots  range is set from 100 to 600.Monte-Carlo trials of 500 are used.The average RMSE of azimuth and elevation angles estimation versus the number of snapshots for the two sources is shown in Figures 7 and 8. From these figures, we observe that RMSE of joint azimuth and elevation angles for source 1 and source 2 decrease with increasing number of snapshots.We can also clearly see that the proposed method has higher estimation accuracy compared to the methods in [17,30].
In the fourth experiment, the effect of the number of antennas on the performance of the proposed method is evaluated.We consider a single source located at (, ) = (65 ∘ , 72 ∘ ), SNR set to 6 dB, and the number of snapshots set to 300.Monte-Carlo trials of 400 are used.The average RMSE of azimuth and elevation angles estimation versus the number of antennas is shown in Figure 9. From the figure, we observe International Journal of Antennas and Propagation that RMSE of joint azimuth and elevation angles decrease with increasing number of antennas.We can also clearly see that the proposed method has better performance compared to the methods in [17,30].

Conclusions
We have proposed a new method for 2D DOA azimuth and elevation angles estimation using L-shaped array.The proposed method has lower complexity and better performance compared with existing methods since constructed data matrices from cross-correlation are almost free of noise.PARAFAC model is derived for automatic pair matching of azimuth and elevation angles for multiple incident sources.In addition, the proposed method does not require spectral peak searching.
[D A  C] are the output identification matrices.Now the diagonal matrices Φ  () in (3) and Φ  (, ) are estimated from the identification matrix D = [ d1 d2 d3 ]  as follows: