Open Access 2D DOA estimation with sparse uniform circular arrays in the presence of mutual coupling

In this article, we consider the uniform circular arrays (UCAs) with the number of antenna elements insufficient to apply the traditional beamspace-based algorithms, which are labeled as sparse UCAs. For such UCAs, we propose a new hybrid approach for 2D direction-of-arrival (DOA) estimation in the presence of mutual coupling. Using the manifold decomposition technique, we present two new formulations of the steering vector in the presence of mutual coupling for sparse UCAs. Then, we introduce the adaptations to a modified uniform circular array rank reduction algorithm. This leads to an algorithm that is able to estimate the azimuth angle without the exact knowledge of mutual coupling. Next, we use a search-free rooting algorithm which expands the steering into a double Fourier series for each estimated azimuth to obtain the elevation angle estimates. The manifold decomposition technique introduces truncation errors. However, the accuracy of the DOA estimates is strongly affected by these errors when the array has a small number of elements. Therefore, expressions describing the truncation errors in the DOA estimates are derived. This allows us to choose an appropriate truncated degree in the manifold separation transformation to enhance the DOA estimate accuracy. Numerical examples are presented to demonstrate the effectiveness of the proposed method.


Introduction
The problem of two-dimensional (2D) directions-of-arrival (DOAs) [1][2][3][4][5][6][7] estimation (i.e., azimuth and elevation angles) has received increasing attention in a variety of applications, such as radar, mobile communications, sonar, and seismology. In general, a planar array is needed when estimates of source azimuth and elevation are required. Such well-known planar arrays include the two-orthogonal uniform linear array (the L-shaped array) [1], the rectangular array [2], and the uniform circular array (UCA) [3][4][5][6][7]. The UCA is able to provide 360°of coverage in the azimuth plane. Moreover, the UCA has uniform performance regardless of angle of arrival. Thus, UCA attracts more attention than other planar arrays recently. Due to the circular symmetry, the beamspace transformation, based on the phase-mode excitation principle, is usually applied to obtain the desired Vandermode structure for the steering vector in the mode space. This transformation results in the development of several DOA estimation algorithms with low computational cost, such as UCA-RB-MUSIC [3], UCA-ESPRIT [3], and uniform circular array rank reduction (UCA-RARE) [4]. However, all these algorithms for UCAs ignore the mutual coupling effect, which ultimately destroys the underlying model assumptions needed for their efficient implementations. Moreover, all these algorithms, based on the traditional beamspace transformation, require a sufficiently large number of elements to avoid aliasing in the steering vector of the mode space.
In this article, we focus on UCAs with the number of antenna elements insufficient to apply the traditional beamspace-based algorithms. In [8], such UCAs are labeled as sparse UCAs and are allowed to adopt an efficient search-free and robust 1D DOA estimation algorithm. This algorithm is based on a modified beamspace transformation and is called as sparse UCA Root-MUSIC. In this algorithm, all relevant phased modes are able to be incorporated in a polynomial rooting procedure leading to biased free estimates when the number of elements of UCA is small. However, the algorithm in [8] only estimates the azimuth angle without considering the mutual coupling, given a fixed elevation angle. In this article, we estimate the 2D DOAs with such sparse UCAs in the presence of mutual coupling, acquiring both the azimuth and elevation angles estimates via a manifold decomposition technique. The straightforward extension of MUSIC to 2D DOA estimations brings on a 2D search over the MUSIC spectrum and has a high computational cost. In [5,7], both the proposed algorithms take the mutual coupling into account and employ the UCA-RARE algorithm to estimate the azimuth angle first. With the open-circuit voltages of the antenna elements expanded in spherical mode, a Root-MUSIC algorithm is able to be performed in the elevation space to obtain the elevation estimates in [5]. In [7], a 1D parameter search replaces the implementation of Root-MUSIC algorithm in the elevation space. In the 1D parameter search for elevation estimates, the elevation-dependent mutual coupling effect can efficiently be compensated by the elevation-dependent receiving mutual impedances. However, this step results in higher computational load. Although these two algorithms are applied for the compact UCAs, they could inspire us for sparse UCAs.
In this article, we use the method proposed in [9] to calculate the mutual coupling. In [6,7,9], computer simulations have shown that this method can produce more accurate DOA estimation results than the open-circuit voltage method. The experiments in [6,7] show that the mutual coupling matrix (MCM) depends on the elevation angle for UCAs. Moreover, the simulation results in [6,7] have shown that the way compensating the mutual coupling with single-elevation-angle receiving-mutual-impedance, computed according to the method in [9], still produces better DOA estimation results than the opencircuit voltage method. In [7], it shows that the variation of the receiving mutual impedances with elevation angle is a process of gradual change. Hence, it is feasible to estimate the elevation angle using mutual coupling compensated with single-elevation-angle receiving-mutualimpedance. In [10], it is shown that any array steering vector can be expanded on a spherical surface to generate an expression containing spherical harmonics, which can be mapped to 2D Fourier basis [11]. In this article, we will extend this expansion to the steering vector in the presence of mutual coupling.
In this article, we propose a new hybrid algorithm for 2D DOA estimation in the presence of mutual coupling for sparse UCAs. Based on the manifold decomposition technique, we will present two new formulations of the steering vector in the presence of mutual coupling for sparse UCAs. One formulation, corresponding to Jacobi-Anger expansion [12], allows applying a modified UCA-RARE algorithm to estimate the azimuth angle without the exact knowledge of mutual coupling and elevation angle. The other formulation, corresponding to Bauer's formula [13], allows executing a Root-MUSIC algorithm in the elevation direction to estimate the elevation angle for each estimated azimuth angle. For sparse UCAs, compared with the original UCA-RARE, the modified UCA-RARE is able to avoid obtaining spurious estimates which only arise from the sparseness of the array elements. Note that the steering vector expansion for estimating the elevation angle in this article differs from that in [5] and has a more universal application [10,14,15]. In fact, these two kinds of decomposition techniques applied in this article can be considered as manifold decomposition transformations [11,[16][17][18]. It is shown that the DOA estimate accuracy usually depends on the truncation error introduced by the transformation [16][17][18][19]. Hence, we analyze the truncation errors for sparse UCAs and derive expressions describing the truncation errors in the DOA estimates. We find that the impact of the truncation error on the estimate accuracy of azimuth angle is weaker than it for the elevation angle estimate. Therefore, a method to choose an appropriate truncated degree for the elevation estimates is presented to enhance the estimate accuracy.
The rest of the article is organized as follows. First, the array signal model is presented in Section 2, followed by the description of the manifold decomposition technique in Section 3. Then, the proposed algorithm for sparse UCA is presented in Section 4. The impact of the truncation errors are analyzed in Section 5. Section 6 shows the simulation results. Finally, Section 7 concludes the article.

Array signal model
Consider a sparse UCA consisting of N identical elements uniformly distributed over the circumference of a circle of radius r. Assume that D narrowband sources, centered on wavelength λ, impinge on the array from directions θ i (i = 1,..., D) and j i (i = 1,..., D), respectively, where θ i [0, π/2] is the elevation angle measured from the Z-axis and j i [0, 2π) is the azimuth angle measured from the X-axis counter-clockwise (see Figure 1). The N× 1 vector received by the array is expressed as where A(θ , φ) = a (θ 1 , φ 1 ) . . . a(θ D , φ D ) is the N × D matrix of the steering vectors, s(t) = [s 1 (t)... s D (t)] T is the D× 1 signal vector, n(t) = [n 1 (t)... n N (t)] T is the N× 1 noise vector. The signal vector s(t) and the vector n(t) of the additive and spatially white noise are assumed to be statistically independent and zero-mean.
If the sparse UCAs are composed of omni-directional antenna elements, its steering vector in the presence of mutual coupling [20] is given by where the N × N matrix c m is the elevation-dependent MCM. Due to the circular symmetry, a model for the MCM of UCAs [5,7,20] can be a complex symmetric circulant matrix. In this article, we only consider N = 2K+1. Therefore, C m can be expressed as The N × 1 vector a(θ, j) is the ideal steering vector and the expression of its nth (n = 1,..., N) component is where g n = 2π(n-1)/N is the angular position of the nth element.
If the sparse UCAs are composed of directional antenna elements, its steering vector in the presence of mutual coupling [21] can be described as where a d (θ, j) is the steering vector having the directional pattern of the form g d (θ, j) and the expression of its nth (n = 1,... N) component is

The manifold decomposition technique
Here, the concepts of the Wavefield modeling for scalarfields are given. Two expressions for decomposing the steering vector (manifold) of array elements with sparse UCAs are presented. These form the theoretical basis of our algorithm.
In [10], it is shown that the array steering vector is able to be decomposed as where Γ s represents the so-called sampling matrix and b(θ, j) is the basis functions of the decomposition. In general, the dimension, i.e., the number of basis functions is infinite in order to hold the equality exactly. Therefore, the sampling matrix can be considered as an operator defined as s : H → C N×1 . The coefficients of the expansion (sampling matrix) map functions defined on H into the Nth-dimensional complex space N×1 . Hence, the sampling matrix is a characteristic of the array only.
Usually there are two kinds of choices for the basis functions b(θ, j). For sparse UCA, one is given by where m =..., -1, 0, 1,... and J m (•) is the Bessel function of the first kind with order m. Its corresponding sampling matrix is s n,m = e −jmγ m (9) where n = 1, 2,..., N. This kind of decomposition can also be considered as the Jacobi-Anger expansion [12], and Equation 4 can be expressed as The other choice is the spherical harmonics and the where t = l 2 +l+m+1 and P m l (cos θ) represents the associated Legendre functions of the lth degree (or level) and mth order (or mode). Its corresponding sampling matrix is s n,t = 4π j l j l (kr) Y m * l π 2 , γ n (12) where j l (•) denotes the spherical Bessel function of the first kind. The Bauer's formula [13] corresponds to this kind of decomposition and is given by , γ n Y m l (θ , φ) (13) Since both J m (kr) and j l (kr) decay exponentially, we can assume that, for m ≫ kr and l ≫ kr, the higher-  order Bessel functions and spherical Bessel functions are negligible. Therefore, the sampling matrix can be truncated by considering a finite number of modes or degrees. Ideally, the resulting truncation error can be made arbitrarily small just by increasing the number of modes or degrees. We assume that the truncated order is M and truncated degree is L for the first and second kinds of decomposition, respectively. The rule to select the truncated order or degree will be discussed in Section 5. In order to distinguish the sampling matrices and basis functions for two kinds of decomposition, let s 1 and b 1 (θ, j) denote the sampling matrix and basis functions for the first kind of decomposition and s 2 and b 2 (θ, j) be the sampling matrix and basis functions for the second kind of decomposition, respectively.
In fact, an alternative expression [10] of b 2 (θ, j) with limited degree L is given by where the (L+1 is the combination of the selection matrix and the coefficients vectors and is given by can be obtained for an arbitrarily l and m using two recurrence expressions. The (2L+1) × 1 vector is described as d(θ) = [e -jLθ ··· 1 ··· e jLθ ] T . More details about Equation 14 can be found in [11]. Apparently, Equation 14 is an expansion of 2D Fourier series.

The hybrid algorithm to DOA estimation
For the ideal UCAs composed of omni-directional antenna elements, exciting the array with the weight For the traditional beamspace transform, there is N > 2 kr and the first term in (15) becomes dominant. However, for the beamspace transform applied to the UCA with N = 2K+1, where K < kr and N ≤ 2 kr , the value of the second term may be significant and cannot be neglected because of the contribution of J k±qN (kr sin θ) with orders K < |k| ≤ kr .
We also label such UCAs as sparse UCAs [8].
Obviously, the algorithms proposed in [3][4][5]7], which are based on the traditional beamspace transform, cannot be employed directly for such UCAs. In this section, we will present a new hybrid algorithm applied to such UCAs. In order to avoid 2D search in MUSIC spectrum, we estimate the DOAs in two steps. Based on the beamspace transformation corresponding to the Jacobi-Anger expansion, first we estimate the azimuth angle using the modified UCA-RARE algorithm, which stems from the original UCA-RARE applied for compact UCA. This algorithm is attractive since it decouples azimuth estimation from elevation estimation and relaxes the assumption of omni-directional element patterns. Then, we perform a Root-MUSIC algorithm to estimate the elevation angle for every estimate azimuth angle using the expansion based on the Bauer's formula.
We start from the signal model of Section 2. Recalling Equation 1, the beamspace array signal model is where the (2K+1) × M weight matrix W k is defined as The corresponding beamspace steering vector is The covariance matrix R of the beamspace data is constructed and an eigendecomposition of R results in a signal and noise subspace where E s and E n denote the signal and noise subspace eigenvectors and the diagonal matrices s and n contain the signal subspace and noise subspace eigenvalues, respectively. The beamspace MUSIC algorithm estimates the DOAs from the D deepest nulls of the MUSIC function

The azimuth angle estimation
When a(θ, j) is expanded as the first kind of decomposition, the beamspace steering vector a b (θ, j) can be expressed as The vector a 1 b (θ , φ) represents the truncation errors term and its element is the summation of j m J m (kr sin θ) e jmj with M < |m| ≤ ∞. This term can be arbitrarily small just by increasing M. The matrix J l is a (2K+1) × (M-K) matrix consisting of the M-K last columns of the unity Note that the expression of H in Equation 22 is restricted to M < 3K+1 [8]. The expression of Equation 21 is equivalent to the one for beamspace manifold in [8]. Now we will extend this transformation to the case considering the mutual coupling.
We should note that [22] where V 1 = I and V n (n > 1) is the (n-1)th power of the cyclic permutation operator given by Rewrite Equation 2 as Using V n a(θ, j) = a(θ, j-2π(n-1)/N) and Equation 21 yields where M n = diag e j2π K(n−1) / N · · · 1 · · · e −j2π K(n−1) / N . Hence, Equation 18 can be re-expressed as where Because of the symmetry of the mutual coupling coefficients and the periodicity of e j2π K(n−1) / N , it is easy to find that the diagonal elements of M s is centro-symmetry. If we neglect the truncation errors term, Equation 27 can be modeled as where g (θ ) = m g (θ ) and m is the first M+1 elements of the diagonal elements of M s . "⊙" denotes the Hadamard product of vectors. The beamspace steering vector in the presence of mutual coupling for compact UCAs [5,7] is a special case of the one for sparse UCAs with H = I. Note that the components of g (θ ) in Equation 29 have the same expression form with the ones in [7] and different expression form from the ones in [5].
Please also note that this transformation is able to be extended to the case that the sparse UCAs are composed of directional antenna elements. Recalling Equation 6, there is In general, the element directional pattern g d (θ, j) [21] can be expressed as Then, Equation 30 can be expressed as It is easy to find that So that we finally get where Similarly, we can get the beamspace manifold in the presence of mutual coupling for arrays composed of directional elements where g (θ ) = m g (θ ) and Observing Equations 29 and 36, the beamspace manifolds for omni-directional elements and directional elements have the same expansion form with different components of g(θ).
Replacing the beamspace steering vector by its factorization (29) or (36), the MUSIC function becomes When N-D ≥ M+1 [4], this structure allows using a rank reduction algorithm, named UCA-RARE. Therefore, we can root the sample polynomial and then find the signal azimuth angle from roots of (39), which are located closest to the unit circle. Note that det{•} is the determinant of a matrix. Making use of a well-known identity for block matrices which holds true for arbitrary matrices B, C, and nonsingular matrices A, D.
where Ψ = T T (1/z)H H HT(z). It is obvious that the roots of det{Ψ}, which are also the roots of (39), are only the spurious roots arising from the sparseness of the array elements and independent of the received data. Hence, a spatial spectrum function for the azimuth estimation with sparse UCAs can be constructed as However, there may be common roots for det{Ψ} and Similar to the Root-MUSIC roots, RARE roots enjoy the so-called conjugate reciprocity property, i.e., if z 0 is a root of P 2 (z), thenz 0 = 1 z * 0 is also a root of P 2 (z). Therefore, there are spurious estimates j i +π for j i <π and j i -π for j i > π. Although there are still spurious estimates (j i + j j )/2 for the case of impinging sources with the same elevation angle (θ i = θ j ), we do not plan to eliminate them in order to avoid cancelling the real root at (j i + j j )/2 when there is a source exactly at (j i + j j )/2. Besides all these spurious estimates, there may be other spurious estimates introduced by the sparseness of the array elements. However, all spurious estimates can be eliminated in the final paired 2D DOA estimation by the elevation estimate in the next step.

The elevation angle estimation
A specifically designed closed-form algorithm similar to UCA-ESPRIT is proposed in the original UCA-RARE algorithm [4] to obtain the elevation estimates. Although it is a search-free implementation, there are some shortcomings that make it somewhat unsuitable for practical application, which are presented in detail in [7]. Hence, we apply for the Root-MUSIC algorithm via decomposing the steering vector into the double Fourier series to estimate the elevation angle. Note that the steering vector expansion in the presence of mutual coupling for estimating the elevation angle in this article differs from that in [5]. The method in [5], which estimates the DOAs for compact UCAs, is based on the open-circuit voltages of the antenna elements expanded in spherical mode, whereas our method, in which the mutual coupling is calculated by the proposed approach in [9], grounds on the manifold decomposition. In [5], the steering vector is expanded into a limited Fourier series of phase modes by considering a general multiport antenna, carrying a current distribution C(r, , z) on the surface S of a cylinder with radius r and height z max . In our proposed method, the steering vector is expanded into a limited Fourier series of phase modes by considering an element on the surface of a unit sphere. Moreover, the truncation degree in [5] is determined by the radius r and height z max together, while the one in our proposed algorithm is only relative to r.
In [7], it shows that the variation of the receiving mutual impedances with elevation angle is a process of gradual change. It means that the receiving mutual impedances do not vary with elevation angle significantly. The simulation results in [6,7] have shown that estimating the elevation angle with single-elevation-angle receiving-mutual-impedance could achieve error accuracy around 1°. Therefore, it is feasible to estimate the elevation angle using mutual coupling compensated with single-elevation-angle receiving-mutual-impedance. We will estimate an initial elevation angle using the MCM obtained at θ = 45°first. Then we can get a more accurate result with the MCM obtained at the initial estimate.
When a(θ, j) is expanded as the second kind of decomposition and b 2 (θ, j) is decomposed as Equation 14, the beamspace steering vector a b (θ, j) can be expressed as where B = W H K s 2 and a 2 b (θ , φ) represent the truncation errors term of the second kind of decomposition. The mth elements of a 2 b (θ , φ) is where t = K+1-n. In the presence of mutual coupling, the corresponding beamspace steering vector a b (θ , φ) , defined in Equation 18, is where B = W H K C m s 2 and a 2 b (θ , φ) represent the corresponding truncation errors term. The nth elements of a 2 b (θ , φ) is If neglect the truncation errors term a 2 b (θ , φ) , the beamspace MUSIC function becomes Apparently, this equation is a polynomial in w = e jθ for each estimate azimuth angle j i . The Root-MUSIC algorithm can be performed to obtain the elevation angle estimates. A method, which extends the steering vector in the elevation field from [0, π] to [0, 2π], is presented in [5] to decrease the implementation times of the Root-MUSIC algorithm. In fact, this method only considers the case that there are only spurious estimate j i + π. This will result in errors for j i >π whose spurious estimate is j i -π. Hence, we prefer to perform Root-MUSIC algorithm for each azimuth estimates separately rather than perform Root-MUSIC algorithm with combining j i and j i + π. Note that there are usually two optimal solutions θ i and π θ i for one azimuth angle j i due to sin θ i = sin (π-θ i ). This characteristic is determined by the symmetry of the circular array's manifold in the elevation range. The real root should be located at (0, π/2] for a circular array. It is clear that the following equation holds true for (θ, j): That shows that the solutions to j i + π or j i -π are 2π-θ [3π/2, 2π] and 2π-(π-θ) = π+θ [π, 3π/2]. So, this algorithm allows eliminating the spurious estimate j i +π or j i -π automatically. Again the spurious estimates (j i +j j )/2 and (j i +j j )/2+π can only keep one result if there are sources with the same elevation angle (θ i = θ j ). Although all spurious azimuth estimates are considered, we only reserve the paired DOAs (θ i , j i ) whose elevation estimates locate at [0, π/2]. The number of such paired estimate (θ i , j i ) may be more than D. Hereby, it is necessary to calculate the MUSIC function for every paired estimate (θ i , j i ). Only the D smallest values of the MUSIC function are considered as the final estimates for the DOAs (θ, j).
For the sparse UCAs composed of directional antenna elements, we could still execute the Root-MUSIC algorithm. Recalling Equations 6 and 13, there is Then, the components of the sampling matrix becomes s n,t = 4π j l j l (kr) Since the value of the azimuth angle has been estimated, g d (θ, j-g n ) is only a function of elevation angle θ and can be labeled as g n d (θ ) . Usually the direction pattern g d (θ) is able to be expressed as a function of cos θ and sin θ. We define w = e jθ . It is easy to get cos θ = (w +w -1 )/2 and sin θ = -j(w-w -1 )/2. In such case, the beamspace steering vector a b (θ , φ) has the same expansion form as it is in Equation 46 but with different components of the sampling matrix. It could sill be written as a polynomial in w. Hence, the Root-MUSIC algorithm is still able to be performed to estimate the elevation angle.
The steps involved in the proposed hybrid algorithm can be summarized below: (1) Compute the sample covariance matrix R = 1 P P p=1 x p x H p by averaging over P data snapshots. Compute the beamspace covariance Form the matrix E s and E n , which spans the estimated signal subspace and the noise subspace, respectively.
(3) Obtain the azimuth angle estimates with Equation 43. All spurious azimuth estimates are reserved. (4) For every reserved azimuth estimate, perform the Root-MUSIC with the MCM obtained at θ = 45°to find an initial elevation estimates. Since the location of the real elevation angle is θ (0, π/2], except for the case that there is a source exactly at j i + π or j i -π, all spurious estimates j i +π or j i -π can be eliminated automatically (see Equation 49). Then, we perform the Root-MUSIC with the MCM obtained at initial estimate to get a more accurate estimate. (5) Calculate the MUSIC function for every paired estimate (θ i , j i ). Take the paired estimates (θ, j) corresponding to the D smallest values of the MUSIC function as the final estimate.

The impact of the truncation errors on the estimation accuracy
As discussed in Sections 3 and 4, the manifold decomposition will introduce truncation errors. Here, a firstorder approximation of the bias based on manifold decomposition is derived for sparse UCAs, and a thumb rule to choose the truncate degree is presented.
The truncation errors terms a i b (θ , φ) (i = 1, 2) in Equations 28 and 47 are dependent on the DOAs (θ, j) and truncation degree L or order M. Let's consider a sparse UCA with N = 11 monopoles tuned to f 0 = 2.4 GHz. The radius is r = λ. The monopole elements are of equal length 3.13 cm and radius is 0.3 mm. All monopole elements are loaded with a terminal load Z 0 = 50Ω. The receiving mutual impedances shown in Table 1 are calculated with receiving-mutual-impedance method [9] for different elevation angles of an impinging source. Figure 2 shows the value of a i b (θ , φ) (i = 1, 2) as a function of both azimuth angle j and elevation angle θ for this UCA with L = M = kr = 7. Notice that it needs more than 15 elements for the traditional beamspace transform. It is easy to see that a 1 b (θ , φ) is always smaller than a 2 b (θ , φ) when these two kinds of decomposition are performed with the same DOAs (θ, j) and truncation degree L or order M. a 1 b (θ , φ) decreases as the elevation angle changes from 90°, while a 2 b (θ , φ) varies over the elevation angle and hold bigger for the angle near 0°and 90°. Both a 1 b (θ , φ) and a 2 b (θ , φ) nearly have the same value over the azimuth angle for a fixed elevation angle.

Analysis of the bias in the azimuth angle estimation
As we know, we estimate the azimuth angle based on the rank reduction theory. For the true j i , there is n HT (φ i ) and the eigendecomposition of F(j i ) is given by If there are m i sources with the same azimuth angle j i and n i sources with the azimuth angle j i +π or j i -π, then Σ s denotes the diagonal matrix containing the D-(m i +n i ) non-zero eigenvalues and Σ n contains the remaining m i +n i zero eigenvalues [4,7]. The matrices U s and U n in turn contain the corresponding eigenvectors, respectively. Define a function where u q represents the qth columns of U n . Due to the truncation errors, y(j i ) ≠ 0 but y(j i ) ≈ 0. An expression for the basis can be found by expanding the first derivative of Equation 54 with respect to j i and evaluating atφ i . For small enough errors, we have [16,17,23]: whereφ i and j i are the estimated and true azimuth angle and y φ i ∂y (φ) ∂φ φ=φ i . Let T'(j) ≜ ∂T (j)/∂j and, similarly, T"(j) ≜ ∂T'(j)/∂j. The first derivative of y(j) is where Re{•} stands for the real part of the argument within the brackets. Similarly, the second derivative of y (j) is yielded as By combining (55)-(57), the bias for the azimuth angle estimates at angle j i can finally be computed from

Analysis of the bias in the elevation angle estimation
Similarly, an expression for the elevation estimation basis can be found by expanding the first derivative of Equation 48 with respect to θ i and evaluating atθ i . For small enough errors, there is [16,17,23]: whereθ i and θ i are the estimated and true elevation angle and f θ i ∂f (θ ) ∂θ θ =θ i . We define the vectors d'(θ) = ∂d(θ)/∂θ and d"(θ) = ∂d'(θ)/∂θ. The first derivative of f(θ) in Equation 48 with respect to θ is And the second derivative of f(θ) is given by Combining (59)-(61), an expression of the bias for the elevation angle estimates at angle θ i can be finally written aŝ In order to verify the approximation (58) and (62), we perform N e = 500 independent experiments using a sparse UCA in Figure 2. In the simulation we used one source with SNR = 25 dB and j = 200°moving in the elevation range θ (0, π/2) since both a 1 b (θ , φ) and a 2 b (θ , φ) are highly dependent on the elevation angle. Figure 3a shows the estimated MSE and the approximated MSE computed from Equation 58 with different M for the azimuth estimation. The approximated bias follows the trend of a 1 b (θ , φ) . However, for the estimation, it works badly when the elevation angle is near 0°and 90°. That is, because the estimate azimuth angles can be any value when the elevation angle is 0°and a 1 b (θ , φ) is large when the elevation angle is near 90°. In fact, the bias computed from Equation 58 can be approximately expressed as a function of the truncation errors when there is only one source (see Appendix). Because the truncation errors are very small when the elevation angle is near 0°(see Figure 2), the   estimate accuracy is enhanced when L increases. Hence, in order to obtain good estimate accuracy, a reasonable criterion to determine the truncated degree L could be where j max = max j l (kr) , l ≤ 2 kr and the predetermined ε is related to the truncated accuracy.

Simulations
A sparse UCA of radius r = λ with N = 11 is employed in all experiments. The distance between each elements is 0.5635λ. The signals and noise in our simulations are assumed to be stationary, zero mean, and uncorrelated Gaussian random processes. Noise is both spatially and temporally white. The truncation order M is 7. In order to reduce the impact of the truncation errors on the estimate accuracy, the truncation degree here is L = 13 for ε = 0.001. The first example shows the spurious azimuth estimates of the original UCA-RARE. We consider the case of two impinging uncorrelated sources at (θ 1 , j 1 ) = (10°, 40°) and (θ 2 , j 2 ) = (30°, 150°). The SNR = 20 dB is quoted per source per array element. Figure 4 depicts the root distributions nearby the unit circle for Equations 39 and 43. As shown in Figure 4a, except for the true azimuth roots 40°and 150°and their corresponding spurious roots 220°and 330°, many other spurious roots introduced by the sparseness of the array element will appear if we employ the original UCA-RARE algorithm. Although the final estimates (θ, j) may be obtained by calculating the MUSIC function, the computational cost due to solving the Root-MUSIC for every azimuth estimate will increase. Therefore, it is better to use the modified UCA-RARE (see Equation 43).
In the second example, one source with SNR = 25 dB and j = 140°is used to move in the elevation range θ (0, π/2). Two kinds of MCM, one obtained at θ = 45°a nd the other obtained at the accurate elevation angle, are employed to estimate the elevation angle. Sample statistics are computed from 500 independent trials. The root-mean-square-errors (RMSEs) and Cramer Rao Bounds (CRB) of the elevation estimate are shown in Figure 5. Although the result for the method with the MCM obtained at θ = 45°is a little inferior to the one with elevation-dependent MCM, the estimate difference for them is below 1°. It shows that it is feasible to get the elevation-dependent MCM by using the MCM obtained at θ = 45°first to obtain an initial estimate.
In the third example, the case of three impinging uncorrelated sources at (θ 1 , j 1 ) = (25°, 60°), (θ 2 , j 2 ) = (25°, 120°), and (θ 3 , j 3 ) = (50°, 330°) are considered. The SNR = 20 dB is quoted per source per array element.  Except the  spurious azimuth estimates 60°+180°= 240°, 120°+180°= 300°, 330°-180°= 150°, (60°+120°)/2 = 90°, and 90°+ 180°= 270°, we can find other spurious azimuth estimates caused by the sparseness of the array element in Figure 6a. Figure 6b depicts the corresponding elevation estimates distributions nearby the unit circle. It is easy to find that only the roots for the true azimuth estimates j i are nearest by the circle and locate at [0, π/2]. The elevation estimates roots for the spurious azimuth estimates 240°, 300°, and 150°are also very close to the circle, but locate at [π, 2π]. The elevation estimates roots corresponding to the spurious azimuth estimates 90°and 270°are quite close to the unit circle and other spurious elevation estimates are all outside the unit circle. Hence, the spurious azimuth estimates 240°, 300°, and 150°can be eliminated automatically in the implementation of Root-MUSIC via judging their corresponding elevation estimates roots location. For the reserved paired DOAs (θ i , j i ) whose elevation roots locate at [0, 2π], besides to determine the final DOAs via computing the MUSIC function, here another approach is provided via judging the distance of elevation roots from the unit circle. We can consider the estimates corresponding to D roots located closest to the unit circle as the final estimates.
The fourth example shows the performance of the proposed algorithm and MUSIC algorithm for different SNR levels. Sample statistics are computed from 500 independent trials. Two signals arrive at the array from two equal power sources from directions (θ 1 , j 1 ) = (10°, 20°) and (θ 2 , j 2 ) = (30°, 60°), respectively. The SNR level is varied from 0 to 35 dB and is quoted per source per array element. For MUSIC algorithm, the DOA estimates are obtained by searching the highest local peaks over the 2D MUSIC spectrum. While for the proposed algorithm, the DOA estimates are found by implementing two computationally efficient rooting algorithms (the modified UCA-RARE and the Root-MUSIC algorithm). Here, the search step for MUSIC algorithm in both elevation and azimuth direction is 0.01°. The RMSEs plots of estimated DOAs and CRBs are shown in Figure 7. It is observed that the RMSEs of the estimates for the proposed algorithm and MUSIC algorithm follow the trend of the corresponding CRBs. The average Matlab-runtime to estimate the DOAs by the proposed algorithm is 0.78 s (simulated on a 2.8-GHz Intel Core i5 CPU and 2.99 G Ram), whereas the average Matlab-runtime to estimate the DOAs by the traditional search of the 2D MUSIC spectrum is 7.43 s. However, the estimated results for the proposed algorithm are comparable to the ones for MUSIC algorithm.
The last example is used to examine the capability of the proposed algorithm to estimate the DOAs of the closely spaced signals. Two uncorrelated sources with DOAs (θ 1 , j 1 ) = (20°, 100°) and (θ 2 , j 2 ) = (20°+δ, 100°+ δ) are considered. The SNR = 25 dB is quoted per source per array element. Again the results are based on 500 independent trials. The DOA of the second source is varied as δ increases from 2°to 30°. For each angle  separation the proposed and the MUSIC algorithms are applied to obtain the DOA estimates of the two impinging signals. The CRBs and RMSEs of the estimate values are shown in Figure 8. For smaller separation angles the accuracy decreases dramatically. The mean bias is smaller than 1°for both two sources when separation angle is larger than 2°. Here, L = 13 for ε = 0.001 and the truncation error of the second kind of manifold decomposition can be nearly neglected for the proposed algorithm. The elevation angle estimation results for both algorithms are comparable. The truncation order is M = 7 for the first kind of manifold decomposition. Its corresponding truncation error for the second source becomes larger along with the increase of the separation angle δ (see Section 5). Such truncation error could impact the DOA estimation accuracy for the proposed algorithm. Moreover, the azimuth estimates for the proposed algorithm are obtained without the exact knowledge of the mutual coupling. For MUSIC algorithm, there is no truncation error and the mutual coupling is compensated accurately. Hence, the azimuth angle estimation results for the proposed algorithm is a little inferior to the ones for MUSIC here.

Conclusions
Several algorithms for DOA estimation with UCAs are based on the traditional beamspace transform, which requires a sufficiently large number of elements to avoid aliasing in the steering vector of the mode space. Sometimes there may be a smaller number of antenna elements for application. We propose a new approach to estimate 2D DOAs for such UCAs. Two kinds of manifold decompositions are applied as the foundation of the proposed algorithm. In the first step, a modified sparse UCA-RARE is performed for the azimuth estimates. This step can be realized without the exact knowledge of elevation angle. It is proved by means of the Jacobi-Anger expansion (a decomposition of the element manifold into phase modes) that the sparse UCA-RARE is still applicable with a limited number of array elements. In the second step, the Root-MUSIC algorithm is used to obtain the elevation estimates via decomposing the manifold with Bauer's formula (an expansion of the array manifold into a double Fourier series). The influence of the truncation errors on the DOA estimate accuracy is analyzed and a method to choose the truncated degree for the elevation estimates is presented. Simulation results show that the proposed algorithm for sparse UCA can obtain good azimuth angle and elevation angle estimate results. The next challenge is to find a computational efficient method to handle the sparse UCAs with much wider inter-element spacing.
(a) (b) Figure 6 The estimated roots distribution for the proposed algorithm: (a) roots of azimuth estimates and (b) roots of elevation estimates.