Two-Dimensional AOA Estimation Based on a Constant Modulus Algorithm

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.


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 as -shaped or parallel ULAs).Then, another method was applied to pair up these estimated elevation and azimuth angles [1][2][3][4][5][6][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][9][10][11].As depicted by Figure 1, for the -shaped ULAs, one ULA will be placed on one axis (-axis) and the other will be placed on a perpendicular axis (-axis).
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.

Problem Formulation
In this section, we present the narrowband received signal model that will be utilised for joint 2D AOA estimation.The L-shaped antenna array has two arrays in two perpendicular axes (Figure 1).The radiation pattern for one array is illustrated by Figure 2. Assume  sources signals, { š  ()}  =1 , impinging upon the array.The signal š  () is represented as š  () =   () cos(2  ) = R{  () exp(2  )}, where   is the carrier frequency, and   () = ∑    ()( − ), where  is an integer which represents the time index, and   () is the low-pass equivalent of š  (),  is the symbol period,   () =     (), where   () ∈ (−1, +1) which represents the symbol parity for a BPSK signal, and   is a positive constant which represents the amplitude of   (), and () is the sampled (raised cosine) pulse shaping function with where  is an integer so that   () =   ().
Since we are looking at 2D AOA estimation, the th source signal arriving at the array will have an elevation angle,   , and an azimuth angle,   .Because of the -shaped antenna orientation shown in Figure 1, the angle   is taken between the -axis and the source signal arrival direction, and the angle   is taken between the -axis and the source signal arrival direction.We use the superscript notation {̂} to denote the estimated value of a variable (for example, ĥ is the estimated value of ℎ).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  = 1 to  = , that are located in the -axis direction, and let us call it r  ().This is represented as where with  being the signal wavelength and  being the distance between consecutive antenna elements.In addition, n  () is the noise vector added to the received signal at the antenna elements (from  = 1 → ) located along the -axis and it is additive white Gaussian noise (AWGN) with a covariance matrix of  2 I × , where I × is the  ×  identity matrix.Finally, (⋅)  represents the transpose operator.Similarly, the sampled received signal at the antenna elements, from  = 1 to  = , that are located in the -axis direction, is given the notation r  () and is represented as where b( Likewise, n  () is the AWGN noise vector added to the received signal at the antenna elements (from  = 1 → ) located along the -axis, and it also has a covariance matrix of  2 I × .So, both ×1 received signal vectors can be written in matrix form as where From now on we will drop the index  from r  () and r  () for simplicity.
Now, we formulate a cross-correlation matrix between the signal received at the array on -axis and the signal received at the array on -axis; that is, we consider the following crosscorrelation matrix: where with   being the power of the th 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 data X = [x 1 x 2 ⋅ ⋅ ⋅ x  ] of size  × , where X = QS (10) with 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 Ŝ or alternatively find a beamformer matrix W, where Ŝ = W  X 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 = (1/) ∑  x  x   .Then, take its eigenvalue decomposition (EVD); that is, R = ÛΣ 2 Û .Set the prewhitening matrix as F = Û Σ−1  , where Σ is a diagonal matrix which contains the square root of the  dominant eigenvalues and Û are their corresponding eigenvectors.Also,  represents the number of sources.Now, prefilter the data, x  = F  x  .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  ×  to  × .Recall that  =  is the number of sources which is assumed to be smaller than  (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 = F  W.
Next, form the matrix where ⊗ is the Kronecker product and {⋅} is the complex conjugate operator, so, for example, ℎ is the complex conjugate of ℎ and vec (H) is the process of stacking the columns of some matrix H into a vector.It is proven in [12] that the null space of Ĉ spans the same subspace as the columns of the beamformer matrix T. So, we compute the null space of Ĉ 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 (‖t  ‖) 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 where ∘ is the Katri-Rao product.
To satisfy both conditions (the basis of the null space of Ĉ and the structured matrix) a joint diagonalisation process is performed which is transformed into the following iterative process [12].
Start with initializing the matrix T by setting it to T = I, where I is the identity matrix.Then, apply the following steps: Then, recondition the matrix M by where where   , u  , and k  are the th singular value, left singular vector, and right singular vector of M. Then project V on M  , as follows: Now, consider the th column of Y and name it y  and apply the inverse of the vectorisation process; that is, Y  = vec −1 (y  ).Then, take the SVD of Y  and consider the left singular vector corresponding to the largest singular value; that is, let and set where t  is the th 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 estimating t  , for all , thus estimating beamformer matrix in the whitened domain (T), we can set the matrix W to be which is the beamformer matrix in the original (nonwhitened) domain; then S is estimated as After presenting the ACMA method, we can apply it to estimate the 2D AOA by applying it to R  .

2D
ACMA Based AOA Estimator.Looking at the definition of the cross-correlation matrix, R  , in (8), we can see that B() is of constant modulus elements.So, we can apply the ACMA method to factorise R  (similar to X in (10) but with size × instead of ×) into the constant modulus matrix U  B()  (similar to S in (10)) and Ψ = Â() P U (similar to Q in (10)), where U is a diagonal matrix which causes the phase ambiguity in the factorised matrices.
Notice that the elements of B() are all of constant modulus with its size of  ×  with U a diagonal matrix which will not affect the relative phases of each column of B() or Â().Also, notice that the matrix P will not affect the complex elements angles in Â() (because P is a real diagonal matrix).Thus, the elevation angles information in Â() is reserved in Ψ and can be deduced by with where ψ (  ) is the th element of the th 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 with Thus, both elevation and azimuth angles are estimated for each source using this 2D ACMA based algorithm.It is worth mentioning that since matrix A() is also constant modulus, then we can apply the ACMA algorithm on the cross-correlation matrix R  (instead of R  ) which is given by Next, we will provide simulation results for the proposed algorithm.

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 4 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 5 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.
Another simulation is carried out to show the span of AOAs within which the 2D ACMA based algorithm can work.So, Figure 6 illustrates the RMSE for the elevation AOA versus different elevation AOAs for a single source. and

Figure 1 :
Figure 1: Layout for the L-shaped antenna array.

Figure 2 :
Figure 2: Radiation pattern for 1D linear antenna array presenting 10 nulls plotted in dB versus AOA in degrees, with 20 elements and interspacing  = 0.12, with a side lobe level of −13 dB.

Figure 3 :
Figure 3: Root mean square error in degrees using the 2D ACMA based method versus SNR in dB.