Bayesian In-Situ Calibration of Multiport Antennas for DoA Estimation: Theory and Measurements

The DoA of radio waves is used for many applications, e.g. the localization of autonomous robots and smart vehicles. Estimating the DoA is possible with a multiport antenna, e.g. an antenna array or a multi-mode antenna (MMA). In practice, DoA estimation performance decisively depends on accurate knowledge of the antenna response, which makes antenna calibration vital. As the antenna surroundings influence its response, it is necessary to measure the entire device with installed antenna to obtain the installed antenna response. Antenna calibration is often done in a dedicated measurement chamber, which can be inconvenient and costly, especially for larger devices. Thus, auto- and in-situ calibration methods aim at making antenna calibration in a measurement chamber redundant. However, existing auto- and in-situ calibration methods are restricted to certain antenna types and certain calibrations. In this paper, we propose a Bayesian in-situ calibration algorithm based on a maximum a posteriori (MAP) estimator, which is suitable for arbitrary multiport antennas. The algorithm uses received signals from a transmitter, noisy external DoA observations, takes multipath propagation into account and does not require synchronization. Furthermore, we take an estimation theoretic perspective and provide an in-depth theoretical discussion of in-situ antenna calibration in unknown propagation conditions based on Bayesian information and the Bayesian Cramér-Rao bound (BCRB). Extensive simulations show that the proposed algorithm operates close to the BCRB and the achieved DoA estimation performance asymptotically approaches the case of a perfectly known antenna response. Finally, we provide an experimental validation, where we calibrate the antenna on a robotic rover and evaluate the DoA estimation performance using measurement data. With the proposed in-situ antenna calibration algorithm, DoA estimation performance is considerably improved compared to using an antenna response obtained by simulation or in a measurement chamber.


I. INTRODUCTION
Direction-of-arrival estimation of radio signals, also called radio direction finding, has a long history in both, research and application [1], [2]. For instance, the DoA is considered for localization of user equipment in cellular networks [3], [4], for cooperative localization in wireless networks [5], [6] and for autonomous driving [7]. Another application The associate editor coordinating the review of this manuscript and approving it for publication was Chinmoy Saha .
where DoA proves useful is robotic exploration in scenarios where global navigation satellite system (GNSS) signals are not available, e.g. extraterrestrial exploration. In order to control the robots, their positions and orientations are required [8], [9]. Positions and orientations can be estimated by using the time-of-arrival (ToA) and DoA of radio signals [10].
On the antenna side, usually antenna arrays are considered for DoA estimation [11]. More recently, DoA estimation with a single MMA has been proposed [12], [13]. An MMA consists of a single antenna element, where VOLUME 10, 2022 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ multiple characteristic modes [14], [15] are excited independently [16], [17]. Thus, an MMA is a multiport antenna, where each port has a distinct far-field antenna response. The different far-field antenna responses of the ports enable DoA estimation of impinging signals. Examples of MMAs can be found in [17]- [19]. MMAs are especially promising for applications with stringent size or shape constraints. For instance, the structure of a vehicle can be used for radiation using characteristic modes [20]. Further examples are found in [21], [22], where the characteristic modes of an airplane structure are used for DoA estimation.
On the signal processing side, a wide range of DoA estimation algorithms exist, see [1], [2], [23] and the references therein. DoA estimation algorithms for MMAs are treated in [12]. However, all of these algorithms have in common, that they rely on exact knowledge of the antenna response. Any mismatch between the assumed antenna model and the real one will impair DoA estimation performance [24].
For that reason, antennas are calibrated, meaning their antenna response is determined. Usually, antenna measurements are performed in a controlled environment without multipath propagation, e.g. a dedicated measurement chamber. The measured antenna response can then be used for DoA estimation using wavefield modeling and manifold separation [25], [26], which provides a mathematical representation of arbitrary antenna responses. However, as the antenna surroundings influence the antenna response, it is necessary to measure the entire device with installed antenna in order to obtain the installed antenna response. This permits the usage of compact and cost-effective near-field measurement chambers. For antennas integrated in e.g. robots, vehicles or airplanes, antenna calibration in a measurement chamber can become costly and impractical.
An attempt to avoid antenna calibration in a measurement chamber are auto-or self-calibration methods, which aim at estimating the antenna parameters together with the wavefield or DoA. However, they suffer from a severe limitation. In general, both wavefield and antenna parameters are not simultaneously identifiable from a collection of received signal snapshots, unless strong assumptions are made [27], [28]. For instance, assuming an ideal antenna array, it is possible to calibrate the gain-phase error of the individual antenna elements [29], [30]. An ideal antenna array consists of omnidirectional or isotropic antenna elements and its antenna response is defined by the steering vector, i.e. the geometry of the array. Restricting the antenna to an ideal uniform linear array (ULA) or uniform circular array (UCA), the mutual coupling matrix, i.e. a linear transformation of the antenna response, can be determined [31], [32]. Unifying gain-phase, mutual coupling and antenna element position calibration, a MAP formulation for auto-calibration can be found in [33], and low complexity variants in [33], [34]. They assume, however, that the perturbations of the antenna response are small or at least not larger than the finite sample errors. The method from [35] allows large antenna element position errors, but is still limited to an ideal antenna array.
More recently, sparsity-based approaches have been proposed [36], [37], which outperform conventional approaches. However, they are limited to ideal, linear arrays. In practice, errors of the antenna response are not necessarily limited to gain-phase errors, mutual coupling and antenna element position errors. Instead, arbitrary nonlinear transformations of the ideal antenna response can occur, e.g. due to manufacturing imperfections and surrounding structure of the installation [24]. Moreover, the known methods are restricted to ideal antenna arrays and cannot be applied to other multiport antennas like an MMA.
In contrast to auto-calibration, in-situ calibration methods use transmitters in (approximately) known directions to calibrate the antenna. In-situ calibration is performed in unknown propagation conditions and without external synchronization. In [38], in-situ calibration of mutual coupling and antenna element positions of an array is introduced. Additionally, also the power patterns of the antenna elements of an array can be calibrated [39]. Mutual coupling and antenna element positions of multiple antenna arrays, which are part of a localization system, are calibrated in [40]. The in-situ calibration methods [38]- [40] are restricted to antenna arrays and specific error models.
Thus, an MMA cannot be calibrated by the known auto-calibration and in-situ calibration methods from literature. Furthermore, also for antenna arrays, arbitrary nonlinear errors of the antenna response can occur, e.g. due to manufacturing imperfections or the influence of the antenna surroundings. The antenna response errors are not restricted to gain-phase errors, mutual coupling or antenna element position errors.
The aim of this paper is to investigate in-situ calibration of arbitrary multiport antennas. Applying wavefield modeling and manifold separation, the in-situ calibration of antenna arrays, collocated antennas, MMAs and other antenna types is treated in a common framework. By this paper, we considerably extend our early work [41]. The main contributions are as follows: • We analyze in-situ calibration of multiport antennas in an estimation theoretic approach. We derive the recursive BCRB, where we distinguish the cases of known and unknown propagation channel. The derivation of the BCRB for the unknown propagation channel is more complicated, as the Bayesian information matrix (BIM) without prior is singular. We also introduce a transformed mean squared error (MSE), as the regular MSE is not a meaningful error metric in this case.
• Based on the derived BCRB, we discuss the antenna response observability and the qualitative behavior of the BCRB. We take an estimation theoretic perspective to investigate the difference between antenna calibration when the propagation channel is known, like in a measurement chamber, and in-situ calibration in unknown propagation conditions.
• We introduce an in-situ antenna calibration algorithm based on MAP estimation, which takes multipath propagation into account. The algorithm can estimate arbitrary antenna responses and thus capture real-world antenna nonidealities including gain-phase offsets, mutual coupling, etc.
• Then, we show the effectiveness of the proposed algorithm by simulations, where the BCRB serves as a theoretical lower bound on the estimation MSE. Asymptotically, the performance approaches the case where the antenna response is perfectly known.
• Finally, we demonstrate the practical applicability of the algorithm by an experimental validation. For that, we perform in-situ calibration of an MMA installed on a robotic rover. We evaluate a measurement data set to compare the DoA estimation performance using the antenna response determined by electromagnetic (EM) simulation, in a near-field measurement chamber, and with the proposed in-situ calibration algorithm. By in-situ calibration, the DoA estimation performance is improved considerably. The paper is organized in eight sections. In Section II, we introduce the signal model and briefly introduce wavefield modeling and manifold separation. In the following Section III, we show how the MAP estimator for in-situ calibration using wavefield modelling is obtained by Bayesian inference. We further present an algorithm for iteratively solving the resulting nonlinear optimization problem. In Section IV, we derive the recursive BCRB. The case where the propagation channel is unknown, which is common in in-situ calibration, receives special attention. The behavior of the BCRB is discussed qualitatively, providing valuable intuition. The performance of the proposed algorithm is evaluated in Section V for the cases of one and multiple impinging signals and compared to the theoretical limits. In Section VI, we introduce the measurement setup and evaluate the algorithm with measurement data. Finally, in Section VII, we briefly discuss alternative approaches in terms of Bayesian filtering and global optimization and the case of using a gyroscope as external sensor, before we conclude the paper with Section VIII.
Notation: Vectors are written in bold lowercase letters and matrices in bold capital letters. (·) T and (·) H stand for vector or matrix transpose and conjugate transpose, respectively. 1 N and 0 N are vectors of ones and zeros with length N , I N is an N × N identity matrix and 0 N an N × N matrix with all zeros. ||a|| is the Euclidean norm. Square brackets refer to an element in a vector [

II. MULTIPORT ANTENNA SIGNAL MODEL
We assume that P superposed signals with index p ∈ {1, . . . , P} with DoAs φ = [φ s 1 , . . . , φ s P ] T and ToAs τ = [τ s 1 , . . . , τ s P ] T arrive at the multiport antenna with M ports. They originate from the same source, but have different DoAs, ToAs, amplitudes and phases, which corresponds to multipath propagation. The antenna is connected to a coherent multichannel receiver with M channels, which operates snapshot-wise and performs a discrete Fourier transform (DFT) on the snapshots. The received signal with snapshot index s ∈ {1, . . . , S} in discrete frequency domain r s (n) = [r 1 (n), . . . , r M (n)] T can then be written as where n ∈ {1, . . . , N } is the DFT bin and s(n, τ ) = s(n)e −j2π τ n−1 N is a delayed version of the transmitted signal s(n). We call α s p and ϕ s p absolute amplitude and phase of the p-th signal to distinguish from the relative amplitudes and phases of the signals at the respective antenna ports. The model follows the assumption that the signal bandwidth is small relative to the carrier frequency [1], [42]. We consider internal receiver noise, such that w s r (n) ∼ CN (0, σ 2 r I M ) is independent and identically distributed (i.i.d.) circular symmetric Gaussian noise.
An arbitrary multiport antenna can be described by its gain pattern g m (φ) and phase pattern m (φ) [42]. The antenna response a m (φ) for ports m ∈ {1, . . . , M } and DoA φ is then defined as Stacking the antenna responses for all ports into a vector yields the antenna response vector Wavefield modeling and manifold separation allows to decompose the antenna response vector (3) into a product of the sampling matrix G ∈ C M ×U , which describes the antenna characteristics, and a basis vector of order U , b(φ) ∈ C U , which is a function of the DoA [25], [26], i.e.
It is required that the antenna response is square integrable and the basis functions of order U are orthonormal on the manifold φ ∈ [−π, π). Considering azimuth only, a suitable basis is given by the Fourier functions, which we arrange in the basis vector The order U of the basis functions can be inferred from the electrical size of the antenna [25], or from the measurement noise floor [43]. Wavefield modeling and manifold separation VOLUME 10, 2022 can be extended to joint estimation of azimuth and elevation [12], [25], [26]. In that case, spherical harmonics or 2D Fourier functions can be used as basis. For this paper, we restrict ourselves to azimuth only, noting that the proposed methods can be extended to joint azimuth and elevation estimation.
During the process of antenna design, usually an EM simulation is performed. The output of this simulation are Q spatial samples of the antenna response e q = e q,1 . . . e q,M T (6) for M ports, DoAs φ q and q = 1, .., Q. In matrix form this yields which is proportional to the electric field strength and can be interpreted as a sampled version of (3). The spatial sampling grid must be dense enough, such that the spatial sample-rate criterion according to the Nyquist theorem is fulfilled. A sampling matrix based on the EM simulation data E sim can then be obtained by the least squares method, However, EM simulation does not account for manufacturing imperfections. The obtained sampling matrixĜ sim will thus deviate from the true G.
To overcome that, it is common to measure the antenna in a dedicated measurement chamber. This yields another set of spatial samples E meas . Similar to (8), the measurement-based sampling matrix is obtained bŷ In order to reflect the true antenna response of the installed antenna, the whole surrounding structure must be in place for the measurement, i.e. the antenna must be measured in its final installation location. The calibration procedure can become laborious and costly for large devices. Thus, we instead propose to calibrate the whole device containing the antenna in-situ, in order to capture the true antenna response of the installed antenna.

III. IN-SITU ANTENNA CALIBRATION ALGORITHM
The log-likelihood function for the received signals, see (1), is given by In [44] it is shown that the log-likelihood function can be concentrated toL with the vectors r s = (r s (1)) T , . . . , (r s (N )) T T and s(τ Concentrated means that the maximum of the log-likelihood function w.r.t. the unknown absolute amplitudes, absolute phases and noise variance has been plugged into the log-likelihood function, such that the unknown variables are reduced to the ones of interest. The joint DoA-ToA maximum likelihood (ML) estimator is then given by and requires knowledge of the true sampling matrix G, i.e. the true antenna response a(φ), see (4). In practice, the true sampling matrix G is unknown and only estimated sampling matrices from EM simulationĜ 0 sim or antenna measurement G 0 meas are available. Deviations of the estimated sampling matrixĜ 0 from the true sampling matrix G result in a model mismatch, which impairs DoA estimation performance [24].
To cope with that, we propose a MAP estimator for in-situ calibration, building on the principles of wavefield modeling and manifold separation.
An exemplary in-situ calibration scenario is depicted in Figure 1. We assume that a DoA observation for the first impinging signal p = 1, is available from an external sensor, which could be e.g. GNSS real-time kinematic (RTK) in an open-sky scenario. We model φ s obs as the true DoA φ s 1 with additive errors following a von Mises distribution w s φ obs ∼ M(0, κ φ obs ) with concentration κ φ obs [45]. The concentration κ φ obs is expected to be high, so the distribution can be approximated by a normal distribution w s φ obs ∼ N (0, σ 2 φ obs = 1/κ obs ), see [45].
The observed DoA log-likelihood function is then given by The observation vector for snapshot s is formed as In order to account for the sequence of snapshots and consider prior information in a Bayesian context, from now on we interpret the sampling matrix G as a random variable. For notational convenience, we vectorize the sampling matrix and split it into real and imaginary parts, The variables g 0 and g 0 RI are defined analogously. We use the complex-valued matrix, complex-valued vector or real-valued vector depending on the context. Empirically, we have found that typical deviations of the installed antenna response compared to the antenna response in free-space can be represented by circular symmetric Gaussian noise on the sampling matrix elements. Following the Bayesian approach, the prior probability density function (pdf) of the sampling matrix elements in logarithm domain is thus which is a circular symmetric Gaussian distribution with meanĝ 0 = E g 0 and variance σ 2 prior sampling matrix. We further define the state vector for snapshot s as and the full state vector as with the DoAs for all snapshots φ 1: We assume here that the model order P is known. In practice, it needs to be estimated, see e.g. [46]. In order to focus on the antenna calibration aspect, we also assume that the ToAs τ s are estimated separately. Under reasonable assumptions detailed in the following and medium to high signal-to-noise ratio (SNR), independent estimation of DoAs and ToAs is possible. In [47], the baseband-carrier correlation is defined in the context of a Cramér-Rao bound (CRB). When the power spectral density (PSD) of the signal is symmetric, the baseband-carrier correlation becomes zero, allowing independent estimation of ToA and DoA. Further support is provided by [8] and the intuition that under the narrowband assumption, ToA is estimated from the observation of a delayed baseband signal, while DoA is estimated from relative amplitudes and phases between the ports. Another prerequisite for the independence of ToA and DoA is that the paths are resolvable and the correct model order is chosen. Using Bayes theorem and the first order Markov assumption, the posterior pdf is written as Inserting (16) and (21) and assuming a stationary sampling matrix, non-informative transition and prior pdfs for φ s and independent noise in (1) and (14) yields We define the function which is proportional to the negative posterior pdf in logarithm domain. The prior pdf of the sampling matrix in logarithm domain, ln p(g 0 ), is given by (19). The concentrated log-likelihood function for the received signalsL r s (φ s , g) is given by (11), where the sampling matrix elements g are considered as unknown and the ToAs τ s are omitted. The log-likelihood function for the observed DoAs L φ s obs (φ s 1 ) is defined by (15). Thus, the maximum a posteriori (MAP) estimator is given bŷ Solving (25) is a challenging nonlinear optimization problem with 2MU + PS unknowns. For example, assuming the parameters considered for the simulations in Section V, this would be 1052 unknowns for one impinging signal and 3052 unknowns for three impinging signals.
To solve (25), we use the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm, which is a quasi-Newton method to solve unconstrained nonlinear optimization problems [48]. BFGS VOLUME 10, 2022 using gradient (26) 6: end for requires the gradient of (28), which has the structure The respective partial derivatives in (26) are derived in Appendix A. For convex functions, global convergence of Broyden-Fletcher-Goldfarb-Shanno (BFGS) is proven [48]. Problem (25) is nonconvex, but ensuring initialization close to the solution, local convergence is sufficient. To ensure close initialization, we propose the sequential procedure outlined in Algorithm 1. For s = 1, the sampling matrix elements are initialized with the prior and the ML estimator (13) is applied to obtain an initial estimate of the DoAsφ 1 . For s ∈ {2, . . . , S}, the (BFGS) algorithm is applied to obtain the estimatex For each snapshot s ∈ {2, . . . , S}, the estimate of the previous snapshot is used to initialize the Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm for the current snapshot. By gradually adding more snapshots, the convergence of the algorithm is improved. To show the convergence of the proposed method for different antenna responses, extensive simulation results are presented in Section V.

IV. ESTIMATION THEORETIC ANALYSIS A. INFORMATION FROM OBSERVATIONS
In order to assess the performance of the proposed algorithm, we derive the BCRB as a theoretical lower bound of the achievable MSE. It also allows us to gain insights into the structure of the estimation problem. We start by calculating the information contained in the observation of a single snapshot s about the parameters of interest, defined by the state vector (21). To consider all unknown parameters, we augment (21) with the nuisance parameters absolute amplitude α s = [α s 1 , . . . , α s P ] T and absolute phase ϕ s = [ϕ s 1 , . . . , ϕ s P ] T for P multipath signals, see (1), to obtaiñ The calculation of the recursive BCRB in Section IV-B requires to quantify the information contained in z s , the observation vector of snapshot s in a Bayesian context [49]. Taking the expectation of the snapshot Fisher information matrix (FIM) I s w.r.t. the augmented state vectorx s , we obtain The snapshot FIM is defined by see [50], with the covariance matrix With (16) we obtain We continue by partitioningĨ s defined by (29) intõ The first main diagonal matrix block refers to the sampling matrix elements which are split into real part I s g R and imaginary part I s g I . The second main diagonal matrix block I s φ corresponds to the DoAs. The third main diagonal matrix block refers to the nuisance parameters, which are split into real-valued absolute amplitudes I s α and phases I s ϕ . The offdiagonal matrix blocks correspond to the relationship between sampling matrix and DoAs I s gφ , the relationship between sampling matrix and nuisance parameters I s gn and the relationship between sampling DoAs and nuisance parameters I s φn , respectively. An extensive derivation of the elements of (34) to (38) is provided in Appendix B.

B. RECURSIVE BAYESIAN Cramér-RAO BOUND 1) KNOWN PROPAGATION CHANNEL
First we examine the case where the propagation channel, expressed by the variables α s and ϕ s , is perfectly known. This is not the main case of interest for this paper, but it can serve as comparison to gain valuable insights. We derive the equivalent Bayesian information matrix (EBIM) [49], adopting the notion of equivalent Fisher information matrix (EFIM) from [51] and extending it to a Bayesian context. The EBIM contains the Bayesian information for a subset of the states, while considering the impact of the other states. It is calculated by applying the Schur complement. In that fashion, the EBIM for the sampling matrix elements is calculated recursively by adding up the information provided by the individual snapshots. The recursion of (39) is initialized with J 0 g,sync , where the subscript sync refers to the perfect propagation channel knowledge. With the prior pdf from (19), J 0 g,sync = 2 prior, J 0 g,sync = 0 2MU is a zero matrix. We define the MSE of the sampling matrix elements as the trace of the MSE matrix divided by the number of complex coefficients MU , Now we can obtain the BCRB for the known propagation channel case as a lower bound on the estimation MSE of the sampling matrix elements, In a similar way, we can obtain the EBIM for the DoAs, which leads to the BCRB on the p-th DoA for the known propagation channel case, 2) UNKNOWN PROPAGATION CHANNEL The main case of interest for in-situ calibration is when the propagation channel is unknown. Calculation of the EBIM for the sampling matrix elements is similar to (39), except that we now have to take the nuisance parameters α s and φ s into account. We thus obtain the EBIM for the sampling matrix elements In the unknown propagation channel case, we have to distinguish between the case with prior and without prior.

a: WITH PRIOR
With prior, we have J 0 g = 2 σ 2 g 0 I 2MU according to the prior pdf (19). As J 0 g is positive definite and by (44), the Schur complement of the expectation of the FIM, which is positive semidefinite, is recursively added. Thus, the sampling matrix EBIM J s g is positive definite with full rank and we can obtain the BCRB of the sampling matrix elements by inverting the EBIM, with the MSE defined by (40). In a similar fashion, we can obtain the EBIM for the DoAs, and the respective BCRB b: WITHOUT PRIOR Without prior, i.e. J 0 g = 0 2MU , one can make an important observation when computing (44) recursively. No matter how many snapshots are added, J s g will always be rank deficient by two. The reason is that in the unknown channel case, the absolute amplitude and phase of the antenna response are not observable. DoA estimation relies on relative amplitudes and phases, therefore absolute offsets do no harm [52]. Still, the VOLUME 10, 2022 non-observability of absolute amplitude and phase implies that by in-situ antenna calibration, instead of the true antenna response, only an antenna response which is equivalent for DoA estimation can be estimated. Thus, the power pattern of the equivalent antenna response does not necessarily represent the antenna gain. Furthermore, the non-observability of absolute amplitude and phase has implications on the estimation theoretic analysis. Since J s g is rank deficient, inversion is not possible. To be able to calculate a BCRB, we set explicit constraints on the sampling matrix in the form see [53]- [55]. For instance, amplitude and phase could be constrained by Since the constraints can be written in the form (48), we can apply the Moore-Penrose to obtain a meaningful BCRB for the constrained problem [53]- [55]. The BCRB for the sampling matrix elements is thus calculated as However, the regular MSE (40) is not useful in this case anymore, as absolute amplitude and phase offsets would count as errors. Instead, in order to calculate a meaningful MSE forĝ s , we have to make sure that its absolute amplitude and phase match the true g, see [54]. This can be achieved by a transformation. We multiplyĝ s by the complex coefficient (ĝ s ) † g to obtain the transformed MSE, The EBIM for the DoAs is also calculated with the pseudoinverse as leading to the BCRB for the p-th DoA

C. QUALITATIVE BEHAVIOR OF THE BCRB
Before discussing the BCRB, we want to gain further insight into the behavior of the sampling matrix EBIM J s g . Since J s g is a real-valued symmetric matrix, an eigendecomposition  with an orthonormal basis Q and a diagonal matrix containing the eigenvalues λ 1 , . . . , λ i , . . . λ 2MU can be performed. Figure 2 shows the evolution of the resulting eigenvalues for a single realization in the unknown channel case with prior. The assumed parameters are stated at the beginning of Section V-B, but are not important for a qualitative assessment. As can be seen, all eigenvalues start off at 2/σ 2 g 0 = 22.22, corresponding to the information provided by the prior. Since E φ s obs − φ (s−1) obs = 1 • , see the beginning of Section V-B, one full 360 • turn corresponds to approximately 360 snapshots. During the first full turn, all eigenvalues except two increase quickly, which means a large information gain. After the first full turn, they increase slowly. Two eigenvalues however stay at their initial value 2/σ 2 g 0 = 22.22. As discussed in the last section, these two eigenvalues correspond to the non-observable absolute amplitude and phase. In the case without prior, the two eigenvalues would be zero, hence the pseudoinverse must be applied in (50) and (52). Having discussed the sampling matrix EBIM, we now want to assess the BCRB behavior qualitatively to gain some insight into the problem at hand. For that we have a look at the BCRBs of the sampling matrix elements g for the case of a single impinging signal, plotted in Figure 3. First we focus on the known propagation channel case. The known channel BCRB without prior has a steep decrease until the completion of the first full 360 • turn after approx. 360 snapshots. The known channel BCRB with prior first decreases slowly, then we see a steep drop around the completion of the first full turn. In this region, the known channel BCRB without prior approaches the known channel BCRB with prior. Obviously, it is vital that the full manifold is covered. Once the full manifold has been observed, the prior is not useful anymore. After the first full turn, both known channel BCRBs decrease slowly with increasing number of snapshots. Here we only see an improvement due to averaging over noisy observations. Now we take at look at the unknown propagation channel BCRBs. The unknown channel BCRB without prior shows a similar behavior as the one for the known channel. However, for the known channel case, the BCRB is almost one order of magnitude lower. The unknown channel BCRB with prior starts off next to the known channel one, but only slightly decreases during the first full turn and then reaches a floor. As we have seen in Figure 2, when s → ∞, two eigenvalues of J s g will remain 2/σ 2 g 0 , corresponding to the non-observable absolute amplitude and phase, while all other eigenvalues increase. Correspondingly, the two dominant eigenvalues of (J s g ) −1 will be σ 2 g 0 /2. In the limit, we thus have For the parameters chosen here, see Section V-B, we have σ 2 g 0 /MU = 0.05, which is exactly the floor value of the unknown channel BCRB with prior from Figure 3.
As discussed in Section IV-B2, the BCRB (45) refers to the regular MSE (40), where absolute amplitude and phase offsets defined by the prior count as errors. It is thus not a useful bound to asses the estimation quality of the sampling matrix elements. In contrast, for the case without prior The true antenna response is constant over time and can thus be estimated perfectly with an infinite number of observations. Still, it is not obvious what a certain error in the estimation of the sampling matrix elementsĝ s means for DoA estimation. Therefore, in Figure 4 we have a look at the BCRB of the DoAφ s 1 . All BCRBs start at σ φ obs . Initially, the observed DoA is more accurate than the estimated DoAs from the received signal with the prior sampling matrixĜ 0 . The unknown channel BCRB without prior decreases steadily, until the first turn has been completed after approx. 360 snapshots, then decreases slowly. The known channel BCRB without prior shows a similar behavior, but exists only for s ≥ 70, when the EBIM (42) is full rank. In contrast to that, both BCRBs with prior undergo a quick decrease during the first few snapshots. The model parameters, which are observable by a single snapshot, are adjusted during these first few snapshots, e.g. amplitude and phase offsets of the different ports. These parameters are not contained in the model explicitly, but are covered by the elements of the sampling matrix G. After approximately one full turn, the BCRB without prior approaches the BCRB with prior. For comparison, the snapshot DoA estimation CRB for known antenna response is plotted. It is slightly direction-dependent, thus snapshot-dependent. After one full turn, the BCRB for a known propagation channel is very close to the snapshot CRB for the case where the true antenna response is known. This is expected, as long as the manifold is sampled dense enough, the SNR of the received signal is high enough and σ φ obs is not too large. Note that the known antenna response CRB is not a lower bound for the BCRB, as it does not consider φ obs , i.e. the BCRB could also be lower. Still it is a good indicator to show when performance is close to optimum. The BCRB without channel knowledge does not fully approach the known antenna response CRB, but comes close after 1000 snapshots. For more accurate observable DoA, i.e. smaller σ φ obs , or more snapshots, it can also approach the known antenna response CRB.

A. SIMULATION SETUP
Now we evaluate the proposed in-situ antenna calibration algorithm by simulation and assess its effectiveness by comparison to the theoretical limits derived in Section IV. To demonstrate that the proposed algorithm works for arbitrary antenna responses, we perform 500 Monte Carlo runs, where for each run we randomly generate a new antenna response. Starting point for the random antenna response VOLUME 10, 2022 generation is a UCA with M = 4 antenna elements and radius R = 0.9/(4 sin(π/M )), which is described by see [1]. We evaluate (58) at a regular angular grid φ 1 , . . . , φ Q with Q = 360, to obtain the UCA sampling matrix with U = 9 basis functions analog to (7) and (8). Then, for each Monte Carlo run a random sampling matrix is generated by distorting the UCA sampling matrix with circular symmetric Gaussian noise with σ g = 0.3. The outlined procedure models deviations of the installed antenna response of a real-world UCA compared to the ideal UCA antenna response (58). Such deviations occur e.g. due to manufacturing imperfections, mutual coupling and the influence of the antenna surroundings. G 0 is obtained for each run by sampling the prior pdf (19) with σ g 0 = 0.3. Furthermore, we assume S = 1000 snapshots, SNR =  .) is a uniform distribution.

B. SINGLE IMPINGING SIGNAL
First, we focus on a single impinging signal P = 1. Figure 5 shows the simulated root-mean-square error (RMSE) for estimating the elements of the sampling matrix G. We compare the untransformed RMSE of the MAP estimator defined by (40) to the BCRB with prior (45). The two curves overlap and decrease during the first full turn, before they flatten out. As discussed for Figure 3, their behavior is explained by the absolute amplitude and phase offset defined by the prior, which adds to the RMSE. In contrast to that, a useful error metric is the transformed RMSE curve calculated by (51). At the beginning, its behavior is similar to the other RMSE curve, however towards the end of the first full turn, approaching s = 360 snapshots, it decreases much more and gets close to the BCRB without prior. The BCRB without prior is not a lower bound for the estimator with prior. However, as Figure 5 shows, after one full turn the information gathered from the observations is much more valuable than the prior, such that the two curves are close to each other.
Next, we have a look at the RMSE of the DoA estimatê φ s 1 , which is plotted in Figure 6. The prior sampling matrix G 0 does not reflect the real sampling matrix very well, hence the initial estimation error is constrained by σ φ obs . After few snapshots, the RMSE drops quickly. As stated in Figure 6, certain parts of the sampling matrix G can be observed on a snapshot basis. After that, the error decreases slowly, until the first turn is completed. Around that point, the error decreases quickly again. For increasing s, it slowly approaches the snapshot CRB for the known antenna response. As discussed in Section IV-C, the snapshot CRB is not a lower bound for the BCRB, but nevertheless an indicator of good performance.
For the final application, the observable DoA φ obs will not be present. Of major interest is thus, how good DoA estimation with the ''trained'' sampling matrixĜ s performs. We thus simulate DoA estimation for a discrete set of DoAs φ q with q ∈ {1, . . . , Q} and Q = 72 spanning a regular grid over the manifold. Figure 7 shows the RMSE over this grid, where for each snapshot the estimated sampling matrixĜ s , obtained by Algorithm 1, is used. The snapshot CRB for known sampling matrix is plotted as benchmark. With the original setting of σ φ obs = 2 • , sub-degree accuracy is achieved after a full turn. After that, the RMSE decreases slowly. For comparison, the whole simulation is also performed with σ φ obs = 0.5 • . In this case, the snapshot CRB is almost reached after one turn. To conclude, the accuracy of the observable DoA decisively determines the performance.

C. MULTIPLE IMPINGING SIGNALS
To assess how the in-situ calibration algorithm performs in the face of multipath propagation, we now assume three impinging signals P = 3. The first signal has SNR = 15 dB, the other two 3 dB less. The signals arrive from fixed locations φ s 2 = φ s 1 + 45 • and φ s 3 = φ s 1 + 90 • with ToAs τ 2 = τ 1 + 48 ns and τ 3 = τ 1 + 96 ns. The first DoA is observable  from another sensor with σ φ obs = 0.5 • . Figure 8 shows the RMSE for the DoAsφ s 1 ,φ s 2 andφ s 3 , as well as the respective BCRBs and snapshot CRBs for known antenna response. The behavior ofφ s 1 , which is observable by φ s obs , is similar to Figure 6. The estimation accuracy ofφ s 2 andφ s 3 relies entirely on an accurately estimated sampling matrixĜ s . Bothφ s 2 and φ s 3 start off with an RMSE above 10 • , which quickly goes down to a few degrees. The result is in line with the findings for Figure 6, that some components of the antenna response or the sampling matrix can be observed by single snapshots. After the first full turn, the RMSE ofφ s 2 andφ s 3 is close to the respective snapshot CRB. The snapshot CRBs forφ s 2 andφ s 3 are higher than forφ s 1 , as the SNR is 3 dB lower. The proposed in-situ calibration algorithm performs well for this artificial multipath sceario. An important prerequisite in practice is a correct model order estimation [46] and separability of the signals from the different propagation paths in at least either DoA or ToA domain. Given that these conditions are fulfilled, we conclude that the algorithm is suitable to perform in-situ calibration in multipath scenarios.

A. MEASUREMENT SETUP
We have performed measurements to assess the practical applicability of the in-situ calibration algorithm. The antenna to be calibrated is an MMA, specifically a dielectric resonator antenna with four modes, which are excited independently. It is an in-house development from the German Aerospace Center (DLR), for details please see [56]. We have chosen U = 13 basis functions, see (5), for the manifold separation in (4).
The MMA together with an assembly is mounted onto a robotic rover shown in Figure 9a. The robotic rover is part of a fleet, which is used for navigation and exploration experiments at DLR, see [10] and the references therein. Part of the assembly is also a Universal Software Radio Peripheral (USRP) N310, a host computer, a Wi-Fi transceiver, batteries and two GNSS RTK receivers. The assembly can optionally be mounted on a turntable, see Figure 9b. The commercial multi-sensor RTK system [57] provides a reference for position and orientation. Internally, fusion of GNSS and inertial observations is performed. The output of the commercial RTK system, together with the known anchor node positions, is used to calculate the observable DoA φ s obs for the in-situ calibration algorithm.
For the physical layer signaling, we use our in-house developed DLR Swarm Communication and Navigation system [10], [58]. The system parameters are 1.68 GHz carrier frequency, 31.25 MHz sampling rate and orthogonal frequency-division multiplexing (OFDM) with fast Fourier transform (FFT) length 1024. The transmitted signal is a Zadoff-Chu sequence of length N = 463, which is mapped onto 925 subcarriers by occupying every second subcarrier. The occupied bandwidth is thus ≈ 28.2 MHz. Channel access VOLUME 10, 2022  is realized by a self-organizing time-division multiple access (TDMA) scheme with 100 ms round-trip schedule. The transmit power for this experiment is −15 dBm.
The system is implemented as software-defined radio (SDR). For the single channel nodes, GNU Radio and Ettus Research B200mini devices are used. For MMAs, DoA information is in general contained in both magnitude and phase of the received signal, see (2) and [12]. Thus, a phase-coherent multichannel receiver is required and the phase and amplitude imbalances between the channels need to be corrected. We use the USRP N310 from Ettus Research. The local oscillator (LO) is provided by an external frequency synthesizer, which enables phase coherency for all four channels [58].

B. IN-SITU ANTENNA CALIBRATION
For the in-situ calibration, the assembly is mounted on the turntable, see Figure 9b, and placed in between the three anchor nodes. An aerial picture of the experiment site is shown in Figure 10. Two full turns of the turntable are performed. During this time, a total of S = 2264 snapshots, consisting of 767 from A1, 750 from A2 and 747 from A3 are received. These are used to run the in-situ calibration Algorithm 1 in post-processing. Figure 11 shows three different antenna responses as power and phase patterns of the MMA. The first antenna response is obtained by EM simulation, where the antenna alone is simulated in free space. The second antenna response is obtained by measurement with an MVG Starlab near-field measurement system. For the phase patterns, but less pronounced also for the power patterns, a difference is visible between simulated and measured patterns. In addition to manufacturing tolerances, this is explained by the fact that the mounting structure and the radome of the manufactured antenna were not considered by the simulation. The third antenna response is the result of the in-situ calibration Algorithm 1. It is closer to the measured than to the simulated antenna response. Still, notable differences are apparent.

C. EVALUATION
In the end, the main interest is the DoA estimation performance. To evaluate that in a fair manner, we use a different data set. To record the evaluation data, the MMA assembly was mounted on a robotic rover, see Figure 9a. The robot was driving in-between and around the three anchor nodes with an average velocity of 0.3-0.6 m/s. A map of the track and the anchor nodes is shown in Figure 12. Starting point for the robot is in the lower left corner. During the 7 min 50 s long track, 4450 snapshots were received from A1, 4336 from A2, 4170 from A3, and recorded. Based on this data, the signal DoA is estimated in post-processing by the ML estimator (11) and (13) Figure 13 shows the absolute DoA estimation error |φ ML − φ| for the signals received from A2 over time. When the antenna response from EM simulation is used, error peaks above 20 • are visible. Using the antenna response from the near-field measurements, the error peaks are reduced to around 10 • to 15 • . With the antenna response obtained by the proposed in-situ calibration algorithm, the estimation error is reduced further. In this case, the error peaks top out at 10 • , but are mostly lower. The experiment area was bumpy grassland. While moving, the robot was shaking to some extent. This could explain the spiky behavior of the estimation error.
For a more detailed analysis, Figure 14 shows the empirical cumulative distribution functions (CDFs) of the DoA estimation error w.r.t. anchor nodes A1, A2 and A3. Again we distinguish between the three different antenna responses. When the measured antenna response is used, 90% of the measurement errors w.r.t. A2 and A3 are below 7.9 • and 8.8 • , respectively, compared to 4.6 • and 4.5 • for the antenna response from in-situ calibration. The signals from A1 show higher estimation errors. Here the 90% values are 10.4 • for the measured antenna response compared to 7.2 • with in-situ calibration. Using the proposed in-situ calibration algorithm results in better performance compared to EM simulation or FIGURE 13. Absolute DoA estimation error with respect to anchor node A2 using the antenna response from EM simulation, measured in a near-field chamber and obtained by the in-situ calibration Algorithm 1.

FIGURE 14.
Empirical cumulative probability of the DoA estimation error for the three anchor nodes A1, A2, A3 using the antenna response from EM simulation, measured in a near-field chamber and obtained by the in-situ calibration Algorithm 1.
a near-field measurement of the antenna. Thus, we conclude that the in-situ antenna calibration algorithm proposed in Section III works well in practice.

A. BAYESIAN FILTERING
Instead of estimating the full posterior (23) at once, a common approach is to estimate it recursively over time as in Bayesian filtering. For the in-situ calibration problem studied in Section III, Bayesian filtering is challenging. Due to the large number of unknown variables, particle filtering approaches are not well suited, as they suffer from the curse of dimensionality. Additionally, the problem is highly nonlinear. Empirically, we have found that the extended Kalman filter (EKF) diverges due to the high nonlinearity. Performing Gauss-Newton iterations during the update step of the filter, as in the iterated extended Kalman filter (IEKF) [59], improves the situation. As long as σ φ obs was not too large, the IEKF did not diverge.

B. GLOBAL OPTIMIZATION TECHNIQUES
Algorithm 1 is a straightforward approach and was found to converge for both, simulated and measurement data. However, its complexity grows with the number of basis functions U , with number of impinging signals P and with the number of snapshots S. Often, in-situ calibration algorithms are applied to recorded data in post-processing, as in Section VI. For post-processing, computational complexity or execution times are usually not critical and are thus not a focus of this paper. In the case where complexity does matter, one could instead try solve for a batch of snapshots or even all snapshots at once. Due to the high nonlinearity and many local minima of the cost function, this is challenging. However, instead of using BFGS as in Algorithm 1, one could try to apply a global optimization algorithm [60].

C. SENSOR BIAS
For the measurement results presented in Section VI, φ obs was obtained by GNSS RTK. When no GNSS reception is possible, other sensors need to be considered for φ obs . Often gyroscopes are used to determine turn rates, which can be integrated to angles relative to a known initial direction. Especially for compact and low-cost microelectromechanical system (MEMS) gyroscopes, the measurements are not only corrupted by zero-mean Gaussian noise, but additionally suffer from a sensor bias. This sensor bias is not constant, but is drifting over time [61]. If a MEMS gyroscope is considered to obtain φ obs , a potential approach could be to include an appropriate model for the sensor bias, see e.g. [62], into the state vector.

VIII. CONCLUSION
In this paper, we present an algorithm to perform in-situ calibration of an arbitrary multiport antenna using wavefield modeling and manifold separation. We derive the MAP estimator to jointly estimate the sampling matrix together with the DoAs using the received signals and a noisy, external sensor. We also analyze in-situ antenna calibration from an estimation theoretic perspective. For this purpose, we derive the BCRB, which also serves as a benchmark for the algorithm. For the unknown propagation channel case with prior sampling matrix, the EBIM is full rank, but the resulting BCRB bounds the regular MSE, where absolute amplitude and phase offsets of the antenna response count as errors. The regular MSE is thus not a meaningful error metric for the sampling matrix. In contrast, for the unknown propagation channel case without prior sampling matrix, the EBIM is rank deficient. We show that a meaningful BCRB can be obtained by the pseudoinverse, which bounds the transformed MSE, where absolute amplitude and phase offsets of the antenna response are compensated. Extensive simulations are performed and compared to the BCRB to show the effectiveness of the proposed algorithm. We also show that by in-situ calibration, the DoA estimation performance asymptotically approaches the case where the antenna response is perfectly known. The algorithm can cope with multipath propagation, given the correct model order and resolvability of the paths. To proof the practicability, measurements with an MMA mounted on a robotic rover have been performed. When the proposed in-situ calibration algorithm is applied, DoA estimation accuracy is improved by 30% to 50% compared to calibration in a near-field measurement chamber.

APPENDIX A GRADIENT OF MAP COST FUNCTION
The elements of the MAP cost function gradient (26) are given by Deriving the likelihood function (11) and applying the chain rule we obtain (68) to (70).
where we define