# Robust Broadband Beamforming
using Bilinear Programming

Nakul Singh, Coleman DeLude, Mark A. Davenport, and Justin Romberg Email:
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>This work was supported in part by COGNISENSE, one of
seven centers in JUMP 2.0, a Semiconductor Research Corporation (SRC) program
sponsored by DARPA, and a grant from Lockheed Martin. School of Electrical and
Computer Engineering
Georgia Institute of Technology

###### Abstract

We introduce a new method for robust beamforming, where the goal is to
estimate a signal from array samples when there is uncertainty in the angle of
arrival. Our method offers state-of-the-art performance on narrowband signals
and is naturally applied to broadband signals. Our beamformer operates by
treating the forward model for the array samples as unknown. We show that the
“true” forward model lies in the linear span of a small number of fixed linear
systems. As a result, we can estimate the forward operator and the signal
simultaneously by solving a bilinear inverse problem using least squares. Our
numerical experiments show that if the angle of arrival is known to only be
within an interval of reasonable size, there is very little loss in estimation
performance compared to the case where the angle is known exactly.

##  I Introduction

Beamforming, the process of forming an estimate of a signal coming from a
fixed direction $\theta$ from the outputs of a multi-sensor array, relies
critically on knowing this direction of arrival $\theta$. When the estimator
is “misaligned,” i.e., using a $\theta$ that is different from the true angle
of arrival $\theta_{\mathrm{true}}$, the quality of the estimate can degrade
quickly as $\theta$ moves away from $\theta_{\mathrm{true}}$. This problem is
especially pronounced in modern arrays with a large number of elements; large
apertures allow us to create tapered beams, but they are also more fragile to
misalignment.

While there is a significant amount of literature on making beamformers
“robust” to unknown $\theta$ (a very nice overview of current methods can be
found in [1]), the most effective approaches known to date 1) apply only to
narrowband signals, 2) require solving costly optimization programs to compute
weight vectors. In this paper, we introduce a relatively simple approach based
on bilinear programming that naturally includes the broadband case while not
being significantly more computationally intensive than broadband beamforming
with a known direction of arrival $\theta$.

Recent work in [2] has shown how broadband beamforming can effectively be
treated as a linear inverse problem from measurements
$\boldsymbol{y}\approx\boldsymbol{A}(\theta)\boldsymbol{\alpha}$ that can be
solved with least squares (this is discussed in more detail in Section II-A).
The forward model $\boldsymbol{A}(\theta)$ of course depends on the angle of
arrival (AOA) $\theta$. In this paper, we show how we can make this broadband
beamforming method robust to uncertainty in $\theta$. We show how each
$A(\theta)$ in a region of angles $\theta\in\Theta$ can be written as a
superposition of a relatively small number of fixed matrices that can be pre-
computed. The result is a bilinear least squares problem. This type of problem
is nonconvex, but is prevalent enough in applied mathematics that many
different heuristics exist for efficient solvers. We demonstrate the
effectiveness of one such heuristic, alternating minimization, in the
numerical experiments in Section IV. These results show that the performance
penalty for not knowing the angle of arrival exactly is small, and the
computation required to achieve this performance is not significantly more
than a standard least squares solver.

In addition to broadband beamforming, the method is also applicable to
narrowband beamforming, where there is a well-established literature (some of
which is briefly reviewed in Section IV). We show in Section IV that our
bilinear robust beamformer also offers state-of-the-art performance for
narrowband signals.

##  II Formulation

In this section, we develop the system model for recovering the signal of
interest when the arrival angle is unknown. We first present the system model
developed in [2] for the case when true AOA is known and then propose a
bilinear system model for the cases where we only have a rough estimate of the
AOA.

###  II-A Known AOA

We consider the beamforming problem where a plane wave signal strikes an $M$
element array at an angle111This can either represent a standard angle or
spherical angle depending on the array configuration. The mathematical
formulation is the same in either case. $\theta$, and $N$ snapshots are
recorded, providing a total of $MN$ observations. We model the signal of
interest $s(t)$ as being bandlimited to $\Omega$, e.g., the spectral support
of $s(t)$ is non-zero for $f\in[-\Omega,\Omega]$. Over the temporal extent of
the $N$ snapshots, we can write $s(t)$ as the superposition of $K$ basis
functions $\psi_{1}(t),\ldots,\psi_{K}(t)$

$s(t)=\sum_{k=1}^{K}\alpha_{k}\psi_{k}(t).$

A natural choice for the $\\{\psi_{k}\\}$ are the Slepian basis functions,
which provide a compact representation for bandlimited signals supported over
a finite interval of time (see [2] for further discussion). Prior to
propagation, it is assumed the signal is modulated to a carrier frequency
$f_{c}$ to form $s_{\text{mod}}(t)=e^{j2\pi f_{c}t}s(t)$.

The received $s_{\text{mod}}(t)$ are then demodulated at each array element,
reducing the carrier to only a phase offset. We model the array measurements
as samples of this demodulated signal. The $n^{\text{th}}$ sample at the
$m^{\text{th}}$ array element is written as

$y_{m}[n]=e^{-j2\pi
f_{c}\tau_{m}(\theta)}\sum_{k=1}^{K}\alpha_{k}\psi_{k}(t_{n}-\tau_{m}(\theta))+\textnormal{noise}.\\\
$ (1)

The $\tau_{m}(\theta)$ in (1) is the delay (relative to array center) of
element $m$ for a plane wave arriving from angle $\theta$; $t_{n}$ denotes the
sample time for the $n^{\text{th}}$ snapshot. Therefore, a forward model for a
single snapshot is given by

$\begin{gathered}\boldsymbol{y}[n]=\boldsymbol{A}_{n}(\theta)\boldsymbol{\alpha}+\textnormal{noise},\\\
\end{gathered}$ (2)

where $\boldsymbol{A}_{n}(\theta)\in\mathbb{C}^{M\times K}$ is the matrix with
entries given by $A_{n}(\theta)[m,k]=e^{-j2\pi
f_{c}\tau_{m}(\theta)}\psi_{k}(t_{n}-\tau_{m}(\theta))$. Stacking the model in
(2) for the $N$ snapshots, we arrive at the system model for the array
beamforming when the AOA is known:

$\begin{gathered}\underbrace{\begin{bmatrix}\boldsymbol{y}[1]\\\ \vdots\\\
\boldsymbol{y}[N]\end{bmatrix}}_{\boldsymbol{y}}=\underbrace{\begin{bmatrix}\boldsymbol{A}_{1}(\theta)\\\
\vdots\\\
\boldsymbol{A}_{N}(\theta)\end{bmatrix}}_{\boldsymbol{A}(\theta)}\boldsymbol{\alpha}+\textnormal{noise}.\end{gathered}$
(3)

To recover the Slepian weights ($\boldsymbol{\alpha}$) from (3), we use least
squares which have a closed-form solution:

$\begin{gathered}\hat{\boldsymbol{\alpha}}=\operatorname*{\text{arg~{}min}}_{\boldsymbol{\alpha}}||\boldsymbol{y}-\mathbf{A(\theta)}\boldsymbol{\alpha}||_{2}^{2}\\\
\implies\hat{\boldsymbol{\alpha}}=(\mathbf{A(\theta)}^{H}\mathbf{A(\theta)})^{-1}\mathbf{A(\theta)}^{H}\boldsymbol{y}.\end{gathered}$
(4)

The Slepian beamforming in (4) works well only when the angle $\theta$ is
known; we show how the performance suffers when the estimated angle deviates
from the true angle in our numerical experiments in Section IV below (in
particular, see Figures 2 and 4).

It is noteworthy to mention that the dimension $K$ of the low-dimensional
subspace spanned by the Slepian basis depends on the true angle of arrival as
$K=\lceil 2\Omega T(\theta)\rceil$ where $T(\theta)$ is the array aperture.
However, for reasonable deviations in the angle of arrival, the variation in
$K$ is insignificant; therefore, in subsequent sections, we will assume the
dimension $K$ to be virtually known.

###  II-B Bilinear formulation for unknown AOA

When $\theta$ is unknown, there is uncertainty in the forward model (3) for
the array samples. We approach this problem by finding a linear embedding for
the rows of the $\boldsymbol{A}_{n}(\theta)$ for an interval $\Theta$ of
possible AOAs around a fixed angle $\theta_{0}$:
$\Theta=[\theta_{0}-\Delta,\theta_{0}+\Delta]$.

To illustrate how this works, let
$\boldsymbol{a}_{\ell}(\theta)\in\mathbb{C}^{K}$ be the $\ell^{\text{th}}$ row
of $\boldsymbol{A}(\theta)$. Note that each value of $\ell=1,\dots,MN$ has a
one-to-one correspondence to a pair $(m,n)$. As $\theta$ changes,
$\boldsymbol{a}_{\ell}(\theta)$ also changes, but does so smoothly. The key
realization is that the collection of vectors
$\\{\boldsymbol{a}_{\ell}(\theta),~{}\theta\in\Theta\\}$ can all be very
closely approximated by vectors from a low dimensional subspace. That is,
there exists a $K\times P$ matrix $\boldsymbol{U}_{\ell}$, for a very small
$P$, such that for every fixed $\theta\in\Theta$ we have
$\boldsymbol{a}_{\ell}(\theta)\approx\boldsymbol{U}_{\ell}\boldsymbol{c}$ for
a short vector $\boldsymbol{c}\in\mathbb{C}^{P}$. As such, we can write
$y_{\ell}=\boldsymbol{c}^{\mathrm{H}}\boldsymbol{U}_{\ell}^{\mathrm{H}}\boldsymbol{\alpha}$
where both $\boldsymbol{c}$ and $\boldsymbol{\alpha}$ are unknown —
$\boldsymbol{\alpha}$ depends on the signal $s(t)$ and $\boldsymbol{c}$
depends on its AOA.

We formalize this carefully as follows. Define the functions

$w_{\ell,k}(\theta)=e^{-j2\pi
f_{c}\tau_{m}(\theta)}\psi_{k}(t_{n}-\tau_{m}(\theta)),~{}\theta\in\Theta,$

for an interval $\Theta$ as above, and collect these into

$\boldsymbol{a}_{\ell}(\theta)=\begin{bmatrix}w_{\ell,1}(\theta)\\\ \vdots\\\
w_{\ell,K}(\theta)\end{bmatrix}.$

Our problem is that we observe

$y_{\ell}=\boldsymbol{a}_{\ell}(\theta^{*})^{\mathrm{H}}\boldsymbol{\alpha}+\mathrm{noise},\quad\ell=1,\ldots,MN,$
(5)

for some unknown $\theta^{*}\in\Theta$. We will show how this can be turned
into a well-posed bilinear problem by using some basic facts about reproducing
kernel Hilbert spaces.

Let $\mathcal{H}$ be the Hilbert space of functions on $\theta\in\Theta$
generated from the functions $w_{\ell,k}(\theta)$

$\displaystyle\mathcal{H}=\operatorname{Span}(\\{w_{\ell,k}(\theta),~{}\ell=1,\dots,MN;k=1,\ldots
K\\}).$ (6)

This is a finite-dimensional Hilbert space of continuous functions, and thus
it will have a reproducing kernel for any inner product
$\langle\cdot,\cdot\rangle$ that we endow it with. $\mathcal{H}$ is
essentially a Slepian space, as it is generated from modulated versions of
small segments of the (approximately) bandlimited functions
$\psi_{1}(t),\ldots,\psi_{K}(t)$.

Let $\boldsymbol{\mathcal{W}}_{\ell}:\mathcal{H}\rightarrow\mathbb{C}^{K}$ be
the linear operator defined by

$\boldsymbol{\mathcal{W}}_{\ell}x=\begin{bmatrix}\left\langle
x,w_{\ell,1}\right\rangle\\\ \vdots\\\ \left\langle
x,w_{\ell,K}\right\rangle\end{bmatrix}$ (7)

and $\boldsymbol{\mathcal{W}}:\mathcal{H}\rightarrow\mathbb{C}^{MNK}$ be the
concatenation of these operators across the array elements and snapshots:

$\boldsymbol{\mathcal{W}}x=\begin{bmatrix}\boldsymbol{\mathcal{W}}_{1}x\\\
\vdots\\\ \boldsymbol{\mathcal{W}}_{MN}x\end{bmatrix}.$

This operator is effectively low rank due to the exponential decay of
eigenvalues of the integral operator of which the Slepian basis are the
eigenfunctions; for a more detailed analysis, refer [3]. By the spectral
theorem, we know that there are (orthonormal) functions
$v_{1},v_{2},\ldots\in\mathcal{H}$ such that

$\displaystyle(\boldsymbol{\mathcal{W}}^{H}\boldsymbol{\mathcal{W}})(\theta,\theta^{\prime})$
$\displaystyle=\sum_{\ell,k}w_{\ell,k}(\theta)^{H}w_{\ell,k}(\theta^{\prime})$
$\displaystyle=\sum_{p=1}^{MNK}\sigma_{p}^{2}v_{p}(\theta)^{H}v_{p}(\theta^{\prime}).$

As noted above, the $w_{\ell,k}(\theta)$ are all just segments of bandlimited
signals over small intervals, we can truncate the sum above at some modest
$P$:

$(\boldsymbol{\mathcal{W}}^{H}\boldsymbol{\mathcal{W}})(\theta,\theta^{\prime})\approx\sum_{p=1}^{P}\sigma_{p}^{2}v_{p}(\theta)^{H}v_{p}(\theta^{\prime}).$

Likewise, if we set

$\displaystyle\boldsymbol{w}(\theta)=\begin{bmatrix}\boldsymbol{a}_{1}(\theta)\\\
\vdots\\\ \boldsymbol{a}_{MN}(\theta)\end{bmatrix},$

we know that the eigenvalue decomposition for the $MNK\times MNK$ matrix can
be similarly truncated222A more detailed analysis of the operator is presented
in appendix A.,

$\displaystyle\boldsymbol{\mathcal{W}}\boldsymbol{\mathcal{W}}^{H}$
$\displaystyle=\int_{\theta\in\Theta}\boldsymbol{w}(\theta)\boldsymbol{w}(\theta)^{\mathrm{H}}~{}d\theta$
$\displaystyle=\sum_{p=1}^{MNK}\sigma_{p}^{2}\boldsymbol{u}_{p}\boldsymbol{u}_{p}^{\mathrm{H}}\approx\sum_{p=1}^{P}\sigma_{p}^{2}\boldsymbol{u}_{p}\boldsymbol{u}_{p}^{\mathrm{H}}.$

We can thus write for any $x\in\mathcal{H}$

$\boldsymbol{\mathcal{W}}x\approx\sum_{p=1}^{P}\sigma_{p}\boldsymbol{u}_{p}\langle
x,v_{p}\rangle,$

and

$\boldsymbol{\mathcal{W}}_{\ell}x\approx\sum_{p=1}^{P}\sigma_{p}\boldsymbol{u}_{p,\ell}\langle
x,v_{p}\rangle,$

where $\boldsymbol{u}_{p,\ell}\in\mathbb{C}^{K}$ is just the
“$\ell^{\text{th}}$ part” in $\boldsymbol{u}_{p}\in\mathbb{C}^{MNK}$:

$\boldsymbol{u}_{p}=\begin{bmatrix}\boldsymbol{u}_{p,1}\\\ \vdots\\\
\boldsymbol{u}_{p,MN}\end{bmatrix}.$

Note that the $\boldsymbol{u}_{p,\ell}$ _are not_ the left singular vectors of
the operator $\boldsymbol{\mathcal{W}}_{\ell}$; they are parts of the singular
vectors of the joint operator $\boldsymbol{\mathcal{W}}$.

Returning to our forward model (5), let $\beta^{*}\in\mathcal{H}$ be the slice
of the reproducing kernel that takes samples at location $\theta^{*}$. Then
$\boldsymbol{a}_{\ell}(\theta^{*})=\boldsymbol{\mathcal{W}}_{\ell}\beta^{*}$
and so (ignoring the noise)

$\displaystyle
y_{\ell}=\boldsymbol{a}_{\ell}(\theta^{*})^{\mathrm{H}}\boldsymbol{\alpha}$
$\displaystyle=\left(\boldsymbol{\mathcal{W}}_{\ell}\beta^{*}\right)^{\mathrm{H}}\boldsymbol{\alpha}$
$\displaystyle\approx\sum_{p=1}^{P}\sigma_{p}\langle\beta^{*},v_{p}\rangle\boldsymbol{u}_{p,\ell}^{\mathrm{H}}\boldsymbol{\alpha}$
$\displaystyle=\sum_{p=1}^{P}c_{p}\boldsymbol{u}_{p,\ell}^{\mathrm{H}}\boldsymbol{\alpha}$
$\displaystyle=\left(\boldsymbol{U}_{\ell}\boldsymbol{c}\right)^{\mathrm{H}}\boldsymbol{\alpha},$

where $c_{p}=\sigma_{p}\langle\beta^{*},v_{p}\rangle$ and the $K\times P$
matrix $\boldsymbol{U}_{\ell}$ is

$\displaystyle\boldsymbol{U}_{\ell}=\begin{bmatrix}\boldsymbol{u}_{1,\ell}&\ldots&\boldsymbol{u}_{P,\ell}\\\
\end{bmatrix}.$

_Note that $\boldsymbol{c}$ is independent of $l$_. Thus our bilinear problem
is to recover $\boldsymbol{c}\in\mathbb{C}^{P}$ and $\alpha\in\mathbb{C}^{K}$
from measurements $y_{\ell}$ given the forward model

$y_{\ell}=\boldsymbol{c}^{\mathrm{H}}(\boldsymbol{U}_{\ell})^{\mathrm{H}}\boldsymbol{\alpha}~{}+~{}\mathrm{noise},$
(8)

where the $\boldsymbol{U}_{\ell}$ are known $K\times P$ matrices that can be
computed beforehand. Ideally, $\boldsymbol{U}_{\ell}$ are calculated after
taking the singular value decomposition of the matrix
$\boldsymbol{\mathcal{W}}\boldsymbol{\mathcal{W}}^{H}$. However, for a large
array with many snapshots, performing the SVD of this matrix can be
computationally expensive. To circumvent this issue, we approximate the
operator $\boldsymbol{\mathcal{W}}$ as a matrix that allows faster computation
of the SVD. Instead of defining the functions $w_{\ell,k}(\theta)$ over a
continuous interval $\Theta$, we discretize the interval by uniformly sampling
$Q$ points $\\{\theta_{1},\ldots,\theta_{Q}\\}\in\Theta$. Corresponding to
these sampled points, we define
$\boldsymbol{W}_{\ell}:\mathbb{C}^{Q}\rightarrow\mathbb{C}^{K}$ as

$\displaystyle\boldsymbol{W}_{\ell}=\begin{bmatrix}w_{\ell,1}(\theta_{1})&&\ldots&&w_{\ell,1}(\theta_{Q})\\\
\vdots&&\ddots&&\vdots\\\
w_{\ell,K}(\theta_{1})&&\ldots&&w_{\ell,K}(\theta_{Q})\\\ \end{bmatrix},$

and $\boldsymbol{W}:\mathbb{C}^{Q}\rightarrow\mathbb{C}^{MNK}$ becomes the
concatenation of these approximate operators. We can now evaluate the SVD of
the $Q\times Q$ matrix $\boldsymbol{W}^{H}\boldsymbol{W}$ much faster since
$Q\ll MNK$ for large arrays. Using the sampled orthonormal functions
$\boldsymbol{v}_{p}$ we get
$\boldsymbol{u}_{p}=\boldsymbol{W}\boldsymbol{v}_{p}/\sigma_{p}$ which are
then used to compute $\boldsymbol{U}_{\ell}$.

The next section describes the alternating minimization method used to solve
the bilinear least squares resulting from the model devised in (8).

##  III Alternating minimization for bilinear least squares

To perform alternating minimization, we first stack the bilinear formulation
derived in (8) across array elements and snapshots to arrive at the following
forward model,

$\begin{gathered}\underbrace{\begin{bmatrix}y_{1}\\\ \vdots\\\
y_{MN}\end{bmatrix}}_{\boldsymbol{y}}=\boldsymbol{C}\underbrace{\begin{bmatrix}\boldsymbol{U}_{1}^{H}\\\
\vdots\\\
\boldsymbol{U}^{H}_{MN}\end{bmatrix}}_{\bar{\boldsymbol{U}}}\boldsymbol{\alpha}+\text{noise},\\\
\boldsymbol{C}\in\mathbb{C}^{MN\times PMN}\text{ s.t.
}C[i,j]=\begin{cases}[c^{H}]_{j}\text{ }(i-1)P\leq j\leq iP\\\ 0\text{
otherwise}.\end{cases}\end{gathered}$

A similar model with block structure over $\boldsymbol{\alpha}$ instead of
$\boldsymbol{c}$ is given by,

$\begin{gathered}\underbrace{\begin{bmatrix}(y_{1})^{*}\\\ \vdots\\\
(y_{MN})^{*}\end{bmatrix}}_{\bar{\boldsymbol{y}}}=\boldsymbol{V}\underbrace{\begin{bmatrix}\boldsymbol{U}_{1}\\\
\vdots\\\
\boldsymbol{U}_{MN}\end{bmatrix}}_{\boldsymbol{U}}\boldsymbol{c}+\text{noise},\\\
\boldsymbol{V}\in\mathbb{C}^{MN\times KMN}\text{ s.t.
}V[i,j]=\begin{cases}[\alpha^{H}]_{j}\text{ }(i-1)K\leq j\leq iK\\\ 0\text{
otherwise}.\end{cases}\end{gathered}$

Using these models, we have two equivalent bilinear least squares problems
that can be solved alternately while fixing one of the variables. Keeping
$\boldsymbol{c}^{t}$ fixed at the $t^{\text{th}}$ iteration, we obtain a
linear least squares problem in $\boldsymbol{\alpha}$:

$\begin{gathered}\hat{\boldsymbol{\alpha}}^{t+1}=\operatorname*{\text{arg~{}min}}_{\boldsymbol{\alpha}}||\boldsymbol{y}-\underbrace{\boldsymbol{C}^{t}\bar{\boldsymbol{U}}}_{\boldsymbol{A}^{t}_{\boldsymbol{\alpha}}}\boldsymbol{\alpha}||_{2}^{2}\\\
\implies\hat{\boldsymbol{\alpha}}^{t+1}=((\boldsymbol{A}^{t}_{\boldsymbol{\alpha}})^{H}\boldsymbol{A}^{t}_{\boldsymbol{\alpha}})^{-1}(\boldsymbol{A}^{t}_{\boldsymbol{\alpha}})^{H}\boldsymbol{y}.\end{gathered}$

Next fixing $\boldsymbol{\alpha}^{t+1}$ we obtain a linear least squares
problem in $\boldsymbol{c}$:

$\displaystyle\begin{gathered}\hat{\boldsymbol{c}}^{t+2}=\operatorname*{\text{arg~{}min}}_{\boldsymbol{c}}||\bar{\boldsymbol{y}}-\underbrace{\boldsymbol{V}^{t+1}\boldsymbol{U}}_{\boldsymbol{A}^{t+1}_{\boldsymbol{c}}}\boldsymbol{c}||_{2}^{2}\\\
\implies\hat{\boldsymbol{c}}^{t+2}=((\boldsymbol{A}^{t+1}_{\boldsymbol{c}})^{H}\boldsymbol{A}^{t+1}_{\boldsymbol{c}})^{-1}(\boldsymbol{A}^{t+1}_{\boldsymbol{c}})^{H}\bar{\boldsymbol{y}}.\end{gathered}$

It is important to note that we solve an unconstrained bilinear least squares
problem. Therefore, the optimizers of the program are not unique, i.e., if the
algorithm recovers weights $\hat{\boldsymbol{\alpha}}$, $\hat{\boldsymbol{c}}$
then $\hat{\boldsymbol{\alpha}}/\gamma$ and $\gamma\hat{\boldsymbol{c}}$ are
also feasible solutions of the problem for $\gamma\neq 0$. To find $\gamma$
associated with the true angle, we can create a sequence of weights
$\boldsymbol{c}_{q}$ for each $\theta_{q}$ sampled from $\Theta$. For a
sufficiently fine-grained sampling, the true $\boldsymbol{c}$ will be well-
approximated by some
$\boldsymbol{c}\in\\{\boldsymbol{c}_{1},\ldots,\boldsymbol{c}_{Q}\\}$. Thus,
if we let

$q^{*}=\operatorname*{\text{arg~{}min}}_{q}\Big{|}\Big{|}\boldsymbol{c}_{q}-\gamma_{q}\hat{\boldsymbol{c}}\Big{|}\Big{|}^{2}_{2}$

where
$\gamma_{q}=\frac{\boldsymbol{c}_{q}^{T}\hat{\boldsymbol{c}}}{\hat{\boldsymbol{c}}^{T}\hat{\boldsymbol{c}}}$,
then $\gamma_{q^{*}}$ is the offset constant linked with the true angle and
$\hat{\boldsymbol{\alpha}}/\gamma_{q*}$ provides the estimate of the true
Slepian weights. Other approaches, like line spectral estimation [4], can also
be adapted for fine-tuning the constant $\gamma$. We empirically observe that
using a sufficiently good initialization for $\boldsymbol{c}^{0}$ results in
$\gamma\approx 1$. Thus, we start the alternating minimization scheme
described above with an initial value of $\boldsymbol{c}^{0}$ corresponding to
the mid-point of the interval $\Theta$, allowing the algorithm to converge in
$10$ to $15$ iterations for the experiments illustrated in the following
section.

##  IV Numerical results

We use the linear Slepian model in (4) to perform two kinds of beamforming to
benchmark the performance of the proposed model. In the first scenario, we
assume that the angle of arrival is known and true Slepian weights
($\boldsymbol{\alpha}_{true}$) are evaluated using (4). In the second case, we
use the midpoint of the uncertainty interval ($\Theta$) as an estimate for the
angle to recover a crude estimate of the Slepian weights. These are then used
for a quantitative comparison against the estimates obtained from the proposed
method. In addition, in the broadband regime, we also compare our algorithm
against conventional approaches like delay and sum beamforming [5]. In the
narrowband regime, we compare the performance with several robust adaptive
beamforming methods.

Figure 1: NMSE of Slepian coefficients ($\boldsymbol{\alpha}$) v/s nominal SNR
for 64-element ULA with 32 snapshots in the broadband regime with bandwidth
$5\text{\,}\mathrm{GHz}$.

We perform this experiment on two array geometries, a one-dimensional uniform
linear array (ULA) and a two-dimensional uniform planar array (UPA), for the
broadband system to validate our findings. In both array designs, the incoming
signal follows a sum of sinusoid model corrupted by complex Gaussian noise
[6]. We first consider a $64$ element ULA with $N=32$ snapshots and three kind
of uncertainty interval $\Theta$, i.e. the interval
$\Theta\equiv[\frac{\pi}{6}-\Delta,\frac{\pi}{6}+\Delta]$ where
$\Delta\in\\{\frac{\pi}{60},\frac{\pi}{40},\frac{\pi}{20}\\}$. Fig.1 plots the
NMSE
$\frac{||\boldsymbol{\hat{\alpha}}-\boldsymbol{\alpha}_{true}||_{2}}{||\boldsymbol{\alpha_{true}}||_{2}}$
for the crude linear Slepian and the proposed bilinear Slepian, we can observe
from Fig. 1 that the proposed method can recover the true Slepian coefficients
as the nominal SNR increases for all levels of uncertainty interval. The
marginal loss in beamforming performance compared to true linear Slepian
confirms this observation as seen in Fig. 2. On the other hand, beamforming
done using rough angle estimates suffers badly at higher nominal SNRs due to
large deviations from the true coefficients.

Next, we use a 16$\times$16 element UPA with $N=32$ snapshots. Two different
uncertainty intervals are defined for the elevation $\theta$ and the azimuth
angle $\phi$;
$\Theta_{\theta}\equiv[\frac{\pi}{4}-\frac{\pi}{40},\frac{\pi}{4}+\frac{\pi}{40}]$
and
$\Theta_{\phi}\equiv[\frac{\pi}{3}-\frac{\pi}{40},\frac{\pi}{3}+\frac{\pi}{40}]$.
Like the ULA experiment, the bilinear method is better at recovering the true
Slepian weights than the crude linear model. The delay and sum beamforming
employs fractional delay filters consisting of sinc interpolators truncated to
$R$ taps. The ideal beamformed SNR is achieved when $R\rightarrow\infty$, for
a truncated/finite $R$, the filter bias increases with an increase in SNR.
This causes deterioration in beamformed SNR as observed in Fig. 2 and Fig. 4.
Further degradation is observed in delay and sum when using a crude angle due
to an inaccurate estimate of the phase delays.

Figure 2: Nominal SNR compared to Beamformed SNR for a 64-element ULA with 32
snapshots in the broadband regime with bandwidth $5\text{\,}\mathrm{GHz}$.
Figure 3: NMSE of Slepian coefficients ($\boldsymbol{\alpha}$) v/s nominal SNR
for 16$\times$16-element UPA with 32 snapshots in the broadband regime with
bandwidth $5\text{\,}\mathrm{GHz}$.

In the narrowband regime, we compare our algorithm against robust adaptive
approaches like RCB [7], DCRCB [8], and MVDR-RAB [9], and the beamformed SNR
is plotted in Fig. 5. RCB uses a spherical uncertainty set to recover a
steering vector estimate. DCRCB adds a norm constraint to the constraints
defined by RCB, slightly improving the beamformed SNR. The spherical
uncertainty constraint in RCB, DCRCB is sensitive to the choice of parameter
$\epsilon$. We observe a performance drop in these algorithms for $\epsilon$
values further away from the true $\epsilon_{0}$. They are also sensitive to
array calibration errors and require perfect knowledge of the array geometry
and noise characteristics. MVDR-RAB, on the other hand, uses minimum prior
information and assumes the angle lies in an angular sector similar to
$\Theta$ defined in our case. The constraint added in MVDR-RAB ensures that
the estimated steering vector does not converge to any interfering signal
outside the interval $\Theta$. The optimization problem solved in MVDR-RAB is
given by

$\begin{gathered}\operatorname*{\text{arg~{}min}}_{\boldsymbol{a}}~{}\boldsymbol{a}^{H}\boldsymbol{R}^{-1}\boldsymbol{a}~{}\text{s.t.
}~{}\boldsymbol{a}^{H}\tilde{\boldsymbol{C}}\boldsymbol{a}\leq\Delta_{0},~{}||\boldsymbol{a}||^{2}_{2}=M,\end{gathered}$
(9)

where
$\tilde{\boldsymbol{C}}=\int_{\tilde{\Theta}}\boldsymbol{a}(\theta)\boldsymbol{a}^{H}(\theta)d\theta$.
In the narrowband system, MVDR-RAB performs as well as the proposed method
with only a marginal drop in performance; however, extending it to broadband
is not trivial due to several drawbacks inherent in the optimization problem
setup. As seen from (9), the norm constraint makes the problem non-convex,
requiring SDP-relaxation to find an optimal solution. In the case of
broadband, it is difficult to implement SDP relaxation due to the matrices
involved. Also, the objective function might not remain convex when extended
to broadband since there is no straightforward analog of the steering vector
for broadband. In contrast, the proposed bilinear algorithm is simple and
works regardless of our operating regime.

Figure 4: Nominal SNR compared to Beamformed SNR for a 16$\times$16 element
UPA with 32 snapshots in the broadband regime with bandwidth
$5\text{\,}\mathrm{GHz}$. Figure 5: Nominal SNR compared to Beamformed SNR for
a 64-element ULA with 32 snapshots in the narrowband regime with bandwidth
$5\text{\,}\mathrm{MHz}$.

## References

  * [1] A. Elbir, K. Mishra, S. Vorobyov, and R. Heath, “Twenty-Five Years of Advances in Beamforming: From convex and nonconvex optimization to learning techniques,” _IEEE Signal Process. Mag._ , vol. 40, no. 4, p. 118–131, 2023.
  * [2] C. DeLude, M. Davenport, and J. Romberg, “Slepian Beamforming: Broadband Beamforming using Streaming Least Squares,” arXiv:2312.03922, December 2023.
  * [3] D. Slepian, “On bandwidth,” _Proceedings of the IEEE_ , vol. 64, no. 3, pp. 292–300, 1976.
  * [4] C. DeLude, R. Sharma, S. Karnik, C. Hood, M. Davenport, and J. Romberg, “Iterative Broadband Source Localization,” _IEEE J. Sel. Areas Inf. Theory_ , vol. 4, pp. 453–469, 2023.
  * [5] H. V. Trees, _Optimum Array Processing: Part IV of Detection, Estimation, and Modulation Theory_. Wiley Interscience, 2002.
  * [6] M. Davenport and M. Wakin, “Compressive sensing of analog signals using Discrete Prolate Spheroidal Sequences,” _Appl. Comput. Harmon. Anal._ , vol. 33, no. 3, pp. 438–472, 2012.
  * [7] J. Li, P. Stoica, and Z. Wang, “On robust Capon beamforming and diagonal loading,” _IEEE Trans. Signal Process._ , vol. 51, no. 7, pp. 1702–1715, 2003.
  * [8] ——, “Doubly constrained robust Capon beamformer,” _IEEE Trans. Signal Process._ , vol. 52, no. 9, pp. 2407–2423, 2004.
  * [9] A. Khabbazibasmenj, S. Vorobyov, and A. Hassanien, “Robust Adaptive Beamforming Based on Steering Vector Estimation With as Little as Possible Prior Information,” _IEEE Trans. Signal Process._ , vol. 60, no. 6, pp. 2974–2987, 2012.
  * [10] J. Weidmann, _Linear Operators in Hilbert Spaces_. Springer, 2012.
  * [11] D. Crane, “The singular value expansion for compact and non-compact operators,” Doctoral Dissertation, Michigan Technological University, 2020.
  * [12] V. Hutson, J. Pym, and M. Cloud, _Applications of Functional Analysis and Operator Theory_ , 2nd ed. Elsevier, 2005.

##  Appendix A Singular value expansion for operator
$\boldsymbol{\mathcal{W}}$

To begin, recall that the finite-dimensional Hilbert space we are examining is

$\displaystyle\mathcal{H}=\text{Span}\left(\\{w_{\ell}(\theta)\\}_{\ell=1}^{L}\right)$

where each basis function $w_{\ell}(\theta)$ operates over the range
$\theta\in\Theta$. We use a slightly different indexing for ease of notation.
The operator we are examining,
$\boldsymbol{\mathcal{W}}:\mathcal{H}\to\mathbb{C}^{L}$, is given by

$\displaystyle\boldsymbol{\mathcal{W}}x=\begin{bmatrix}\langle
x,w_{1}\rangle\\\ \langle x,w_{2}\rangle\\\ \vdots\\\ \langle
x,w_{L}\rangle\end{bmatrix}.$

It is an operator between two finite-dimensional Hilbert spaces and,
therefore, is compact [10, Chapt. 4]. The first step in determining the
Singular value expansion (SVE) is determining the adjoint of the operator.
Note that the inner product for $\mathcal{H}$ is given by $\langle
w,x\rangle_{\mathcal{H}}=\int_{\Theta}x(\theta)\bar{w}(\theta)d\theta$,a and
the innner product for $\mathbb{C}^{L}$ is the standard vector inner product
$\langle\cdot,\cdot\rangle_{\mathbb{C}}$. Then we have

$\displaystyle\left\langle\boldsymbol{\mathcal{W}}x,y\right\rangle_{\mathbb{C}}$
$\displaystyle=\sum_{\ell=1}^{L}y_{\ell}\overline{\langle
x,w_{\ell}\rangle}_{\mathcal{H}}=\sum_{\ell=1}^{L}y_{\ell}\langle
w_{\ell},x\rangle_{\mathcal{H}}$
$\displaystyle=\sum_{\ell=1}^{L}y_{\ell}\int_{\Theta}w_{\ell}(\theta)\overline{x}(\theta)d\theta$
$\displaystyle=\int_{\Theta}\overline{x}(\theta)\sum_{\ell=1}^{L}y_{\ell}w_{\ell}(\theta)d\theta$
$\displaystyle=\left\langle
x,\sum_{\ell=1}^{L}y_{\ell}w_{\ell}\right\rangle_{\mathcal{H}},$

Hence the adjoint $\boldsymbol{\mathcal{W}}^{H}:\mathbb{C}^{L}\to\mathcal{H}$
is given by

$\displaystyle\boldsymbol{\mathcal{W}}^{H}y=\sum_{\ell=1}^{L}y_{\ell}w_{\ell}(\theta).$

Since the adjoint of a compact operator is also compact, therefore the
composite operators $\boldsymbol{\mathcal{W}}\boldsymbol{\mathcal{W}}^{H}$ and
$\boldsymbol{\mathcal{W}}^{H}\boldsymbol{\mathcal{W}}$ are also compact [10,
Chapt. 4]. It is easy to verify they are also self-adjoint and, therefore,
admit an eigendecomposition [11, Chapt. 1]. It is the case that the
eigenvectors of $\boldsymbol{\mathcal{W}}^{H}\boldsymbol{\mathcal{W}}$ are the
sequence of functions $v_{\ell}$, and similarly, the eigenvectors of
$\boldsymbol{\mathcal{W}}\boldsymbol{\mathcal{W}}^{H}$ are the sequence of
vectors $\boldsymbol{u}_{\ell}$ [12, Chapt. 7 ]. Through similar reasoning, it
is also the case that the non-zero eigenvalues of these self-adjoint operators
are the same and equal to $\\{\sigma^{2}_{\ell}\\}$. From this point, it is
not too difficult to see that

$\displaystyle\boldsymbol{\mathcal{W}}^{H}\boldsymbol{\mathcal{W}}x$
$\displaystyle=\sum_{\ell=1}^{L}w_{\ell}(\theta)\langle
w_{\ell},x\rangle_{\mathcal{H}}$
$\displaystyle=\sum_{\ell=1}^{L}\sigma_{\ell}^{2}v_{\ell}(\theta)\langle
v_{\ell},x\rangle_{\mathcal{H}}$

where the second line follows from the spectral theorem [11]. We can state
this more concisely as

$\displaystyle\boldsymbol{\mathcal{W}}^{H}\boldsymbol{\mathcal{W}}=\sum_{\ell=1}^{L}w_{\ell}(\theta)\otimes
w_{\ell}(\theta^{\prime})=\sum_{\ell=1}^{L}\sigma_{\ell}^{2}v_{\ell}(\theta)\otimes
v_{\ell}(\theta^{\prime})$

where the addition of a $\theta^{\prime}$ is meant to emphasize that the inner
product against the rightmost eigenfunction is taken independently from the
leftmost. In a similar manner, we can solve for the left eigenfunctions by
examining

$\displaystyle\boldsymbol{\mathcal{W}}\boldsymbol{\mathcal{W}}^{H}x$
$\displaystyle=\begin{bmatrix}\langle
w_{1},\sum_{\ell=1}^{L}x_{\ell}w_{\ell}(\theta)\rangle\\\ \langle
w_{2},\sum_{\ell=1}^{L}x_{\ell}w_{\ell}(\theta)\rangle\\\ \vdots\\\ \langle
w_{L},\sum_{\ell=1}^{L}y_{\ell}w_{\ell}(\theta)\rangle\end{bmatrix}.$
$\displaystyle=\int_{\Theta}\begin{bmatrix}w_{1}(\theta)\overline{w_{1}}(\theta)&\dots&w_{L}(\theta)\overline{w_{1}}(\theta)\\\
\vdots&\ddots&\vdots\\\
w_{1}(\theta)\overline{w_{L}}(\theta)&\dots&w_{L}(\theta)\overline{w_{L}}(\theta)\end{bmatrix}\begin{bmatrix}x_{1}\\\
x_{2}\\\ \vdots\\\ x_{L}\end{bmatrix}d\theta.$

which is just an $L\times L$ symmetric matrix operator, and as with all such
matrices, we can from an eigendecomposition on it to get a set of eigenvectors
$\\{\boldsymbol{u}_{\ell}\\}_{\ell=1}^{L}$. Putting it all together, we have

$\displaystyle\boldsymbol{\mathcal{W}}=\sum_{\ell=1}^{L}\sigma_{\ell}\boldsymbol{u}_{\ell}\otimes
v_{\ell}$

such that

$\displaystyle\boldsymbol{\mathcal{W}}x=\sum_{\ell=1}^{L}\sigma_{\ell}\boldsymbol{u}_{\ell}\langle
x,v_{\ell}\rangle_{\mathcal{H}}.$

Thus, the calculation and justification of the factorization present in the
robust beamforming method are concluded.