# Compressed Sensing Constant Modulus Constrained Projection Matrix Design and
High-Resolution DoA Estimation Methods

Khaled Ardah and Martin Haardt The authors gratefully acknowledge the support
of the German Research Foundation (DFG) under contract no. HA 2239/6-2
(EXPRESS II).K. Ardah and M. Haardt are with Communications Research
Laboratory (CRL), TU Ilmenau, Ilmenau, Germany (e-mail: {khaled.ardah,
martin.haardt}@tu-ilmenau.de).

###### Abstract

This paper proposes a compressed sensing-based high-resolution direction-of-
arrival estimation method called gradient orthogonal matching pursuit (GOMP).
It contains two main steps: a sparse coding approximation step using the well-
known OMP method and a sequential iterative refinement step using a newly
proposed gradient-descent method. To enhance the recoverability, we further
propose an efficient projection matrix design method, which considers the
constant modulus constraints imposed by the projection matrix hardware
components. Simulation results show the effectiveness of the proposed methods
as compared to benchmark algorithms.

###### Index Terms:

Off-grid DoA estimation, projection matrix design, constant modulus
constraints, compressed sensing.

##  I Introduction

Direction-of-arrival (DoA) estimation is a fundamental problem in a large
variety of applications including radar and wireless communications. Over the
years, various DoA estimation methods have been proposed [1], including the
subspace-based methods MUSIC and ESPRIT [2, 3]. While subspace-based methods
have been shown to perform asymptotically optimal, these methods suffer from a
performance degradation in difficult scenarios such as high noise power and
low number of measurements. Recently, compressed sensing (CS) techniques [4]
have been introduced as an attractive alternative to subspace-based methods.
By exploiting the sparsity property, it has been shown that CS-based
techniques exhibit a good estimation accuracy, even in difficult scenarios and
are able to recover the signal of interest from far fewer samples than
required by traditional acquisition systems [5, 6].

Generally, CS-based methods can be categorized into on-grid, off-grid, and
gridless [7]. The former approximates the DoA parameters using a given
dictionary, which is formed by sampling the desired domain into a number of
grid points. A well-known on-grid CS-based technique is the orthogonal
matching pursuit (OMP) [8]. However, if the true DoAs belong to a continuous
domain, which is the case in general, on-grid methods suffer from the sampling
errors, since the number of grid-points are practically limited. To overcome
this issue, off-grid methods extend the on-grid methods by introducing a
refinement step to reduce or even eliminate the sampling errors [9]. On the
other hand, gridless methods [10] do not require a dictionary, but rather
operate directly on the continuous domain, thus, avoiding sampling errors.
However, gridless methods are, generally, very complex especially with large-
scale arrays, which limits their practical implementation.

In this paper, as our first contribution, we propose a high-resolution CS-
based DoA estimation method called Gradient OMP (GOMP). Specifically, GOMP
contains two steps. In Step I, the OMP method is used to approximate the true
DoA parameters to their best on grid points. After that, the approximated DoA
parameters are refined sequentially and iteratively in Step II using a novel
proposed method to eliminate the sampling errors. Simulation results show that
GOMP exhibits a high-resolution DoA estimation accuracy with low computational
complexity.

In CS-based methods, the sensing matrix should be designed carefully so that
it satisfies a certain property, e.g., the mutual coherence to enhance
recoverability of desired signal [4, 6, 5, 11]. In general, the sensing matrix
comprises a tunable $N\times M$ projection matrix and a fixed $M\times P$
dictionary matrix, where $N$, $M$, and $P$, with $N\leq M\leq P$, denote the
number of dedicated radio-frequency (RF) chains, the number of array elements,
and the number of sampling (grid) points, respectively. Considering such a
structure, several sensing matrix design methods via mutual coherence
minimization have been proposed, e.g., in [12, 13, 14, 15, 16, 6]. Assuming
that $M=P$, we have proposed in [12] a sensing matrix design method called
sequential mutual coherence minimization (SMCM), where the original nonconvex
problem is relaxed and divided into several convex sub-problems that are
updated sequentially using alternating optimization. On the other hand,
gradient-descent (GD)-based sensing matrix design methods are proposed in [13]
and [14]. In contrast to [13], the method in [14] imposes the unit-norm
constraints directly into the objective function, leading to a better
solution.

In many applications, however, the projection matrix is realized using, e.g.,
a network of tunable phase-shifters [17]. This implies that the entries of the
projection matrix are subject to constant modulus (CM) constraints so that
only their phase angles are adjustable. In [6], we have proposed a GD-based
projection matrix design method, which takes into account, for the first time,
the CM constraints imposed by the projection matrix hardware components. In
this paper, as our second contribution, we propose a new GD-based projection
matrix design method by extending and enhancing the proposed methods in [6]
and [14] while taking into account the CM constraints. Using computer
simulations, we show that the new proposed GD-based method outperforms the
baseline methods leading to a higher DoA estimation accuracy.

##  II System model

We assume that $K$ narrow-band far-field source signals impinge on an array
composed of $M$ omni-directional sensors. The output vector at the $\ell$th
time sample is given by

$\displaystyle\bar{{\bf y}}_{[\ell]}={\bf A}(\bm{\nu}){\bf
x}_{[\ell]}+{\bf\bar{n}}_{[\ell]}\in\mathbb{C}^{M},\quad\ell\in\\{1,\dots,L\\},$
(1)

where ${\bf x}_{[\ell]}\in\mathbb{C}^{K}$ is the vector of waveforms emitted
by the $K$ sources, ${\bf\bar{n}}_{[\ell]}\in\mathbb{C}^{M}$ contains the
spatially and temporally white circular Gaussian sensor noise, and $L$ denotes
the number of available time samples. The matrix ${\bf A}(\bm{\nu})=[{\bf
a}(\nu_{1}),\dots,{\bf a}(\nu_{K})]\in\mathbb{C}^{M\times K}$ denotes the true
array steering matrix, whose $k$th column is the array response vector ${\bf
a}(\nu_{k})$ corresponding to the DoA of the $k$th source $\nu_{k}$,
$k\in\\{1,\dots,K\\}$. Here, $\nu_{k}=\frac{2\pi d}{\lambda}\sin(\theta_{k})$
is the DoA parameter in the spatial frequency domain,
$\theta_{k}\in[-90^{\circ},90^{\circ}]$ is the DoA parameter in the angular
domain, $d$ is the inter-sensor spacing, $\lambda$ is the wavelength. We
denote $\bm{\nu}=[\nu_{1},\dots,\nu_{K}]^{\mathsf{T}}\in\mathbb{R}^{K}$ as the
true DoA parameter vector. The array response vector ${\bf a}(\nu)$ describes
a manifold denoted as $\mathbb{M}^{M}$. In this paper, we assume a uniform
linear array (ULA), where ${\bf a}(\nu)$ is given by ${\bf
a}(\nu)=[1,e^{j\nu},\dots,e^{j(M-1)\nu}]^{\mathsf{T}}\in\mathbb{C}^{M}$.

Let $\bm{\Phi}\in\mathbb{C}^{N\times M}$ denote the analog projection matrix,
which compresses the output of $M$ antennas to $N$ active RF channels. We
assume that the analog domain is realized in hardware using a network of phase
shifters. Therefore, the entries of $\bm{\Phi}$ must satisfy the CM
constraints as $|\phi_{n,m}|=1,\forall n,m$, where $\phi_{n,m}$ is the
$(n,m)$th entry of $\bm{\Phi}$. Then, the complex baseband array output in
(1), after projection, can be expressed as

$\displaystyle{{\bf y}}_{[\ell]}=\bm{\Phi}\bar{{\bf y}}_{[\ell]}=\bm{\Phi}{\bf
A}(\bm{\nu}){\bf x}_{[\ell]}+{\bf n}_{[\ell]}\in\mathbb{C}^{N},$ (2)

where ${\bf n}_{[\ell]}=\bm{\Phi}{\bf\bar{n}}_{[\ell]}\in\mathbb{C}^{N}$. Let
${\bf X}=[{\bf x}_{[1]},\dots,{\bf x}_{[L]}]\in\mathbb{C}^{K\times L}$ and
${\bf N}=[{\bf n}_{[1]},\dots,{\bf n}_{[L]}]\in\mathbb{C}^{N\times L}$. Then,
the measurement matrix ${\bf Y}=[{\bf y}_{[1]},\dots,{\bf
y}_{[L]}]\in\mathbb{C}^{N\times L}$ can be written as

$\displaystyle{\bf Y}=\bm{\Phi}{\bf A}(\bm{\nu}){\bf X}+{\bf
N}\in\mathbb{C}^{N\times L}.$ (3)

In Section III, we propose the Gradient OMP (GOMP) method for estimating the
DoA parameter vector $\bm{\nu}$ and the signal matrix ${\bf X}$ from ${\bf Y}$
with a high-resolution. Next, we propose in Section IV a new GD-based
projection matrix design method while taking into account the CM constraints
imposed by the projection matrix hardware components.

##  III Proposed GOMP For DoA estimation

###  III-A The single-source case

To simplify the exposition, we first assume a single-source case, i.e., $K=1$
and that $\bm{\Phi}$ is given and fixed. Therefore, after removing the source
index $k$, (3) simplifies to

$\displaystyle{\bf Y}=\bm{\Phi}{\bf a}(\nu){\bf x}^{\mathsf{T}}+{\bf
N}\in\mathbb{C}^{N\times L},$ (4)

where ${\bf x}\in\mathbb{C}^{L}$. To estimate the spatial frequency ${\nu}$
and the source signal ${\bf x}$, the problem can be written as

$\displaystyle\underset{{\bf x},\nu}{\arg\min}\|{\bf Y}-\bm{\Phi}{{\bf
a}}({\nu}){\bf x}^{\mathsf{T}}\|^{2}_{\text{F}},$ (5)

which is non-convex due to its joint optimization. Therefore, we decouple the
recovery procedure between ${\nu}$ and ${\bf x}$, where we update for one
variable assuming the other is fixed.

First, its not hard to see that the optimal solution to ${\bf x}$, in the
least squares sense, for a fixed $\nu$ is given as

$\displaystyle{\bf x}=(\\{\bm{\Phi}{{\bf a}}({\nu})\\}^{+}{\bf
Y})^{\mathsf{T}}\in\mathbb{C}^{L}.$ (6)

On the other hand, to estimate $\nu$ for a fixed ${\bf x}$, we start by
writing the steering vector ${\bf a}(\nu)$ using its first order Taylor series
expansion as

$\displaystyle{\bf a}(\nu)\approx{\bf a}(\mathring{\nu})+{\bf
g}(\mathring{\nu})\cdot\delta\in\mathbb{C}^{M},$ (7)

where $\mathring{\nu}$ is a known close approximation to $\nu$ (e.g., the best
on-grid sampling frequency), $\delta=\nu-\mathring{\nu}$ is the sampling error
between the true $\nu$ and the approximated $\mathring{\nu}$, and ${\bf
g}(\mathring{\nu})$ is the gradient vector of ${\bf a}(\mathring{\nu})$ with
respect to $\mathring{\nu}$ which is given as

$\displaystyle{\bf g}(\mathring{\nu})=j\text{diag}\\{0,1,\dots,M-1\\}{\bf
a}(\mathring{\nu})\in\mathbb{C}^{M}.$ (8)

By substituting (7) into (4), the matrix ${\bf Y}$ can be written as

$\displaystyle{\bf Y}\approx\bm{\Phi}{\bf a}(\mathring{\nu}){\bf
x}^{\mathsf{T}}+\bm{\Phi}{\bf g}(\mathring{\nu}){\bf
x}^{\mathsf{T}}\cdot\delta+{\bf N},$ (9)

where the only unknown variable is the sampling error $\delta$. To this end,
we note from (7) that for any given $\mathring{\nu}$, the true DoA spatial
frequency $\nu$ can also be obtained as

$\displaystyle\nu=\mathring{\nu}+\delta,$ (10)

where $\delta$ here can be seen as the step-size that indicates the shift
value (i.e., the sampling error correction) that should be added to or
subtracted from $\mathring{\nu}$, depending on its sign, so that we have a
perfect estimation of $\nu$. Therefore, for a given approximation
$\mathring{\nu}$, finding the true DoA parameter $\nu$ requires finding the
best step-size $\delta$. The problem can be written as

$\displaystyle\underset{\delta}{\arg\min}\|{\bf Y}-\bm{\Phi}{\bf
a}(\mathring{\nu}){\bf x}^{\mathsf{T}}-\bm{\Phi}{\bf g}(\mathring{\nu}){\bf
x}^{\mathsf{T}}\cdot\delta\|_{\text{F}}^{2}.$ (11)

By applying the $\text{vec}\\{{\bf A}\text{diag}\\{{\bf b}\\}{\bf C}\\}=({\bf
C}^{\mathsf{T}}\diamond{\bf A}){\bf b}$ property, where $\diamond$ denotes the
Khatri-Rao product, the vectorized form of ${\bf Y}$ can be expressed as

$\displaystyle{\bf y}=\text{vec}\\{{\bf Y}\\}=({\bf x}\diamond\bm{\Phi}{\bf
a}(\mathring{\nu}))+({\bf x}\diamond\bm{\Phi}{\bf
g}(\mathring{\nu}))\cdot\delta+{{\bf n}},$ (12)

where ${{\bf n}}=\text{vec}\\{{\bf N}\\}$. Therefore, (11) can be rewritten as

$\displaystyle\underset{\delta}{\arg\min}\|{\bf y}-({\bf
x}\diamond\bm{\Phi}{\bf a}(\mathring{\nu}))-({\bf x}\diamond\bm{\Phi}{\bf
g}(\mathring{\nu}))\cdot\delta\|_{2}^{2},$ (13)

so that a real-valued solution to $\delta$, for fixed $\mathring{\nu}$ and
${\bf x}$, can be obtained as

$\displaystyle\delta=\Re\big{\\{}\\{{\bf x}\diamond\bm{\Phi}{\bf
g}(\mathring{\nu})\\}^{+}({\bf y}-({\bf x}\diamond\bm{\Phi}{\bf
a}(\mathring{\nu})))\big{\\}}.$ (14)

Using (6), (10) and (14), a solution to (5), i.e., to ${\bf x}$ and to $\nu$
can be obtained as follows: Let us denote $\mathring{\nu}$ by ${\nu}^{(0)}$,
i.e., the first estimate of $\nu$ obtained using a CS-based on-grid method,
e.g. OMP. Then, the first estimate of ${\bf x}$, i.e., ${\bf x}^{(0)}$ can be
obtained using (6). After that, we refine $\delta^{(i)}$, $\nu^{(i)}$, and
${\bf x}^{(i)}$ iteratively, as summarized in Algorithm 1, where $i$ denotes
the inner iteration index. The proposed iterative steps are repeated until a
maximum number of iterations $I_{\text{max}}$ is reached. Note that, in Step
9, we have included an update acceptance condition, which implies that GOMP
accepts the new updates ${\nu^{(i+1)}}$ and ${{\bf x}^{(i+1)}}$ only if they
reduce the current cost-function $\epsilon^{(i)}$. This means that Algorithm 1
has a monotonic and a guaranteed convergence, since we always have
$\epsilon_{(i+1)}\leq\epsilon_{(i)}$ and the cost-function $\epsilon$ is
bounded as $0\leq\epsilon_{(i)}\leq\epsilon_{(0)}$.

Algorithm 1 GOMP: The single-source case.

1:Input: Measurement matrix ${\bf Y}$ and initials ${\nu}^{(0)}$ and ${\bf
x}^{(0)}$

2:Compute $\epsilon^{(0)}=\|{\bf Y}-\bm{\Phi}{\bf a}(\nu^{(0)})({\bf
x}^{(0)})^{\mathsf{T}}\|^{2}_{\text{F}}$

3:Choose $I_{\text{max}}$ and set $i=1$

4:while $i\leq I_{\text{max}}$ do

5: For given ${\bf x}^{(i-1)}$ and ${\nu}^{(i-1)}$, compute ${\delta}^{(i)}$
using (14)

6: For given ${\nu}^{(i-1)}$ and $\delta^{(i)}$, compute $\nu^{(i)}$ using
(10)

7: For given $\nu^{(i)}$, compute ${\bf x}^{(i)}$ using (6)

8: Compute $\epsilon^{(i)}=\|{\bf Y}-\bm{\Phi}{\bf a}(\nu^{(i)})({\bf
x}^{(i)})^{\mathsf{T}}\|^{2}_{\text{F}}$

9: if $\epsilon^{(i)}>\epsilon^{(i-1)}$ then

10: Return ${\hat{\nu}}={\nu}^{(i-1)}$ and ${\bf\hat{x}}={\bf x}^{(i-1)}$ and
terminate

11: end if

12:end while

13:Output: The estimated parameters ${\hat{\nu}}$ and ${\bf\hat{x}}$

###  III-B The general case

For the general case, i.e., $K\geq 1$, the signal model in (3), i.e., ${\bf
Y}=\bm{\Phi}{\bf A}(\bm{\nu}){\bf X}+{\bf N}$ can also be written as

$\displaystyle{\bf Y}=\bm{\Phi}{\bf a}({\nu_{k}}){\bf
x}^{\mathsf{T}}_{k}+\sum_{k^{\prime}\neq k}^{K}\bm{\Phi}{\bf
a}({\nu_{k^{\prime}}}){\bf x}^{\mathsf{T}}_{k^{\prime}}+{\bf
N}\in\mathbb{C}^{N\times L},$ (15)

where ${\bf x}_{k}\in\mathbb{C}^{L}$ is the $k$th column-vector of ${\bf
X}^{\mathsf{T}}\in\mathbb{C}^{L\times K}$ and ${\bf
a}({\nu_{k}})\in\mathbb{C}^{M}$ is the $k$th column vector of ${\bf
A}(\bm{\nu})$. Therefore, Algorithm 1 can be utilized to estimate ${\nu}_{k}$
and ${\bf x}_{k},\forall k$, sequentially as summarized in Algorithm 2. In
Step 2, for a given measurement matrix ${\bf Y}$ and a dictionary matrix
${\bf\mathring{A}}(\bm{\mathring{\nu}})\in\mathbb{C}^{M\times P}$ with $P\gg
K$, we use the CS-based on-grid method OMP to obtain the initial estimates[1]
[1] [1] Note that (3) can be written in a sparse-form as [11]: ${\bf
Y}\approx\bm{\Phi}{\bf\mathring{A}}(\bm{\mathring{\nu}}){\bf\mathring{X}}+{\bf
N}\in\mathbb{C}^{N\times L}$, where
${\bf\mathring{A}}(\bm{\mathring{\nu}})=[{\bf
a}(\mathring{\nu}_{1}),\dots,{\bf
a}(\mathring{\nu}_{P})]\in\mathbb{C}^{M\times P}$ denotes a dictionary matrix
formed using the grid points
$\bm{\mathring{\nu}}=[\mathring{\nu}_{1},\dots,\mathring{\nu}_{P}]^{\mathsf{T}}\in\mathbb{R}^{P}$,
with $P\gg K$ denoting the number of grid points, and
${\bf\mathring{X}}\in\mathbb{C}^{P\times N}$ is an ${K}$ row-sparse matrix.
Therefore, the on-grid CS-based techniques, e.g., OMP [8] can readily be
applied to estimate $\bm{\nu}^{(0)}$ and ${\bf X}^{(0)}$. of the spatial
frequency vector $\bm{\nu}^{(0)}\in\mathbb{C}^{K}$ and the source signal
matrix ${\bf X}^{(0)}\in\mathbb{C}^{K\times L}$. After that, Algorithm 2 runs
for a maximum of $J_{\text{max}}$ outer iterations to refine the initial
parameters ${\nu}^{(0)}_{k}$ and ${\bf x}^{(0)}_{k},\forall k$, sequentially.
Specifically, to refine the $k$th source parameters, i.e.,
${\nu}^{(j-1)}_{k}$, and ${{\bf x}}^{(j-1)}_{k}$ at the $j$th outer iteration,
we first calculate ${\bf Y}^{(j)}_{k}$ as shown in Step 9, where ${\bf
Y}^{(j)}_{\text{new},k}$ is calculated using the updated parameters in the
current $j$th iteration, i.e., $\nu^{(j)}_{1},\dots,\nu^{(j)}_{k-1}$ and ${\bf
x}^{(j)}_{1},\dots,{\bf x}^{(j)}_{k-1}$ while ${\bf Y}^{(j)}_{\text{old},k}$
is calculated using the updated parameters in the previous $(j-1)$th
iteration, i.e., $\nu^{(j-1)}_{k+1},\dots,\nu^{(j-1)}_{K}$ and ${\bf
x}^{(j-1)}_{k+1},\dots,{\bf x}^{(j-1)}_{K}$. Given ${\bf Y}^{(j)}_{k}$,
${\nu}^{(j-1)}_{k}$, and ${{\bf x}}^{(j-1)}_{k}$, Algorithm 1 is then used in
Step 10 to obtain ${\nu}^{(j)}_{k}$ and ${{\bf x}}^{(j)}_{k}$.

Algorithm 2 GOMP: The general case.

1:Input: Measurement matrix ${\bf Y}$, dictionary
${\bf\mathring{A}}\in\mathbb{C}^{M\times P}$

2:For given ${\bf Y}$ and ${\bf\mathring{A}}$, obtain $\bm{\nu}^{(0)}$ and
${\bf X}^{(0)}$ using OMP

3:Choose $J_{\text{max}}$ and set $j=1$

4:while $j\leq J_{\text{max}}$ do

5: Set $\bm{\nu}^{(j)}=\bm{0}_{K}$ and ${\bf X}^{(j)}=\bm{0}_{L\times K}$

6: for $k=1$ to $K$ do

7: Get ${\bf Y}^{(j)}_{\text{new},k}=\sum_{k^{\prime}=1}^{k-1}\bm{\Phi}{\bf
a}({\nu^{(j)}_{k^{\prime}}})({\bf x}^{(j)}_{k^{\prime}})^{\mathsf{T}}$

8: Get ${\bf Y}^{(j)}_{\text{old},k}=\sum_{k^{\prime}=k+1}^{K}\bm{\Phi}{\bf
a}({\nu^{(j-1)}_{k^{\prime}}})({\bf x}^{(j-1)}_{k^{\prime}})^{\mathsf{T}}$

9: Get ${\bf Y}^{(j)}_{k}={\bf Y}-{\bf Y}^{(j)}_{\text{new},k}-{\bf
Y}^{(j)}_{\text{old},k}$

10: Get $[{\nu}^{(j)}_{k},{{\bf x}}^{(j)}_{k}]\leftarrow\text{GOMP}\big{(}{\bf
Y}^{(j)}_{k},\nu^{(j-1)}_{k},{\bf x}^{(j-1)}_{k}\big{)}$

11: end for

12:end while

13:Output: The estimated ${\bm{\hat{\nu}}}$ and ${{\bf\hat{X}}}$

##  IV Proposed EGD for projection matrix design

From the above, it is clear that the estimation accuracy and the complexity of
the proposed GOMP method (i.e., Algorithm 2) is highly dependent on the
initial on-grid estimates $\bm{\nu}^{(0)}$ and ${\bf X}^{(0)}$. Specifically,
we have observed that if the OMP method correctly approximates the true
frequencies $\bm{\nu}$ to its nearest on-grid points $\bm{\nu}^{(0)}$ (i.e.,
$\nu^{(0)}_{k}=\underset{\mathring{\nu}_{p}\in\bm{\mathring{\nu}}}{\arg\min}|\nu_{k}-\mathring{\nu}_{p}|,\forall
k$), the refinement steps of Algorithm 1 always lead to a high accuracy
estimation. Like any CS-based method, the performance of the OMP is dependent
on the employed sensing matrix that is given by
${\bf\Psi}={\bf\Phi}{\bf\mathring{A}}(\bm{\mathring{\nu}})\in\mathbb{C}^{N\times
P}$, which should be designed carefully so that it satisfies a certain
property, e.g., by minimizing its mutual coherence $\mu_{\max}({\bf\Psi})$
defined as

$\displaystyle\mu_{\max}({\bf\Psi})=\max_{i\neq
j}\frac{|{\bm{\psi}}_{i}^{\text{H}}{\bm{\psi}}_{j}|}{\|{\bm{\psi}}_{i}\|_{2}\|{\bm{\psi}}_{j}\|_{2}},$
(16)

with columns ${\bm{\psi}}_{p}\in\mathbb{C}^{N},p\in\\{1,\dots,P\\}$. Clearly,
a large coherence $\mu_{\max}({\bf\Psi})$ means that there exist, at least,
two highly correlated columns in ${\bf\Psi}$, which may confuse any sparse
coding technique like OMP. Nonetheless, if
$K<\frac{1}{2}\Big{(}1+\frac{1}{\mu_{\max}({\bf\Psi})}\Big{)}$, the OMP is
guaranteed to correctly approximates $\bm{\nu}$ to the nearest on-grid points
$\bm{\nu}^{(0)}$ with overwhelming probability [6].

Formally, by assuming that the dictionary
${\bf\mathring{A}}(\bm{\mathring{\nu}})\in\mathbb{C}^{M\times P}$ is fixed,
the problem of sensing matrix design reduces to finding the projection matrix
${\bf\Phi}\in\mathbb{C}^{N\times M}$ with CM entries so that the coherence
$\mu_{\max}({\bf\Psi})$ is minimized, which can be expressed as

$\displaystyle\underset{{\bf\Phi}}{\min}$
$\displaystyle\quad\mu_{\max}({\bf\Psi})$ (17) s.t.
$\displaystyle\quad\|{{\bm{\psi}}_{p}}\|_{2}=1,\forall p,\text{ and
}|\phi_{n,m}|=1,\forall n,m.$

Note that problem (17) is non non-convex and NP-hard [12]. Therefore, we
propose to solve (17) indirectly as [6]

$\displaystyle\underset{{\bf\Phi}}{\min}\text{ }$
$\displaystyle\quad\|{\bf\Psi}^{\text{H}}{\bf\Psi}-{\bf
I}_{P}\|^{2}_{\text{F}}$ (18) s.t.
$\displaystyle\quad\|{{\bm{\psi}}_{p}}\|_{2}=1,\forall p,\text{ and
}|\phi_{n,m}|=1,\forall n,m.$

Similarly to [14], we reformulate (18) so that the unit-norm constraints are
directly embedded into the objective function. To this end, let
$\bm{\Psi}={\bf Q}{\bf D}$, where ${\bf
Q}={\bf\Phi}{\bf\mathring{A}}(\bm{\mathring{\nu}})=[{\bf q}_{1},\dots,{\bf
q}_{P}]\in\mathbb{C}^{N\times P}$ and ${\bf D}=\text{
diag}\Big{\\{}\frac{1}{\|{\bf q}_{1}\|_{2}},\dots,\frac{1}{\|{\bf
q}_{P}\|_{2}}\Big{\\}}\in\mathbb{R}^{P\times P}$. Since the columns of the
matrix $\bm{\Psi}$ have unit-norm, (18) can be rewritten as

$\displaystyle\underset{{\bf\Phi}}{\min}\text{ }$
$\displaystyle\quad\eta({\bf\Phi})$ (19) s.t.
$\displaystyle\quad{\bf\Psi}={\bf Q}{\bf D}\text{ and }{\bf
Q}={\bf\Phi}{\bf\mathring{A}}(\bm{\mathring{\nu}}),$ $\displaystyle\quad{\bf
D}=\text{ diag}\Bigg{\\{}\frac{1}{\|{\bf q}_{1}\|_{2}},\dots,\frac{1}{\|{\bf
q}_{P}\|_{2}}\Bigg{\\}},$ $\displaystyle\quad|\phi_{n,m}|=1,\forall n,m,$

where $\eta({\bf\Phi})=\|{\bf D}{\bf Q}^{\text{H}}{\bf Q}{\bf D}-{\bf
I}_{P}\|^{2}_{\text{F}}$ and ${\bf D}={\bf D}^{\mathsf{H}}$ due to its
diagonal and real-valued structure. To obtain a solution for (19), we propose
a CM constrained GD-based method, which updates the projection matrix
${\bf\Phi}$ iteratively as

$\displaystyle{\bf\Phi}^{(t)}={\bf\Pi}\bigg{(}{\bf\Phi}^{(t-1)}-\zeta\cdot\frac{\partial\eta({\bf\Phi})}{\partial{\bf\Phi}}\bigg{|}_{{\bf\Phi}={\bf\Phi}^{(t-1)}}\bigg{)},$
(20)

where $t$ is the iteration index, $\zeta$ is the step-size, ${\bf\Pi}(\cdot)$
is a projection function that imposes the CM constraints on the entries of the
input matrix entry-wise, i.e., ${\bf\Pi}(z)=\frac{z}{|z|}$, and
$\frac{\partial\eta({\bf\Phi})}{\partial{\bf\Phi}}$ is the gradient of
$\eta({\bf\Phi})$ w.r.t ${\bf\Phi}$, which is given as [14]

$\displaystyle\frac{\partial\eta({\bf\Phi})}{\partial{\bf\Phi}}=$
$\displaystyle\text{ }4{\bf Q}{\bf D}{\bf E}{\bf
D}{\bf\mathring{A}}(\bm{\mathring{\nu}})^{\text{H}}-2{\bf\Phi}{\bf\mathring{A}}(\bm{\mathring{\nu}}){\bf
R}{\bf\mathring{A}}(\bm{\mathring{\nu}})^{\text{H}},$ (21)

where ${\bf E}={\bf D}{\bf Q}^{\text{H}}{\bf Q}{\bf D}-{\bf I}_{P}$, ${\bf
R}=\text{diag}\\{[{\bf C}]_{[1,1]},\dots,[{\bf C}]_{[P,P]}\\}$, and ${\bf
C}=2{\bf E}{\bf D}{\bf Q}^{\text{H}}{\bf Q}{\bf D}^{3}$. Note that, the update
step in (20) is a direct extension of the proposed unconstrained GD method in
[14] to accounts for the CM constraints. Our results show that both the
unconstrained method of [14] and its direct extension in (20) achieve a mutual
coherence that is far from the known theoretical Welch lower-bound
$\beta=\sqrt{\frac{P-N}{N(P-1)}}$ [4]. To enhance their performance, we first
note that ${\bf E}$ in (21) represents the error matrix of the objective
function $\eta({\bf\Phi})$. Here, we propose to apply a shrinking operator on
${\bf E}$ element-wise to get $\underline{{\bf E}}$, where the $(k,j)$th entry
of $\underline{{\bf E}}$, i.e., $\underline{e}_{k,j}$ is obtained as [6]

$\displaystyle\underline{e}_{k,j}=\begin{cases}0,&\big{|}{e}_{k,j}\big{|}<\alpha\cdot{\beta}\\\
\frac{{e}_{k,j}}{|{e}_{k,j}|}\cdot\Big{(}\big{|}{e}_{k,j}\big{|}-\alpha\cdot{\beta}\Big{)},&\text{otherwise},\end{cases}$
(22)

where ${e}_{k,j}$ is the $(k,j)$th entry of ${{\bf E}}$, $\alpha\geq 1$ is a
relaxation parameter, and $\beta$ is the Welch lower-bound defined above. In
[6], we have shown that the choice of $\alpha$ highly impacts the mutual
coherence $\mu_{\max}(\bm{\Psi})$ of the obtained sensing matrix, where a
very-large or a very-small value leads to a high $\mu_{\max}(\bm{\Psi})$.
Therefore, $\alpha$ should be chosen carefully depending on the given
dictionary matrix.

After a closer look at (22), we can observe that for a very tight threshold
$\bar{\beta}=\alpha\cdot{\beta}$, the resulting error matrix $\underline{{\bf
E}}$ becomes a sparse matrix, where some of its entries that are smaller than
$\bar{\beta}$ will be set to zero. The direct implication of such a shrinking
operator is that the new projection matrix ${\bf\Phi}^{(t)}$ will be updated
so that it mainly reduces the entries that are larger than $\bar{\beta}$. In
summary, the proposed GD method for sensing/projection matrix design is given
by Algorithm 3.

Algorithm 3 Proposed GD method for sensing matrix design.

1:Inputs: Initial ${\bm{\Phi}}^{(0)}\in\mathbb{C}^{N\times M}$ and dictionary
${\bf\mathring{A}}(\bm{\mathring{\nu}})$.

2:Select $T_{\max}$ and step-size $\zeta$. Set $t=1$

3:while $t\leq T_{\max}$ do

4: For given ${\bf\Phi}^{(t-1)}$, compute ${\bf E}^{(t-1)}$

5: Apply the shrinking (22) on ${\bf E}^{(t-1)}$ to get $\underline{{\bf
E}}^{(t-1)}$

6: For given $\underline{{\bf E}}^{(t-1)}$ and ${\bf\Phi}^{(t-1)}$, compute
${\bf\Phi}^{(t)}$ using (20)

7:end while

8:Output: Sensing matrix
${\bf\Psi}^{\star}={\bf\Phi}^{\star}{\bf\mathring{A}}(\bm{\mathring{\nu}})$

##  V Numerical Results

Example 1: Evaluation of Algorithm 3. In this example, we set $N=16$, $M=64$,
and design the dictionary matrix as
${\bf\mathring{A}}(\bm{\mathring{\nu}})=[{\bf
a}(\mathring{\nu}_{1}),\dots,{\bf
a}(\mathring{\nu}_{P})]\in\mathbb{C}^{M\times P}$, where
$\mathring{\nu}_{p}=\frac{2\pi(p-1)}{P},p\in\\{1,\dots,P\\}$, so that
${\bf\mathring{A}}(\bm{\mathring{\nu}}){\bf\mathring{A}}(\bm{\mathring{\nu}})^{\mathsf{H}}=P{\bf
I}_{M}$. For comparison, we include the results when the projection matrix
${\bm{\Phi}}$ is designed from a DFT matrix[2] [2] [2] Let ${\bf W}_{64}$
denotes the $64\times 64$ DFT matrix. Then, the projection matrix is given as
${\bm{\Phi}}=[{\bf W}_{64}]_{[{\bf n},:]}\in\mathbb{C}^{16\times 64}$, where
${\bf n}=[1,5,9,13,17,\dots,64]^{\mathsf{T}}$., randomly[3] [3] [3] The
$(i,j)$th entry of ${\bm{\Phi}}$ is given as
$[{\bm{\Phi}}]_{[i,j]}=e^{j\theta_{i,j}},\theta_{i,j}\in[0,2\pi]$., using the
GD-based method of [6], or using the extended GD-based method of [14] to
accounts for the CM constraints[4] [4] [4] Note that the main difference
between our proposed GD-based method in Algorithm 3 and the GD-based method in
[6] is that the former imposes the unit-norm constraints directly into the
objective function, while the latter imposes the unit-norm constraints after
updating the projection matrix on every iteration. On the other hand, the main
difference between Algorithm 3 and the extended GD-based method of [14] is
that the former uses the error matrix $\underline{{\bf E}}$ obtained using
(22) when updating ${\bm{\Phi}}$ on every iteration, while the latter uses the
error matrix ${\bf E}$ when updating ${\bm{\Phi}}$ on every iteration.. For
the GD-based methods, we assume that the initial projection matrix, i.e.,
${\bm{\Phi}}^{(0)}$ is obtained using the closed-form approach discussed in
[6], since we have observed that it always leads to a better result than a
random initialization approach. Moreover, the relaxation parameter $\alpha$ is
chosen as the one leading to the lowest coherence among several tested
candidates.

Fig. 1 shows the mutual coherence $\mu_{\text{max}}(\bm{\Psi})$ versus the
iteration index $t$ results for different values of $P$. It can be seen that
the proposed GD-based method has the best performance in all of the considered
scenarios. This shows the effectiveness of the proposed shrinking operator
(22) on the projection matrix design[5] [5] [5] It is worth mentioning that
the performance gain of the proposed GD-based method over the other benchmark
methods holds true even in the unconstrained projection matrix design, i.e.,
when we ignore the CM constraints. The simulation results are not included
here due to space limitations., when compared to the GD-based method of [14]
and the effectiveness of embedding the unit-norm constraints into the
objective function, when compared to the GD-based method of [6]. To
investigate the impact of the projection matrix deign on the DoA estimation
accuracy, Fig. 2 shows the mean-squared errors (MSE) versus the signal-to-
noise-ratio (SNR) results, which are defined as
$\text{MSE}(\bm{\mathring{\nu}})=\mathbb{E}\\{\sum_{k}|\nu_{k}-\mathring{\nu}_{k}|^{2}\\}$
and $\text{SNR}=\mathbb{E}\\{{\|{\bf Y}-{\bf
N}\|_{\text{F}}^{2}}\\}/\mathbb{E}\\{{\|{\bf N}\|^{2}_{\text{F}}}\\}$ (in dB).
Here, we assume that the OMP method [8] uses the dictionary matrix
${\bf\mathring{A}}\in\mathbb{C}^{M\times P}$ to approximate
$\bm{\nu}=[\nu_{1},\dots,\nu_{K}]$ with $\nu_{k}\in[0,2\pi],\forall k$, to the
best on-grid points
$\bm{\mathring{\nu}}=[\mathring{\nu}_{1},\dots,\mathring{\nu}_{K}]$. The
figure clearly shows that the accuracy of the OMP can be significantly
improved by designing a projection matrix with a low mutual coherence
$\mu_{\max}\\{{\bf\Psi}\\}$, especially with smaller number of time samples,
$L$, and dictionary grid points, $P$.

Figure 1: Mutual coherence $\mu_{\text{max}}(\bm{\Psi})$ versus iteration
index $t$. Figure 2: MSE versus SNR [$K=5$, $N=16$, $M=64$].

Example 2: Evaluation of Algorithm 2. Here, we show simulation results
comparing our proposed GOMP method with two high-resolution DoA estimation
methods: CS-based Newtonized OMP (NOMP) [9] and DFT beamspace ESPRIT (BS-
ESPRIT) [3]. To comply with the BS-ESPRIT requirements as discussed in [3,
Lemma 1], we assume that $\nu_{k}\in[0,\nu_{\max}],\forall k$, where
${\nu_{\max}}=\frac{2\pi(N-1)}{M}$. For CS-based OMP, NOMP, and GOMP methods,
we assume that the proposed GD-based method is used to design the projection
matrix ${\bf\Phi}$ for a given dictionary matrix that is formed as
${\bf\mathring{A}}(\bm{\mathring{\nu}})=[{\bf
a}(\mathring{\nu}_{1}),\dots,{\bf
a}(\mathring{\nu}_{P})]\in\mathbb{C}^{M\times P}$ with
$\mathring{\nu}_{p}=\frac{\nu_{\max}(p-1)}{P},p\in\\{1,\dots,P\\}$. Note that
in this case,
${\bf\mathring{A}}(\bm{\mathring{\nu}}){\bf\mathring{A}}(\bm{\mathring{\nu}})^{\mathsf{H}}\neq
P{\bf I}_{M}$. Finally, for NOMP and GOMP, we set the maximum number of
iterations as $I_{\text{max}}=10$ and $J_{\text{max}}=5$.

Fig. 3 shows MSE($\bm{\hat{\nu}}$) and $\text{MSE}(\bm{\mathring{\nu}})$
versus SNR results for various system settings, where $\bm{\hat{\nu}}$ is the
high-resolution estimated parameters vector and MSE($\bm{\hat{\nu}}$) is
defined similarly to $\text{MSE}(\bm{\mathring{\nu}})$. From the figure, we
can see that the CS-based methods have better DoA estimation accuracy than the
subspace-based BS-ESPRIT method, especially in the low SNR regime. Note that,
BS-ESPRIT requires a sufficient number of time samples, i.e., $L\geq K$ to
have an accurate DoA estimation. On the other hand, this condition is not
required by the CS-based methods, although in the expense of higher
computational complexity. Moreover, Fig. 3 shows that GOMP, in average, has a
similar estimation accuracy compared to NOMP. Nonetheless, in our simulations,
we noted that in some realizations, GOMP has higher accuracy than NOMP and
vice-versa on the other realizations.

Figure 3: MSE versus SNR assuming $N=16$, $M=64$, $P=64$.

##  VI Conclusions

In this paper, we have proposed the GOMP method for high-resolution DoA
estimation and a new GD-based method for projection/sensing matrix design
while taking into account the constant modulus constraints imposed by the
projection matrix hardware components. Simulation results are provided showing
the effectiveness of the proposed methods.

## References

  * [1] P. Stoica and K. Sharman, “Maximum likelihood methods for direction-of-arrival estimation,” _IEEE Transactions on Acoustics, Speech, and Signal Processing_ , vol. 38, no. 7, pp. 1132–1143, 1990.
  * [2] P. Gupta and S. Kar, “MUSIC and improved MUSIC algorithm to estimate direction of arrival,” in _Proc. International Conference on Communications and Signal Processing (ICCSP)_ , 2015, pp. 0757–0761.
  * [3] J. Zhang and M. Haardt, “Channel estimation and training design for hybrid multi-carrier mmwave massive MIMO systems: The beamspace ESPRIT approach,” in _Proc. 25th European Signal Processing Conference (EUSIPCO)_ , Sep. 2017, pp. 385–389.
  * [4] D. L. Donoho, “Compressed sensing,” _IEEE Transactions on Information Theory_ , vol. 52, no. 4, pp. 1289–1306, Apr. 2006.
  * [5] H. Reboredo, F. Renna, R. Calderbank, and M. R. D. Rodrigues, “Bounds on the number of measurements for reliable compressive classification,” _IEEE Trans. Signal Process._ , vol. 64, no. 22, pp. 5778–5793, Aug. 2016\.
  * [6] K. Ardah, M. Haardt, T. Liu, F. Matter, M. Pesavento, and M. E. Pfetsch, “Recovery under side constraints,” in _Compressed Sensing in Information Processing_ , G. Kutyniok, H. Rauhut, and R. J. Kunsch, Eds. Switzerland: Springer Nature, 2021, published in the book series Applied and Numerical Harmonic Analysis.
  * [7] H. Boche, R. Calderbank, G. Kutyniok, and J. Vybíral, _A Survey of Compressed Sensing_. Cham: Springer International Publishing, 2015, pp. 1–39.
  * [8] B. L. Sturm and M. G. Christensen, “Comparison of orthogonal matching pursuit implementations,” in _Proc. of the 20th European Signal Processing Conference (EUSIPCO)_ , Aug. 2012, pp. 220–224.
  * [9] B. Mamandipoor, D. Ramasamy, and U. Madhow, “Newtonized orthogonal matching pursuit: Frequency estimation over the continuum,” _IEEE Transactions on Signal Processing_ , vol. 64, no. 19, pp. 5066–5081, Jun. 2016\.
  * [10] K. Ardah, A. L. F. de Almeida, and M. Haardt, “A gridless CS approach for channel estimation in hybrid massive MIMO systems,” in _Proc. IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , May 2019, pp. 4160–4164.
  * [11] K. Ardah, B. Sokal, A. L. F. de Almeida, and M. Haardt, “Compressed sensing based channel estimation and open-loop training design for hybrid analog-digital massive MIMO systems,” in _Proc. IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , May 2020, pp. 4597–4601.
  * [12] K. Ardah, M. Pesavento, and M. Haardt, “A novel sensing matrix design for compressed sensing via mutual coherence minimization,” in _Proc. IEEE 8th International Workshop on Computational Advances in Multi-Sensor Adaptive Processing (CAMSAP)_ , Dec. 2019, pp. 66–70.
  * [13] V. Abolghasemi, S. Ferdowsi, B. Makkiabadi, and S. Sanei, “On optimization of the measurement matrix for compressive sensing,” in _Proc. 18th European Signal Processing Conference_ , Aug. 2010, pp. 427–431.
  * [14] H. Bai, C. Hong, and X. Li, “Construction of unit-norm tight frame based preconditioner for sparse coding,” in _Proc. IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , 2021, pp. 5420–5424.
  * [15] C. Lu, H. Li, and Z. Lin, “Optimized projections for compressed sensing via direct mutual coherence minimization,” _Signal Processing_ , vol. 151, pp. 45 – 55, 2018.
  * [16] L. Yu, G. Li, and L. Chang, “Optimizing projection matrix for compressed sensing systems,” in _Proc. 8th International Conference on Information, Communications Signal Processing (ICICS)_ , 2011, pp. 1–5.
  * [17] K. Ardah, G. Fodor, Y. C. B. Silva, W. C. Freitas, and F. R. P. Cavalcanti, “A unifying design of hybrid beamforming architectures employing phase shifters or switches,” _IEEE Trans. Veh. Technol._ , vol. 67, no. 11, pp. 11 243–11 247, Nov. 2018.