# Two-Dimensional AOA Estimation Based on a Constant Modulus Algorithm

Authors: Saleh O. Al-Jazzar; Zoubir Hamici; Sami Aldalahmeh
Journal: International Journal of Antennas and Propagation
 (2017)

---

## Abstract

We propose a two-dimensional (2D) angle of arrival (AOA) estimator using the algebraic constant modulus algorithm (ACMA). This algorithm was originally introduced to estimate the one-dimensional (1D) AOA. An extension to estimate and automatically pair the elevation and azimuth angles for different sources is derived and proved in this paper. The ACMA method factorises a matrix into two different matrices; one is of constant modulus and contains the azimuth AOA information; however the second was previously ignored. In this paper we will prove that this second matrix contains the elevation AOA information. Thus, 2D AOA estimation is proved possible using the ACMA method. Simulation results are presented to illustrate the proposed method’s performances under different conditions.


---

## Body

## 1. Introduction

Direction of arrival denotes the direction from which a propagating wave is impinging on an antenna array, while the angle of arrival estimates this direction by measuring the difference of phases or time delays at the individual elements of the array using the wavenumber. Angle of arrival (AOA) estimation has gained great attention due to its importance in enhancing the received signal quality and/or locating wireless transmitters. Many estimators were developed in the literature to estimate the AOAs. Initially, these methods were developed to estimate one-dimensional AOAs. Different antenna array structures were assumed for the AOA estimation such as the uniform linear antenna array (ULA) and circular antenna array structures.In many practical applications it is required to estimate both the elevation and azimuth AOAs. So, many methods were developed initially to estimate the elevation and azimuth angles separately using two ULAs (such asL-shaped or parallel ULAs). Then, another method was applied to pair up these estimated elevation and azimuth angles [1–7]. Other methods were also developed in which the elevation and azimuth angles were estimated jointly without the need for any pairing technique; these methods require the use of two ULAs [1, 8–11]. As depicted by Figure 1, for the L-shaped ULAs, one ULA will be placed on one axis (x-axis) and the other will be placed on a perpendicular axis (z-axis).Figure 1
Layout for theL-shaped antenna array.In this paper we present an elevation and azimuth AOA estimator based on the algebraic constant modulus algorithm (ACMA) presented in [12]. The process is performed blindly at the receiver without the need for any training. The ACMA method, upon which the 2D AOA estimator is based, performs factorisation on a received signal block to retrieve the constant modulus (CM) part contained in this block. This is valid under the assumption that the original transmitted signal has a CM feature.A cross-correlation matrix is initially formed from the antenna model; then the ACMA algorithm is applied on it to get two matrices; one is CM and contains the azimuth AOAs information for different sources and the second, previously ignored, holds their corresponding elevation AOAs. Thus, a 2D AOA estimator can be achieved by utilising both matrices of the ACMA factorisation. Also, automatic pairing between different elevation and azimuth AOAs is achieved for different sources. The proposed estimator does not require any training, so it is considered a blind estimator. The proposed estimator is named the 2D ACMA based algorithm.


## 2. Problem Formulation

In this section, we present the narrowband received signal model that will be utilised for joint 2D AOA estimation. TheL-shaped antenna array has two arrays in two perpendicular axes (Figure 1). The radiation pattern for one array is illustrated by Figure 2. Assume K sources signals, {sˇk(t)}k=1K, impinging upon the array. The signal sˇk(t) is represented as sˇk(t)=sktcos⁡2πfct=R{sktexp⁡j2πfct}, where fc is the carrier frequency, and sk(t)=∑iαk(i)g(t-iT), where i is an integer which represents the time index, and sk(t) is the low-pass equivalent of sˇk(t), T is the symbol period, αk(i)=Akβk(i), where βk(i)∈(-1,+1) which represents the symbol parity for a BPSK signal, and Ak is a positive constant which represents the amplitude of sk(t), and g(t) is the sampled (raised cosine) pulse shaping function with(1)gnT=1,ifn=00,otherwise,where n is an integer so that sk(nT)=αk(n).Figure 2
Radiation pattern for 1D linear antenna array presenting 10 nulls plotted in dB versus AOA in degrees, with 20 elements and interspacingd=0.12, with a side lobe level of −13 dB.Since we are looking at 2D AOA estimation, thekth source signal arriving at the array will have an elevation angle, θk, and an azimuth angle, ϕk. Because of the L-shaped antenna orientation shown in Figure 1, the angle ϕk is taken between the x-axis and the source signal arrival direction, and the angle θk is taken between the z-axis and the source signal arrival direction. We use the superscript notation ^ to denote the estimated value of a variable (for example, h^ is the estimated value of h). So, to start developing the received signal model, let us first consider the sampled received signal at the antenna elements (after the matched filter stage), from m=1 to m=M, that are located in the z-axis direction, and let us call it rz(i). This is represented as(2)rzi=∑k=1Kαkiaθk+nzi,where a(θk)=a1(θk)⋯aM(θk)T and (3)amθk=exp⁡j2πλd×m-1cos⁡θkwith λ being the signal wavelength and d being the distance between consecutive antenna elements. In addition, nz(i) is the noise vector added to the received signal at the antenna elements (from m=1→M) located along the z-axis and it is additive white Gaussian noise (AWGN) with a covariance matrix of σ2IM×M, where IM×M is the M×M identity matrix. Finally, ·T represents the transpose operator.Similarly, the sampled received signal at the antenna elements, fromm=1 to m=M, that are located in the x-axis direction, is given the notation rx(i) and is represented as(4)rxi=∑k=1Kαkbϕk+nxi,where bϕk=b1(ϕk)⋯bM(ϕk)T and (5)bmϕk=exp⁡j2πλd×m-1cos⁡ϕk.Likewise, nx(i) is the AWGN noise vector added to the received signal at the antenna elements (from m=1→M) located along the x-axis, and it also has a covariance matrix of σ2IM×M. So, both M×1 received signal vectors can be written in matrix form as (6)rzi=Aθαi+nzi,rxi=Bϕαi+nxi,where (7)Aθ=aθ1⋯aθK,Bϕ=bϕ1⋯bϕK,αi=α1i⋯αKi.From now on we will drop the index i from rz(i) and rx(i) for simplicity.Now, we formulate a cross-correlation matrix between the signal received at the array onz-axis and the signal received at the array on x-axis; that is, we consider the following cross-correlation matrix:(8)Rzx=ErzrxH,Rzx=AθPsBϕH,where (9)Ps≔EααH=diag⁡p1⋯pk⋯pKwith pk being the power of the kth source. This formation of the cross-correlation matrix is the matrix upon which the ACMA algorithm will be applied. Next, we give a description on the ACMA method which will be used in the 2D AOA estimation.


### 2.1. The ACMA Algorithm

Given a block of dataX=x1x2⋯xN of size M×N, where(10)X=QSwith S being the matrix composed of constant modulus elements and Q being some full rank matrix. Then, the aim of the ACMA method is to factorise X to get back Q^ and S^ or alternatively find a beamformer matrix W, where S^=WHX with some phase ambiguity.The ACMA described in [12] starts by a prewhitening step (which is not considered in the original ACMA method proposed in [13]). To do so, we compute R^x=1/N∑nxnxnH. Then, take its eigenvalue decomposition (EVD); that is, R^x=U^Σ^2U^H. Set the prewhitening matrix as F=U^sΣ^s-1, where Σ^s is a diagonal matrix which contains the square root of the D dominant eigenvalues and U^s are their corresponding eigenvectors. Also, D represents the number of sources. Now, prefilter the data, x_n=FHxn. A whitening transformation or Principal Components Analysis (PCA) is a dimensionality reduction algorithm that transforms a vector of random variables with a known correlation matrix into a set of new variables whose correlation is the identity matrix which means that such signals are uncorrelated. This prewhitening step will reduce the matrix size over which the ACMA algorithm is performed from M×N to D×N. Recall that D=K is the number of sources which is assumed to be smaller than M (the number of antenna elements). Now, we will define the matrix T which is similar to the beamformer matrix W but in the whitened domain; that is, T=FHW.Next, form the matrix(11)C_^=1N∑nx̅-n⊗x_nx̅-n⊗x_nH-vec⁡R_^xR_^xH,where ⊗ is the Kronecker product and {·¯} is the complex conjugate operator, so, for example, h¯ is the complex conjugate of h and vecH is the process of stacking the columns of some matrix H into a vector. It is proven in [12] that the null space of C_^ spans the same subspace as the columns of the beamformer matrix T. So, we compute the null space of C_^ which we will denote as V to help find T. The source separation beamformer matrix T, in the whitened domain constraint, will be found by an iterative process of joint diagonalisation such that the solutions of the beamformer columns elements are scaled to unit-norm vectors (ti) as a decorrelation step defined by the properties of the constant modulus algorithm in the presence of noise [12]. So the solution can be written in the form (12)T=argminti,Mis  invertableVM-T¯∘T,where ∘ is the Katri-Rao product.To satisfy both conditions (the basis of the null space ofC_^ and the structured matrix) a joint diagonalisation process is performed which is transformed into the following iterative process [12].Start with initializing the matrixT by setting it to T=I, where I is the identity matrix. Then, apply the following steps: (13)M=VHT¯∘T.Then, recondition the matrix M by (14)M′=recondM≔∑jσj′ujvjH,whereσj′=1,σj<0.5σj,σj≥0.5,where σj, uj, and vj are the jth singular value, left singular vector, and right singular vector of M. Then project V on M′, as follows: (15)Y=VM′.Now, consider theith column of Y and name it yi and apply the inverse of the vectorisation process; that is, Yi=vec-1yi. Then, take the SVD of Yi and consider the left singular vector corresponding to the largest singular value; that is, let (16)Yi=ΣjγjpjqjHand set (17)ti=p1,where ti is the ith column vector of T. The process of joint diagonalisation is repeated until convergence. The described joint diagonalisation process is closely related to the (canonical) polyadic decomposition which is used for factor retrieval, dimensionality reduction, and signal denoising to name but a few. This multiway array (tensor) decomposition is able to take into account, for example, spatial, temporal, and spectral information simultaneously and retrieve factors under mild conditions [14, 15].After estimatingti, for all i, thus estimating beamformer matrix in the whitened domain (T), we can set the matrix W to be (18)W=FTwhich is the beamformer matrix in the original (nonwhitened) domain; then S is estimated as (19)S^=WHX.After presenting the ACMA method, we can apply it to estimate the 2D AOA by applying it to Rzx.


### 2.2. 2D ACMA Based AOA Estimator

Looking at the definition of the cross-correlation matrix,Rzx, in (8), we can see that B(ϕ) is of constant modulus elements. So, we can apply the ACMA method to factorise Rzx (similar to X in (10) but with size M×M instead of M×N) into the constant modulus matrix UHB^(ϕ)H (similar to S in (10)) and Ψ^=A^(θ)P^sU (similar to Q in (10)), where U is a diagonal matrix which causes the phase ambiguity in the factorised matrices.Notice that the elements ofB^(ϕ) are all of constant modulus with its size of K×M with U a diagonal matrix which will not affect the relative phases of each column of B^(ϕ) or A^(θ). Also, notice that the matrix P^s will not affect the complex elements angles in A^(θ) (because P^s is a real diagonal matrix). Thus, the elevation angles information in A^(θ) is reserved in Ψ^ and can be deduced by (20)θ^k=cos-1⁡arg⁡ψ^avg2πd/λwith (21)ψ^avg=∑m=1M-1ψ^m+1θk/ψ^mθkM-1,where ψ^m(θk) is the kth element of the mth column of Ψ^ and ψ^avg is the average of all relative phases, not affected by phase ambiguity. Similarly, the azimuth angles can be deduced from B^(ϕ) by (22)ϕ^k=cos-1⁡arg⁡B^avg2πd/λwith (23)B^avg=∑m=1M-1b^m+1ϕk/b^mϕkM-1.Thus, both elevation and azimuth angles are estimated for each source using this 2D ACMA based algorithm. It is worth mentioning that since matrixA(θ) is also constant modulus, then we can apply the ACMA algorithm on the cross-correlation matrix Rxz (instead of Rzx) which is given by (24)Rxz=RzxH=BϕPsAθH.Next, we will provide simulation results for the proposed algorithm.


## 2.1. The ACMA Algorithm

Given a block of dataX=x1x2⋯xN of size M×N, where(10)X=QSwith S being the matrix composed of constant modulus elements and Q being some full rank matrix. Then, the aim of the ACMA method is to factorise X to get back Q^ and S^ or alternatively find a beamformer matrix W, where S^=WHX with some phase ambiguity.The ACMA described in [12] starts by a prewhitening step (which is not considered in the original ACMA method proposed in [13]). To do so, we compute R^x=1/N∑nxnxnH. Then, take its eigenvalue decomposition (EVD); that is, R^x=U^Σ^2U^H. Set the prewhitening matrix as F=U^sΣ^s-1, where Σ^s is a diagonal matrix which contains the square root of the D dominant eigenvalues and U^s are their corresponding eigenvectors. Also, D represents the number of sources. Now, prefilter the data, x_n=FHxn. A whitening transformation or Principal Components Analysis (PCA) is a dimensionality reduction algorithm that transforms a vector of random variables with a known correlation matrix into a set of new variables whose correlation is the identity matrix which means that such signals are uncorrelated. This prewhitening step will reduce the matrix size over which the ACMA algorithm is performed from M×N to D×N. Recall that D=K is the number of sources which is assumed to be smaller than M (the number of antenna elements). Now, we will define the matrix T which is similar to the beamformer matrix W but in the whitened domain; that is, T=FHW.Next, form the matrix(11)C_^=1N∑nx̅-n⊗x_nx̅-n⊗x_nH-vec⁡R_^xR_^xH,where ⊗ is the Kronecker product and {·¯} is the complex conjugate operator, so, for example, h¯ is the complex conjugate of h and vecH is the process of stacking the columns of some matrix H into a vector. It is proven in [12] that the null space of C_^ spans the same subspace as the columns of the beamformer matrix T. So, we compute the null space of C_^ which we will denote as V to help find T. The source separation beamformer matrix T, in the whitened domain constraint, will be found by an iterative process of joint diagonalisation such that the solutions of the beamformer columns elements are scaled to unit-norm vectors (ti) as a decorrelation step defined by the properties of the constant modulus algorithm in the presence of noise [12]. So the solution can be written in the form (12)T=argminti,Mis  invertableVM-T¯∘T,where ∘ is the Katri-Rao product.To satisfy both conditions (the basis of the null space ofC_^ and the structured matrix) a joint diagonalisation process is performed which is transformed into the following iterative process [12].Start with initializing the matrixT by setting it to T=I, where I is the identity matrix. Then, apply the following steps: (13)M=VHT¯∘T.Then, recondition the matrix M by (14)M′=recondM≔∑jσj′ujvjH,whereσj′=1,σj<0.5σj,σj≥0.5,where σj, uj, and vj are the jth singular value, left singular vector, and right singular vector of M. Then project V on M′, as follows: (15)Y=VM′.Now, consider theith column of Y and name it yi and apply the inverse of the vectorisation process; that is, Yi=vec-1yi. Then, take the SVD of Yi and consider the left singular vector corresponding to the largest singular value; that is, let (16)Yi=ΣjγjpjqjHand set (17)ti=p1,where ti is the ith column vector of T. The process of joint diagonalisation is repeated until convergence. The described joint diagonalisation process is closely related to the (canonical) polyadic decomposition which is used for factor retrieval, dimensionality reduction, and signal denoising to name but a few. This multiway array (tensor) decomposition is able to take into account, for example, spatial, temporal, and spectral information simultaneously and retrieve factors under mild conditions [14, 15].After estimatingti, for all i, thus estimating beamformer matrix in the whitened domain (T), we can set the matrix W to be (18)W=FTwhich is the beamformer matrix in the original (nonwhitened) domain; then S is estimated as (19)S^=WHX.After presenting the ACMA method, we can apply it to estimate the 2D AOA by applying it to Rzx.


## 2.2. 2D ACMA Based AOA Estimator

Looking at the definition of the cross-correlation matrix,Rzx, in (8), we can see that B(ϕ) is of constant modulus elements. So, we can apply the ACMA method to factorise Rzx (similar to X in (10) but with size M×M instead of M×N) into the constant modulus matrix UHB^(ϕ)H (similar to S in (10)) and Ψ^=A^(θ)P^sU (similar to Q in (10)), where U is a diagonal matrix which causes the phase ambiguity in the factorised matrices.Notice that the elements ofB^(ϕ) are all of constant modulus with its size of K×M with U a diagonal matrix which will not affect the relative phases of each column of B^(ϕ) or A^(θ). Also, notice that the matrix P^s will not affect the complex elements angles in A^(θ) (because P^s is a real diagonal matrix). Thus, the elevation angles information in A^(θ) is reserved in Ψ^ and can be deduced by (20)θ^k=cos-1⁡arg⁡ψ^avg2πd/λwith (21)ψ^avg=∑m=1M-1ψ^m+1θk/ψ^mθkM-1,where ψ^m(θk) is the kth element of the mth column of Ψ^ and ψ^avg is the average of all relative phases, not affected by phase ambiguity. Similarly, the azimuth angles can be deduced from B^(ϕ) by (22)ϕ^k=cos-1⁡arg⁡B^avg2πd/λwith (23)B^avg=∑m=1M-1b^m+1ϕk/b^mϕkM-1.Thus, both elevation and azimuth angles are estimated for each source using this 2D ACMA based algorithm. It is worth mentioning that since matrixA(θ) is also constant modulus, then we can apply the ACMA algorithm on the cross-correlation matrix Rxz (instead of Rzx) which is given by (24)Rxz=RzxH=BϕPsAθH.Next, we will provide simulation results for the proposed algorithm.


## 3. Simulation Results

Simulations are performed with MATLAB on a dual core computer, i5 processor clocked at 2.5 GHz with a RAM of 6 GB. Results are presented for 2D AOA estimation using a 2D ACMA based algorithm. The number of sources is two with elevation anglesθ1,2=75°,45° and azimuth angles ϕ1,2=80°,50°. The number of ensembles over which the simulation results are averaged is 1000 ensembles. The near far ratio (NFR) which is the ratio of power between both sources is set at −20 dB. Figure 3 shows the root mean square error (RMSE) for the elevation and azimuth AOAs estimation using the 2D ACMA based method versus different signal to noise ratios (SNRs). The number of antenna elements and number of ensembles over which the cross-correlation matrix is averaged are variable as shown in the figure. The results show that the performance is enhanced as the number of antenna elements, ensembles number, and SNR increase.Figure 3
Root mean square error in degrees using the 2D ACMA based method versus SNR in dB.Figure4 shows the RMSE for the elevation and azimuth AOAs estimation using the 2D ACMA based method versus different elevation angle difference between the two sources. The elevation angle for the second source is θ2=45° and the first source elevation angle is computed by adding the difference to θ2. The results show that the performance is enhanced as the elevation angle difference increases.Figure 4
Root mean square error in degrees using the 2D ACMA based method versus AOA difference between two sources in degrees.Figure5 shows the RMSE for the elevation and azimuth AOAs estimation using the 2D ACMA based method versus different NFR ratios between both sources. The results show that the performance is degraded as the NFR ratio increases.Figure 5
Root mean square error in degrees using the 2D ACMA based method versus NFR in dB.Another simulation is carried out to show the span of AOAs within which the 2D ACMA based algorithm can work. So, Figure6 illustrates the RMSE for the elevation AOA versus different elevation AOAs for a single source. M and SNR are set to 10 and 10 dB, respectively. The 2D ACMA based algorithm has a high accuracy in the range 20° to 160° with an error less than 0.1° for 1000 ensembles (lower curve), while it achieves an error of 0.3° within a range between 35° and 145° for 100 ensembles (upper curve). For both results, the azimuth is fixed at 80°. Meanwhile, Figure 7 shows the span behaviour for two values of the antenna size M. Lower curve is for M=20, while upper curve is for M=10. Clearly, the 2D ACMA based performance is enhanced as the size of antenna M increases but the span over which the AOAs are estimated decreases. This latter behaviour illustrates the compromise between range and accuracy which is expected. The oscillations present in both figures, Figures 6 and 7, are due to the radiation pattern of antenna arrays shown in Figure 2 which is oscillatory by nature.Figure 6
Root mean square error in degrees using the 2D ACMA based method versus AOA span in degrees.Figure 7
Root mean square error in degrees using the 2D ACMA based method versus AOA span in degrees.Figure8 shows the locus of 10 runs of the 2D ACMA based algorithm where the elevation and azimuth AOAs true values and their corresponding estimations are plotted in the same figure. M, NFR, and SNR are set to 10, −5 dB, and 10 dB, respectively, which shows the high accuracy of the algorithm.Figure 8
Locus of the true and estimated elevation and azimuth AOAs in degrees for 10 different runs.In terms of processing time performance, the time that takes the 2D ACMA based to estimate two pairs of the elevation and azimuth AOAs for two sources, evaluated with MATLAB, is about 8 milliseconds for a number of ensembles equal to 100.


## 4. Conclusion

In this paper a 2D AOA estimator based on the ACMA method is presented. The 2D extended CM estimator utilises theL-shaped structure to form the 2D antenna array. The model of the antenna array is represented by a cross-correlation matrix through which the elevation and azimuth angles are estimated simultaneously for different sources. The estimator does not require any training nor it requires any pairing process to combine the elevation and azimuth angles for each source. Simulation results show the effectiveness of the new 2D ACMA based algorithm applied on an L-shaped uniform linear arrays. A high accuracy multisource AOAs estimation of 0.1° is achieved within a wide range of 140° for both perpendicular antenna arrays.



---

*Source: 3214021-2017-01-18.xml*