# Azimuth/Elevation Directional Finding with Automatic Pair Matching

Authors: Nizar Tayem
Journal: International Journal of Antennas and Propagation
 (2016)

---

## Abstract

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 proposed method 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.


---

## Body

## 1. Introduction

The problem of estimating the two-dimensional direction-of-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 multiple-input-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 x-axis and the other on the z-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–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–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–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–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 Section2. Analysis of computational complexity is presented in Section 3. Simulation results are presented in Section 4. Conclusions are drawn in Section 5.


## 2. Signal Model and Proposed Method

The geometry of the proposed L-shaped array is shown in Figure1. The distance between adjacent antenna elements is d where d=λ/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 N-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 z-axis and x-axis is (2N-1). Consider K narrowband noncoherent sources in far-field of the antenna array.Figure 1
Configuration of array geometry for 2D DOA azimuth and elevation angles.( N - 1 ) × 1 collected signal vectors received at x, y, z, and w subarrays are defined as follows:(1)xt=x1tx2t⋯xN-1tT,yt=y2ty3t⋯yNtT,zt=z1tz2t⋯zN-1tT,wt=w2tw3t⋯wNtT,where t represents the snapshot index and superscript T represents the transpose operation. The received (N-1)×1 signal vectors in (1) can be represented as follows:(2)xt=Azθst+nxt,(3)yt=AzθΦzθst+nyt,(4)zt=Axθ,ϕst+nzt,(5)wt=Axθ,ϕΦxθ,ϕst+nwt,(6)Azθ=azθ1azθ2⋯azθK,(7)Axθ,ϕ=axθ1,ϕ1axθ2,ϕ2⋯axθK,ϕK,(8)st=s1ts2t⋯sKtT,where the array manifold vector for kth source azθk=1e-jαk⋯e-jN-1αk, αk=2πdcosθk/λ, and axθk,ϕk=1e-jβk⋯e-jN-1βkT, βk=2πdsinθkcosϕk/λ. s(t) is the complex envelope vector of k incident sources and at snapshot t, nx(t),ny(t),nz(t), and nw(t) are the Gaussian white noise vectors of zero mean and variance σ2. Azθ and Axθ,ϕ are array manifold matrices.The matricesΦzθ in (3) and Φxθ,ϕ in (5) are K×K diagonal matrices containing information about the azimuth angle ϕk and elevation angle θk which can be presented as follows:(9)Φzθ=diage-iα1e-iα2⋯e-iαK,Φxθ,ϕ=diage-iβ1e-iβ2⋯e-iβK.In the proposed method, we first construct three cross-correlation matrices (z subarray, x subarray), (w subarray, x subarray), and (z subarray, y subarray) as follows:(10)Rzx=EztxtH=Axθ,ϕEstsHtAzHθ+EnztnxHt=Axθ,ϕRsAzHθ+Nzx,whereNzx=σ20⋯000⋯0⋮⋮⋱⋮00⋯0∈CN-1×N-1,Rwx=EwtxtH=Axθ,ϕΦxθ,ϕEstsHtAzHθ+EnwtnxHt=Axθ,ϕΦxθ,ϕRsAzHθ+Nwx,Rzy=EztytH=Axθ,ϕEstsHtΦzHθAzHθ+EnztnyHt=Axθ,ϕRsΦzHθAzHθ=Axθ,ϕΦzHθRsAzHθ+Nzy,where the superscript H represents the conjugate and transpose operations. Note that nwt,nx(t) and nwt,nx(t) are spatially independent of each other. Therefore, Nwx=EnwtnxHt=0 and Nzy=EnztnyHt=0, where 0 matrix has a dimension of N-1×N-1 with all entries zero. Additive noises are not correlated with incident signals. For uncorrelated sources st=s1(t)s2(t)⋯sK(t)T, z(t),x(t),w(t),x(t) and z(t),y(t) are wide-sense stationary sequences. As a consequence, the correlation matrix of the signal sources Rs=Es(t)sH(t) is a diagonal matrix where its entries represent the power of signal sources. The matrices RsΦzHθ=ΦzHθRs since both matrices RsandΦzHθ are diagonal.The cross-correlation matrices in (10) are concatenated to form a new data matrix R as follows:(11)R=RzxRzyRwx=Axθ,ϕRsAzHθAxθ,ϕΦzθRsAzHθAxθ,ϕΦxHθ,ϕRsAzHθ+NzxNzyNwx.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 ϕk and elevation θk angles for each signal source.Using the definition of PARAFAC model, the outer producta∘b∘c of the three vectors a∈CI×1, b∈CJ×1 and c∈CK×1 can be expressed in a third-order tensor form as Q∈CI×J×K with typical elements defined as qijk=aibjck. Q can be expressed as a sum of tensor product:(12)Q=∑u=1Uau∘bu∘cu,where au, bu, and cu are uth columns of the following load matrices A, B, and C. These matrices for a given PARAFAC model can be defined as follows:(13)A=a1a2⋯aU,A∈CI×U,B=b1b2⋯bU,B∈CJ×U,C=c1c2⋯cU,C∈CK×U.The PARAFAC decomposition in (12) can also be represented in another matrix form of 3D tensor matrix Q∈CI×J×K and can be represented using three slice matrices Q1, Q2, Q3 as follows:(14)Q1=A⊙BCT,Q2=B⊙CAT,Q3=C⊙ABT,where A⊙B=a1⊗b1a2⊗b2⋯aU⊗bU 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]):(15)κA+κB+κC≥2K+1,where κA, κB, κC denote the maximum number of arbitrary linearly independent columns of matrices A, B, and C, respectively. The arbitrary permutation and scaling implies that there exists a triplet matrix (A~,B~,C~) related to (A,B,C) as follows:(16)A=A~ΠΔ1,B=B~ΠΔ2,C=C~ΠΔ3,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)N-1×N-1×3 can be constructed using (11) as follows: (17)R:,:,1R:,:,2R:,:,3=RzxRzyRwx=Axθ,ϕRsAzHθAxθ,ϕΦzθRsAzHθAxθ,ϕΦxHθ,ϕRsAzHθ+NzxNzyNwx.LetC=AzHθ; according to the definition of Khatri-Rao product the matrix in (17) can be transformed as follows:(18)R1=Ax⊙CTDT,(19)R2=CT⊙DAxT,(20)R3=D⊙AxC,(21)D=d1d2d3=Λ-1RsΛ-1ΦzθRsΛ-1Φxθ,ΦRs,where Λ-1Rs represents the row vector data built from the diagonal elements of the diagonal matrix Rs. The uniqueness of (18), (19), and (20) will be guaranteed if the following inequality holds:(22)rankD+rankA+rankC≥2K+2.For different DOAs and independent sourcesC and Ax 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)=3K; for multiple incident sources K≥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–33]. TALS method can be applied to solve the matrices D,C, and Ax 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 Ax 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 Ax using TALS is as follows.Define the cost functions for finding the matricesD,C, and Ax as(23)F1D,C,Ax;R1=R1-Ax⊙CTDTF,(24)F2D,C,Ax;R2=R2-CT⊙DAxTF,(25)F3D,C,Ax;R3=R3-D⊙AxCF,where ·F stands for Frobenius norm. Given the estimation of matrices C^ and A^x, the matrix D can be found from (23) as follows:(26)DT=argminDR3-Ax⊙CTDTF2,(27)D^=A^x⊙C^†R3T,where J† represents the pseudoinverse of matrix J. Ax can be also obtained by minimizing the cost function in (24) and keeping C^ and D^ fixed.(28)A^x=CT⊙D†R2T.Similarly, an estimation of matrixC^ can be obtained as(29)C^=D^⊙A^x†R1.The process in (27), (28), and (29) will continue until matrices D^, C^, and A^x 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 matrixD^ 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^, C^, and A^x (as described in [24]) as follows:(30)DAxCi=comfacR,K,where R is the input data, K is the number of sources, i represents the number of iterations, and DAxC are the output identification matrices.Now the diagonal matricesΦzθ in (3) and Φxθ,ϕ are estimated from the identification matrix D^=d^1d^2d^3Tas follows:(31)φ^zk=d^2kd^1k,where φ^zk is the estimated kth entry of the diagonal matrix Φ^zθ. Similarly, kth diagonal entry of Φ^xθ,ϕ can be obtained as follows: (32)φ^xk=d^3kd^1k.The azimuth angleϕ^k and elevation angle θ^k for kth source can be estimated as follows: (33)θ^k=cos-1∠φ^zk×λ2πd,ϕ^k=cos-1∠φ^xk×λ2πdsin⁡θ^k.The estimated azimuth and elevation angles in (33) are for kth source. In case of multiple sources, the following pairs are obtained: θ^1,ϕ^1,θ^2,ϕ^2,…,θ^K,ϕ^K. The estimated matrices D^, C^, and A^x have the same column permutation matrix. This implies automatic pair matching since kth column of the steering matrix Ax corresponds to kth column of matrix D^.The procedure of the 2D DOA proposed method is summarized as follows.Step 1.
Constructx,y,z, and w subarrays according to (2), (3), (4), and (5) from multiple snapshots of the array data x(t)y(t)z(t)w(t),  fort=1,…,L.Step 2.
Estimate the cross-correlation matricesR^zx, R^wx, and R^zy from multiple snapshots in (10) as follows:(34)R^zx=1L∑t=1LztxHt,R^wx=1L∑t=1LwtxHt,R^zy=1L∑t=1LztyHt.Step 3.
Concatenate the estimated cross-correlation matricesR^zx,R^wx,R^zy according to (34).Step 4.
Construct the three-way array (TWA)N-1×N-1×3 according to (17).Step 5.
Apply alternative least squares method to estimateD^, C^, and A^x which minimize the cost function in (23), (24), and (25).Step 6.
Repeat Step5 until convergence.Step 7.
Obtain the estimated diagonal matricesφ^zk and φ^xk from the identification matrix D^.Step 8.
Estimate the 2D DOA azimuth and elevation angles according to (33).


## 3. 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 L total snapshots, N number of antenna elements, n number of iterations, and K 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 O3N-1L+n3K3+9N-12K. The complexity of the novel L-shaped array method is in the order of O4N-1L+n3K3+12N-12K and the complexity of the fourth-order cumulant method is in the order of O212N+12L+n3K3+122N+12K. 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].


## 4. 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 (K=1, K=2 and K=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:(35)RMSE=Eθ^i-θi+ϕ^i-ϕi,where i represents the source index, EQ′ represents the expectation value of a random variable Q′, and (θ^i,ϕ^i) 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 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 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.Figure 2
Comparison of RMSE of joint azimuth-elevation versus SNR for source 1 located at(75∘,60∘) using proposed method, novel L-shaped, and cumulant method.Figure 3
Comparison of RMSE of joint azimuth-elevation versus SNR for source 2 located at(80∘,70∘) using proposed method, novel L-shaped, and cumulant method.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.Figure 4
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 proposed method.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].Figure 6
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 cumulant method [17].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 L 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].Figure 7
RMSE of azimuth-elevation versus number of snapshots at SNR = 10 dB for source 1 located at(45∘,50∘) using proposed method, novel L-shaped, and cumulant method.Figure 8
RMSE of azimuth-elevation versus number of snapshots at SNR = 10 dB for source 2 located at(65∘,80∘) using proposed method, novel L-shaped, and cumulant method.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 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].Figure 9
RMSE of azimuth-elevation versus number of antennas at SNR = 6 dB for single source located at(65∘,72∘) using proposed method, novel L-shaped, and cumulant method.


## 5. 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.



---

*Source: 5063450-2016-12-12.xml*