A New Conformal Map for Polynomial Chaos Applied to Direction-of-Arrival Estimation via UCA Root-MUSIC

The effects of random array deformations on Direction-of-Arrival (DOA) estimation with root-Multiple Signal Classification for uniform circular arrays (UCA root-MUSIC) are characterized by a conformally mapped generalized Polynomial Chaos (gPC) algorithm. The studied random deformations of the array are elliptical and are described by different Beta distributions. To successfully capture the erratic deviations in DOA estimates that occur at larger deformations, specifically at the edges of the distributions, a novel conformal map is introduced, based on the hyperbolic tangent function. The application of this new map is compared to regular gPC and Monte Carlo sampling as a reference. A significant increase in convergence rate is observed. The numerical experiments show that the UCA root-MUSIC algorithm is robust to the considered array deformations, since the resulting errors on the DOA estimates are limited to only 2 to 3 degrees in most cases.


Introduction
Direction-of-arrival (DOA) estimation is a fundamental enabler of current and nextgeneration wireless communication systems. With the arrival of 5G and the development of 6G, it is of great importance to understand how DOA techniques perform under imperfect conditions, especially as the effects of hardware impairments become more important when the operating frequency increases [1]. Simulation times can be large and Monte Carlo (MC) analysis is often too time-consuming due to the excessive number of realizations that must be evaluated. Therefore, a fitting stochastic framework for the characterization of uncertainty propagation within these systems is Generalized Polynomial Chaos (gPC) [2,3]. As a versatile technique, gPC has been applied extensively to study the effects of randomness on antenna performance and radio wave propagation [4][5][6][7][8][9][10][11]. With the appropriate approach, gPC can even compete with MC when a high number of random variables are used [12][13][14][15][16][17]. In the context of localization, gPC has also been used-for example, in [18,19], where the effects of element displacements on DOA estimation were investigated; in [20], where the effects of random element gain and phase variations were studied, and in [21], where the uncertainty in multiple angle-of-arrival measurements was translated into an uncertainty in position estimation.
Being based on polynomials, however, gPC has its own drawbacks, one of the most important ones being slow convergence in the presence of function singularities [22,23]. Conformal maps can alleviate these problems, as demonstrated on Maxwell's source problem in [24]. Conformal maps were used previously in [25] to compute more accurate quadrature rules for integrands with similar non-polynomial behavior, inspired by [26,27].
In this work, we extend the techniques described in [24] by introducing a novel map that compensates for (apparent) singularities on the real axis. We illustrate the effectiveness of this new conformal map by characterizing the effects of elliptical array deformations on root-Multiple Signal Classification for uniform circular arrays (UCA root-MUSIC) [28], as it is a well-established and efficient DOA estimation algorithm. The UCA topology is a relatively simple array geometry that enables DOA estimation of the azimuth angle over a full 360 • field of view, besides some more limited capability to provide estimates of the elevation angle. However, the dedicated UCA-root-MUSIC algorithm heavily relies on the inherent circular symmetry of the antenna array. Therefore, this algorithm and array topology are the ideal candidates to study the effects of random deformations on DOA estimation techniques. To be concise and keep the focus on the presence of singularities, we do not include white noise, in contrast to [18][19][20].
In the next section, key concepts are introduced and the novel conformal map is revealed in the appropriate mathematical framework. In Sections 3 and 4, the computational results are presented and discussed, and in the final section, Section 5, a conclusion is given.

gPC Approximation
Generalized Polynomial Chaos (gPC) [2,3] approximates a function f of a random variable x by an expansion in a well-chosen orthogonal polynomial basis {Φ k }. The orthogonal polynomials in question are associated with the weight function w, representing the probability density function of the random variable being used, typically according the Askey scheme [2]. If f is a function that is (computationally) expensive to evaluate, a gPC approximation of f can provide a computationally cheap substitute.
Assuming that w is a weight function with support [−1, 1], the gPC orthogonal projection of degree N is defined as Once the expansion coefficients t k are known, the first two moments of f can be calculated and estimated by while the standard deviation on f can be estimated by The error in the L 2 -norm · L 2 ,w can be written in terms of the standard deviation σ and its estimator: The quality of the gPC approximation will depend upon the analyticity of f in the complex plane [3,22,23]. The analyticity of f can be described by a Bernstein ellipse E ρ , which is defined as the ellipse with foci −1 and 1 and ρ as the sum of its semiminor and semimajor axis; see Figure 1. From Bernstein's theorem for the convergence of the Chebyshev projection P Ch N [22,23] and the fact that the gPC projection P N minimizes the L 2 -error [3], it follows that, if f is analytically continuable to the open Bernstein ellipse E ρ ⊆ C, the error on the Nth degree gPC approximation, according to the L 2 -norm, is limited to with · ∞ the supremum norm and | f (x)| ≤ M for x ∈ E ρ . It is clear from Equation (5) that it is preferable to have a high ρ, i.e., a large Bernstein ellipse, for optimal convergence. Singularities in the complex plane, close to the interval of interest [−1, 1], can substantially lower ρ and will, therefore, be detrimental to the gPC approximation.  Figure 1. The E 1. 6 Bernstein ellipse. Its size ρ is equal to the sum of its semiminor and semimajor axes-in this case, 1.6. Its foci −1 and 1 are shown by the black dots.
In most practical applications, the integral in Equation (1) is not analytically computable. In 1D cases, the coefficients t k are generally approximated numerically using a Gauss quadrature rule [3]: with (x i , w i ) being the N + 1 quadrature nodes and weights of polynomial accuracy 2N + 1 associated with w. According to [3], Section 3.6, the aliasing error on the polynomial expansion of f is given by In essence, the aliasing error follows from the fact that the lower-order polynomials {Φ k } cannot be distinguished from the higher-order polynomials {Φ j } on a finite grid. As seen in Equation (7), it can be interpreted as the error that is introduced by using the lower-order discrete expansion of the higher-order polynomials instead of the higher-order polynomials themselves. In [3], it is mentioned that the aliasing error induced by using Formula (6) is usually of the same order as the projection error in Equation (5). Hence, the aliasing error will also benefit from a higher ρ.

Conformally Mapped gPC
In [24], a framework was established to incorporate conformal maps into the gPC algorithm for enlarging the Bernstein ellipse, based on earlier research from [25]. Consider a map g that is conformal in an open region Ω ⊆ C with subdomain [−1, 1] ⊆ R, with g([−1, 1]) = [−1, 1] and g(±1) = ±1. This map can be used to define a new variablex: The symbol • is used to denote function composition, i.e., By choosing g > 0 in Ω, the above equation becomes Assuming that the orthogonal polynomials {Φ k } associated with w are known, the orthogonal polynomials {Φ k } associated withw can be constructed by the Modified Chebyshev Algorithm [29]. Next, the quadrature pointsx i and associated weightsw i ofw can be calculated using the Golub-Welsch Algorithm [30]. The Modified Chebyshev Algorithm is based upon a set of integrals called the modified moments, defined by Only for select combinations of weight function and conformal map can these integrals be calculated analytically. In other cases, one has to make use of quadrature rules (or other computational methods) to evaluate these integrals numerically. As the quadrature nodes and weights belonging tow are not yet known at this point in the algorithm, this integral has to be reformulated as an integral with weight function w, whose quadrature nodes and weights are known. One means of achieving this is by rewriting the integral as with (x i , w i ) the quadrature nodes and weights of polynomial accuracy 2N m − 1 associated with w. Within this research, the value of N m was set to 50, a value for which sufficient accuracy was reached.
Instead of a polynomial expansion of f , in the conformally mapped gPC algorithm introduced by [24], an expansion of f • g is performed, using the newly constructed Φ k polynomials: By substituting the inverse mapx = g −1 (x) at both sides of Equation (12), a mapped approximation of f is obtained: A schematic overview of this algorithm is shown in Scheme 1.
x-spacex -space Scheme 1. An overview of the conformally mapped gPC algorithm. A new random variablex is defined by applying a conformal map x = g(x) on the random variable x. An adapted quadrature rule (x i ,w i ) is constructed for the new weight functionw(x) via the modified moments m k , along with its own orthogonal polynomials {Φ k }. Generalized Polynomial Chaos (gPC) is applied to the composite function [ f • g] and a mapped gPC approximation of f is found after performing the inverse conformal mapx = g −1 (x).
Analogously to classic gPC, the coefficientst k can be approximated by means of discrete projection, yielding the quadrature nodes and weights of polynomial accuracy 2N + 1 associated withw.
An upper bound similar to the one in Equation (5) can be derived for the conformally mapped gPC expansion [3,[22][23][24]: In this case,ρ is the size of the open Bernstein ellipse Eρ, in which f • g is analytically continuable. It is apparent from Equations (5) and (15) that the conformal map g needs to be chosen in such a way thatρ exceeds ρ, thus achieving a faster convergence. In other words, f • g needs to be analytically continuable in a larger open Bernstein ellipse than f .
To the best of the authors' knowledge, only one class of conformal map has been applied in the context of polynomial-based methods. This set of maps shifts singularities directly above or below the [−1, 1] interval, away from the origin, in order to achieve a larger Bernstein ellipse. The Ellipse-to-Strip map [25], the Sausage map [25], the Kosloff Tal-Ezer (KTE) map [27] and the Ellipse-to-Slit map [31] all fall into this category. In Figure 2, examples are shown of these established maps. However, it is possible to encounter function singularities in other locations of the complex plane, such as on the real axis, close to the [−1, 1] interval. In these situations, there is a need for a new conformal map, which is proposed further in Section 2.3.

Tanh Map
Assume that f has a singularity on the real axis at location p with |p| > 1. Due to this singularity, the size of the Bernstein ellipse E ρ associated with f (in the x-space) is limited to ρ = |p| + p 2 − 1. Ideally, one would use a map that gives rise to a Bernstein ellipse Eρ, associated with f • g in thex-space, that is larger than E ρ . It is clear that a suitable map for this purpose should shift the singularity away from the [−1, 1] interval. We propose the tanh map for this purpose: An illustration of this map for a real x andx, and for different values of κ, is shown in Figure 3. The parameter κ is added to make the map adaptable to different values of the singularity's location p.
Trivial map κ = 1.5 κ = 2 κ = 2.5 The tanh(z) function is periodic with period πj, i.e., tanh(z) = tanh(z ± πj), and has simple poles at π 2 j ± πjl, with l ∈ Z. As the map needs to be bijective in order to define the inverse map, the domain of the above map is restricted to the strip around the real axis between its closest poles: The inverse of Equation (16) can now be defined as Using the tanh map, the singularity will shift to a position |p| = |g −1 (p; κ)| > |p|. As g −1 (p; κ) needs to be defined and g −1 (x; κ) has branch cuts −∞, − 1 tanh(κ) and 1 tanh(κ) , +∞ , κ is fundamentally limited to Two Bernstein ellipses in thex-space can be defined: one that only takes into account the shifted singularityp with size (20) and one ellipse that only takes into account the singularities ± π 2κ j introduced by the tanh map, with sizeρ Depending on the values of p and κ, either Eρ re or Eρ im will be the largest ellipse. As the presence of singularities is restrictive for the gPC algorithm,ρ is equal to the size of the smallest of both ellipses, i.e.,ρ = min(ρ re ,ρ im ). The largest and optimal value of ρ, namedρ eq , is reached when both ellipses overlap (ρ re =ρ im ) at a certain κ eq , which can be found by solving Equations (20) and (21). Figure 4 illustrates this procedure for a singularity located at p = ±1.1125. The value of κ eq as a function of the singularity position |p| is shown in Figure 5, along with the two regions in (|p|, κ)-space, one whereρ re ≤ρ im and one whereρ re ≥ρ im .

Re (d)
Eρ re Eρ im Eρ eq Figure 4. The Bernstein ellipse with size ρ = 1.6 (singularity located at p = ±1.1125) in the x-space (a) and the corresponding Bernstein ellipses Eρ re and Eρ im when using κ = 1.05κ eq (b), κ = 0.95κ eq (c) and κ = κ eq (d). The singularities are depicted with hollow dots. Since Eρ eq has the singularities of g, being ± π 2κ j, on its border, it will map to an infinitely large region g(Eρ eq ; κ eq ) in the x-space, as can be seen in Figure 6. This is similar to the Ellipse-to-Strip [25] and the Ellipse-to-Slit maps [31]. The condition that f is analytically continuable within this infinite region is rather strict, and if this is not the case, the effective Bernstein ellipse will be smaller than Eρ eq . Luckily, the region of analyticity for f will shrink very quickly in comparison to the corresponding Bernstein ellipse Eρ. Additionally, when comparingρ with ρ in Figure 7, one can conclude that, in practice, as long as singularities above and below the [−1, 1] interval are reasonably far away, the convergence rate gain does not suffer much when applying the tanh map. ρ eq 0.95ρ eq , 0.9ρ eq , 0.85ρ eq , 0.8ρ eq Figure 6. Starting with a Bernstein ellipse E 1.6 in the x-space, the Bernstein ellipses in thex-space with sizesρ eq = 2.72, 0.95ρ eq = 2.59, 0.9ρ eq = 2.45, 0.85ρ eq = 2.31 and 0.8ρ eq = 2.18 are shown (b). The corresponding regions g(Eρ; κ eq ) in which analytic continuability of f is assumed are shown in (a). The original Bernstein ellipse E 1.6 is shown with a dashed line as a reference and the singularities are depicted with hollow dots. Size of Bernstein ellipse ρ Size of Bernstein ellipseρ˜ρ eq 0.95ρ eq , 0.9ρ eq , 0.85ρ eq , 0.8ρ eq Figure 7. The size of Eρ as a function of the size of E ρ when using the tanh map.ρ eq corresponds to the maximally achievable value ofρ.

Setup
Consider a nine-element uniform circular array (UCA) deployed in the azimuth plane. As the effects of the displacement of individual antenna elements on root-MUSIC have already been studied in earlier research-for example, in [18,19]-to showcase the proposed method, we assume that the array deforms in the shape of an ellipse with constant circumference and constant distance between the elements along the elliptical arc. This deformation is characterized by a single (random) variable: an "extended" eccentricity e defined by with 2a being the width of the array (along φ = 0 • ) and 2b the height (along φ = 90 • ). Different array deformations are illustrated in Figure 8. We limit the eccentricity to the interval [−0.9, 0.9], since values outside of this interval were deemed too unrealistic and the deformation would become too large for the DOA algorithm to provide a DOA estimation. However, it is standard practice [24,25] to rescale random variables to [−1, 1]. Therefore, the random variable x is introduced as e = 0.9 · x.
Since the UCA root-MUSIC algorithm was calibrated with a circular array in mind, as the array deforms, the estimated DOAs will deviate from their correct positions. This error propagation is described by a functionφ = f (x), withφ being the DOA estimation corresponding to one specific source. From a practical viewpoint, it is non-intuitive to consider the analytic continuation of f (x) in the complex plane beyond the [−1, 1] interval. However, one can state that f , the relation between deformation and (erroneous) DOA estimation, behaves in a highly erratic manner when x approaches the edges of the interval [−1, 1]. This behavior is equivalent to the presence of nearby singularities on the real axis. Therefore, the principles of Sections 2.1-2.3 apply to this case, even though f (x) has no closed form and the analytic continuation is not known.

Results
The operating frequency of the considered antenna array is set to 3.5 GHz (wavelength λ ≈ 8.57 cm), being the center of a 5G band [34]. The dipole elements are all of length λ/2, as is the radius of the UCA. The array is excited by six plane waves along directions φ = 50 • , 70 • , 165 • , 220 • , 305 • and 350 • in the azimuth plane. DOA estimation is performed with the UCA root-MUSIC algorithm [28]. The full-wave NEC2++ simulator [35] is applied to rigorously simulate the complete antenna array, including all mutual coupling effects. For conciseness, we only discuss the behavior of one of the six DOA estimates, being the one corresponding to the source at φ = 50 • .
In Figures 10-12, the absolute and relative errors of µ and σ are presented as a function of N for the different Beta distributions. These values are calculated with classic gPC and mapped gPC for varying values of κ. The expansion coefficients are approximated with discrete projection according to Equations (6) and (14). As µ is equal to t 0 andt 0 , the error shown in subfigures (c) is only the aliasing error introduced by the discrete projection. The error on the estimation of σ in subfigures (d) can, in addition to the aliasing error, be directly linked to the L 2 -error according to Equation (4). As there are no analytical solutions to compare the results with, reference values are computed numerically by means of Monte Carlo simulation, using Latin Hypercube Sampling (LHS) with 10 6 samples [36]. These values are displayed in Table 1. The implemented UCA root-MUSIC algorithm has a precision of around 10 −6 degrees, which is why, in some graphs, convergence halts at a relative error of around 2 × 10 −8 and 10 −6 for µ and σ, respectively.    In Figure 13, a comparison between the resulting classic and mapped gPC approximations of f , using the Beta(x; 2, 2) distribution and N = 15, is shown. In Figure 14, the comparison of the resulting empirical cumulative distribution functions (CDFs) is plotted.

Discussion
The discussion is presented in two parts. First, the advantages of the use of the tanh map are analyzed, based on the deformed UCA application. Afterwards, the effects of the random array deformations on the UCA root-MUSIC algorithm are discussed.

Comparison of Classic and Mapped gPC
In Figures 10-12, we see a general improvement when using mapped gPC over classic gPC, as, in all cases, mapped gPC reaches the precision bound the fastest. Comparing the results for the estimation of µ, in subfigures (a) and (c), in which only the aliasing error is present, we can confirm that the aliasing error does indeed benefit from an increase in the size of the Bernstein ellipse. As expected, this increase in convergence rate is also present in the results for the estimation of σ in subfigures (b) and (d), which affirms the principles from Sections 2.2 and 2.3. Note that the error on σ is linked to the L 2 -error in Equations (5) and (15) via Equation (4).
One aspect that should be mentioned is that, according to Equations (5) and (15), maximizing the size of the Bernstein ellipse only maximizes the rate of convergence, i.e., the incline of the convergence curves in subfigures (c) and (d). The vertical position of the convergence curves will, however, depend on other factors besides the size of the Bernstein ellipse. Therefore, it is possible that the fastest-converging method is not the one with the smallest error, especially in the cases with Beta(x; 2, 2) and Beta(x; 3, 3) as weight functions, where the precision floor is reached relatively quickly.
Two factors influence this phenomenon: first, the M/M parameter in Equations (5) and (15), which is an upper bound of | f | in its supposed region of analyticity, being either E ρ for classic gPC or g(Eρ; κ) for mapped gPC. Although it is difficult to establish closed-form mathematical relations for the value of M/M, for the mapped gPC case, one can state that an increase in κ will causeM to either increase or stay the same, as clarified in Equation (24). In other words, an increase in convergence rate due to an increase in κ can be paired with an upward vertical shift in the convergence curve.
Another factor is the aliasing error, defined by Equation (7), which is a function of the used weight function w/w and the expansion polynomials {Φ k }/{Φ k }, adding a dependency on w and κ. An exact evaluation of Equation (7) is difficult. Upper bounds to the aliasing error are available for, among others, the Legendre and Chebychev polynomial expansions [37]; however, these are not readily applicable to this mapped gPC context. The dependency of both these factors on κ and w explains why we see different performance for the different κ values as the weight function changes, even though f and its singularities, and therefore also κ eq , stay the same. Unfortunately, it is difficult to establish closed-form mathematical relations for these dependencies. Figure 13 shows that a better fit is achieved at the extremities of the interval when using mapped instead of classic gPC, resulting in a lower supremum error and L 2 -error. As this erratic behavior of f in these regions has a large influence on the statistical moments of the function, a better fit at the edges will have a significant impact on the accuracy of the estimation of µ and σ. Another benefit of the mapped approach is a better approximation of the CDF at the far left and far right sides, as seen in Figure 14.

Consequences for the UCA Root-MUSIC Algorithm
In the situations studied in this paper, the introduction of a random deformation of the array causes a bias in the DOA estimator of 0.05 to 0.5 degrees and a standard deviation of the DOA estimation of 1 to 2.5 degrees (depending on the distribution shape; see Table 1). All things considered, this makes the UCA root-MUSIC algorithm rather robust against the array deformations studied in this work.

Conclusions and Future Work
The non-polynomial behavior of the root-MUSIC DOA estimation as a function of the elliptical deformation of the UCA can be compared to the presence of singularities on the real axis, close to the interval of interest, which have a detrimental effect on the convergence of the classic gPC algorithm. Luckily, using the newly defined tanh conformal map, these singularities are moved further away from the domain of the random variable, which makes for a better characterization of the erratic behavior of the DOA estimation and, as a result, causes a considerable increase in the convergence rate of the first-and second-order statistics in comparison to classic gPC.
We conclude from the simulations that the errors induced in the DOA estimation due to the elliptical deformation of the UCA are limited to only a few degrees in most cases. However, when the eccentricity reaches an absolute value of around 0.7, the DOA estimations become very volatile, with much larger errors of up to 10 degrees.
In future work, it can be of interest to develop and research even more specific conformal maps so that each type of function singularity can be dealt with in an efficient manner. As for the DOA estimation in 5G and 6G wireless communication networks, it might be advisable to look at other types of antenna arrays, such as rectangular arrays, which are currently integrated into 5G base stations and handsets [38]. Additionally, it could be interesting to look at other, modified types of MUSIC, such as or spatial/backward/time smoothing MUSIC, which are better equipped to deal with highly correlated signals, as encountered due to multipath propagation in indoor environments [39]. The technique could also be extended to hybrid techniques, such as SpotFi, that rely on a combination of timeof-flight and angle-of-arrival with MUSIC to perform accurate localization with common WiFi infrastructures [40].