Super-Resolution DOA Estimation for Arbitrary Array Geometries Using a Single Noisy Snapshot

TL;DR

This paper introduces a gridless super-resolution method for direction-of-arrival estimation from a single noisy snapshot, applicable to arbitrary sensor array geometries, and demonstrates high-resolution results in simulations.

Contribution

It extends gridless super-resolution DOA estimation to arbitrary array geometries with noisy data, using a dual polynomial approach and polynomial rooting techniques.

Findings

High-resolution DOA estimation achieved in simulations.

Effective in white and colored noise scenarios.

Applicable to circular and random planar arrays.

Abstract

We address the problem of search-free DOA estimation from a single noisy snapshot for sensor arrays of arbitrary geometry, by extending a method of gridless super-resolution beamforming to arbitrary arrays with noisy measurements. The primal atomic norm minimization problem is converted to a dual problem in which the periodic dual function is represented with a trigonometric polynomial using truncated Fourier series. The number of terms required for accurate representation depends linearly on the distance of the farthest sensor from a reference. The dual problem is then expressed as a semidefinite program and solved in polynomial time. DOA estimates are obtained via polynomial rooting followed by a LASSO based approach to remove extraneous roots arising in root finding from noisy data, and then source amplitudes are recovered by least squares. Simulations using circular and random planar arrays show high resolution DOA estimation in white and colored noise scenarios.

Figures

Equations

$$\mathbf{y}=\mathcal{S}x+\mathbf{n},\text{where}{y}_m=n_m+\underset{-\pi }{\overset{\pi }{\int }}{a}_m(\theta )x(\theta )d\theta ,$$

$$a_m(\theta )=e^{-j2\pi f{\tau }_m(\theta )},$$

$$2\pi f{\tau }_m(\theta )=2\pi (|{\mathbf{p}}_m|/\lambda )\cos (\theta -\angle {\mathbf{p}}_m),$$

$$\underset{x}{\min }\Vert x{\Vert }_{\text{orgmathcal}A}^{},s.t.\Vert \mathbf{y}-\mathcal{S}x{\Vert }_2^{}\le \delta ,$$

$$\underset{\mathbf{c}\in {\mathbb{C}}^M}{\max }\Re \{{\mathbf{c}}^H\mathbf{y}\}-\delta {\Vert \mathbf{c}\Vert }_2,\text{s.t. }\Vert \mathcal{S}(\theta {)}^H\mathbf{c}{\Vert }_{\infty }\le 1,$$

$$b(\theta )=\sum _{k=-N}^N\sum _{m=1}^M({\alpha }_m[k]c_m)e^{jk\theta },$$

$${\hat{\alpha }}_m^{}[k]\simeq (1/P)\sum _{l=-N}^N{a}_m^{\ast }(l\Delta \theta )e^{-j(2\pi /P)lk},$$

$$P=15.9|\mathbf{p}|/\lambda +\mathrm{27.03.}$$

$$b(\theta )$$

$$\mathbf{h}={\left[{\hat{B}}_{-N}{\hat{B}}_{-(N-1)}\dots {\hat{B}}_N\right]}^T={\mathbf{G}}^H\mathbf{c},$$

$$\underset{\mathbf{c},\mathbf{H}}{\min }\Re \{{\mathbf{c}}^H\mathbf{y}\}-$$

$$\displaystyle\ \sum_{i=1}^{P-j}\bm{H}_{i,i+j}=\left\{\begin{array}[]{cl}1,&\mbox{ $j=0$}\\ 0&\mbox{ $j=1,\ldots,P-1$},\end{array}\right.$$

$$\underset{\mathbf{x}}{\min }\Vert \mathbf{x}{\Vert }_1^{},s.t.\Vert \mathbf{y}-{\mathbf{A}}_{\text{aug}}\mathbf{x}{\Vert }_2^{}\le ϵ,$$

$$\underset{\mathbf{x}}{\min }\frac{1}{2}\Vert \mathbf{y}-{\mathbf{A}}_{\text{aug}}\mathbf{x}{\Vert }_2^{}+\beta \Vert x{\Vert }_1^{},$$

Full text

**Index Terms— ** Super-resolution, off-grid problem, sparse DOA estimation, arbitrary array geometry, single snapshot.

1 Introduction

Direction-of-arrival (DOA) estimation can be very challenging when snapshots are limited and sources are coherent as in the case of fast moving sources and multipath arrivals. Under these conditions, high resolution DOA methods such as MVDR and MUSIC [1, 2] fail due to inaccurate estimation of spatial covariance matrix and self signal cancellation.

Sparsity based methods for DOA estimation inspired by compressed sensing (CS) [3, 4, 5, 6] can tackle coherent sources and single snapshot. However, the CS based approaches are limited by the finite discrete grid of angles used to form the basis, leading to the off-grid problem [7] when the source directions do not lie on the grid. To improve performance, greedy algorithms with a highly coherent dictionary (finer search grids) are used in [8, 9], but they are computationally demanding. The off-grid DOA approaches [10, 11, 12, 13] applicable for arbitrary arrays use a Taylor series approximation of array steering vectors on fixed grids, or iterative methods with dynamic grids to tackle the grid mismatch. However, their performance and accuracy depends on the grid density or they require noncovex optimization. Recent gridless super-resolution approaches using convex optimization [14, 15, 16, 17] eliminate the off-grid problem by forming the basis in the continuous angle domain and provide high accuracy, but they are not applicable to arbitrary geometries.

In this paper, we develop a search-free DOA estimation method for arrays of arbitrary geometry under the challenging conditions of coherent sources and a single noisy snapshot. This extends our earlier work [18] on super-resolution DOA estimation for arbitrary geometry, to noisy measurements. The DOA estimation problem for arbitrary geometry is solved as a dual maximization problem. By exploiting the periodicity and band-limited nature of the dual function, we can represent it with a finite trigonometric polynomial using Fourier series (FS). The proposed approach is motivated by [19, 20, 21, 22], where root-MUSIC is extended to arbitrary arrays. The modified dual problem can then be expressed as a finite semidefinite program (SDP), and solved. Finally, the search-free DOA estimates are obtained through polynomial rooting of a nonnegative polynomial formed from the dual polynomial. To remove the extraneous roots arising in the noisy case, we use a LASSO-like approach related to [23, 24].

2 Data Model

Consider an $M$-element array of arbitrary geometry, which receives signals from $L$ narrowband far-field sources with complex amplitude $s_{l}$ and azimuth DOA $\theta_{l}$, $l=1,\dots,L$. We define the sparse source function $x(\theta)$ in the continuous angle domain $\theta\in(-\pi,\pi]$ with impulses as $x(\theta)=\sum_{l=1}^{L}s_{l}\delta(\theta-\theta_{l})$. Then the $M\times 1$ observed array snapshot vector $\bm{y}$ is

[TABLE]

$m=1,\dots,M$ and $\bm{n}\in\mathbb{C}^{M}$ is the received additive noise across the array. The linear measurement operator $\mathcal{S}$ represents the array manifold over $\theta$, whose $m$-th component $a_{m}(\theta)$ is the response of the $m$-th sensor for a source at direction $\theta$.

[TABLE]

where $\tau_{m}(\theta)$ is the propagation delay with respect to a reference.111We prefer to study $a_{m}(\theta)$ as a function of $\theta$. On the other hand, at a specific angle $\theta_{1}$, $[a_{m}(\theta_{1})]\in\mathbb{C}^{M}$ is the steering vector for direction $\theta_{1}$. For narrowband sources of frequency $f$ and propagation speed $v$, the wavelength is $\lambda=v/f$. Using $\tau_{m}(\theta)={\langle\bm{p}_{m},\bm{u}_{\theta}\rangle}/{v}$, we simplify the exponent in (2) as

[TABLE]

where $\bm{p}_{m}$ is the position vector of the $m$-th sensor with respect to a reference, and $\bm{u}_{\theta}$ is a unit vector in direction $\theta$.

3 Proposed Method

Assuming the sources are sparse in angle, $x(\theta)$ could be recovered from noisy measurements $\bm{y}=\mathcal{S}x+\bm{n}$ via [16]

[TABLE]

where $\delta$ satisfies the condition that $\lVert\bm{n}\rVert^{\phantom{.}}_{2}\leq\delta$. $\lVert.\rVert^{\phantom{.}}_{\orgmathcal{A}}$ denotes the atomic norm [25] which is a continuous analogue of the $l_{1}$ norm, i.e., $\lVert x\rVert^{\phantom{.}}_{\orgmathcal{A}}=\sum_{l=1}^{L}\lvert s_{l}\rvert.$ Here $\mathcal{S}$ does not represent Fourier measurements, unlike [14, 15, 16, 17]. The primal problem (4) is infinite dimensional and difficult to solve. Therefore, we work with the corresponding dual maximization problem

[TABLE]

where $\bm{c}$ is the dual variable (see details in [17, 16]). The dual function defined by $\mathcal{S}(\theta)^{H}\bm{c}$ has unit magnitude in the direction of actual sources, irrespective of geometry. For a uniform linear array (ULA), $\mathcal{S}(\theta)^{H}\bm{c}$ is, in fact, an $(M-1)^{\text{th}}$ degree polynomial in $e^{j\theta}$, and (5) is then solved using an SDP [15, 16, 17]. The polynomial structure arises from the fact that sensor delays, $\tau_{m}(\theta)$ in (2), for a ULA are integer multiples of a constant. For arbitrary arrays, $\mathcal{S}(\theta)^{H}\bm{c}$ cannot be directly expressed as a polynomial, but we overcome this difficulty with a Fourier domain (FD) representation of the dual function that provides a polynomial form for the SDP.

3.1 Fourier Domain Representation of the Dual Function

We review the Fourier series representation of the dual function [18] here for completeness. The function $b(\theta)=\mathcal{S}(\theta)^{H}\bm{c}$ is periodic in $\theta$ with period $2\pi$ as it is a linear combination of smooth (band-limited) periodic functions, $a^{*}_{m}(\theta)$, $m=1,\dots,M$. Thus, $b(\theta)$ has a Fourier series (FS) which can be truncated if its Fourier coefficients ${B}_{k}\!\approx\!0$ for $|k|>N$. Each $a^{*}_{m}(\theta)$, being periodic, has a FS with coefficients $\alpha_{m}[k]$, related to $B_{k}$ via $B_{k}=\sum_{m}\alpha_{m}[k]c_{m}$. So we have

[TABLE]

which is a finite degree polynomial in $z=e^{j\theta}$. As a result, we can determine $N$ for FS truncation by examining the FS coefficients at each sensor, $\alpha^{\phantom{.}}_{m}[k]$, which depend solely on the array geometry and not on the measured signals.

Assuming a sufficiently large number of DFT points $(P=2N+1)$ for dense sampling in $\theta$, the FS coefficients $\alpha^{\phantom{.}}_{m}[k]$ can be estimated from samples of $a^{*}_{m}(\theta)$ using the DFT [19, 26] as,

[TABLE]

where $\Delta\theta=2\pi/P$, and $k=-N,\ldots,0,1,\ldots,N$. Note that circular indexing of the DFT is exploited in (7).

Next, we conduct a numerical study of the FS for the continuous function $a^{*}_{m}(\theta)$ defined in (2, 3) to determine the value of $P$ needed for various array geometries. The FS coefficients of $a^{*}_{m}(\theta)$ can be approximated numerically by a very long DFT to get $\hat{\alpha}[k]$. From (2, 3), the magnitude $|\hat{\alpha}_{m}[k]|$ depends only on $|\bm{p}|/\lambda$, the normalized distance of the sensor from the origin. This is because $(\theta-\angle\bm{p})$ is a shift in the argument of $a_{m}(\theta)$ which changes only the phase of its FS coefficients. We use a long DFT to obtain FS coefficients for many different values of $|\bm{p}|/\lambda$, and display the magnitude $|\hat{\alpha}[k]|^{2}$ as an image in Fig. 1, which confirms that $\hat{\alpha}[k]$ is bandlimited.

We observe that as $|\bm{p}|/\lambda$ increases, the bandwidth of the FS grows and hence the distance of the farthest sensor from origin in an array controls the minimum $P$ needed to get an accurate DFT representation. The index $N$ where $|\hat{\alpha}[k]|\approx 0$ for $|k|>N$ depends on choosing a threshold $\gamma$ for the squared magnitude of the FS. Figure 1 shows the case for $\gamma=-160$ dB below the maximum. A linear approximation derived for $|\bm{p}|/\lambda\geq 2$ gives an excellent estimate for $P=2N+1$. This minimum value of $P$ is important for reducing the computational complexity of the SDP. For $\gamma=-160$ dB, the linear estimate is,

[TABLE]

An example verifying the match between predicted and observed $P$ needed to ensure success is presented in [18].

Using the DFT representation in (7), the dual function $b(\theta)$ can be related to a dual polynomial $\hat{b}(z)$ as

[TABLE]

Combining (6) and (9), we recognize that the coefficients $\hat{B}_{k}$ can be written in matrix-vector form with $\bm{h}\in\mathbb{C}^{P}$ being

[TABLE]

where $\bm{G}^{H}=\begin{bmatrix}\hat{\alpha}_{m}[k]\end{bmatrix}_{P\times M}$ is a matrix whose $m$-th column contains the FS coefficients of $a^{*}_{m}(\theta)$, and $\bm{c}$ is the dual vector.

3.2 Semidefinite Programming

We convert the infinite number of constraints in the dual problem (5) into finite-dimensional matrix constraints as in [17, 16], by using the uniform boundedness of the function $\mathcal{S}(\theta)^{H}\bm{c}$ in (5) and hence that of its FD representation given by the dual polynomial $\hat{b}(z)$, to obtain the following SDP.

[TABLE]

where $\operatorname{\Re}\{\cdot\}$ denotes the real part. $\bm{H}$ is a positive semidefinite matrix satisfying the constraints in ($11$). The SDP ($11$) has $n=P^{2}/2+M$ optimization variables and is solvable in polynomial time by interior-point methods [27]. The observed time complexity was found to be much less than the worst case $\mathcal{O}(n^{3})$. The dual polynomial $\hat{b}(z)$ is the desired output after the SDP, so its coefficient vector is constructed from the optimal $\bm{c}_{*}$ via $\bm{h}_{*}=\bm{G}^{H}\bm{c}_{*}$.

3.3 Recovery of Source DOAs and Amplitudes

For sufficiently large $P$, the approximation of the dual function $\mathcal{S}(\theta)^{H}\bm{c}_{*}$ by the dual polynomial $\hat{b}(e^{j\theta})$ is highly accurate. Based on the constraint (5), $|\hat{b}(e^{j\theta})|$ would be equal to one for true DOAs, and less than one elsewhere [15]. To locate the angles $\theta$ where the magnitude of the dual polynomial is one, we form a nonnegative polynomial $p(z)=1-|\hat{b}(z)|^{2}$ from the dual polynomial coefficients $\bm{h}_{*}$. The coefficients of $|\hat{b}(z)|^{2}$, denoted by $r_{k}$, are the autocorrelation coefficients of $\bm{h}_{*}$, i.e., $r_{k}=\sum_{j}h_{j}h_{j-k}^{*}$. The angles of the zeros of $p(z)$ on the unit circle include the DOAs of the sources.

Due to numerical issues of polynomial rooting at low SNRs, the SDP might provide extraneous unit-circle zeros that do not correspond to true sources. Therefore, the DOAs are finally recovered by a sparsity-promoting $\ell_{1}$ problem.

[TABLE]

where $\bm{A}_{\text{aug}}$ is a dictionary of steering vectors that has steering vectors for a discrete set of angles as its columns. This discrete set includes the angles of the unit-circle roots from the SDP, as well as additional angles drawn from a uniform distribution in $(-180^{\circ},180^{\circ}]$. Then ($12$) is written as the following LASSO-like problem and solved using convex optimization.

[TABLE]

The support of the solution $\bm{x}_{*}$ yields the DOAs of interest.

Once we estimate the DOAs, the amplitudes of the sources are recovered by least squares $\hat{\bm{s}}=\bm{A}(\hat{\bm{\theta}})^{\dagger}\bm{y}$ where † denotes the pseudo-inverse. The columns of the matrix $\bm{A}(\hat{\bm{\theta}})$ are the steering vectors for the estimated DOAs $\hat{\bm{\theta}}$.

To summarize, the steps in the proposed method are:

1. For the geometry, compute $\bm{G}^{H}=\begin{bmatrix}\hat{\alpha}_{m}[k]\end{bmatrix}$ via (7).

2. Estimate the noise level, and set $\delta$.

3. Using $\bm{G}$, $\bm{y}$ and $\delta$, solve the SDP in ($11$) to find the optimal $\bm{c}_{*}$.

4. Get the optimal dual polynomial coefficients via $\bm{h}_{*}=\bm{G}^{H}\bm{c}_{*}$.

5. Estimate DOAs $\hat{\bm{\theta}}$ by finding angles of unit circle roots of $p(z)$.

6. Eliminate extraneous zeros via the $\ell_{1}$ sparsity optimization ($13$).

7. Recover the source amplitudes $\hat{\bm{s}}$ by least squares.

4 Simulations

Results for the uniform circular array (UCA) and random planar array (RPA) geometries are presented in Sections 4.1 and 4.2. Performance is compared with the conventional delay-sum beamformer (CBF). All simulations consider a single snapshot and multiple coherent sources [28], which are complex sinusoids of the same frequency with constant phase difference. We implemented the SDP ($11$) using CVX [29]. For DOA estimation, we use only those roots of $p(z)$ that lie within a distance of $0.02$ from the unit circle.

4.1 Simulations for Uniform Circular Array (UCA)

Two examples using a 40-element UCA are presented here. The array radius is $r=2\lambda$, and the uniform sensor separation is $d=(\pi/10)\lambda$. With the reference point at the center of the array, $|\bm{p}_{m}|=r$ for all sensors.

In the first example in Fig. 2(b), we study the angular resolution of the proposed method by considering two equal magnitude sources of SNR $20$ dB separated by $10^{\circ}$. Since the noise in practice is often colored, we simulate noise with $1/f$ spectral decay along frequency for this example. As seen in Fig. 2(b), the CBF is not able to resolve the two closely located sources, whereas estimates from the unit-circle zeros in Fig. 2(a) are very accurate. This reinforces that the proposed approach offers higher resolution than existing methods for single snapshot DOA estimation. The approach only assumes additive noise and this example also verifies its applicability to colored noise scenarios. Additive white Gaussian noise is used in the rest of the examples.

In Fig. 3, we consider five equal magnitude sources at $5$ dB SNR. Due to the lower SNR, the set of unit-circle zeros of $p(z)$ in Fig. 3(a) includes three extraneous zeros in addition to the five zeros that correspond to the true DOAs. Using the $\ell_{1}$ norm based DOA recovery in ($13$), we eliminate those unwanted roots as shown in Fig. 3(b). The nonzero elements in the $\ell_{1}$ recovery result are the final estimated DOAs. The amplitudes from the $\ell_{1}$ recovery are expected to be inaccurate due to shrinkage operation. Once we estimate the DOAs, the amplitudes can be recovered via least squares. The CBF is unable to resolve two among the five sources, and shows high side lobes as well inaccurate source amplitude estimates. On the other hand, the proposed approach in combination with the $\ell_{1}$ recovery accurately estimates the DOAs of all three sources. Both UCA examples validate the ability of the proposed method to estimate DOAs accurately for an arbitrary 2-D array.

4.2 Simulation for Random Planar Array (RPA)

In Fig. 4(a), we consider an RPA with $30$ sensors. The minimum sensor spacing is $d=\lambda/4$, and the distance of the farthest sensor from origin is around $2\lambda$. The proposed method resolves both sources as shown in Fig. 4(b), whereas, the CBF results in a single peak at $65^{\circ}$ (CBF result not shown).

4.3 Performance Evaluation Vs. SNRs

We now evaluate the performance of the approach for various SNRs, and the sensitivity of the method to the value of noise norm upper-bound $\delta$. The RMSE in DOA estimation for different SNRs and $\delta$ is provided in Figs. 5(a) and 5(b) for two sources of separation $10^{\circ}$ and $30^{\circ}$, respectively. $e_{n}$ is the expected value of noise norm $e_{n}=\mathbb{E}[\lVert\bm{n}\rVert^{\phantom{.}}_{2}]$. For i.i.d noise $\mathcal{N}(0,\sigma_{n})$ across the sensors, $e_{n}=\sigma_{n}\sqrt{M}$. The simulation considers 50 Monte Carlo trials with random source DOAs for each SNR. The performance depends on the minimum separation between sources. For larger separations, a smaller DOA error was observed. Note that at $10^{\circ}$ separation, the CBF is unable to resolve the sources at all SNRs (see Fig. 2). Regarding the choice of $\delta$, an underestimation of $\delta$ was observed to cause many extraneous unit circle roots, but the $\ell_{1}$ recovery could remove those additional roots. The overestimation of $\delta$, on the other hand, resulted in fewer roots on the unit circle, but they were slightly less accurate. In general, with $\ell_{1}$ recovery processing, an underestimated $\delta$ provided better results than the overestimated one. As the SNR improves, the performance becomes less dependent on the choice of $\delta.$ For SNR above $30$ dB in Fig. 5(b), the estimates are nearly perfect. The parameters involved in the approach are: $\delta$, $\beta$ and two thresholds, one for unit-circle zero detection, another for discarding low magnitude coefficients in the $\ell_{1}$ recovery.

5 Discussion

We have presented a search-free super-resolution DOA estimation and beamforming method for arbitrary geometry arrays, which is applicable for a single noisy snapshot, and correlated or uncorrelated sources. Further SNR improvement should be possible using multiple snapshots. The upper bound of the noise norm $\delta$ in (4) needs to be estimated in practice. However, unlike traditional high resolution approaches, the proposed method does not require knowledge of the number of sources. We made comparisons with the CBF, but not with traditional high resolution DOA approaches such as MUSIC and MVDR as they fail in the single snapshot case and coherent signal conditions, though they are applicable for arbitrary arrays. Moreover, existing sparsity based gridless super-resolution approaches are applicable only for ULAs. Simulation results prove that the new method can perform high resolution search-free DOA estimation for arbitrary geometries, using a single noisy snapshot.