2D-DOA Estimation and Auto-Calibration in UCAs via an Integrated Wideband Dictionary

Zavareh Bozorgasl, Hao Chen, Mohammad J. Dehghani Z. Bozorgasl .H. Chen is with the Department of Electrical and Computer Engineering, Boise State University, Boise, ID, 83712 (email: [email protected]) M. J. Dehghani is with the Department of Electrical and Electronics Engineering, Shiraz University of Technology, Shiraz, Iran (email: [email protected]).

Abstract

In this paper, we present a novel auto-calibration scheme for the joint estimation of the two-dimensional (2-D) direction-of-arrival (DOA) and the mutual coupling matrix (MCM) for a signal measured using uniform circular arrays. The method employs an integrated wideband dictionary to mitigate the detrimental effects of the discretization of the continuous parameter space over the considered azimuth and elevation angles. This leads to a reduction of the computational complexity and obtaining of more accurate DOA estimates. Given the more reliable DOA estimates, the method also allows for the estimation of more accurate mutual coupling coefficients. The method utilizes an integrated dictionary in order to iteratively refine the active parameter space, thereby reducing the required computational complexity without reducing the overall performance. The complexity is further reduced by employing only the dominant subspace of the measured signal. Furthermore, the proposed method does not require a constraint on the prior knowledge of the number of nonzero coupling coefficients nor suffer from ambiguity problems. Moreover, a simple formulation for 2-D non-numerical integration is presented. Simulation results show the effectiveness of the proposed method.

Index Terms—Direction-of-Arrival, Sensor Array Processing, Sparse Representation, Integrated Dictionary, Auto-Calibration, Uniform Circular Array, Mutual Coupling Matrix.

I Introduction

Estimating the direction-of-arrival (DOA) of a plane wave impinging on a sensor array, a classical problem in array signal processing, has been widely investigated in recent decades [1, 2, 3]. Although the literature mostly deals with DOA estimation techniques which use the assumption of a perfectly known measurement array, the resulting steering vector will generally only be approximately known. Indeed, the manifold will suffer from effects such as mutual coupling, which may cause a substantial degradation of the performance of applied algorithms [4, 5, 6]. Two of the most commonly used array configurations are the uniform linear array (ULA) and the uniform circular array (UCA). As the mutual coupling effect will be more pronounced for a circular than for a linear array, the detrimental effects of such coupling will affect the resulting performance more, suggesting the current study.

In order to alleviate the mutual coupling effects, two main categories of algorithms have been proposed in the literature. The first of these employ offline calibration algorithms[7, 8], which require calibration sources, i.e., some sources with perfectly known locations, to determine the array manifold matrix. In contrast, the second category uses auto-calibration or online calibration algorithms[9, 10, 11, 12], which do not require such source calibration, instead aiming to compensate for the calibration automatically. As finding additional calibration sources is not possible in many practical applications, the auto-calibration algorithms are often preferable.

Moreover, as the discretization of the parameter space leads to undesirable off-grid effects, there has recently been increasing interest in continuous parameter space techniques[13, 14]. One such alternative is to use a continuous dictionary together with an atomic norm penalty[15, 16, 17]. While such an approach often yields an accurate signal reconstruction, this form of problem formulation often leads to optimization problems with a high computational complexity. In this work, instead, we propose the use of an integrated wideband dictionary, as introduced [18], which is formed over subsets of the continuous parameter space. The essence of using this kind of dictionary is the ability to discard the non-activated subsets and preserving, then refining the activated ones for further zooming without the risk of missing components. These steps are then iterated until one achieves the desired resolution.

Common problems of earlier auto-calibration, including ambiguity problems [11, 12] and/or the requirement of a priori knowledge of the number of nonzero mutual coupling coefficients [19, 20], both typically being infeasible in practice. The here proposed method, which iteratively estimates both the DOAs and the mutual coupling coefficients, suffers from neither of these shortcomings. For the DOA estimation stage, instead of using a conventional grid-based dictionary which would necessitate a high computational complexity to allow for a sufficiently fine grid, a wideband dictionary framework is employed. By forming the used dictionary consisting of subsets of the continuous parameter space, and applying the noted screening process, the non-activated sets may be discarded in the first stage, whereas the activated subsets are retained and refined for further processing. In this sense, the initial DOA estimation stage will constitute an initial coarse estimate that will not suffer from the usual off-grid effects, such as missing active components. Via the further steps, i.e., refining the remaining activated bands, the closely spaced DOAs are separated. Next, the MCM estimation is formed utilizing the complex symmetric circular Toeplitz structure of the MCM, as was proposed for 1-D DOA estimation in [21, 22]. Here, this mutual coefficient estimation stage is further extended to the 2-D DOA estimation problem for the UCA configuration.

The remainder of the paper is organized as follows. Section II outlines the data model for 2-D UCAs. Section III describes the details of the new method for joint iterative 2-D DOA and MCM estimation. Section IV provides some simulations to demonstrate the effectiveness of the proposed method. Finally, Section V concludes the work.

Notation: Bold lower-case (upper-case) letters are used throughout to denote vectors (matrices). $\mathbb{C}$ denotes the set of complex numbers. $(\cdot)^{T}$ and $(\cdot)^{H}$ represent the transpose and the complex conjugate transpose, respectively. The notations $\|\cdot\|_{k}$ and $\|\cdot\|_{F}$ stand for the k-norm of a vector and the Frobenius norm of a matrix, respectively. $\hat{(\cdot)}$ is an estimate of $(\cdot)$ , and $\text{diag}(\mathbf{x})$ is a diagonal matrix with $\mathbf{x}$ being its diagonal elements. $\mathbf{X}(:,i)$ and $\mathbf{X}(j,:)$ show the $i$ -th column and the $j$ -th row of $\mathbf{X}$ , respectively. $\text{vec}(\mathbf{X})$ stands for the vector with the elements consists of columns of $\mathbf{X}$ stacked together.

II SIGNAL MODEL

Suppose that $K$ narrowband far-field sources simultaneously impinge onto a uniform circular array (UCA) of omnidirectional sensors from angles of arrival $(\phi,\theta)$ , where $\phi$ denotes the azimuth angle, i.e., $\phi\in[0,360)$ degrees, measuring counterclockwise from the $x$ -axis, and $\theta$ denotes the elevation angle, i.e., $\theta\in[0,90]$ degrees, measuring down from the $z$ -axis. The UCA geometry, with $N$ equispaced sensors on the circumference of the array, is shown in Fig. 1 in the $xy$ -plane. The marked arrow represents an example of a received signal geometrical description.

Refer to caption — Figure 1: The considered UCA receiver model for DOA and MCM estimation.

Therefore, by considering an array of $N$ sensors which are placed uniformly on a circle of radius $r$ , the array output $\mathbf{x}\in\mathbb{C}^{N\times T}$ may be expressed as

\mathbf{x}(t)=\mathbf{A}(\boldsymbol{\phi},\boldsymbol{\theta})+\mathbf{e}(t),\quad t=t_{1},t_{2},\ldots,t_{T},

(1)

where $\mathbf{A}(\boldsymbol{\phi},\boldsymbol{\theta})=[\mathbf{a}(\phi_{1},\theta_{1}),\mathbf{a}(\phi_{2},\theta_{2}),\ldots,\mathbf{a}(\phi_{K},\theta_{K})]$ is the array manifold/steering matrix, with

\mathbf{a}(\phi_{i},\theta_{i})=\exp(j\xi_{i}\cos(\phi_{i}-\gamma_{n})),\quad n=0,1,\ldots,N-1,

denoting the steering vector of the source from $(\phi_{i},\theta_{i})$ . Here, $\gamma_{n}=\frac{2\pi n}{N}$ is the displacement of the $n$ th element of the array from the $x$ axis, $\xi_{i}=k_{0}r\sin\theta_{i}$ , with the wave number $k_{0}=\frac{2\pi}{\lambda}$ , whereas $\mathbf{s}(t)=[s_{1}(t),s_{2}(t),\ldots,s_{K}(t)]^{T}$ and $\mathbf{e}(t)=[e_{1}(t),e_{2}(t),\ldots,e_{N}(t)]^{T}$ denote the source data and an additive white Gaussian noise vector, respectively. Finally, $t$ indexes each snapshot and $T$ is the number of snapshots. The signal vector $\mathbf{s}(t)$ and the noise vector $\mathbf{e}(t)$ are assumed to be zero mean circularly symmetric and statistically independent processes.

III THE PROPOSED AUTOCALIBRATION METHOD

By taking mutual coupling into consideration, the observation model in (1) can be reformulated as

\mathbf{X}=\mathbf{C}\mathbf{A}(\phi,\theta)\mathbf{S}+\mathbf{E},

(2)

where $\mathbf{X}=[\mathbf{x}(t_{1}),\mathbf{x}(t_{2}),\ldots,\mathbf{x}(t_{T})]$ , $\mathbf{S}=[\mathbf{s}(t_{1}),\mathbf{s}(t_{2}),\ldots,\mathbf{s}(t_{T})]$ , and $\mathbf{E}=[\mathbf{e}(t_{1}),\mathbf{e}(t_{2}),\ldots,\mathbf{e}(t_{T})]$ are the array output, as well as the source and noise data matrices, respectively. Furthermore, $\mathbf{C}$ denotes the mutual coupling matrix which is assumed to have a complex symmetric circular Toeplitz structure such that

\mathbf{C}=\text{toeplitz}(\tilde{\mathbf{c}},\tilde{\mathbf{c}})\in\mathbb{C}^{N,N},

(3)

where

\mathbf{\tilde{c}}=[c_{1},c_{2},\ldots,c_{L},c_{L-1},\ldots,c_{3},c_{2}]\in\mathbb{C}^{1,N},\quad L=\frac{N}{2}+1,\text{ for }N\text{ even,}

and

\tilde{\mathbf{c}}=[c_{1},c_{2},\ldots,c_{L},\ldots,c_{3},c_{2}]\in\mathbb{C}^{1,N},\quad L=\left(\frac{N+1}{2}\right),\text{ for }N\text{ odd.}

Here, $L$ denotes the degrees of freedom of the MCM, i.e., the number of the unknown elements in the first row of the matrix $\mathbf{C}$ . We note that because of the inverse relationship of mutual coupling and distance between each pairs of sensors $0<|c_{L}|<\ldots<|c_{2}|<|c_{1}|=1$ .

Algorithm 1 summarizes the proposed method for estimating the DOA and the MCM. Each substep will be elaborated upon in the following subsections.

Algorithm 1 The proposed SVD-based sparse WB method

1:Assign an initial MCM such that

\mathbf{C}_{0}=\mathbf{I}_{N\times N}

2:Choose the number of zooming steps and number of bands to use for both azimuth and elevation; select a suitable

\alpha

for Eq. (11).

3:Compute the number of sources and take the SVD of Eq. (2).

4:Compose the DOA integrated wideband dictionary with elements defined in Eq. (6).

5:Estimate the azimuth and elevation (as stated in Subsection III-B) using Eq. (10).

6:Estimate the MCM as detailed in Subsection III-C.

7:Form the DOA estimates on the active bands (composing of the zoomed dictionary given in Eq. (6)), and then using Eq. (10) include the MCM in step 6. The iterative zooming procedure enables discarding the non-activated parts of each dictionary.

8:Estimate the MCM as stated in Subsection III-C

9:Repeat steps 7 and 8 until obtaining the desired resolution.

III-A Sparse SVD-based Representation

The overall complexity of the proposed method may be notably reduced by projecting the measured data onto the $K$ -dimensional signal subspace. This can be formed by the $K$ dominant singular vectors of $\mathbf{X}$ . Form the singular value decomposition (SVD) as:

\mathbf{X}=\mathbf{U}_{s}\Sigma_{s}\mathbf{V}_{s}^{H}+\mathbf{U}_{c}\Sigma_{c}\mathbf{V}_{c}^{H},

(4)

where $\mathbf{U}_{s}\in\mathbb{C}^{N\times K}$ and $\mathbf{V}_{s}\in\mathbb{C}^{T\times K}$ are the singular vectors corresponding to the $K$ largest singular values, and $K$ represents the number of signals (sources). To do this, the number of sources, $K$ , must also be estimated. Here, this is done by forming the resulting solution over a set of potential $K\in\{1,2,\ldots\}$ and then selecting the estimated $K$ as the one resulting in the best signal modeling when combined with a BIC penalty to compensate for the measuring model order [23]. Let $\mathbf{X}_{\text{SV}}\triangleq\mathbf{X}V_{s}$ , $\mathbf{S}_{\text{SV}}\triangleq\mathbf{S}V_{s}$ , and $\mathbf{E}_{\text{SV}}\triangleq\mathbf{E}V_{s}$ , such that

\mathbf{X}_{\text{SV}}=\mathbf{C}\mathbf{A}(\boldsymbol{\phi},\boldsymbol{\theta})\mathbf{S}_{\text{SV}}+\mathbf{E}_{\text{SV}},

(5)

Indeed, by using the $N\times K$ dimensional matrix $\mathbf{X}_{\text{SV}}$ in place of $\mathbf{X}$ , one may significantly reduce the resulting computational complexity. The reduction will be significant in case of a small number of sources and/or a large number of time samples.

III-B DOA Estimation

Proceeding to determine the DOAs, let $(\phi_{i},\theta_{i})_{i=1}^{Q}$ for $i=1,2,...,Q$ , stands for the set of potential DOA candidates, with $Q$ indicating the number of candidates. Denote the dictionary over the set of candidates $\mathbf{A}(\hat{\phi},\hat{\theta})=[\mathbf{a}(\hat{\phi}_{1},\hat{\theta}_{1}),\mathbf{a}(\hat{\phi}_{2},\hat{\theta}_{2}),...,\mathbf{a}(\hat{\phi}_{Q},\hat{\theta}_{Q})]$ , where $\mathbf{a}(\hat{\phi}_{i},\hat{\theta}_{i})=[\mathbf{a}_{0}(\hat{\phi}_{i},\hat{\theta}_{i}),\mathbf{a}_{1}(\hat{\phi}_{i},\hat{\theta}_{i}),...,\mathbf{a}_{N-1}(\hat{\phi}_{i},\hat{\theta}_{i})]^{T}$ and each candidate is formed using the integrated wideband element, $\mathbf{a}_{n_{r}}(\phi,\theta)$ , such that

\mathbf{a}_{n_{r}}(\phi_{i},\theta_{i})=\int_{\phi_{i}}^{\phi_{i+1}}\int_{\theta_{i}}^{\theta_{i+1}}\exp\left(jk_{0}r\sin\theta\cos\left(\phi-\frac{2\pi n_{r}}{N}\right)\right)d\theta d\phi,

(6)

with $n_{r}=0,1,...,N-1$ , and where $[\phi_{i},\phi_{i+1})$ and $[\theta_{i},\theta_{i+1})$ are the intervals for the $i$ th azimuth and elevation candidate. Indeed, instead of using a finely spaced narrowband dictionary of the candidate DOAs, the 2-D angles domain, i.e., azimuth $\phi:[0^{\circ},360^{\circ})$ and elevation $\theta:[0^{\circ},90^{\circ})$ , is here divided into $B_{1}$ and $B_{2}$ continuous bands, respectively. As long as $B_{1}$ and $B_{2}$ are selected to be sufficiently large, this process would reduce the risk of any off-grid components and provide an initial coarse estimation of the angle of arrivals. Since all the bands are covered in the dictionary, no power is off-grid, which in turn avoids a non-sparse solution due to dictionary mismatch. In its initial step, this will yield a coarse estimation of the regions of interest, whereafter one may discard the non-activated bands and form the refined dictionary over the activated bands. This process may then be repeated until one reaches the desired resolution. In other words, instead of making a dictionary over a fine grid to limit the performance degradation due to off-grid effects, thereby necessitating a substantial computational burden, we use the integrated dictionary to determine and refine the estimate of the regions where sources are present. The process of zooming in some steps (here, 3 steps) is illustrated in Fig. 2, for the case of a single angle $\theta$ . As shown in [18], the average computation time of employing the wideband dictionary is substantially less than what would be required by the corresponding narrowband dictionary. The main difference of this zooming with respect to the earlier zooming methods is in utilizing an integrated dictionary for bands which significantly improves the performance.

Suppose the arrow in Fig. 2 shows the location of one of the sources. Then, one may describe the process of zooming in three stages as follows:

1.

Initially, one forms an integrated dictionary over the coarse grid ( $\mathbf{\theta}_{0}$ ).
2.

Then, one identifies the regions that may contain a source (as an example, the right side of $\mathbf{\theta}_{1}$ is a possible region).
3.

Finally, one composes an integrated dictionary with finer grids over the possible regions identified in step 2 (the third step is shown in the $\mathbf{\theta}_{2}$ region).

Remark: In Appendix A, we introduce an approach to form the 2-D integral in equation (6).

Proceeding, the azimuth and elevation estimation are formulated using the LASSO [24], which solves the resulting penalized regression problem. To this end, let $\mathbf{\hat{A}}=\mathbf{\hat{C}}\mathbf{A}(\hat{\phi},\hat{\theta})\in\mathbb{C}^{M\times Q}$ and then

\mathbf{D}=\begin{bmatrix}\mathbf{\hat{A}}\\ \vdots\\ \mathbf{\hat{A}}\end{bmatrix}_{MK\times Q}

(7)

Then,

\begin{bmatrix}\mathbf{X}_{\text{SV}}(:,1)\\ \vdots\\ \mathbf{X}_{\text{SV}}(:,K)\end{bmatrix}=\mathbf{D}\mathbf{\hat{s}}+\begin{bmatrix}\mathbf{N}_{\text{SV}}(:,1)\\ \vdots\\ \mathbf{N}_{\text{SV}}(:,K)\end{bmatrix}

(8)

Let $\mathbf{x}=\text{vec}(\mathbf{X}_{\text{sv}})\in\mathbb{C}^{MK}$ and $\mathbf{n}=\text{vec}(\mathbf{N}_{\text{sv}})\in\mathbb{C}^{MK}$ . Then,

\mathbf{x}=\mathbf{D}_{MK\times Q}\mathbf{\hat{s}}+\mathbf{n},

(9)

allowing the LASSO formulation¹, i.e., $\mathbf{\hat{s}}$ , results in

\mathbf{\hat{s}}=\arg\min_{\mathbf{\hat{s}}}\left\{\left\|\mathbf{x}-\mathbf{D}\mathbf{\hat{s}}\right\|_{2}^{2}+\gamma\left\|\mathbf{\hat{s}}\right\|_{1}\right\},

(10)

where $\mathbf{\hat{s}}\in\mathbb{C}^{Q\times 1}$ and $\gamma$ is a user chosen parameter controlling the sparsity of the solution, here being selected as

\gamma=\alpha\max_{i}|\mathbf{a}_{i}^{H}\mathbf{x}|,

(11)

where $\mathbf{a}_{i}$ , $i=1,2,\ldots,Q$ are columns of the $\mathbf{D}$ . The constant $\alpha$ is a user-chosen parameter allowing the performance of each dictionary to be evaluated (as in the simulation results) by varying values of this parameter within the range $0<\alpha\leq 1$ , with $\gamma_{max}=\alpha\max_{i}|\mathbf{a}_{i}^{H}\mathbf{x}|$ denoting the smallest tuning parameter value that gives a solution with the coefficients of zero [18, 24].

III-C MCM Estimation

Proceeding, we estimate the mutual coupling matrix coefficients by utilizing the complex symmetric circular Toeplitz structure of the MCM in a UCA. Suppose $\mathbf{R}_{x}$ denotes the array output covariance matrix and let $\lambda_{n}$ and $\mathbf{e}_{n}$ denote its eigenvalues and the corresponding eigenvectors, respectively. Then,

\mathbf{R}_{x}=\sum_{n=1}^{N}\lambda_{n}\mathbf{e}_{n}\mathbf{e}_{n}^{H}=\mathbf{E}_{s}\Lambda_{s}\mathbf{E}_{s}^{H}+\mathbf{E}_{n}\Lambda_{n}\mathbf{E}_{n}^{H},

(12)

where $\mathbf{E}_{s}=[\mathbf{e}_{1},\mathbf{e}_{2},\ldots,\mathbf{e}_{K}]\in\mathbb{C}^{N\times K}$ denotes the signal subspace containing the principal eigenvectors corresponding to the $K$ maximum eigenvalues, with the noise subspace $\mathbf{E}_{n}$ contains the remaining $N-K$ eigenvectors.

The signal subspace spans the same space as the array manifold matrix, i.e., $\text{span}\{\mathbf{E}_{s}\}\perp\text{span}\{\mathbf{CA(\phi,\theta)}\}$ . Moreover, the signal subspace and noise subspace are orthogonal such that $\text{span}\{\mathbf{E}_{s}\}\perp\text{span}\{\mathbf{E}_{n}\}$ , implying that $\text{span}\{\mathbf{E}_{s}\}\perp\text{span}\{\mathbf{E}_{n}\}$ , Implying that $\text{span}\{\mathbf{E}_{n}\}\perp\text{span}\{\mathbf{CA(\phi,\theta)}\}$ .
By supposing that we have DOAs (as obtained in Subsection III-B), the MCM could be estimated as the matrix $\mathbf{C}$ minimizing:

J=\left\|\mathbf{E}_{n}^{H}\mathbf{C}\mathbf{A}(\phi,\theta)\right\|_{F}^{2}=\sum_{k=1}^{K}\left\|\mathbf{E}_{n}^{H}\mathbf{C}\mathbf{a}(\phi_{k},\theta_{k})\right\|^{2}

(13)

where $\|\cdot\|_{F}$ denotes the Frobenius norm.

By using the circular symmetry of the mutual coupling matrix, one may instead estimate the vector which constructs the mutual coupling matrix, i.e., by rewriting $\mathbf{C}\mathbf{a}(\phi,\theta)=\mathbf{F}\{\mathbf{a}(\phi,\theta)\}\mathbf{c}\in\mathbb{C}^{L\times 1}$ where $\mathbf{F}$ denotes the resulting transform matrix, and $\mathbf{c}=[c_{1},c_{2},\ldots,c_{L}]^{T}\in\mathbb{C}^{L}$ . By extending the relation of the mutual coupling introduced in [21, 22] to 2-D and exploiting the symmetric circulant property of $\mathbf{C}$ , the matrix $\mathbf{F}\{\mathbf{a}(\phi,\theta)\}$ may be defined as:

\mathbf{F}\{\mathbf{a}(\phi,\theta)\}=\mathbf{F}_{1}\{\mathbf{a}(\phi,\theta)\}+\mathbf{F}_{2}\{\mathbf{a}(\phi,\theta)\}+\mathbf{F}_{3}\{\mathbf{a}(\phi,\theta)\}+\mathbf{F}_{4}\{\mathbf{a}(\phi,\theta)\}

(14)

$\displaystyle\left[\mathbf{F}_{1}\{\mathbf{a}(\phi,\theta)\}\right]_{i,j}$	$\displaystyle=\begin{cases}\left[\mathbf{a}(\phi,\theta)\right]_{i+j-1,1}&\text{if }i+j\leq N+1,\\ 0&\text{otherwise},\end{cases}$	(15)
$\displaystyle\left[\mathbf{F}_{2}\{\mathbf{a}(\phi,\theta)\}\right]_{i,j}$	$\displaystyle=\begin{cases}\left[\mathbf{a}(\phi,\theta)\right]_{j-i+1,1}&\text{if }i\geq 2,\\ 0&\text{otherwise},\end{cases}$
$\displaystyle\left[\mathbf{F}_{3}\{\mathbf{a}(\phi,\theta)\}\right]_{i,j}$	$\displaystyle=\begin{cases}\left[\mathbf{a}(\phi,\theta)\right]_{N+i+j-1,1}&\text{if }i<j\leq p,\\ 0&\text{otherwise},\end{cases}$
$\displaystyle\left[\mathbf{F}_{4}\{\mathbf{a}(\phi,\theta)\}\right]_{i,j}$	$\displaystyle=\begin{cases}\left[\mathbf{a}(\phi,\theta)\right]_{i+j-N-1,1}&\text{if }2\leq j\leq p,i+j\geq N+2,\\ 0&\text{otherwise}.\end{cases}$

where $i=1,2,\ldots,N$ , $j=1,2,\ldots,L$ , and $p=\left\lfloor\frac{N+1}{2}\right\rfloor$ . Applying this transformation to (13) yields

J=\sum_{k=1}^{K}\mathbf{c}^{H}\mathbf{F}^{H}\{\mathbf{a}(\phi_{k},\theta_{k})\}\mathbf{E}_{n}\mathbf{E}_{n}^{H}\mathbf{a}(\phi_{k},\theta_{k})\mathbf{F}\{a(\phi_{k},\theta_{k})\}\mathbf{c}=\mathbf{c}^{H}\mathbf{U}(\phi,\theta)\mathbf{c},

(16)

where

\mathbf{U}(\phi,\theta)=\sum_{k=1}^{K}\mathbf{F}^{H}\{\mathbf{a}(\phi_{k},\theta_{k})\}\mathbf{E}_{n}\mathbf{E}_{n}^{H}\mathbf{F}\{\mathbf{a}(\phi_{k},\theta_{k})\}^{H}\in\mathbb{C}^{L\times L},

(17)

By imposing the constraint that $c_{1}=1$ , (16) may be used as the cost function to be minimized in order to obtain the desired mutual coupling vector, i.e.,

\mathbf{c}=\mathbf{e}_{\min}\left[\mathbf{U}(\phi,\theta)\right],

(18)

where $\mathbf{e}_{\min}[\cdot]$ denotes the eigenvector corresponding to the smallest eigenvalue.

IV Simulation Results

In this section, some simulations will be carried out to investigate the performance of the proposed method. We consider a UCA with $N=15$ sensors, the radius of the UCA being $r=\lambda$ , where $\lambda$ denotes the center wavelength of the narrow-band signals. We further assume three far-field sources from angle of arrivals $(18.3^{\circ},243.4^{\circ})$ , $(83.6^{\circ},60^{\circ})$ , and $(73.9^{\circ},357.8^{\circ})$ for which $T=200$ snapshots are measured. The signals and the additive sensor noises are stationary, mutually uncorrelated, and zero mean. The signals are modeled as circularly symmetric Gaussian processes with variance $\sigma_{i}^{2}$ for the $i$ -th process. The corrupting noise is assumed to be i.i.d. white circularly symmetric Gaussian process with variance $\sigma^{2}$ , such that the input SNR of the $i$ -th signal is $10\log(\sigma_{i}^{2}/\sigma^{2})$ . The coupling coefficient is set to be $c_{2}=0.79+j0.432$ and $c_{3}=0.35+j0.16$ .

Fig. 3 presents the probability of correctly estimating the number of angle of arrivals using an integrated wideband dictionary. Suppose $B^{\theta}_{i}$ and $B^{\phi}_{i}$ denote the number of bands over elevation and azimuth, respectively. Also, we have $i=1,2,3$ for each stage of the dictionary. We here consider for the zooming steps of the elevation $B^{\theta}_{1}=30$ , $B^{\theta}_{2}=10$ , $B^{\theta}_{3}=3$ , and for the azimuth $B^{\phi}_{1}=120$ , $B^{\phi}_{2}=10$ , and $B^{\phi}_{3}=3$ (the grey curve), and $B^{\theta}_{1}=30$ , $B^{\theta}_{2}=5$ , $B^{\theta}_{3}=3$ , $B^{\phi}_{1}=120$ , $B^{\phi}_{2}=5$ , and $B^{\phi}_{3}=3$ (the brown curve). In this figure, the correct model order is given for varying values of $\alpha$ for (11). As is clear from Fig. 3, using $\alpha\leq 0.4$ for the three-stage dictionary (the grey curve) is a suitable region for utilizing the proposed method.

To evaluate the precision of DOA estimation methods, the average root-mean-square errors (RMSEs) of the elevation and the azimuth estimation are used as follows:

\text{RMSE}=\sqrt{\frac{1}{KM}\sum_{k=1}^{K}\sum_{m=1}^{M}\left[(\theta_{k}^{m}-\theta_{k})^{2}+(\phi_{k}^{m}-\phi_{k})^{2}\right]},

(19)

for $M=500$ number of Monte Carlo simulations, where $(\theta_{k}^{m},\phi_{k}^{m})$ denotes the estimate of $(\theta_{k},\phi_{k})$ in the $m$ -th Monte Carlo run (out of $M$ runs). Fig. 4 shows the resulting RMSE when the SNR varies from 0 dB to 20 dB. The results are compared with the LASSO estimates resulting when the dictionary is composed of discrete prolate spheroidal sequences (DPSS) [25, 26], and to the 2-D extension of the method in [20], here termed the Iterative-MUSIC estimator. This Iterative-MUSIC is found to have the highest computational complexity while offering the worst performance of the discussed methods.

Fig. 5 shows the RMSE of coupling coefficients, being defined as

\text{RMSE}_{\text{coupling coeff.}}=\frac{\sqrt{\sum_{m=1}^{M}\|\mathbf{\hat{c}}^{m}-\mathbf{c}\|^{2}}}{\|\mathbf{c}\|}\times 100\%,

(20)

where $\mathbf{\hat{c}}^{m}$ is the estimation of coupling coefficient vector, $\mathbf{c}$ , in the $m$ -th Monte Carlo experiment. As shown in Fig. 5, the coupling coefficients estimation of the proposed method is very effective and also has a very high precision.

In order to compare the performance in terms of the number of snapshots, we fix the SNR at $5$ dB, $M=200$ and vary the snapshot number. The results are shown in Fig. 6. The other parameters are set as before. As seen in Fig. LABEL:fig:fig3 to Fig. 6, the proposed method offers preferable performance of the three methods.

V Conclusion

In this work, we have introduced a self-calibration algorithm which is able to jointly estimate the 2-D DOAs and the MCM for a UCA configuration. The DOA estimation step is formed using a sparse representation framework in combination with an integrated dictionary elements, which span bands of the desired parameter space. The MCM estimation step exploits the circular symmetry of the mutual coupling matrix resulting for UCAs. This method has the advantages of allowing for an initial coarse gridding of the dictionary, without the risk of suffering from off-grid effects, as well as not posing any constraints on the coupling matrix, nor resulting in any ambiguity in the resulting DOA estimates. Several simulation results confirm the superior performance of the proposed method.

Appendix A Simplifying the 2-D integral of the dictionary elements in order to allow for non-numerical computations

To simplify the 2-D integration in (6) and attain a formula for computing the integration, one may use the Taylor series expansion for $e^{z}$ about $z=0$ , i.e.,

e^{z}=\sum_{n=0}^{\infty}\frac{z^{n}}{n!}

(21)

$\displaystyle\nu\triangleq$	$\displaystyle\int_{\phi_{i}}^{\phi_{i+1}}\int_{\theta_{i}}^{\theta_{i+1}}e^{j\alpha r\sin\theta\cos(\phi-\lambda)}d\theta d\phi$
$\displaystyle=$	$\displaystyle\int_{\phi_{i}}^{\phi_{i+1}}\int_{\theta_{i}}^{\theta_{i+1}}\sum_{n=0}^{\infty}\frac{\left(\alpha\cos(\phi-\lambda)\right)^{n}}{n!}\sin^{n}\theta d\theta d\phi$
$\displaystyle=$	$\displaystyle\int_{\phi_{i}}^{\phi_{i+1}}\sum_{n=0}^{\infty}\left(\frac{\alpha^{n}}{n!}\cos^{n}(\phi-\lambda)\int_{\theta_{i}}^{\theta_{i+1}}\sin^{n}\theta d\theta\right)d\phi$
$\displaystyle=$	$\displaystyle\sum_{n=0}^{\infty}\frac{\alpha^{n}}{n!}\int_{\phi_{i}}^{\phi_{i+1}}\cos^{n}(\phi-\lambda)d\phi\int_{\theta_{i}}^{\theta_{i+1}}\sin^{n}\theta d\theta$
$\displaystyle=$	$\displaystyle\sum_{n=0}^{\infty}\frac{\alpha^{n}}{n!}\int_{-\lambda_{n}}^{\phi_{i+1}-\lambda_{n}}\cos^{n}\phi d\phi\int_{\theta_{i}}^{\theta_{i+1}}\sin^{n}\theta d\theta$	(22)

Computing the first and the second terms yields

\nu=\left(\int_{\phi_{i}}^{\phi_{i+1}-\lambda_{n}}\phi d\phi\right)\left(\int_{\theta_{i}}^{\theta_{i+1}}d\theta\right)+\alpha\left(\sin\int_{\phi_{i}}^{\phi_{i+1}-\lambda_{n}}d\phi\right)\left(-\cos\int_{\theta_{i}}^{\theta_{i+1}}\theta d\theta\right)+\sum_{n=2}^{\infty}\frac{\alpha^{n}}{n!}\int_{\phi_{i}}^{\phi_{i+1}-\lambda_{n}}\cos^{n}\phi d\phi\int_{\theta_{i}}^{\theta_{i+1}}\sin^{n}\theta d\theta

(23)

Integration by parts for the third term of the above equation yields

\int\sin^{n}\theta d\theta=-\frac{1}{n}\sin^{n-1}\theta\cos\theta+\frac{n-1}{n}\int\sin^{n-2}\theta d\theta

(24)

and

\int\cos^{n}\phi d\phi=\frac{1}{n}\cos^{n-1}\phi\sin\phi+\frac{n-1}{n}\int\cos^{n-2}\phi d\phi

(25)

One may thus avoid a numerical integration by computing the recursive factors in the terms. To do this, the number of terms (for this work, computing about $\sim 10$ terms suffices and gives the same precision as the built-in integral2 function in Matlab) in this series could be computed to get the required precision.

Acknowledgment

The first author would like to thank Dr. Andreas Jakobsson (Dept. of Mathematical Statistics, Lund University) because of his insightful comments and discussions.

References

[1] H. Krim and M. Viberg, “Two decades of array signal processing research: the parametric approach,” IEEE Signal Process Mag, vol. 13, no. 4, pp. 67–94, 1996.
[2] H. C. So, Source localization: Algorithms and analysis. New York: Wiley-IEEE Press, 2011.
[3] Y. Liu, J. Chai, Y. Zhang, Z. Liu, M. Jin, and T. Qiu, “Low-complexity neural network based doa estimation for wideband signals in massive mimo systems,” AEU- Int J Electron Commun, vol. 138, p. 153853, 2021.
[4] Z. Zheng and C. Yang, “Direction-of-arrival estimation of coherent signals under direction-dependent mutual coupling,” IEEE Commun Lett, vol. 25, no. 1, pp. 147–151, 2020.
[5] F. Mei, H. Xu, W. Cui, B. Ba, and Y. Wang, “A transformed coprime array with reduced mutual coupling for doa estimation of non-circular signals,” IEEE Access, vol. 9, pp. 125 984–125 998, 2021.
[6] Z. Zheng, C. Yang, W. Q. Wang, and H. C. So, “Robust doa estimation against mutual coupling with nested array,” IEEE Signal Process Lett, vol. 27, pp. 1360–1364, 2020.
[7] C. M. See, “Sensor array calibration in the presence of mutual coupling and unknown sensor gains and phases,” Electronics Letters, vol. 30, no. 5, pp. 373–374, 1994.
[8] S. Liu, L. Yang, and S. Yang, “Robust joint calibration of mutual coupling and channel gain/phase inconsistency for uniform circular array,” IEEE Anten Wirel Propag, vol. 15, pp. 1191–1195, 2015.
[9] W. Hu and Q. Wang, “Doa estimation for uca in the presence of mutual coupling via error model equivalence,” IEEE Wirel Comm Letter, vol. 9, no. 1, pp. 121–124, 2019.
[10] R. Goossens and H. Rogier, “Direction-of-arrival and polarization estimation with uniform circular arrays in the presence of mutual coupling,” AEU- Int J Electron Commun, vol. 62, no. 3, pp. 199–206, 2008.
[11] K. Wang, J. Yi, F. Cheng, Y. Rao, and X. Wan, “Array errors and antenna element patterns calibration based on uniform circular array,” IEEE Anten Wirel Propag, vol. 20, no. 6, pp. 1063–1067, 2021.
[12] G. J. Jiang, Y. N. Dong, X. P. Mao, and Y. T. Liu, “Improved 2d direction of arrival estimation with a small number of elements in uca in the presence of mutual coupling,” AEU- Int J Electron Commun, vol. 71, pp. 131–138, 2017.
[13] J. Zhao, J. Liu, F. Gao, W. Jia, and W. Zhang, “Gridless compressed sensing based channel estimation for uav wideband communications with beam squint,” IEEE Trans Veh Technol, vol. 70, no. 10, pp. 10 265–10 277, 2021.
[14] P. Stoica and P. Babu, “Sparse estimation of spectral lines: Grid selection problems and their solutions,” IEEE Trans Signal Process, vol. 60, no. 2, pp. 962–967, 2011.
[15] M. Wagner, Y. Park, and P. Gerstoft, “Gridless doa estimation and root-music for non-uniform linear arrays,” IEEE Trans Signal Process, vol. 69, pp. 2144–2157, 2021.
[16] Y. Chi and Y. Chen, “Compressive two-dimensional harmonic retrieval via atomic norm minimization,” IEEE Trans. Signal Process., vol. 63, no. 4, pp. 1030–1042, 2014.
[17] Z. Yang and L. Xie, “Enhancing sparsity and resolution via reweighted atomic norm minimization,” IEEE Trans. Signal Process., vol. 64, no. 4, pp. 995–1006, 2015.
[18] M. Butsenko, J. Swärd, and A. Jakobsson, “Estimating sparse signals using integrated wideband dictionaries,” IEEE Trans. Signal Process., vol. 66, no. 16, pp. 4170–4181, 2018.
[19] D. Y. Gao, B. H. Wang, and Y. Guo, “Comments on "blind calibration and doa estimation with uniform circular arrays in the presence of mutual coupling",” IEEE Anten. Wirel. Propag. Lett., vol. 5, pp. 566–568, 2006.
[20] J. Xie, Z. S. He, and H. Y. Li, “A fast doa estimation algorithm for uniform circular arrays in the presence of unknown mutual coupling,” Progress In Electromagnetics Research C, vol. 21, pp. 257–271, 2011.
[21] M. Wang, X. Ma, S. Yan, and C. Hao, “An autocalibration algorithm for uniform circular array with unknown mutual coupling,” IEEE Anten. Wirel. Propag. Lett., vol. 15, pp. 12–15, 2015.
[22] B. Friedlander and A. J. Weiss, “Direction finding in the presence of mutual coupling,” IEEE Trans. Anten. Propag., vol. 39, no. 3, pp. 273–284, 1991.
[23] M. Wax and T. Kailath, “Detection of signals by information theoretic criteria,” IEEE Trans. Acoust. Speech Signal Process., vol. 33, no. 2, pp. 387–392, 1985.
[24] R. Tibshirani, J. Bien, J. Friedman, T. Hastie, N. Simon, J. Taylor, and R. J. Tibshirani, “Strong rules for discarding predictors in lasso-type problems,” Journal of the Royal Statistical Society: Series B (Statistical Methodology), vol. 74, no. 2, pp. 245–266, 2012.
[25] M. A. Davenport and M. B. Wakin, “Compressive sensing of analog signals using discrete prolate spheroidal sequences,” Applied and Computational Harmonic Analysis, vol. 33, no. 3, pp. 438–472, 2012.
[26] D. Slepian, “Prolate spheroidal wave functions, fourier analysis, and uncertainty—v: The discrete case,” Bell System Technical Journal, vol. 57, no. 5, pp. 1371–1430, 1978.