Iterative Calibration and Direction-of-Arrival Estimation for Uniform Circular Arrays Affected by Mutual Coupling

This letter presents an iterative method to estimate directions-of-arrival (DOA) of plane wave signals impinging on a uniform circular array (UCA), accounting for the possible presence of mutual coupling (MC). After a tailored transformation that maps the data gathered by the UCA into a virtual uniform linear array (VULA), the mutual coupling effect is modeled via gain and phase errors affecting the resulting VULA. Therefore, a coordinate descent (CD) iterative procedure is designed to approximate the maximum likelihood estimates of the VULA calibration factors and of the sources power. Leveraging the resulting calibrated VULA array manifold, the DOAs are obtained using the estimation of signal parameters via rotational invariant technique (ESPRIT). Simulation results demonstrate the superior performance of the proposed method with respect to the baseline uncalibrated ESPRIT, in terms of DOA estimation accuracy in the presence of mutual coupling, making it viable to be implemented in actual sensing systems due to its moderate computational complexity.


I. INTRODUCTION
The problem of estimating the directions-of-arrival (DOA) of plane wave signals impinging on an antenna array is a long-standing topic in sensing systems and represents the subject of extensive research [1], [2], [3], [4], [5].Notably, in scenarios where a 360 • field of view is required with possible extraction of elevation information, the uniform circular arrray (UCA) is often the preferred option.
The majority of the signal processing architectures presented in the open literature assume an exact knowledge of the actual antenna array manifold.In practice, this assumption is rarely satisfied, due for instance (and primarily) to mutual coupling, which refers to the alteration of the electromagnetic characteristics of each array element as a result of the leakage from adjacent radiating elements [6], [7], [8].The presence of mutual coupling can have a significant impact on the array sensing performance, degrading the resolution capability, robustness to the interference of adaptive algorithms, and, needless to say, accuracy of DOA estimation [6], [9], [10], [11].This demands the development of adaptive and accurate methods capable of mitigating the performance loss induced by mismatch effects.
Motivated by the above reasons, in this letter, a high-resolution iterative technique for the estimation of DOAs is investigated for a UCA accounting, at the design stage, for the possible presence of mutual coupling.First of all, a phase-mode excitation-based transformation [12], [13], [14] is employed, which allows modeling the collected signal as received by an equivalent virtual uniform linear array (VULA).Then, following [14], the mutual coupling effect is modeled as gain and phase errors affecting the VULA.In this context, the main contribution of this letter is the development of a coordinate descent (CD) iterative procedure to approximate the maximum likelihood (ML) estimates of the calibration factors of the VULA, i.e., the mutual coupling coefficients, as well as the ML estimates of the sources amplitude.Leveraging the resulting calibrated VULA array manifold, accurate DOA estimates are obtained by resorting to the estimation of signal parameters via rotational invariant technique (ESPRIT) [15].Simulation results show that the proposed method is capable of achieving better performance, in terms of DOA estimation accuracy, with respect to the standard UCA-ESPRIT, which neglects the presence of mutual coupling.
Notations: Boldface is used for vectors a (lower case), and matrices A (upper case).I denotes the identity matrix, whose size is determined from the context.Moreover, diag(x) indicates the diagonal matrix whose ith diagonal element is the ith element of x.The transpose and the conjugate transpose operators are denoted by the symbols (•) T , and (•) † , respectively.C N is the set of N-dimensional column vectors of complex numbers.The letter j represents the imaginary unit (i.e., j = √ −1).For any complex number x, |x| indicates the modulus of x.For any x ∈ C N , x denotes its Euclidean norm.

II. PROBLEM FORMULATION
Let us consider a narrowband UCA sensing system operating at wavelength λ.Specifically, the array is composed of M omnidirectional antennas arranged on a circle of radius r, with the array elements lying on the plane z = 0 at angular positions γ m = 2π m/M, m = 1, . . ., M, with respect to the x-axis (see Fig. 1).Let us assume that the UCA collects L snapshots from an environmental scenario comprising K < L narrow-band, far field, uncorrelated sources s k impinging on the array from azimuth directions θ k , with k = 1, . . ., K, whose elevations, in this study, are set to 0 degrees.Accounting for the presence of mutual coupling, after radio frequency chain processing and digital sampling, the received signal collected at the lth snapshot can be expressed as [12] where 1) C ∈ C M×M is the mutual coupling matrix (MCM); 2) A(θ) = [a(θ 1 ), . . ., a(θ k ), . . ., a(θ K )] ∈ C M×K is the ideal antenna array manifold matrix, associated the K impinging sources; 3) a(θ k ) = [e j(2π r/λ) cos(γ 1 −θ k ) , . . ., e j(2π r/λ) cos(γm−θ k ) , . . . . . ., e j(2π r/λ) cos(γ M −θ k ) ] T ∈ C M , k = 1, . . ., K is the spatial steering vector for the kth source; 4) s l = [s (1)  l , . . ., s (k) l , . . ., s (K ) l ] T ∈ C K , l = 1, . . ., L are the vectors of (possibly random) complex amplitude signals of the K sources at the lth snapshot, modeled at the calibration procedure design stage as unknown deterministic parameters; 5) n l ∈ C M , l = 1, . . ., L are independent and identically distributed (IID) zero-mean circularly symmetric Gaussian random vectors with mean square values σ 2 n , assumed statistically independent from the sources.In practice, the mutual coupling phenomenon induces a mismatch between the actual and the ideal steering matrix, described by the matrix C. Formally, the actual steering matrix Ã(θ) is given by Ã Interestingly, for an UCA, in absence of additional impairments and asymmetries, the MCM can be well described as a complex symmetric circulant Toeplitz matrix [14], [16].Specifically, due to the UCA geometric symmetries, the MCM C presents only ρ = (M + 1)/2 degrees of freedom [16].Therefore, it can be completely described by the complex elements of its first column In practical scenarios, the MCM is unknown and demands its ad-hoc estimate.Furthermore, even in conditions where an offline calibration could be accomplished, bespoke signal processing architectures encompassing an array autocalibration are actually required, since the mutual coupling effect is not necessarily stationary (for instance, due to changes in the environmental conditions) [17].

A. VULA Array Manifold Description
Let us now consider the well-known phase-mode excitation-based transformation [12], [13], which allows to describe (with a good approximation) the UCA antenna array manifold matrix A(θ) in terms of a beamspace manifold B(θ) with a significant structure modeled according to a typical ULA signature.In the ideal case with no mutual coupling effects, [13] where is the beamformer matrix of the mode domain transformation [14], with Q ≤ (M − 1)/2 the maximum mode order of the transformation of interest, J m the Bessel function of the first kind of order m, evaluated at (2π r/λ), = 2Q + 1 the number of elements in the equivalent VULA and w h = [e − jh2π/M , . . ., e − jh2 mπ/M , . . ., e − jh2π ] T .
Following [14], after the appropriate transformation of the UCA manifold into a VULA via (2), it is possible to consider the MC effect as gain and phase errors affecting the VULA observations.Specifically, the actual beamspace manifold of the VULA is given by B(θ) = JWCA(θ) ∈ C ×K , which is tantamount to considering the VULA useful returns as affected by gain and phase errors [14], i.e., with H = ∈ C × being the diagonal complex matrix of the coupling induced error terms (see [14, Sec.III.A] for the expressions of the corresponding gains and phases matrices).

III. DOA ESTIMATION IN THE PRESENCE OF MUTUAL COUPLING
Leveraging the aforementioned transformation, the lth received snapshot x l can be modeled as received by the VULA as [13], [14] with H = diag(h) and ñl the resulting noise term pertaining to the lth snapshot after the beamspace transformation.Therefore, after standard argumentation [18], the ML estimates of H (tantamount to estimate and ), S, and θ can be obtained solving min where S = [s 1 , . . ., s L ] T ∈ C K×L .Since the objective function in ( 5) is not convex, the computation of the ML estimates represents a quite challenging task for which an analytic closed-form solution is not available.In addition, due to the high number of unknowns, methods based on multidimensional grid search are computationally infeasible.Therefore, to ensure analytic tractability with a reasonable computational cost, it becomes mandatory the design of approximated solution methods.In this respect, the idea pursued in this letter is to start from an initial estimate of the unknowns and then to refine it via an alternating optimization procedure over different blocks of variables.Specifically, the output of the standard UCA-ESPRIT method [13], i.e., the vector of the DOAs estimates, denoted by θ(0) , is used to initialize the proposed procedure along with the complex gain matrix set to Ĥ(0) = I.Hence, leveraging θ = θ(0) , a refined estimate of the error terms can be achieved by solving the following least square problem: min This formulation is suitable for an iterative optimization based on the CD framework, whereby it alternates between the minimization over Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

Algorithm 1: MC-UCA-ESPRIT.
Input: {x l } l=1,...,L , J, W , 1 , 2 .1) Compute θ (0) via standard UCA-ESPRIT method [13]; 2) Transform the received snapshots as z l = JW x l , l = 1, . . ., L; Compute ŝ(n+1) l via (7) using Ĥ = Ĥ(n) , l = 1, . . ., L; 8) Compute ĥ(n+1) via ( 9) using Sl = diag(B( θ)ŝ (n+1) l ), l = 1, . . ., L; 9) Set n = n + 1; 10) Compute h and S, namely, optimizing one block variable at a time while keeping the other fixed.Given the starting point of the iterations, since each marginal problem admits a unique solution, it follows that the final point (after any finite number of iterations) of the overall CD procedure is uniquely determined.Although the nonconvex nature of the problem, it is possible to claim the stationarity of any cluster point provided by the indefinitely long sequence of points [19], [20].First, by fixing h = ĥ, (6) reduces to min S L l=1 z l − ĤB( θ)s l 2 , with Ĥ = diag( ĥ), whose optimal solution is given by ŝl = B( θ) † Ĥ † ĤB( θ) Then, leveraging (7), the CD method demands the solution to the following minimization problem: whose optimal solution is given by ĥ where Sl = diag(B( θ)ŝ l ), l = 1, . . ., L. Thus, ( 7) and ( 9) are repeated until convergence (guaranteed by Wright [21]) is achieved.Now, leveraging the estimate of h achieved by the CD-based procedure, respectively, an updated and refined estimate of the DOAs can be actually computed using the standard UCA-ESPRIT method [13] applied to the calibrated sample covariance matrix of the snapshots where δ is an estimate of the noise mean square value obtained averaging the smallest (M − K ) eigenvalues of 1 L L l=1 x l x † l .Algorithm 1 summarizes the devised CD-ESPRIT-based solution technique, referred to in the following as MC-UCA-ESPRIT, whereby 1 > 0 and 2 > 0 are the user-defined exit thresholds for the inner and the outer cycles, respectively.

IV. PERFORMANCE ANALYSIS
In this section, the estimation performance of the proposed method is analyzed assuming an UCA with M = 7 elements, λ/R = 5, and two equal-power zero-mean circularly symmetric Gaussian random sources with common mean square value |s| 2 , located at θ 1 = 175 • and θ 2 = 45 • with respect to the reference x-axis, respectively.To assess the capabilities of the devised procedure, the root-mean-square error (RMSE) of the estimates is used as figure of merit, computed via standard Monte Carlo counting techniques over 1000 independent trials.Moreover, the results are compared with the standard UCA-ESPRIT [13] procedure, which does not account for the mutual coupling effect.For the execution of Algorithm 1, the exit conditions are set as 1 = 10 −12 and 2 = 10 −5 , yielding a typical number of iterations equal to 40 for the innermost cycle and 3 for the outermost one.Two scenarios are analyzed.In the former, the system operates in an ideal environment, i.e., without the presence of mutual coupling.In the latter scenario, the presence of mutual coupling effect is accounted for.
Ideal case without mutual coupling: Without mutual coupling, the data covariance matrix is given by M = M S + σ 2 n I, where σ 2 n is the white noise power level, assumed (without loss of generality) equal to 0 dB, and M S = |s| 2 K k=1 a(θ k )a(θ k ) † , is the signal covariance contribution.
The results, in terms of RMSE for θ 1 versus the signal-to-noise ratio (SNR), defined as SNR = |s| 2 /σ 2 n , are depicted in Fig. 2(a), for three different values for the number of snapshots, i.e., L = 100, L = 300, and L = 1000.
An inspection of the figure reveals that the standard UCA-ESPRIT outperforms the proposed method, with a gap between the curves in the order of 6 dB at RMSE = 1 • for L = 100 and L = 300.This is not surprising, as the actual signal model matches the (ideal) conditions assumed at the UCA-ESPRIT design stage, i.e., absence of mutual coupling effect.Nevertheless, the proposed method is still capable of providing adequate performance, assuming that there is a sufficient number of snapshots, which allows to accurately estimate the data covariance matrix.In addition, in the high SNR regime (SNR ≥ 15 dB), the devised method provides accurate DOA estimates in all the analyzed cases, with RMSE values less than 1 • .Finally, as expected, as the number of observed data L increases, both methods yield better estimation performance.
Actual case with mutual coupling: Assuming the same system parameters as in the previous case, the effectiveness of the proposed estimation strategy is now studied considering the presence of mutual coupling, characterized with c 1 = [1, 0.6 − 0.56 j, 0.42 + 0.43 j, 0.23 − 0.27 j, 0.23 − 0.27 j, 0.42 + 0.43 j, 0.6 − 0.56 j] T .To this end, the signal covariance matrix is given by M S = |s| 2 K  k=1 Ca(θ k )(Ca(θ k )) † .The resulting RMSE curves, displayed in Fig. 2(b), show that the presence of mutual coupling induces a performance degradation of the standard UCA-ESPRIT, resulting in an estimation accuracy loss, at SNR = 15 dB and for L = 1000, of 1.3 • .Moreover, the UCA-ESPRIT exhibits a performance saturation in the high SNR regime, reaching RMSE values, for L = 1000, of approximately 1.5 • .Conversely, the MC-UCA-ESPRIT method is able to yield adequate estimation performance, achieving lower RMSE values than its counterpart in all the analyzed cases when the SNR is sufficiently high, with an accuracy improvement of 1 • at SNR = 15 dB for L = 300 and L = 1000.Moreover, for a given number of snapshots, the proposed estimation algorithm provides RMSE values that become smaller and smaller as the SNR increases, corroborating the effectiveness of this strategy to cope with the mutual coupling effect.Finally, results similar to those depicted in Fig. 2 are attained for the estimates related to θ 2 .

V. CONCLUSION
This letter has addressed the estimation of the source DOAs in the presence of mutual coupling among the array elements of a UCA.Leveraging a bespoke transformation that allows modeling the collected data as gathered by a VULA, the mutual coupling effect has been described in terms of gain and phase error factors impairing the "equivalent" ULA elements.To effectively calibrate the VULA, an approximate ML estimator of the unknowns has been designed.The involved optimization problem has then been tackled via a CD-based procedure, which is a fast-converging technique able to provide monotonic decreasing objective function values.Finally, using the calibrated VULA array manifold, the sources DOA have been estimated using the ESPRIT method.Simulation results have demonstrated the capabilities of the proposed method to yield accurate DOA estimates in the presence of mutual coupling.
Possible future research studies might include the extension of the framework to also provide sources elevation estimation, the design of alternative estimation strategies via compressed sensing techniques [22], as well as the extension of the proposed approach to the case of spherical antenna arrays [11].

Fig. 2 .
Fig. 2. Estimation performance w.r.t.θ 1 .RMSE curves for the cases without and with mutual coupling are reported in (a) and (b), respectively.