Block- and Rank-Sparse Recovery for Direction Finding in Partly Calibrated Arrays

TL;DR

This paper introduces a novel sparse recovery method for direction finding in partly calibrated arrays, leveraging block-sparsity and low-rank structures, with improved performance over existing methods in challenging scenarios.

Contribution

The paper presents a new technique combining nuclear norm and 1 norm minimization for direction finding in partly calibrated arrays, including a reformulation for efficient implementation and a gridless extension.

Findings

Outperforms the RARE method in low SNR scenarios

Effective with low sample numbers and correlated signals

Applicable to arbitrary subarray topologies

Abstract

A sparse recovery approach for direction finding in partly calibrated arrays composed of subarrays with unknown displacements is introduced. The proposed method is based on mixed nuclear norm and 1 norm minimization and exploits block-sparsity and low-rank structure in the signal model. For efficient implementation a compact equivalent problem reformulation is presented. The new technique is applicable to subarrays of arbitrary topologies and grid-based sampling of the subarray manifolds. In the special case of subarrays with a common baseline our new technique admits extension to a gridless implementation. As shown by simulations, our new block- and rank-sparse direction finding technique for partly calibrated arrays outperforms the state of the art method RARE in difficult scenarios of low sample numbers, low signal-to-noise ratio or correlated signals.

Equations

$$r_m^{(p)}={\rho }_m^{(p)}+{\eta }^{(p)},$$

$$\mathbf{Y}=\mathbf{A}\left(\mathbf{\mu },\mathbf{\eta }\right)\Psi +\mathbf{N},$$

$$\mathbf{A}(\mathbf{\mu },\mathbf{\eta })=\left[\mathbf{a}\left({\mu }_1,\mathbf{\eta }\right),\dots ,\mathbf{a}\left({\mu }_L,\mathbf{\eta }\right)\right],$$

$$\mathbf{a}(\mu ,\mathbf{\eta })$$

$$\displaystyle\boldsymbol{B}(\mu)={\rm blkdiag}\big{(}\boldsymbol{b}^{(1)}(\mu),\ldots,\boldsymbol{b}^{(P)}(\mu)\big{)}$$

$$\boldsymbol{b}^{(p)}\left(\mu\right)=\big{[}1,e^{\text{\rm j}\mu\,\rho_{2}^{(p)}},\ldots,e^{\text{\rm j}\mu\,\rho_{M_{p}}^{(p)}}\big{]}^{\text{\sf T}},$$

$$\mathbf{\phi }(\mu ,\mathbf{\eta })=[1,{\alpha }^{(2)}{\mathrm{e}}^{\text{j}\pi \mu {\eta }^{(2)}},\dots ,{\alpha }^{(P)}{\mathrm{e}}^{\text{j}\pi \mu {\eta }^{(P)}}{]}^{\text{T}}$$

$$\displaystyle\boldsymbol{B}(\boldsymbol{\mu})=\big{[}\boldsymbol{B}(\mu_{1}),\ldots,\boldsymbol{B}(\mu_{L})\big{]}$$

$$\Phi (\mathbf{\mu },\mathbf{\eta })$$

$$\mathbf{A}\left(\mathbf{\mu },\mathbf{\eta }\right)=\mathbf{B}(\mathbf{\mu })\Phi (\mathbf{\mu },\mathbf{\eta })$$

$$\mathbf{Y}=\mathbf{B}(\mathbf{\mu })\Phi (\mathbf{\mu },\mathbf{\eta })\Psi +\mathbf{N},$$

$$\mathbf{Y}=\mathbf{A}(\mathbf{\nu },\mathbf{\eta })\mathbf{X}+\mathbf{N}.$$

$$[\mathbf{X}{]}_{k,n}=$$

$$\underset{\mathbf{X}}{\min }\frac{1}{2}{\Vert \mathbf{A}(\mathbf{\nu },\mathbf{\eta })\mathbf{X}-\mathbf{Y}\Vert }_{\text{F}}^2+\lambda \sqrt{N}\Vert \mathbf{X}{\Vert }_{2,0},$$

$$\displaystyle\|\boldsymbol{X}\|_{2,0}=\big{|}\{k\,:\,\|\boldsymbol{x}_{k}\|_{2}>0\}\big{|}.$$

$$\underset{\mathbf{X}}{\min }\frac{1}{2}{\Vert \mathbf{A}(\mathbf{\nu },\mathbf{\eta })\mathbf{X}-\mathbf{Y}\Vert }_{\text{F}}^2+\lambda \sqrt{N}\Vert \mathbf{X}{\Vert }_{2,1}.$$

$$\Vert \mathbf{X}{\Vert }_{2,1}=\sum _{k=1}^K{\Vert {\mathbf{x}}_k\Vert }_2$$

$$\mathbf{Y}=\mathbf{B}(\mathbf{\nu })\Phi (\mathbf{\nu },\mathbf{\eta })\mathbf{X}+\mathbf{N},$$

$$\mathbf{Q}$$

$${\mathbf{Q}}_k=\mathbf{\phi }({\nu }_k,\mathbf{\eta }){\mathbf{x}}_k^{\text{T}},\text{ for }k=1,\dots ,K.$$

$$\mathbf{Y}=\mathbf{B}\mathbf{Q}+\mathbf{N},$$

$$\underset{\mathbf{Q}}{\min }\frac{1}{2}{\Vert \mathbf{B}\mathbf{Q}-\mathbf{Y}\Vert }_{\text{F}}^2+\lambda \sqrt{N}\sum _{k=1}^K\mathrm{rank}({\mathbf{Q}}_k).$$

$$\left\|\boldsymbol{Q}_{k}\right\|_{*}={\rm Tr}\big{(}(\boldsymbol{Q}_{k}^{\text{\sf H}}\boldsymbol{Q}_{k})^{1/2}\big{)}=\sum_{i=1}^{r}\sigma_{k,i},$$

$$\underset{\mathbf{Q}}{\min }\frac{1}{2}{\Vert \mathbf{B}\mathbf{Q}-\mathbf{Y}\Vert }_{\text{F}}^2+\lambda \sqrt{N}{\Vert \mathbf{Q}\Vert }_{\ast ,1},$$

$${\Vert \mathbf{Q}\Vert }_{\ast ,1}=\sum _{k=1}^K{\Vert {\mathbf{Q}}_k\Vert }_{\ast }.$$

$${\hat{\mathbf{Q}}}_k={\hat{\mathbf{U}}}_k{\hat{\mathbf{\Sigma }}}_k{\hat{\mathbf{V}}}_k^{\text{T}}\text{for }k=1,\dots ,K,$$

$${\hat{\mathbf{x}}}_k={\hat{\sigma }}_{k,1}[{\hat{\mathbf{u}}}_{k,1}{]}_1{\hat{\mathbf{v}}}_{k,1}\text{and}\hat{\mathbf{\phi }}({\nu }_k,\mathbf{\eta })=\frac{{\hat{\mathbf{u}}}_{k,1}}{[{\hat{\mathbf{u}}}_{k,1}{]}_1}.$$

$$\underset{\mathbf{Q}}{\min }\frac{1}{2}{\Vert \mathbf{B}\mathbf{Q}-\mathbf{Y}\Vert }_{\text{F}}^2+\lambda \sqrt{N}{\Vert \mathbf{Q}\Vert }_{\ast ,1}$$

$$\displaystyle\min_{\boldsymbol{S}\in\mathcal{B}^{P\,+}_{K}}{\rm Tr}\big{(}(\boldsymbol{B}\boldsymbol{S}\boldsymbol{B}^{\text{\sf H}}+\lambda\boldsymbol{I})^{-1}\hat{\boldsymbol{R}}\big{)}+{\rm Tr}(\boldsymbol{S}),$$

$$\hat{\mathbf{Q}}=\hat{\mathbf{S}}{\mathbf{B}}^{\text{H}}(\mathbf{B}\hat{\mathbf{S}}{\mathbf{B}}^{\text{H}}+\lambda \mathbf{I}{)}^{-1}\mathbf{Y}$$

Full text

1 Introduction

Direction finding with sensor arrays has applications in various fields of signal processing such as wireless communications, radar, sonar or astronomy. In these types of applications it is desired to achieve a high angular resolution and to identify a large number of sources. This can be achieved by sensor arrays with a large aperture and a large number of sensors [1]. However, a large aperture size makes it difficult to achieve and maintain precise array calibration. Possible reasons for imperfect calibration are inaccuracies in the sensor positions, timing synchronization errors, or other unknown gain and phase offsets among sensors [2, 3, 4, 5, 6, 7, 8]. Standard approaches to this problem usually rely on either offline or online calibration. Offline calibration of the overall array is performed using reference sources at known positions and can easily become a challenging and time consuming task [2]. Alternatively, several online calibration techniques have been proposed which use calibration sources at unknown positions [3, 4, 5, 6, 7], but the complexity of these techniques is prohibitively high, and performance can be severely limited in the case of large sensor position errors [7]. Moreover, these techniques cannot be employed in scenarios with imperfect time synchronization of sensors or other unknown sensor gain and phase offsets.

One way to overcome the calibration problem is to partition the overall array into smaller subarrays which are themselves comparably easy to calibrate. This type of array is referred to as partly calibrated array (PCA) and has received considerable interest in recent years. Generally, direction finding approaches for this type of arrays can be classified into non-coherent and coherent methods. In the non-coherent case the subarrays independently perform estimation of the directions of arrival (DoAs) or the signal covariance matrix to communicate these estimates to a central processor, where further processing is performed to achieve an improved joint estimate [9, 10, 11, 12, 13]. In the coherent approach parameter estimation is performed based on joint coherent processing of all available sensor measurements, e.g., by computing a global sample covariance matrix, and imperfect calibration among the different subarrays is taken account of in the estimation process [14, 15, 16, 17, 18, 19, 8, 20, 21]. In this work we consider the latter of the two approaches.

A prominent class of DoA estimation methods is based on subspace separation. In [14, 15, 16, 17, 18] the authors consider PCAs composed of multiple identical subarrays. Such types of array exhibit multiple shift invariances and methods such as the multiple invariant MUSIC and MODE [14] or multiple invariance ESPRIT [16, 17, 18] can be used to provide DoA estimates in a search-free fashion. In [19] it is assumed that the PCA is composed of identically oriented linear arrays that can be transformed to a uniform linear array by linear translations of the subarrays. The authors present the root-RARE algorithm which admits search-free DoA estimation. The root-RARE method in [19] was modified to the spectral RARE in [8] which admits application to arbitrary array topologies at the cost of increased computational complexity. Subspace-based methods are well investigated and are shown to asymptotically achieve an estimation performance close to the Cramér-Rao bound at low computational complexity. However, these subspace-based methods often have difficulties in certain practical scenarios. First, correlated source signals, e.g., in multipath environments, can significantly reduce the estimation performance. Second, subspace-based methods yield poor performance in the case of low number of snapshots and low signal-to-noise ratio, e.g., in fast changing environments.

Recently, sparse recovery (SR) methods came into focus of DoA estimation studies. As reported in [22, 23], SR methods provide high-resolution parameter estimation performance without the aforementioned shortcomings of subspace-based methods. Moreover, SR methods are computationally tractable since they can be formulated as convex optimization problems. While classical SR methods aim at recovering sparse signal vectors and admit grid-based parameter estimation [24, 25], the special case of fully calibrated arrays (FCAs) of uniform linear topology with possibly missing sensors allows for gridless SR methods as proposed in [26, 27, 28]. In [12, 13] the authors propose grid-based and gridless SR methods applicable for non-coherent processing in PCAs, where joint sparsity in the subarray signal representations is exploited. SR methods for coherent processing in PCAs have been presented in  citesteffens2017shiftinvariance, 6882328. The method in [20] is based on the recently proposed SPARROW formulation [23] and exploits multiple shift-invariances in PCAs composed of identical subarrays to provide gridless parameter estimation. In [21] the well-known $\ell_{2,1}$ mixed-norm minimization approach [22] for FCAs is generalized to grid-based SR in PCAs of arbitrary topology by means of a mixed nuclear norm [29, 30] and $\ell_{1}$ norm, termed here as $\ell_{*,1}$ mixed-norm. As shown by numerical experiments [21], $\ell_{*,1}$ mixed-norm minimization clearly outperforms the spectral RARE [8] in frequency resolution performance for low signal-to-noise ratio and low number of snapshots.

In this paper we consider the $\ell_{*,1}$ mixed-norm minimization problem proposed in [21] and derive an equivalent compact reformulation, termed as COmpact Block- and RAnk-Sparse recovery (COBRAS). The COBRAS formulation has a reduced number of optimization parameters as compared to the original $\ell_{*,1}$ mixed-norm minimization problem and we provide efficient implementations of the COBRAS formulation by means of semidefinite programming (SDP). While the SDP implementation is based on grid-based sampling of the subarray manifolds and applicable to arbitrary array topologies, we furthermore present a search-free implementation of our COBRAS formulation for the special case of linear subarrays with a common baseline. We show by extensive numerical experiments that the COBRAS approach outperforms the state of the art methods in difficult scenarios. In summary, our main contributions are given as:

• •

We introduce a sparse recovery approach for coherent processing in PCAs using $\ell_{*,1}$ mixed-norm minimization.

• •

We derive a compact reformulation of the $\ell_{*,1}$ mixed-norm minimization problem, termed as COBRAS.

• •

We develop a computationally efficient grid-based SDP implementation of the COBRAS formulation for arbitrary array topologies, and

• •

an efficient gridless SDP implementation of the COBRAS formulation for PCAs composed of subarrays with a common baseline.

The paper is organized as follows: Section 2 introduces the PCA signal model. The $\ell_{2,1}$ and $\ell_{*,1}$ mixed-norm minimization problems for FCAs and PCAs are discussed in Section 3. The COBRAS formulation is derived in Section 4 while grid-based and gridless SDP implementations are provided in Sections 5 and 6. Numerical results are presented in Section 7 before the paper is concluded in Section 8.

Notation: Boldface uppercase letters $\boldsymbol{X}$ denote matrices, boldface lowercase letters $\boldsymbol{x}$ denote column vectors, and regular letters $x,N$ denote scalars, with j denoting the imaginary unit. Superscripts $\boldsymbol{X}^{\text{\sf T}}$ and $\boldsymbol{X}^{\text{\sf H}}$ denote transpose and conjugate transpose of a matrix $\boldsymbol{X}$, respectively. The term $\mathcal{B}^{P\,+}_{K}$ denotes the set of positive semidefinite block-diagonal matrices composed of $K$ blocks of size $P\times P$ on the main diagonal . We write $[\boldsymbol{X}]_{m,n}$ to indicate the element in the $m$th row and $n$th column of matrix $\boldsymbol{X}$. The statistical expectation of a random variable $x$ is denoted as $\text{E}\{x\}$, and the trace of a matrix $\boldsymbol{X}$ is referred to as ${\rm Tr}(\boldsymbol{X})$. The Frobenius norm and the $\ell_{p,q}$ mixed-norm of a matrix $\boldsymbol{X}$ are referred to as $\|\boldsymbol{X}\|_{\text{\sf F}}$ and $\|\boldsymbol{X}\|_{p,q}$, respectively, while the $\ell_{p}$ norm of a vector $\boldsymbol{x}$ is denoted as $\|\boldsymbol{x}\|_{p}$. The term ${\rm diag}(x_{1},\ldots,x_{K})$ denotes a diagonal matrix with the elements $x_{1},\ldots,x_{K}$ on its main diagonal while ${\rm blkdiag}(\boldsymbol{X}_{1},\ldots,\boldsymbol{X}_{K})$ denotes a block-diagonal matrix composed of submatrices $\boldsymbol{X}_{1},\ldots,\boldsymbol{X}_{K}$ on its main block-diagonal.

2 Signal Model

Consider a linear array of arbitrary topology, composed of $M$ omnidirectional sensors, as depicted in Figure 1. Assume the overall array is partitioned into $P$ subarrays with $M_{p}$ sensors in subarray $p$, for $p=1,\ldots,P$, such that $M=\sum_{p=1}^{P}M_{p}$. We define $\boldsymbol{\eta}=[\eta^{(2)},\ldots,\eta^{(P)}]^{\text{\sf T}}$ as the vector containing the $P-1$ unknown inter-subarray displacements $\eta^{(2)},\ldots,\eta^{(K)}$ expressed in half signal wavelength and relative to the first subarray, i.e., $\eta^{(1)}=0$. Furthermore, let $\rho_{m}^{(p)}$, for $m=1,\ldots,M_{p}$, $p=1,\ldots,P$, denote the perfectly known intra-subarray position of the $m$th sensor of subarray $p$ relative to the first sensor in the subarray, hence $\rho_{1}^{(p)}=0$, and expressed in half signal wavelength. Consequently, the position of sensor $m$ in subarray $p$, relative to the first sensor in the first subarray, can be expressed as

[TABLE]

for $m=1,\ldots,M_{p}$ and $p=1,\ldots,P$.

Moreover, assume a number of $L$ narrowband and far-field sources are illuminating the sensor array from angular directions $\theta_{1},\ldots,\theta_{L}$, as illustrated in Figure 1. The corresponding spatial frequencies are defined as $\mu_{l}=\cos\theta_{l}\in[-1,1)$, for $l=1,\ldots,L$, and comprise the vector $\boldsymbol{\mu}=[\mu_{1},\ldots,\mu_{L}]^{\text{\sf T}}$. A total of $N$ signal snapshots are obtained at the output of each subarray $p$ and collected in the $M_{p}\times N$ subarray measurement matrix $\smash{\boldsymbol{Y}^{(p)}}$, for $p=1,\ldots,P$, where $\smash{[\boldsymbol{Y}^{(p)}]_{m,n}}$ denotes the output of the $m$th sensor in the $p$th subarray at time instant $n$. The subarray measurement matrices are collected in the $M\times N$ array measurement matrix $\boldsymbol{Y}=[\smash{\boldsymbol{Y}^{(1)\text{\sf T}}},\ldots,\smash{\boldsymbol{Y}^{(P)\text{\sf T}}}]^{\text{\sf T}}$, which is modeled as

[TABLE]

where $\boldsymbol{\varPsi}\in\mathbb{C}^{L\times N}$ is the source signal matrix and $\boldsymbol{N}\in\mathbb{C}^{M\times N}$ denotes a spatio-temporal white Gaussian sensor noise matrix. The $M\times L$ array steering matrix $\boldsymbol{A}\left(\boldsymbol{\mu},\boldsymbol{\eta}\right)$ in (2) is given by

[TABLE]

and represents the response of the entire array, where $\boldsymbol{a}(\mu,\boldsymbol{\eta})$ denotes the steering vector for spatial frequency $\mu$ and subarray displacements $\boldsymbol{\eta}$. Based on the sensor position definition in (1), the array steering vectors can be factorized as

[TABLE]

where the $M\times P$ block-diagonal matrix

[TABLE]

contains the perfectly known subarray steering vectors

[TABLE]

for $p=1,\ldots,P$, on its diagonal, and the $L\times 1$ vector

[TABLE]

takes account of the subarray displacement shifts ${\rm e}^{\text{\rm j}\pi\mu\eta^{(p)}}$, for $p=2,\ldots,P$, depending on the spatial frequencies in $\boldsymbol{\mu}$ and the subarray displacements in $\boldsymbol{\eta}$, and further unknown shifts $\alpha^{(p)}$, e.g., gain/phase or timing offsets among the subarrays [8]. In relation to (3), let us define the $M\times PL$ matrix

[TABLE]

containing all subarray responses for the spatial frequencies in $\boldsymbol{\mu}$, and the $PL\times L$ block-diagonal matrix

[TABLE]

composed of the subarray shift vectors in (7). Using (8) and (9), the overall array steering matrix (3) can be factorized as

[TABLE]

such that the overall array measurement matrix in (2) is equivalently modeled as

[TABLE]

which forms the basis for the $\ell_{*,1}$ mixed-norm minimization problem discussed in the following section.

3 State-of-the-Art

In this section we will shortly review the $\ell_{2,1}$ mixed-norm minimization approach for FCAs before turning to the $\ell_{*,1}$ mixed-norm minimization approach for PCAs.

3.1 Fully Calibrated Array

We first consider the case of an FCA where the subarray displacements in $\boldsymbol{\eta}$ are perfectly known. Based on the signal model in (2) we introduce a sparse representation of the measurement matrix as

[TABLE]

The $M\times K$ overcomplete dictionary matrix $\boldsymbol{A}(\boldsymbol{\nu},\boldsymbol{\eta})$ is obtained by sampling the field-of-view in $K\gg L$ spatial frequencies $\boldsymbol{\nu}=[\nu_{1},\ldots,\nu_{K}]^{\text{\sf T}}$. For ease of presentation we assume that the frequency grid is sufficiently fine, such that the true frequencies in $\boldsymbol{\mu}$ are contained in the frequency grid $\boldsymbol{\nu}$, i.e., $\{\mu_{l}\}_{l=1}^{L}\subset\{\nu_{k}\}_{k=1}^{K}$. In Section 6 we present an extension of our proposed formulation for subarrays with a common baseline which does not rely on the on-grid assumption. The $K\times N$ sparse signal matrix $\boldsymbol{X}$ in (12) contains elements

[TABLE]

for $k=1,\ldots,K$, $l=1,\ldots,L$. Thus $\boldsymbol{X}$ exhibits a row-sparse structure, i.e., the elements in a row of $\boldsymbol{X}$ are either jointly zero or primarily non-zero. Based on the sparse representation (12), the frequency estimation problem can be formulated as the mixed-norm minimization problem

[TABLE]

where $\lambda>0$ is a regularization parameter determining the sparsity, i.e., the number of non-zero rows in the minimizer $\smash{\hat{\boldsymbol{X}}}$. Row-sparsity is enforced by minimizing the $\ell_{2,0}$ mixed-norm in (14), which is defined as the number of non-zero rows $\boldsymbol{x}_{k}$ of the matrix $\boldsymbol{X}=[\boldsymbol{x}_{1},\ldots,\boldsymbol{x}_{K}]^{\text{\sf T}}$, i.e., the cardinality of the union support set according to

[TABLE]

Since the problem in (14) is NP-hard, several approximation methods have been proposed in the literature, including convex relaxation to the well-known $\ell_{2,1}$ mixed-norm minimization problem [31, 32, 22, 23]

[TABLE]

The $\ell_{2,1}$ mixed-norm in (16) is defined as

[TABLE]

and induces a non-linear coupling among the elements in each row $\boldsymbol{x}_{k}$, $k=1,\ldots,K$, of the matrix $\boldsymbol{X}$ such that the $\ell_{1}$ norm, i.e., the nonnegative summation, is performed on the $\ell_{2}$ norms of the rows in $\smash{\hat{\boldsymbol{X}}}$. Given a minimizer $\smash{\hat{\boldsymbol{X}}=[\hat{\boldsymbol{x}}_{1},\ldots,\hat{\boldsymbol{x}}_{K}]^{\text{\sf T}}}$ of (16), the frequency estimation problem reduces to finding the local maxima in the vector of the signal $\ell_{2}$ row-norms $\hat{\boldsymbol{x}}^{\ell_{2}}=[\|\hat{\boldsymbol{x}}_{1}\|_{2},\ldots,\|\hat{\boldsymbol{x}}_{K}\|_{2}]^{\text{\sf T}}$ and assigning the corresponding frequency grid points to the set $\{\hat{\mu}\}$ of estimated frequencies.

For the PCA case with uncertain array response $\boldsymbol{A}(\boldsymbol{\nu},\boldsymbol{\eta})$ due to the unknown displacements in $\boldsymbol{\eta}$, the $\ell_{2,1}$ mixed-norm minimization approach in (16) cannot be applied and a more sophisticated approach has to be devised, as discussed in the following subsection.

3.2 Partly Calibrated Array

Analogous to the FCA case in (12), we introduce a sparse representation of the signal model in (11) for the PCA case as

[TABLE]

where the row-sparse matrix $\boldsymbol{X}$ is defined similarly as for the FCA case in (13). Furthermore, the $M\times PK$ overcomplete subarray dictionary matrix $\boldsymbol{B}(\boldsymbol{\nu})$ and the $PK\times K$ overcomplete subarray shift matrix $\boldsymbol{\varPhi}(\boldsymbol{\nu},\boldsymbol{\eta})$ are defined in correspondence to (8) and (9), respectively.

In the PCA case, the inter-subarray displacements in $\boldsymbol{\eta}$ are unknown and thus represent additional estimation variables, hence the subarray shifts in $\boldsymbol{\varPhi}(\boldsymbol{\nu},\boldsymbol{\eta})$, which depend on the spatial frequencies in $\boldsymbol{\nu}$ and the subarray displacements $\boldsymbol{\eta}$, have to be appropriately included in the sparse estimation problem. To this end we introduce a model that couples among the variables $\boldsymbol{x}_{k}$ in the rows of $\boldsymbol{X}=[\boldsymbol{x}_{1},\ldots,\boldsymbol{x}_{K}]^{\text{\sf T}}$ and the subarray shifts in $\boldsymbol{\varphi}(\nu_{k},\boldsymbol{\eta})$, for $k=1,\ldots,K$. We define the $KP\times N$ extended signal matrix $\boldsymbol{Q}$ as

[TABLE]

containing the products of the subarray shifts and the signal waveforms. Note that in this formulation the number of the unknown complex-valued signal variables is increased to $KPN$ elements in the matrix $\boldsymbol{Q}$, as compared to the total $K(N+P-1)$ complex-valued unknowns in both $\boldsymbol{X}$ and $\boldsymbol{\varPhi}(\boldsymbol{\nu},\boldsymbol{\eta})$. On the other hand, due to the block structure of the subarray shift matrix $\boldsymbol{\varPhi}(\boldsymbol{\nu},\boldsymbol{\eta})$ as defined in (9), the matrix $\boldsymbol{Q}=[\boldsymbol{Q}_{1}^{\text{\sf T}},\ldots,\boldsymbol{Q}_{K}^{\text{\sf T}}]^{\text{\sf T}}$ in (19) enjoys a special structure as it is composed of $K$ stacked rank-one matrices

[TABLE]

Using the formulation in (19), the sparse representation for the PCA case in (18) is equivalently described by

[TABLE]

where for ease of presentation we use $\boldsymbol{B}=\boldsymbol{B}(\boldsymbol{\nu})$ to denote the dictionary matrix in (21) and throughout the paper. An SR approach to take account of the special structure of the signal matrix $\boldsymbol{Q}$ in (19) is given as

[TABLE]

The formulation in (22) takes twofold advantage of the sparsity assumption. First, minimization of the rank-terms encourages low-rank blocks $\hat{\boldsymbol{Q}}_{1},\ldots,\hat{\boldsymbol{Q}}_{K}$ in the minimizer $\hat{\boldsymbol{Q}}$. Second, minimizing the sum-of-ranks provides a block-sparse structure of $\hat{\boldsymbol{Q}}$, i.e., the elements in each block $\hat{\boldsymbol{Q}}_{k}$, for $k=1,\ldots,K$, are either jointly zero or primarily non-zero. However, the problem in (22) is NP-hard and computationally intractable.

The nuclear norm represents a tight convex approximation of the rank function and it has been successfully applied in a variety of rank minimization problems [29, 30]. The definition of the nuclear norm is given as

[TABLE]

where $r=\min(P,N)$ and $\sigma_{k,i}$ is the $i$th singular value of $\boldsymbol{Q}_{k}$. Along these lines it has been proposed in [21] to approximate the sparse estimation problem (22) by the following convex minimization problem

[TABLE]

where $\left\|\boldsymbol{Q}\right\|_{*,1}$ denotes the $\ell_{*,1}$ mixed-norm, computed as

[TABLE]

Similar to (22), the problem in (24) motivates low-rank blocks $\hat{\boldsymbol{Q}}_{1},\ldots,\hat{\boldsymbol{Q}}_{K}$ and a block-sparse structure in the minimizer $\hat{\boldsymbol{Q}}=[\hat{\boldsymbol{Q}}_{1}^{\text{\sf T}},\ldots,\hat{\boldsymbol{Q}}_{K}^{\text{\sf T}}]^{\text{\sf T}}$. Note that the PCA formulations in (22) and (24) reduce to the FCA formulations in (14) and (16), respectively, in the case of a single subarray, i.e., $P=1$.

Performing singular value decomposition on the matrix blocks in $\hat{\boldsymbol{Q}}$, i.e.,

[TABLE]

the signal waveform $\hat{\boldsymbol{x}}_{k}$ and subarray shifts $\hat{\boldsymbol{\varphi}}(\nu_{k},\boldsymbol{\eta})$ corresponding to the spatial frequency $\nu_{k}$ can be recovered according to

[TABLE]

The left and right singular vectors $\hat{\boldsymbol{u}}_{k,1}$ and $\hat{\boldsymbol{v}}_{k,1}$ in (27) correspond to the largest singular value $\hat{\sigma}_{k,1}$ of $\hat{\boldsymbol{Q}}_{k}$ and normalization to the first element $[\hat{\boldsymbol{u}}_{k,1}]_{1}$ of $\hat{\boldsymbol{u}}_{k,1}$ in (27) is performed to take account of the structure of the subarray shift vectors $\hat{\boldsymbol{\varphi}}(\nu_{k},\boldsymbol{\eta})$ according to (7).

While the $\ell_{*,1}$ mixed-norm minimization problem in (24) provides a tractable approach for SR in PCAs, it suffers from high computational complexity in the case of a large number of snapshots $N$ and grid points $K$. To overcome this difficulty we provide in the following section a compact reformulation of problem (24).

4 Compact Block- and Rank-Sparse Recovery

One of the main results of this paper is formulated in the following theorem:

Theorem 1** (Problem Equivalence).**

The block- and rank-sparsity inducing $\ell_{*,1}$ mixed-norm minimization problem

[TABLE]

is equivalent to the convex problem

[TABLE]

with $\hat{\boldsymbol{R}}=\boldsymbol{Y}\boldsymbol{Y}^{\text{\sf H}}/N$ and $\mathcal{B}^{P\,+}_{K}$ denoting the sample covariance matrix and the set of positive semidefinite block-diagonal matrices composed of $K$ blocks of size $P\times P$, respectively. The equivalence holds in the sense that a minimizer $\smash{\hat{\boldsymbol{Q}}}$ for problem (28) can be factorized as

[TABLE]

where $\smash{\hat{\boldsymbol{S}}}$ is a minimizer for problem (29).

A proof of the equivalence is provided in Appendix Acknowledgment, while a proof of the convexity of (29) is provided in Section 5.2 by establishing equivalence to a semidefinite program.

In addition to relation (30), it can be shown (see Appendix Acknowledgment) that a minimizer $\smash{\hat{\boldsymbol{S}}}={\rm blkdiag}(\smash{\hat{\boldsymbol{S}}}_{1},\ldots,\smash{\hat{\boldsymbol{S}}}_{K})$ of (29) relates to the signal matrix $\smash{\hat{\boldsymbol{Q}}}=[\smash{\hat{\boldsymbol{Q}}}_{1}^{\text{\sf T}},\ldots,\smash{\hat{\boldsymbol{Q}}}_{K}^{\text{\sf T}}]^{\text{\sf T}}$ according to

[TABLE]

for $k=1,\ldots,K$, such that the block-support of $\smash{\hat{\boldsymbol{Q}}}$ is equivalently represented by the block-support of the matrix $[\smash{\hat{\boldsymbol{S}}}_{1}^{\text{\sf T}},\ldots,\smash{\hat{\boldsymbol{S}}}_{K}^{\text{\sf T}}]^{\text{\sf T}}$. Similarly, the rank of the matrix blocks $\smash{\hat{\boldsymbol{Q}}}_{k}$ is equivalently represented by the matrix blocks $\smash{\hat{\boldsymbol{S}}}_{k}$, i.e., ${\rm rank}(\smash{\hat{\boldsymbol{Q}}}_{k})={\rm rank}(\smash{\hat{\boldsymbol{S}}}_{k})$, for $k=1,\ldots,K$.

We observe that the problem in (29) only relies on the measurement matrix $\boldsymbol{Y}$ through the sample covariance matrix $\hat{\boldsymbol{R}}$, leading to a significantly reduced problem size, especially in the case of large number of snapshots $N$. In this context we term the formulation in (29) as COmpact Block- and RAnk-Sparse recovery (COBRAS). The compact formulation (29) contains $KP^{2}$ real-valued optimization parameters in the positive semidefinite matrix $\boldsymbol{S}$, as opposed to the $2KPN$ real-valued optimization parameters in $\boldsymbol{Q}$ in problem (28). Consequently, in the case of a large number of snapshots $N>P/2$, the reformulation (29) has reduced computational complexity as compared to (28).

5 Implementation of the COBRAS Formulation

The $\ell_{*,1}$ mixed-norm minimization problem (28) has been well investigated in literature and implementations based on the coordinate descent method [21], the STELA algorithm [33] and semidefinite programming (SDP) [34] have been proposed. Here we will shortly revise the SDP implementation of the $\ell_{*,1}$ mixed-norm minimization problem (28) to highlight the reduction in computational complexity obtained by employing the COBRAS formulation in (29).

5.1 SDP Form of the $\ell_{*,1}$ Mixed-Norm Minimization Problem

As discussed in [29], minimization of the nuclear norm

[TABLE]

for some convex set $\mathcal{C}$, can be expressed as the SDP

[TABLE]

where $\boldsymbol{P}_{k,1}=\boldsymbol{P}_{k,1}^{\text{\sf H}}$ and $\boldsymbol{P}_{k,2}=\boldsymbol{P}_{k,2}^{\text{\sf H}}$ are auxiliary variables of size $P\times P$ and $N\times N$, respectively. The SDP formulation (33) admits simple implementation of the nuclear norm minimization problem using standard convex solvers, such as SeDuMi [35].

Based on the equivalence of (32) and (33), the $\ell_{*,1}$ mixed-norm minimization problem (28) can be equivalently formulated as

[TABLE]

Note that with the auxiliary variables in $\boldsymbol{P}_{k,1}$ and $\boldsymbol{P}_{k,2}$, for $k=1,\ldots,K$, the problem (34) has $K(P+N)^{2}$ real-valued optimization variables as opposed to the problem formulation in (28) which has $2KPN$ real-valued optimization variables. For a large number of grid points $K$ or snapshots $N$ the SDP formulation (34) becomes intractable and alternative implementations are required as presented in the next subsection. We remark that the problem of large snapshot number has been addressed in previous literature by matching the signal subspace of the measurements $\boldsymbol{Y}$ instead of the measurements itself, see [22, 34], leading to a reduced number of effective signal snapshots, however at the expense of potential performance degradation, e.g., in the case of correlated source signals.

5.2 SDP Form of the COBRAS Method

In order to solve the COBRAS formulation in (29) by means of a tractable SDP which can be treated by standard convex solvers consider the following corollaries [36]:

Corollary 1**.**

The COBRAS formulation in (29) is equivalent to the convex semidefinite program

[TABLE]

where $\boldsymbol{Z}_{N}$ is a Hermitian matrix of size $N\times N$.

To see the equivalence between the two problems we note that $\boldsymbol{B}\boldsymbol{S}\boldsymbol{B}^{\text{\sf H}}+\lambda\boldsymbol{I}\succ\boldsymbol{0}$ is positive definite for any $\lambda>0$ and consider the Schur complement of the constraint (35b)

[TABLE]

which implies

[TABLE]

Since in problem (35) ${\rm Tr}(\boldsymbol{Z}_{N})$ is minimized, it can be proved by contradiction that the relation in (37) must hold with equality, proving the equivalence of (29) and (35).

Corollary 2**.**

The COBRAS formulation in (29) admits the equivalent problem formulation

[TABLE]

The proof to Corollary 2 follows similar arguments as in the proof of Corollary 1 and is therefore omitted here. In contrast to (35), the size of the semidefinite constraint in (38) is independent of the number of snapshots $N$. It follows that either problem formulation (35) or (38) can be selected to solve (29), depending on the number of snapshots $N$ and the resulting size of the semidefinite constraint. The problems (35) and (38) have $KP^{2}$ real-valued optimization variables in $\boldsymbol{S}$ and additional $N^{2}$ or $M^{2}$ real-valued parameters in $\boldsymbol{Z}_{N}$ and $\boldsymbol{Z}_{M}$, respectively. Thus, in the undersampled case $N<M$ it is preferable to use the SDP formulation in (35), while in the oversampled case $N\geq M$ it is preferable to apply the SDP formulation in (38). We remark that the subspace matching approach discussed in [22, 34] can be applied to formulation (35) as well. A further investigation of the subspace matching approach is, however, beyond the scope this paper. The various equivalent SDP implementations and the corresponding number of variables and constraints are listed in Table 1.

6 Gridless COBRAS Implementation

While the SDP formulations in Section 5 are applicable to arbitrary array topologies, we consider in the following the special case of linear subarrays with a common baseline, where the sensors within each subarray are located at integer multiples of a baseline $\delta$, i.e., $\rho_{m}^{(p)}=\delta d_{m}^{(p)}$ with $d_{m}^{(p)}\in\mathbb{Z}$ for $m=1,\ldots,M_{P}$ and $p=1,\ldots,P$. This type of array topologies admits the extension of the COBRAS formulation to gridless frequency estimation.

We start noting that strong duality holds for problem (38) and consider the Lagrange dual problem, which is given as

[TABLE]

where $\boldsymbol{\varUpsilon}_{0}$ is an $M\times M$ positive semidefinite matrix and $\boldsymbol{\varUpsilon}_{1}$ is of size $M\times M$ and does not exhibit specific structure. Complementary slackness requires that

[TABLE]

for $k=1,\ldots,K$, i.e., if $\boldsymbol{S}_{k}\neq\boldsymbol{0}$ then $\boldsymbol{I}_{P}-\boldsymbol{B}^{\text{\sf H}}(\nu_{k})\,\boldsymbol{\varUpsilon}_{0}\boldsymbol{B}(\nu_{k})$ must be singular, such that

[TABLE]

Condition (41) indicates that instead of solving the primal problem (29) and identifying the block-support from $\boldsymbol{S}$, we can equivalently solve the dual problem (39) and identify the block-support from the roots of (41).

Let us consider the limiting case of an infinitesimal frequency grid spacing, i.e., $\lim_{K\rightarrow\infty}\nu_{k}-\nu_{k-1}=0$ for $k=2,\ldots,K$, such that the frequency becomes a continuous parameter $\nu$. By introducing the variable

[TABLE]

the subarray steering matrices in (5) can be equivalently described as

[TABLE]

where the subarray steering vectors are given as

[TABLE]

with $d_{m}^{(p)}\in\mathbb{Z}$ for $m=1,\ldots,M_{p}$ and $p=1,\ldots,P$. By the definition in (43), the matrix product $\boldsymbol{B}^{\text{\sf H}}(z)\boldsymbol{\varUpsilon}_{0}\boldsymbol{B}(z)$ in constraint (39c) constitutes a trigonometric matrix polynomial of degree $D=\max_{p,m_{p}}d^{(p)}_{m_{p}}$, according to

[TABLE]

with matrix coefficients $\boldsymbol{K}_{i}$ of size $P\times P$ [37]. In the continuous case the constraint (39c) is replaced by the constraint

[TABLE]

which provides an upper bound on the matrix polynomial (45) and can be implemented by semidefinite programming, e.g., by the problem formulation (84) derived in Appendix Appendix A - Proof of Problem Equivalence or by other techniques discussed in [37]. Once the continuous implementation of problem (39) is solved, the spatial frequencies can be recovered by finding the roots for which the left-hand side of (46) becomes singular, e.g., by rooting the continuous counterpart of (41) using the techniques discussed in [38].

6.1 Related Work

In a recent work [20] it has been shown that PCAs composed of identical subarrays admit gridless compressed sensing by means of the SPARROW formulation [23]. Interestingly, for the special case considered in this section of PCAs composed of linear subarrays with a common baseline (and possibly missing sensors in particular subarrays) the dual problem formulation (39) can equivalently be derived by means of the SPARROW formulation.

Let us consider the $\ell_{2,1}$ mixed-norm minimization problem for FCAs given in (16) which can equivalently be formulated as the SPARROW problem [23]

[TABLE]

where $\tilde{\lambda}>0$ is a regularization parameter and $\boldsymbol{S}={\rm diag}(s_{1},\ldots,s_{K})\succeq\boldsymbol{0}$ is of size $K\times K$. Problem (47) can be formulated as the semidefinite program

[TABLE]

The Lagrange dual problem of (48) is given as

[TABLE]

and with strong duality it follows from complementary slackness that

[TABLE]

As previously discussed in the context of condition (41), the support of the vector $\boldsymbol{s}=[s_{1},\ldots,s_{K}]^{\text{\sf T}}$, i.e., the spatial frequency estimates, can equivalently be identified by rooting the function in (50).

Making use of the notation $\boldsymbol{a}(\nu_{k},\boldsymbol{\eta})=\boldsymbol{B}(\nu_{k})\boldsymbol{\varphi}(\nu_{k},\boldsymbol{\eta})$, as introduced in (4), condition (50) can be rewritten as

[TABLE]

where

[TABLE]

Condition (51) is fulfilled if

[TABLE]

which is identical to the constraint (39c) in problem (39). Replacing the constraint (49c) in problem (49) by the condition (54) and further using (53) and the substitutions $\lambda=\tilde{\lambda}/\|\boldsymbol{\varphi}(\nu_{k},\boldsymbol{\eta})\|_{2}^{2}$ and $\boldsymbol{\varUpsilon}_{1}=\tilde{\boldsymbol{\varUpsilon}}_{1}$ shows that for the PCA case the dual problem (49) can be reformulated as the dual problem (39). As demonstrated in the previous section, condition (54) can be extended to an infinitesimal grid spacing, resulting in a matrix polynomial constraint, such that the resulting gridless estimation problem can be implemented by semidefinite programming (see Appendix Appendix A - Proof of Problem Equivalence).

7 Numerical Results

For experimental performance evaluation of our proposed COBRAS method we compare its estimation performance to the state-of-the-art methods spectral RARE [8] and root-RARE [19] as well as the Cramér-Rao bound (CRB) [8].

For all simulations we use circular complex Gaussian source signals $\boldsymbol{\varPsi}$ with covariance matrix $\text{E}(\boldsymbol{\varPsi}\boldsymbol{\varPsi}^{\text{\sf H}})=N\boldsymbol{I}$, if not specified otherwise. We further consider spatio-temporal white circular complex Gaussian sensor noise $\boldsymbol{N}$ with covariance matrix $\text{E}(\boldsymbol{N}\boldsymbol{N}^{\text{\sf H}})=\sigma^{2}N\boldsymbol{I}$ and define the signal-to-noise ratio (SNR) as ${\rm SNR}~{}=~{}1/\sigma^{2}$. The vector $\boldsymbol{r}^{(p)}=[r_{1}^{(p)},\ldots,r_{M_{p}}^{(p)}]^{\text{\sf T}}$ contains the global sensor positions $\smash{r_{m}^{(p)}}$ of subarray $p$, for $m=1,\ldots,M_{p}$, $p=1,\ldots,P$, expressed in half-wavelength, as defined in (1). If not stated otherwise, we perform $T=1000$ Monte Carlo trials for each experimental setup and compute the statistical error.

The estimation performance of the COBRAS method strongly depends on proper selection of a regularization parameter $\lambda$. While regularization parameter selection is a research field of its own, in this paper we follow a heuristic approach and select the regularization parameter as

[TABLE]

which has shown good estimation performance in all investigated scenarios.

We remark that the RARE and COBRAS method make different assumptions on the availability of a-priori knowledge. While the RARE method requires knowledge of the number of source signals, the regularization parameter selection for the COBRAS method according to (55) requires knowledge of the noise power. However, since estimation of these parameters itself might affect the frequency estimation performance of the RARE and COBRAS methods, we apply the standard assumption of perfectly known number of source signals and noise power and investigate the achievable performance under these idealized assumptions.

7.1 Arbitrary Array Topologies and Grid-Based Estimation

In the first scenario we consider a PCA with a large aperture, composed of $M=11$ sensors which are partitioned in $P=4$ linear subarrays with 3,2,3, and 3 sensors, respectively. The sensor positions for each subarray are $\boldsymbol{r}^{(1)}=[0.0,0.6,2.3]^{\text{\sf T}}$, $\boldsymbol{r}^{(2)}=[12.2,13.0]^{\text{\sf T}}$, $\boldsymbol{r}^{(3)}=[21.5,22.8,23.6]^{\text{\sf T}}$, and $\boldsymbol{r}^{(4)}=[37.6,38.5,41.1]^{\text{\sf T}}$, and we assume no additional gain/phase offsets among the subarrays, i.e., $\boldsymbol{\alpha}=[1,1,1,1]^{\text{\sf T}}$ in (7). We further consider $L=3$ uncorrelated Gaussian source signals with spatial frequencies $\boldsymbol{\mu}=[0.5011,0.4672,-0.2007]^{\text{\sf T}}$.

The array topology does not admit a direct implementation of the gridless COBRAS and the root-RARE methods such that we limit the experiments in this subsection to the investigation of the grid-based COBRAS method and the spectral RARE method. For both grid-based methods we use a gird of $K=400$ grid points according to $\boldsymbol{\nu}=[-1.000,-0.995,-0.999,\ldots,0.995]^{\text{\sf T}}$.

To investigate the frequency estimation performance, we compute the root-mean-square error of the frequency estimates in $\hat{\boldsymbol{\mu}}$ as

[TABLE]

where $\hat{\mu}_{l}(t)$ denotes the frequency estimate of signal $l$ in trial $t$ and $|\mu_{1}-\mu_{2}|_{\rm wa}=\min_{i\in\mathbb{Z}}|\mu_{1}-\mu_{2}+2i|$ denotes the wrap-around distance for two frequencies $\mu_{1},\mu_{2}\in[-1,1)$. Since the ${\rm RMSE}$ computation (56) requires the number of estimated source signals $\hat{L}$ to be equal to the true number of source signals $L$, we have to consider two special cases: in the case of overestimation of the model order, $\hat{L}>L$, we select the $L$ frequency estimates with the largest corresponding magnitudes, whereas we select $L-\hat{L}$ additional random spatial frequencies in the case of underestimation $\hat{L}<L$.

In the first experiment the signal-to-noise ratio (${\rm SNR}$) is fixed to ${\rm SNR}=6{\rm\,dB}$, while the number of snapshots $N$ is varied. Figure 3 clearly demonstrates that our proposed grid-based COBRAS technique outperforms the spectral RARE for low number of signal snapshots $N$. While the spectral RARE method is not able to always resolve the two closely spaced signals with spatial frequencies $\mu_{1}=0.5011$ and $\mu_{2}=0.4672$ for $N\leq 500$ signal snapshots, our proposed COBRAS method resolves the signals for any $N\geq 30$ snapshots.

In a second experiment we fix the number of snapshots as $N=20$ and vary the SNR. As can be observed from Figure 3 the grid-based COBRAS method shows superior threshold performance as compared to the spectral RARE. While the spectral RARE can reliably resolve the two closely spaced sources only for ${\rm SNR}\geq 22{\rm\,dB}$, our proposed COBRAS can do so for ${\rm SNR}\geq 8{\rm\,dB}$. For high ${\rm SNR}$, spectral RARE reaches a bias in the RMSE which is caused mainly by the finite grid. A similar bias effect can be observed for the grid-based COBRAS method. However, for the grid-based COBRAS method the bias is larger than for the spectral RARE method and it is not only caused by the finite grid, as discussed in the following subsection.

7.2 Resolution Performance and Estimation Bias

For further investigation of the spatial frequency estimation bias, we consider a uniform linear array of $M=9$ sensors, partitioned into $P=3$ identical, uniform linear subarrays of 3 sensors each, without additional gain/phase offsets, i.e., $\boldsymbol{\alpha}=[1,1,1]^{\text{\sf T}}$ in (7). For the experiment we consider $L=2$ uncorrelated signals and fix the spatial frequency of the first signal as $\mu_{1}=0.505$ while the spatial frequency of the second signal is varied according to $\mu_{2}=\mu_{1}-\Delta\mu$ with $10^{-2}\leq\Delta\mu\leq 1$. For all grid-based estimation methods we make use of a uniform grid of $K=200$ points according to $\boldsymbol{\nu}=[-1,-0.99,-0.98,\ldots,0.99]^{\text{\sf T}}$. The SNR and number of snapshots are fixed as ${\rm SNR}=0{\rm\,dB}$ and $N=20$.

First, we observe from Figure 6 that the spectral RARE performs significantly worse than root-RARE in terms of threshold performance, i.e., the spectral RARE cannot always resolve the two signals for a frequency separation of $\Delta\mu\lessapprox 0.4$ while the root-RARE can resolve the signals for $\Delta\mu\gtrapprox 0.12$. The reason for this difference in resolution performance is that the root-RARE method locates the roots of the corresponding matrix polynomial in the entire complex plane, while the spectral RARE only searches minima on the unit circle (see also [39, 40]). In contrast to that, the grid-based and the gridless COBRAS methods both show rather similar estimation performance, comparable to that of the root-RARE method, and reach the CRB for sufficiently large frequency separation. This observation can be explained by the fact that both dual COBRAS optimization problems provide matrix polynomials with the roots of interest constrained on the unit circle, as discussed in Section 6. This explains the similar performance results of the grid-based and gridless COBRAS methods. The only difference between the grid-based and gridless COBRAS methods is that in the first case the roots are generated on a grid of candidate frequencies on the unit circle, while in the latter case the roots are continuously located on the unit circle.

In a slightly modified experiment we fix the SNR and number of snapshots to ${\rm SNR}=20{\rm\,dB}$ and $N=50$, respectively. While the root-RARE method performs close to the CRB for the region of interest, the spectral RARE can not always resolve the signals for $\Delta\mu\lessapprox 0.06$ and reaches an estimation bias for large source separation, which is caused by the finite frequency grid. Furthermore, it can be observed that the estimation performance of the COBRAS methods deviates from that of the root-RARE method. For large frequency separation $\Delta\mu\gtrapprox 0.2$ the grid-based COBRAS method reaches the grid bias, similar to the spectral RARE. However, also for low frequency separation $\Delta\mu\lessapprox 0.2$ both methods do not reach the CRB. This can be explained by an inherent frequency estimation bias for SR methods (see also [22, 23]). For further investigation we compute the spatial frequency estimation bias as

[TABLE]

where the mean estimate for spatial frequency $\mu_{l}$ is computed as ${\rm Mean}(\hat{\mu}_{l})=1/T\sum_{t=1}^{T}\hat{\mu}_{l}(t)$.

For the given scenario, the estimation bias is displayed in Figure 6. In the case of low frequency separation $\Delta\mu\lessapprox 0.2$ both COBRAS methods show a relatively large estimation bias of $\text{Bias}(\hat{\boldsymbol{\mu}})\approx 0.01$. For larger frequency separation $\Delta\mu\gtrapprox 0.2$, the bias of the grid-based COBRAS method is mainly determined by the finite grid, while the bias of the gridless COBRAS method shows to be periodic in $\Delta\mu$. In difficult scenarios, with low SNR and low number of snapshots as for the previous setup, the estimation bias is below the CRB, such that it is negligible in the RMSE performance. The frequency estimation bias is a well known phenomenon in SR research [22, 23] and bias mitigation techniques have been discussed, e.g., in [41].

7.3 Correlated Signals

As discussed in the previous subsection, gridless COBRAS and root-RARE show approximately equal resolution performance for uncorrelated signals in difficult scenarios with low SNR, low number of snapshots and uncorrelated signals. This situation changes in the case of correlated signals, where preprocessing in form of subspace separation, as required for the RARE method, becomes difficult. For further investigation of this aspect we consider a PCA of $M=9$ sensors partitioned into $P=3$ subarrays of 3,4 and 2 sensors with positions $\boldsymbol{r}^{(1)}=[0,1,3]^{\text{\sf T}}$, $\boldsymbol{r}^{(2)}=[17.4,18.4,19.4,21.4]^{\text{\sf T}}$ and $\boldsymbol{r}^{(3)}=[24.8,25.8]^{\text{\sf T}}$. Furthermore we consider gain/phase offsets among the subarrays according to $\boldsymbol{\alpha}=[1,\,0.7\cdot{\rm e}^{\text{\rm j}\frac{2}{3}\pi},\,1.2\cdot{\rm e}^{\text{\rm j}\frac{1}{4}\pi}]^{\text{\sf T}}$ in (7). The SNR and number of snapshots are selected as ${\rm SNR}=0{\rm\,dB}$ and $N=30$. We consider $L=2$ source signals with spatial frequencies $\boldsymbol{\mu}=[0.505,\,0.105]^{\text{\sf T}}$ and a source covariance matrix given as

[TABLE]

where the correlation coefficient $\rho$ is assumed to be real-valued and varied in the experiment. For the grid-based estimation methods we consider a grid of $K=200$ candidate frequencies, defined as in the previous subsection. As seen from Figure 7, the spectral and root-RARE methods fail to properly estimate the spatial frequencies for high correlation ($\rho>0.6$) while the grid-based and gridless COBRAS methods still show estimation performance close to the CRB, since these methods do not require subspace separation.

7.4 Array Calibration Performance

Besides estimation of the spatial frequencies, the COBRAS method also admits estimation of the subarray shifts in $\boldsymbol{\varphi}$ as defined in (7). Since the RARE methods do not provide direct estimation of the subarray shifts, we use the method presented in [15] in equation (11) 111Without the restriction that the complex phase terms must be of unit magnitude. on the basis of the spatial frequency estimates obtained by the RARE methods.

The setup under investigation consists of a PCA of $M=10$ sensors partitioned into $P=4$ subarrays of 3,2,3 and 2 sensors at positions $\boldsymbol{r}^{(1)}=[0,2,3]^{\text{\sf T}}$, $\boldsymbol{r}^{(2)}=[10.1,\,11.1]^{\text{\sf T}}$, $\boldsymbol{r}^{(3)}=[27.4,\,28.4,\,30.4]^{\text{\sf T}}$ and $\boldsymbol{r}^{(4)}=[54.8,\,56.8]^{\text{\sf T}}$. The subarray gain/phase offsets are set as $\boldsymbol{\alpha}=[1,\,1.3\cdot{\rm e}^{\text{\rm j}\frac{2}{3}\pi},\,0.7\cdot{\rm e}^{-\text{\rm j}\frac{1}{4}\pi},\,0.9\cdot{\rm e}^{-\text{\rm j}\frac{3}{5}\pi}]^{\text{\sf T}}$ in (7). We consider $L=3$ uncorrelated source signals with spatial frequencies $\boldsymbol{\mu}=[0.605,\,0.255,\,-0.305]^{\text{\sf T}}$ and the number of snapshot is set to $N=20$.

Figure 9 displays the frequency estimation error of the different methods for varying SNR, where both COBRAS methods show the best thresholding performance but reach an estimation bias for ${\rm SNR}\geq 15{\rm\,dB}$. Similarly, the spectral RARE algorithm reaches an estimation bias which is caused by the finite grid. On the other hand, the root-RARE performs asymptotically optimal and reaches the CRB for high SNR. The corresponding subarray shift estimation performance is displayed in Figure 9, where the root-mean-square error is computed according to

[TABLE]

with $\hat{\boldsymbol{\varphi}}_{l}(t)$ being the displacement phase vector estimate for signal $l$ in Monte Carlo trial $t$. As can be observed from Figure 9, the subarray shift estimation method in [15], based on the frequency estimates obtained from the RARE methods, achieves a relatively large estimation bias for high SNR. In contrast to that, the grid-based and gridless COBRAS methods show a significantly reduced estimation error, which demonstrates the advantage of joint frequency and displacement phase estimation.

7.5 Computational Complexity

To investigate the computation time of the COBRAS formulation, we perform simulations in Matlab using the SeDuMi solver [35] with the CVX interface [42, 43] on a machine with an Intel Core i5-760 CPU @ $2.80\;{\rm GHz}\times 4$ and $8\,{\rm GByte}$ RAM. We consider a scenario with two independent complex Gaussian sources with static spatial frequencies $\mu_{1}=0.505$ and $\mu_{2}=-0.205$ and a uniform linear PCA of $M=9$ sensors partitioned into $P=3$ identical and uniform linear subarrays of 3 sensors. We neglect subarray gain/phase offsets, i.e., $\boldsymbol{\alpha}=[1,1,1]^{\text{\sf T}}$ in (7).

For the first experiment the SNR is fixed at ${\rm SNR}=0{\rm\,dB}$ while the number of snapshots $N$ is varied. Figure 11 shows the average computation time for $T=100$ Monte Carlo runs of the SDP implementation of the $\ell_{*,1}$ mixed-norm minimization (34) and the grid-based COBRAS formulations (35) and (38) with a grid size of $K=100$, as well as the gridless (GL-) COBRAS formulation in (84). The computation time is measured only for solving the corresponding optimization problem in CVX. Pre-processing steps, such as computation of the sample covariance matrix, or post-processing steps, such as peak-search or polynomial rooting, are not included into this consideration. As can be observed from Figure 11, for a number of $N<5$ snapshots all grid-based methods exhibit approximately equal computation time. For $5\leq N<40$ the $\ell_{*,1}$ mixed-norm minimization problem has largest computation time while the COBRAS formulation (35) requires longest computation time for $N>40$, due to the large dimension of the semidefinite constraint (35b). Regarding the computation time of the grid-based COBRAS formulation using the sample covariance matrix (38) we observe that it is relatively constant for any number of snapshots $N$ and lower than for the other implementations especially for large number of snapshots $N>10$. The lowest computation time is required for solving the GL-COBRAS implementation (84).

Figure 11 shows the average computation time for $T=100$ Monte Carlo runs for a varying number of grid points $K$ and a fixed number of $N=9$ signal snapshots, corresponding to the case where the signal subspace matching techniques are applied (compare [22, 34]). For all grid-based methods the number of SDP constraints grows nearly linear with the number of grid points $K$. Clearly, the $\ell_{*,1}$ mixed-norm minimization approach has the largest computation time for any investigated number of grid points $K$. The grid-based COBRAS formulations (35) and (38) show approximately equal computation time, since both formulations have SDP constraints of identical dimension. Since the GL-COBRAS formulation is independent of the grid size $K$, it is constant for all grid size numbers $K$ in Figure 11 and provides the fastest computation of all methods under investigation.

8 Conclusion

Partly calibrated arrays are attractive setups for direction of arrival estimation. Computationally efficient subspace-based methods such as RARE and ESPRIT show asymptotically optimal performance, but have problems in difficult scenarios such as low sample size, low signal-to-noise ratio or correlated source signals. Sparse recovery in the form of $\ell_{*,1}$ mixed-norm minimization has been shown to be an attractive alternative to subspace-based methods in these difficult scenarios. In this paper we have derived a compact, equivalent formulation of the $\ell_{*,1}$ mixed-norm minimization, referred to as COmpact Block- and RAnk-Sparse recovery (COBRAS). The COBRAS formulation is attractive especially in the case of a large number of signal snapshots. For the special case of subarrays with common baseline we have presented an extension to gridless estimation, referred to as gridless COBRAS (GL-COBRAS).

As shown by numerical results, the grid-based COBRAS significantly outperforms the spectral RARE method in terms of thresholding performance for closely spaced source signals. Furthermore, the COBRAS method outperforms the RARE method in the case of strongly correlated source signals and in the calibration (subarray shift estimation) performance. A drawback of the $\ell_{*,1}$ mixed-norm minimization approach and the COBRAS formulation is the estimation bias which becomes significant in the asymptotic case of large number of snapshots or high signal-to-noise ratio. However, if higher estimation accuracy is required, the COBRAS estimates can be used to provide initial estimates, e.g., for a subsequent maximum likelihood estimator.

Acknowledgment

This work was supported by the German Research Foundation (DFG) within the DFG priority program on Compressed Sensing in Information Processing (CoSIP DFG-SPP 1798).

Appendix A - Proof of Problem Equivalence

The proof of Theorem 1 relies on the following lemma, see [44, 30]:

Lemma 1**.**

The nuclear norm of the $P\times N$ matrix $\boldsymbol{Q}_{k}$ is equivalently computed by the minimization problem

[TABLE]

where $\boldsymbol{\varGamma}_{k}$ and $\boldsymbol{G}_{k}$ are complex matrices of dimensions $P\times r$ and $r\times N$, respectively, with $r=\min(N,P)$.

Proof of Lemma 1.

Let us define the compact singular value decomposition of matrix $\boldsymbol{Q}_{k}$ as

[TABLE]

such that the factorization terms of $\boldsymbol{Q}_{k}=\boldsymbol{\varGamma}_{k}\boldsymbol{G}_{k}$ can be expressed as

[TABLE]

for $\boldsymbol{\varSigma}_{k}=\boldsymbol{\varPi}_{k,1}\boldsymbol{\varPi}_{k,2}$ of size $r\times r$ and some $r\times r$ arbitrary unitary matrix $\boldsymbol{W}_{k}$, i.e., $\boldsymbol{W}_{k}^{\text{\sf H}}\boldsymbol{W}_{k}=\boldsymbol{I}_{r}$. Based on (62) it holds that

[TABLE]

where the inequality stems from the Cauchy-Schwartz inequality and is fulfilled with equality if and only if $\boldsymbol{\varPi}_{k,1}=\boldsymbol{\varPi}_{k,2}=\boldsymbol{\varSigma}^{1/2}_{k}$. In this case, the matrix factors in (62) are given as

[TABLE]

Furthermore, by the arithmetic-geometric-mean inequality it follows that

[TABLE]

where equality holds if $\|\boldsymbol{\varGamma}_{k}\|_{\text{\sf F}}=\|\boldsymbol{G}_{k}\|_{\text{\sf F}}$, such that the minimum of (60) is given by \left\|\boldsymbol{Q}_{k}\right\|_{*}=\frac{1}{2}\big{(}\|\boldsymbol{\varGamma}_{k}\|_{\text{\sf F}}^{2}+\|\boldsymbol{G}_{k}\|_{\text{\sf F}}^{2}\big{)} with $\boldsymbol{\varGamma}_{k}$ and $\boldsymbol{G}_{k}$ given by (64).

∎

Proof of Theorem 1.

Based on Lemma 1, the $\ell_{*,1}$ mixed-norm of the source signal matrix $\boldsymbol{Q}=[\boldsymbol{Q}_{1}^{\text{\sf T}},\ldots,\boldsymbol{Q}_{K}^{\text{\sf T}}]^{\text{\sf T}}$, as defined in (25), is equivalently computed by

[TABLE]

where $r=\min(N,P)$, $\boldsymbol{\varGamma}={\rm blkdiag}(\boldsymbol{\varGamma_{1}},\ldots,\boldsymbol{\varGamma}_{K})$ is taken from the set $\mathcal{B}^{P\times r}_{K}$ of block-diagonal matrices composed of $K$ blocks of size $P\times r$ on the main diagonal, and $\boldsymbol{G}=[\boldsymbol{G}_{1}^{\text{\sf T}},\ldots,\boldsymbol{G}_{K}^{\text{\sf T}}]^{\text{\sf T}}$ is a $(Kr)\times N$ complex matrix composed of blocks $\boldsymbol{G}_{k}$, for $k=1,\ldots,K$. Inserting equation (66) into the $\ell_{*,1}$ mixed-norm minimization problem in (28) we formulate the minimization problem

[TABLE]

For a fixed matrix $\boldsymbol{\varGamma}$, the minimizer $\smash{\hat{\boldsymbol{G}}}$ of problem (66) has the closed form expression

[TABLE]

where the last equation is derived from the Woodbury matrix identity [45, p.151]. Reinserting the optimal matrix $\smash{\hat{\boldsymbol{G}}}$ into equation (67) and using basic reformulations of the objective function results in the concentrated minimization problem

[TABLE]

Upon summarizing $\boldsymbol{Y}\boldsymbol{Y}^{\text{\sf H}}/N=\hat{\boldsymbol{R}}$ and defining the positive semidefinite block-diagonal matrix

[TABLE]

we can rewrite (69) as

[TABLE]

Neglecting the factor $\lambda N/2$ in (71), we arrive at formulation (29). Using equations (61), (64) and the definition of $\smash{\hat{\boldsymbol{S}}}={\rm blkdiag}(\smash{\hat{\boldsymbol{S}}}_{1},\ldots,\smash{\hat{\boldsymbol{S}}}_{K})$ in (70) we conclude that

[TABLE]

as given in (31). Making further use of (68) and the factorization in (66) we obtain

[TABLE]

which corresponds to relation (30).

∎

Appendix B - SDP Form of the Matrix Polynomial Constraint

As discussed in Section 6, the matrix $\boldsymbol{M}(z)$ represents a matrix polynomial of degree $D$ in the variable $z$ [37]. To see the relation between matrix $\boldsymbol{\varUpsilon}_{0}$ and the matrix coefficients $\boldsymbol{K}_{i}$, let us define the $((D+1)P)\times P$ matrix

[TABLE]

and introduce the $M\times((D+1)P)$ permutation and selection matrix $\boldsymbol{J}$ such that the subarray steering matrix can be expressed as

[TABLE]

Inserting (75) in (45) yields

[TABLE]

where $\boldsymbol{F}=\boldsymbol{J}^{\text{\sf H}}\boldsymbol{\varUpsilon}_{0}\boldsymbol{J}$ is of size $((D+1)P)\times((D+1)P)$ and is composed of the $P\times P$ blocks $\boldsymbol{F}_{i,j}$, for $i,j=1,\ldots D+1$, as

[TABLE]

Equation (76) is also referred to as the Gram matrix representation of the polynomial $\boldsymbol{M}(z)$, and $\boldsymbol{F}$ is referred to as the corresponding Gram matrix [37].

We define the block trace operator for matrix $\boldsymbol{F}$ as

[TABLE]

i.e., the summation of the $P\times P$ submatrices $\boldsymbol{F}_{i,i}$, for $i=1,\ldots,D+1$, on the main diagonal of matrix $\boldsymbol{F}$. Furthermore, let us define the $(D+1)\times(D+1)$ elementary Toeplitz matrix $\boldsymbol{\varTheta}_{i}$, with ones on the $i$th diagonal and zeros elsewhere, as well as the elementary block Toeplitz matrix $\boldsymbol{\varXi}_{i}=\boldsymbol{\varTheta}_{i}\otimes\boldsymbol{I}_{P}$.

Using the block trace operator (78) and the elementary block Toeplitz matrices $\boldsymbol{\varXi}_{i}$, the matrix coefficients $\boldsymbol{K}_{i}$ in (45) can be computed from the Gram matrix $\boldsymbol{F}$ in (76) as

[TABLE]

i.e., the summation of the $P\times P$ submatrices on the $i$th block-diagonal of the Gram matrix $\boldsymbol{F}$. Note that the mapping (45) is unique, i.e., for any PCA steering matrix block $\boldsymbol{B}(z)$ and matrix $\boldsymbol{\varUpsilon}_{0}$ the coefficients $\boldsymbol{K}_{i}$, $i=1,\ldots,D$, of the matrix polynomial $\boldsymbol{M}(z)$ are unique. However, the Gram matrix $\boldsymbol{F}$ in (76) is not unique, i.e., a matrix polynomial $\boldsymbol{M}(z)$ generally admits different Gram matrix representations.

Let us define a second matrix polynomial which has constant value $\boldsymbol{I}_{P}$ as

[TABLE]

such that the corresponding Gram matrix $\boldsymbol{H}$ of size $((D+1)P)\times((D+1)P)$ fulfills

[TABLE]

By using (76), (80) and (81) we can express the constraint (46) as

[TABLE]

which is fulfilled for

[TABLE]

Applying (81) and (83) in problem (39) we can define the gridless frequency estimation problem

[TABLE]

Given a minimizer $\smash{\hat{\boldsymbol{\varUpsilon}}}_{0}$ to problem (84) the frequency estimation problem reduces to finding roots for which the constraint (46) becomes singular, as discussed in Section 6.