# Angle of Arrival Estimation Using Cholesky Decomposition

Authors: Saleh O. Al-Jazzar
Journal: International Journal of Antennas and Propagation
 (2012)

---

## Abstract

An angle of arrival (AOA) estimator is presented. Many applications require accurate
AOA estimates such as wireless positioning and signal enhancement using space-processing techniques.
The proposed AOA estimator depends on the Cholesky decomposition of the received signal autocorrelation
matrix. The resultant decomposed matrices are used to modify the crosscorrelation matrix
of the received signals at the antenna array doublets. The proposed method is named the Cholesky-decomposition-based-AOA (CDBA) estimator. In comparison with the TLS-ESPRIT algorithm which
utilizes the eigenvalue decomposition (EVD) of the received signal autocorrelation and crosscorrelation
matrices, the CDBA method has better performance than the TLS-ESPRIT algorithm especially in low
signal-to-noise-ratio (SNR) cases. Simulations for the proposed CDBA method are shown to assess its
performance.


---

## Body

## 1. Introduction

Accurate estimation of the received signal angle of arrival (AOA) can be very beneficial for signal reception enhancement [1, 2] and wireless positioning [3, 4]. If the AOA estimates were not estimated accurately, then wireless devices would not be located accurately or the received signal bit error rate would be high. Thus, AOA estimators are developed and modified in such cases to increase their accuracy.Many AOA estimators were proposed in the literature [5–9]. The authors in [5] present an iterative method to estimate the AOA. In [6], the authors propose estimating the AOA using the frequency domain of the received signals which requires taking their Fourier transform. The authors in [7] introduce an iterative AOA estimator which depends on the generalized expectation maximization algorithm. In [8], the authors present an AOA estimator for circular arrays. Another very popular AOA estimator is the TLS-ESPRIT algorithm proposed in [9] which utilizes subspace techniques to estimate the AOA.In this paper, an AOA estimator is proposed which enhances the AOA estimation. The proposed method is named the Cholesky-decomposition-based-AOA (CDBA) estimator. The CDBA method utilizes the received signal at the two sides of the antenna array doublets. An autocorrelation matrix of the received signal at one side of the antenna array doublets is calculated. Then, the received signal autocorrelation matrix is decomposed using the Cholesky decomposition. A crosscorrelation matrix between the received signals at both sides of the antenna doublets is calculated. The decomposed matrices of the autocorrelation matrix together with the crosscorrelation matrix are used to form a new matrix from which the AOAs are estimated.The CDBA method does not require taking the Fourier transform of the received signals nor it requires any iterative or searching procedure to estimate the AOA. Also, in comparison with the well-known TLS-ESPRIT algorithm of [9], the accuracy of the CDBA method is higher than that of the TLS-ESPRIT method especially in low signal-to-noise-ratio (SNR) cases.The paper is organized as follows: Section2 introduces the system model that forms the foundation for the proposed estimator. The proposed CDBA estimator is presented in Section 3. Section 4 presents the simulated performance of the CDBA estimator. Finally, conclusions are shown in Section 5.


## 2. System Model

This section presents the narrowband received signal model that will be utilized for the AOA estimation. The antenna array is formed fromM uniform antenna doublets (i.e., 2M total antenna elements). Each antenna doublet is formed of two antenna elements spaced by a distance dd. We assume a K BPSK sources signals, {s~k(t)}k=1K, impinging upon the antenna array. The signal s~k(t) is represented as s~k(t)=sk(t)cos(2πfct)=ℜ{sk(t)exp(j2πfct)}, where fc is the carrier frequency sk(t)=∑i‍αk(i)g(t-iT), i is an integer which represents the time index and sk(t) is the complex low-pass equivalent of s~k(t), T is the symbol period, αk(i)=Akβk(i), where βk(i) is discrete βk(i)∈{-1,+1} which represents the symbol parity and Ak is a positive constant which represents the amplitude of sk(t), and g(t) is the (raised cosine) pulse-shaping function with
(1)g(nT)={1ifn=00otherwise,
where n is an integer so that sk(nT)=αk(n).The AOA is taken between the antenna array axis, and the source signal arrival direction and will be given the notation(θk). In this paper, the superscript notation {^} is used 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 first side of the antenna doublets, and let us call it r1(i), and is represented as
(2)r1(i)=∑k=1K‍αk(i)a(θk)+n1(i),
where
(3)a(θk)=[a1(θk)⋯aM(θk)]T,am(θk)=exp{j2πλd×(m-1)cos(θk)}
where λ is the signal wavelength, and d is the distance between consecutive antenna elements. In addition, n1(i) is the noise vector added to the received signal at the antenna elements (from m=1→M) located along the first side of the antenna doublets, which is additive white Gaussian noise (AWGN) and has a covariance matrix of σ2IM×M, where IM×M is the M×M identity matrix. Finally, (·)T represents the transpose operation.The(M×1) received signal vector at the first side of the antenna doublets set can be written in matrix form as
(4)r1(i)=A(θ)α(i)+n1(i),
where
(5)A(θ)=[a(θ1)⋯a(θK)],θ=[θ1⋯θK]T,α(i)=[α1(i)⋯αK(i)]T.
From now on the term (i) will be dropped from all terms for simplicity.The received signal vector at the second side of the antenna doublets set will be given the notationr2. Thus,
(6)r2=A(θ)Zα+n2,
where n2 is an AWGN vector at the second side of the antenna array doublets, and
(7)Z=diag⁡([z1⋯zk⋯zK]),
with
(8)zk=exp⁡{j2πλdd×cos⁡(θk)}.
Next section will describe the proposed CDBA estimator.


## 3. Proposed Cholesky-Decomposition-Based AOA (CDBA) Estimator

To implement the CDBA method, the crosscorrelation matrix (R21) between r2 and r1 is formulated as follows:
(9)R21=E[r2r1H]=A(θ)PsZA(θ)H,
where
(10)Ps:=E[ααH]=diag⁡([p1⋯pk⋯pK]),
where pk is the power of the kth source.Another matrix considered by the CDBA algorithm is the autocorrelation matrix ofr1, as follows:
(11)R~11=E[r1r1H]=A(θ)PsA(θ)H+σ2IM×M.
Take the noise free version of R~11, and let us call it R11, where R11=R~11-σ2IM×M, that is,
(12)R11=A(θ)PsA(θ)H.
Then, the eigenvalue decomposition (EVD) of R11 is
(13)R11=VΓVH,
where V is the eigenvector matrix of R11, and Γ is an M×M diagonal matrix with its diagonal that contains the eigenvalues of R11, that is,
(14)Γ=diag⁡([γ1⋯γK0⋯0]),
where the K nonzero values of the diagonal elements of Γ correspond to the K sources which will be called γ1→K≡γ1→γK. Each element of γ1→K corresponds to one of the K sources.Comparing (12) with (13), then we have
(15)A(θ)PsA(θ)H=VΓVH.
Comparing both sides of (15), then we have
(16)A(θ)=VT,
where T is the appropriate matrix to change the basis.Substituting (16) into (15), then we have
(17)Γ=TPsTH.
The matrix Γ is an M×M matrix, and Ps is a K×K diagonal matrix. Also, the nonzero diagonal elements of Γ are only the first K diagonal elements. So, (17) can be written as
(18)[γ1⋯00⋯0⋮⋱⋮⋮⋱⋮0⋯γK0⋯00⋯00⋯0⋮⋱⋮⋮⋱⋮0⋯00⋯0]M×M=T[p1⋯0⋮⋱⋮0⋯pK]K×KTH.
Looking at (18), it can be deduced that the matrix T should be M×K. Also, the lower (M-K)×K part of T should be all zero elements, that is,
(19)T=[t1,1⋯t1,K⋮⋮tK,1⋯tK,K0⋯0⋮⋮0⋯0]M×K=[Ω0],
where
(20)Ω=[t1,1⋯t1,K⋮⋮tK,1⋯tK,K]K×K.
But, the left hand side of (18) is a diagonal matrix, and the matrix Ps is diagonal as well. So, from (18) it can deduced that the matrix Ω is diagonal too, that is,
(21)Ω=[t1,1⋯0⋮⋮0⋯tK,K]K×K.
Thus, from (18), (19), and (21), each of the nonzero diagonal elements of Γ is given by
(22)γk=pk|tk,k|2,
which has a positive real value.Consequently, the EVD ofR~11 would be given by
(23)R~11=VΓ~VH,
where Γ~ is an M×M diagonal matrix and is given by Γ~=Γ+σ2IM×M, that is,
(24)Γ~=diag⁡([γ1+σ2⋯γK+σ2σ2⋯σ2]).Substituting (16) into (23), then we have
(25)R~11=A(θ)T+Γ~T+HA(θ)H,
where the notation (·)+ is the pseudoinverse for a given matrix.Now, taking the Cholesky decomposition ofR~11, then we have
(26)R~11=L~L~H,
where L~ is a lower triangular matrix that results from the Cholesky decomposition.Comparing (26) with (25), it can be deduced that
(27)L~=A(θ)T+Γ~1/2F,
where F is an orthogonal rotating matrix.A new matrix(Ψ) is formed by multiplying the left and right hand sides of R21 by L~-1 and L~-H, respectively, as follows:
(28)Ψ=L~-1R21L~-H.
But, R21=A(θ)PsZA(θ)H as shown in (9) and L~=A(θ)T+Γ~1/2F as shown in (27). Thus,
(29)Ψ=F-1Γ~-1/2TPsZTHΓ~-H/2F-H.
Now, let
(30)Φ=Γ~-1/2TPsZTHΓ~-H/2.
Also, let G=F-1. Then, (29) can be written as
(31)Ψ=GΦGH.
Looking at (31), if it is proved that the matrix Φ is a diagonal matrix, then (31) will be the EVD of Ψ, where G and Φ are the corresponding eigenvectors and eigenvalues matrices, respectively.To prove thatΦ is a diagonal matrix, consider the matrix Γ~-1/2 which is a diagonal matrix of size M×M and the matrix T which is an M×K matrix where its upper K×K part is diagonal, and its lower (M-K)×K part contains zero elements, then
(32)Γ~-1/2T=TΓ~s-1/2,
where Γ~s is a diagonal matrix of size K×K where its diagonal elements are the same as the first K diagonal elements of Γ~, that is,
(33)Γ~s=diag⁡([γ1+σ2⋯γK+σ2]).
Thus, (30) can be written as
(34)Φ=TΓ~s-1/2PsZΓ~s-H/2TH.
Since Γ~s-1/2, Ps and Z are diagonal matrices, then we have
(35)Φ=T|Γ~s|-1PsZTH.
But, the matrix Γ~s is diagonal matrix where its elements are positive real. So, |Γ~s|=Γ~s, and
(36)Φ=TΓ~s-1PsZTH.Substituting forT, Γ~s, Ps, and Z in (36), then we have
(37)Φ=diag⁡([ϕ1⋯ϕK0⋯0]),
where
(38)ϕk=pk|tk,k|2zkγk+σ2.
But, γk is defined in (22) as γk=pk|tk,k|2. Thus, (38) can be written as
(39)ϕk=γkzkγk+σ2.
Thus, it is clear that Φ is a diagonal matrix and it is proved now that (31) is the EVD of Ψ.Defining the term(40)yk=γkγk+σ2,
then (39) can be written as
(41)ϕk=ykzk.
Now, because yk has a positive real value, then
(42)∠ϕk=∠zk.
From the definition of zk, it is clear that the AOAs are implicated in ∠zk. Also, from (42) the phase of zk equals the phase of ϕk, and ϕk for k=1→K are the diagonal elements of Γ. Thus, to find θk, consider the kth diagonal element of Φ as follows:
(43)θ^k=cos⁡-1(λ∠(ϕ^k)2πdd),
where ϕ^k is the kth diagonal element of the estimated Φ^.Thus, the required AOAs are estimated.


### 3.1.  Summary of the CDBA Method

(1)
First step: calculateR21 from (9).
(2)
Second step: calculateR~11 from (11).
(3)
Third step: calculateL~ by taking the Cholesky decomposition of R~11.
(4)
Fourth step: calculate theΨ matrix by Ψ=L~-1R21L~-H.
(5)
Fifth step: calculate the eigenvalues ofΨ matrix.
(6)
Sixth step: calculate the AOAs from theK largest eigenvalues of Ψ matrix as shown in (43).


### 3.2. Analysis for the Effect of Noise Variance

The performance of the AOA estimators in high noise power (low SNR) cases is of interest. A measure of the noise power is the noise variance. A study of the noise variance effect on the performance of the CDBA algorithm is derived and compared with that of the TLS-ESPRIT.As shown in (43), the AOAs are estimated from the eigenvalues (ϕk). Thus, to start let us consider the effect of noise variance (power) on ϕk.Consider the magnitude ofϕk as follows:
(44)|ϕk|=|ykzk|.
But yk has a positive real value (recall that yk=(γk/(γk+σ2)), and σ2 and γk have positive real values). Thus, |ykzk|=|yk||zk|=yk|zk|. But, zk=exp{j(2π/λ)dd×cos(θk)}, so it is clear that |zk|=1. So, the magnitude of the eigenvalue (ϕk) is given by
(45)|ϕk|=yk=γkγk+σ2.
The noise variance effect on the ϕk can be measured by defining the variable ρk which is defined as the absolute difference between |ϕk| and its noise-free version (|ϕ¯k|), where ϕ¯k=ϕk|σ2=0=1. Thus, ρk is defined as follows:
(46)ρk=||ϕk|-|ϕ¯k||=σ2γk+σ2.
Clearly, the minimum value for ρk is 0 when σ2=0. Whereas its maximum value is 1 as σ2→∞.As for the TLS-ESPRIT algorithm, the AOA estimation starts by taking the EVD of an autocorrelation matrix of the received signals (call itR~TLS) [9]. Let RTLS be defined as the noise-free version of R~TLS. Also, let μ¯k for k=1→K be one of the matrix RTLS eigenvalues. Then, each eigenvalue of R~TLS (which is given the notation μk) will be given as follows:
(47)μk=μ¯k+σ2.
Defining ϵk as the absolute difference between |μk| and its noise free version (|μ¯k|), then,
(48)ϵk=σ2.
Thus the minimum value for ϵk is 0 when σ2=0. Whereas its maximum value approaches ∞ as σ2→∞.Comparingρk and ϵk in (46) and (48), respectively, it is clear that in low noise variance cases, the deviation between the estimated eigenvalues and the true eigenvalues will be small for both algorithms (CDBA and TLS-ESPRIT); thus, it is expected that both algorithms will have close performance in low noise variance cases (i.e., high SNR).In high noise variance cases, the deviation between the estimated eigenvalues and the true eigenvalues will have a maximum value of 1 in CDBA algorithm, whereas, in the TLS-ESPRIT the deviation will approach∞. Recall that the CDBA algorithm estimates the AOAs from the eigenvalues (ϕk) for k=1→K, and the TLS-ESPRIT estimates the AOAs from the eigenvectors of R~TLS (see [9] for more insight), and any error in estimating the eigenvalues of R~TLS will directly change their corresponding eigenvectors causing an error in estimating the AOAs. Also, it is clear that in high noise variance cases (i.e., low SNR) the EVD process in the CDBA algorithm is much less affected by the noise than the EVD process in the TLS-ESPRIT algorithm. Thus, it is expected that the CDBA algorithm will perform better than the TLS-ESPRIT in high noise variance cases (i.e., low SNR).The performance of the proposed CDBA method will be shown in the simulation section (Section4).


## 3.1.  Summary of the CDBA Method

(1)
First step: calculateR21 from (9).
(2)
Second step: calculateR~11 from (11).
(3)
Third step: calculateL~ by taking the Cholesky decomposition of R~11.
(4)
Fourth step: calculate theΨ matrix by Ψ=L~-1R21L~-H.
(5)
Fifth step: calculate the eigenvalues ofΨ matrix.
(6)
Sixth step: calculate the AOAs from theK largest eigenvalues of Ψ matrix as shown in (43).


## 3.2. Analysis for the Effect of Noise Variance

The performance of the AOA estimators in high noise power (low SNR) cases is of interest. A measure of the noise power is the noise variance. A study of the noise variance effect on the performance of the CDBA algorithm is derived and compared with that of the TLS-ESPRIT.As shown in (43), the AOAs are estimated from the eigenvalues (ϕk). Thus, to start let us consider the effect of noise variance (power) on ϕk.Consider the magnitude ofϕk as follows:
(44)|ϕk|=|ykzk|.
But yk has a positive real value (recall that yk=(γk/(γk+σ2)), and σ2 and γk have positive real values). Thus, |ykzk|=|yk||zk|=yk|zk|. But, zk=exp{j(2π/λ)dd×cos(θk)}, so it is clear that |zk|=1. So, the magnitude of the eigenvalue (ϕk) is given by
(45)|ϕk|=yk=γkγk+σ2.
The noise variance effect on the ϕk can be measured by defining the variable ρk which is defined as the absolute difference between |ϕk| and its noise-free version (|ϕ¯k|), where ϕ¯k=ϕk|σ2=0=1. Thus, ρk is defined as follows:
(46)ρk=||ϕk|-|ϕ¯k||=σ2γk+σ2.
Clearly, the minimum value for ρk is 0 when σ2=0. Whereas its maximum value is 1 as σ2→∞.As for the TLS-ESPRIT algorithm, the AOA estimation starts by taking the EVD of an autocorrelation matrix of the received signals (call itR~TLS) [9]. Let RTLS be defined as the noise-free version of R~TLS. Also, let μ¯k for k=1→K be one of the matrix RTLS eigenvalues. Then, each eigenvalue of R~TLS (which is given the notation μk) will be given as follows:
(47)μk=μ¯k+σ2.
Defining ϵk as the absolute difference between |μk| and its noise free version (|μ¯k|), then,
(48)ϵk=σ2.
Thus the minimum value for ϵk is 0 when σ2=0. Whereas its maximum value approaches ∞ as σ2→∞.Comparingρk and ϵk in (46) and (48), respectively, it is clear that in low noise variance cases, the deviation between the estimated eigenvalues and the true eigenvalues will be small for both algorithms (CDBA and TLS-ESPRIT); thus, it is expected that both algorithms will have close performance in low noise variance cases (i.e., high SNR).In high noise variance cases, the deviation between the estimated eigenvalues and the true eigenvalues will have a maximum value of 1 in CDBA algorithm, whereas, in the TLS-ESPRIT the deviation will approach∞. Recall that the CDBA algorithm estimates the AOAs from the eigenvalues (ϕk) for k=1→K, and the TLS-ESPRIT estimates the AOAs from the eigenvectors of R~TLS (see [9] for more insight), and any error in estimating the eigenvalues of R~TLS will directly change their corresponding eigenvectors causing an error in estimating the AOAs. Also, it is clear that in high noise variance cases (i.e., low SNR) the EVD process in the CDBA algorithm is much less affected by the noise than the EVD process in the TLS-ESPRIT algorithm. Thus, it is expected that the CDBA algorithm will perform better than the TLS-ESPRIT in high noise variance cases (i.e., low SNR).The performance of the proposed CDBA method will be shown in the simulation section (Section4).


## 4. Simulation Results

Simulations for the proposed CDBA estimator were completed to assess its performance. The results in this section were averaged over 1000 ensemble runs. The elements of each antenna array were separated by a half-wavelength (i.e.,d=λ/2), and without any loss of generality the distance between the elements of each antenna doublet (dd) was chosen to be λ/2 as well. The number of sources was set to 2. A1 was set to 1.2, and A2 was set to 1. The proposed CDBA method was compared with the TLS-ESPRIT method of [9].Figure1 shows the root-mean-square error (RMSE) of the AOA estimation in degrees for the proposed CDBA method compared with the TLS-ESPRIT method for different SNRs. Since the power for both users was not equal, then the SNR was taken for the second user; that is, if the SNR was set to 5 dB, then the second user would have a 5 dB SNR and the first user would have SNR of 20log1.2+5=6.58 dB. The angles for the two users were set to θ1=30° and θ2=75°. The number of antenna doublets was set to M=3. The number of snap shots was set to W=200 and W=1000 (where W is the number of snap shots over which the correlation matrices were estimated). The results in Figure 1 show that the proposed CDBA method gave better performance than the TLS-ESPRIT method in low SNR cases. This result agrees with the discussion presented in Section 3.2.Figure 1
Root mean square error (RMSE) of angle estimation in degrees versus signal-to-noise-ratio (SNR) in dB for(θ1,2=30°,75°).Figure2 shows the RMSE of the AOA estimation in degrees for the proposed CDBA method compared with the TLS-ESPRIT method for different number of antenna doublets. The angles for the two users were set to θ1=30° and θ2=75°. The number of snap shots was set to W=200. The results in Figure 2 show that the CDBA method outperformed the TLS-ESPRIT algorithm for different number of antenna doublets.Figure 2
Root mean square error (RMSE) of angle estimation in degrees versus number of antenna doublets for(θ1,2=30°,75°).Another important issue in comparing AOA estimators is their performance when the AOAs of the received signals are close to each other. Table1 shows the RMSE of the AOA estimation in degrees for the proposed CDBA method compared with the TLS-ESPRIT method for different AOA separations for two users. The number of antenna doublets was set to M=6. The number of snap shots was set to W=200. The angular deviation was performed by setting θ1 to be taken from the following equation: θ1=θ2+δ where δ is the angular deviation with θ2=75°. Table 1 shows that the CDBA method outperformed the TLS-ESPRIT algorithm for different AOA separations.Table 1
Root mean square error (RMSE) of angle estimation in degrees for different AOA separations.
AOA separation(δ)
TLS-ESPRIT
CDBA
1°
11.5223°
3.1770°
1.5°
6.8057°
1.6336°
2°
4.2739°
0.7330°Also, to compare both methods computational complexity, two functions in MATLAB (tic.m and toc.m) were used to measure the time it takes each method to estimate the AOA. The CDBA and the TLS-ESPRIT methods took 0.281 and 0.297 seconds, respectively, to perform one run of AOA estimation for two users. Clearly, this indicates that the CDBA method has less computational complexity than the TLS-ESPRIT method.


## 5. Conclusion

In this paper, an AOA estimator is proposed which is named the CDBA method. The CDBA method is applied by taking the Cholesky decomposition of the received signal autocorrelation matrix. The resultant decomposed matrices are used to modify the crosscorrelation matrix of the received signals at the antenna array doublets. The proposed CDBA method has better performance than the TLS-ESPRIT method in estimating the AOAs especially in low SNR cases. The performance of the proposed CDBA method was assessed and compared to the TLS-ESPRIT method.



---

*Source: 803617-2012-06-27.xml*