# DOA Estimation in Nonuniform Sensor Noise

Majdoddin Esfandiari, and Sergiy A. Vorobyov M. Esfandiari and S. A. Vorobyov
are with Dept. Signal Processing and Acoustics, Aalto University, PO Box
15400, 00076 Aalto, Finland. Emails<EMAIL_ADDRESS><EMAIL_ADDRESS>Corresponding author is S. A. Vorobyov.This work was
supported in part by Academy of Finland under Research Grant 319822.

###### Abstract

The problem of direction-of-arrival (DOA) estimation in the presence of
nonuniform sensor noise is considered and a novel algorithm is developed. The
algorithm consists of three phases. First, the diagonal nonuniform sensor
noise covariance matrix is estimated using an iterative procedure that
requires only few iterations to obtain an accurate estimate. The asymptotic
variance of one iteration is derived for the proposed noise covariance
estimator. Second, a forward-only rooting-based DOA estimator as well as its
forward-backward averaging extension are developed for DOA estimation. The DOA
estimators take advantage of using second-order statistics of signal subspace
perturbation in constructing a weight matrix of a properly designed
generalized least squares minimization problem. Despite the fact that these
DOA estimators are iterative, only a few iterations are sufficient to reach
accurate results. The asymptotic performance of these DOA estimators is also
investigated. Third, a newly designed DOA selection strategy with reasonable
computational cost is developed to select $L$ actual sources out of $2L$
candidates generated at the second phase. Numerical simulations are conducted
in order to establish the considerable superiority of the proposed algorithm
compared to the the existing state-of-the-art methods in challenging scenarios
in both cases of uniform and nonuniform sensor noise.

###### Index Terms:

DOA estimation, subspace method, nonuniform noise, generalized least squares
(GLS), small sample size.

##  I Introduction

DIRECTION-of-arrival (DOA) estimation is an active field of research for
decades due to multiple traditional and new important applications and due to
significance of source localization in many practical scenarios. The notable
current applications of interest are, for example, wireless communication,
automotive radar, and sonar where DOA estimation is an essential task [1, 2,
3, 4, 5, 6, 7, 8, 9].

In the context of DOA estimation, several assumptions can be regarded
concerning the structure of the second-order statistics of the observation
noise. Most common assumptions are that the noise is uniform white and that
the noise is nonuniform white. Moreover, a spatially block-correlated noise
assumption may be more accurate in some applications [10].

For the uniform white noise assumption, i.e., when the noise powers are
identical across all array sensors, many well-known subspace-based methods
such as MUSIC [11, 12, 13], root-MUSIC [14, 15, 16], and ESPRIT [17, 18] have
been developed, and are based on decomposing the signal or its sample
covariance matrix (SCM) into two disjoint subspaces of noise and signal. The
popularity of the subspace-based methods over, for example, near optimal
maximum likelihood (ML) estimators [19, 20, 21], comes from the fact that they
achieve good estimation accuracy along with affordable computational
complexity in contrast to the prohibitive computational cost of implementing
ML estimators. In recent years, several works have proposed competitive
algorithms to fill the performance gap between subspace-based methods and ML
estimators including root-swap root-MUSIC [22], enhanced principal-singular-
vector utilization for modal analysis (EPUMA) [23], standard ESPRIT using
generalized least squares (SE GLS) and unitary ESPRIT using generalized least
squares (UE GLS) [24], partial relaxation (PR)-based approaches [25, 26],
root-clustering algorithm and root-certificate algorithm [27]. The aim of the
aforementioned methods is to achieve an adequate performance in challenging
scenarios like scarcity of available data samples, low signal-to-noise ratio
(SNR), and/or presence of some correlated or even coherent sources. In [22], a
new concept called root swap has been introduced and recognized as the main
reason behind the collapse of the root-MUSIC algorithm in finding the correct
roots in the challenging environments. To remedy this phenomena, the root-swap
root-MUSIC was proposed which identifies the correct roots associated with the
sources’ DOAs by exploiting the deterministic ML (DML) or stochastic ML (SML)
[28] objective functions instead of deciding based on the closeness of the
absolute values of the roots to unity. In [27], a new criterion for root
selection based on the algebraic structure of the noise subspace has been
proposed. Using this criterion and the relationship between the source
localization problem and the problem of computing the approximate greatest
common divisor (GCD) for polynomials, algorithms that learn the number of
sources and estimate their locations have been proposed. In [23], the authors
have developed a method called EPUMA via solving a particular generalized
least squares (GLS) problem by taking into account the second-order statistics
of the estimated signal subspace, and also producing more DOA candidates than
the number of sources and then selecting final DOAs using DML or SML cost
functions. Similar to EPUMA, the authors of [24] have extended ESPRIT and
unitary ESPRIT to their generalized versions by exploiting the signal subspace
perturbation as the source of error in the shift invariance equation (SIE).
The aforementioned error has been minimized using the GLS. Adopting the main
idea of [24], enhanced ESPRIT-based methods have been developed in [29], which
first produces $2L$ DOA candidates for $L$ sources and then selects the final
$L$ DOAs. The PR framework has been introduced in [25] to relax the manifold
structure of sources partially, resulting in four methods by applying the PR
concept to four previously presented in the literature DOA estimation
techniques. Note that most of the subspace based methods can be extended by
means of using forward-backward averaging (FBA) [30], and forward-backward
spatial smoothing (FBSS) [31] techniques for the cases of correlated or
coherent signals.

The uniform noise assumption can be, however, violated in some practical
scenarios, which appear increasingly more often, and the case of nonuniform
white noise has drawn considerable attention in the last two decades. The
deterministic ML estimator along with Cramer-Rao bound (CRB) for both
deterministic and stochastic source models have been proposed in [32], while
[33] have developed the stochastic ML algorithm for the case of nonuniform
noise. A simple method has been devised in [34] possessing less computational
cost than ML estimators as well as improving the precision of the DOA
estimates. In [35], the authors have developed two iterative methods referred
to as iterative ML subspace estimation (IMLSE) and iterative least squares
subspace estimation (ILSSE), which estimate the signal subspace and noise
covariance matrices based on the ML and least squares (LS) criteria,
respectively. Moreover, it has been shown in [36] that the signal and noise
subspaces are separable via applying the eigendecomposition (ED) of the so-
called reduced covariance matrix for the case of uncorrelated sources, while a
rank minimization approach has been suggested for the correlated sources case.
The performance degradation caused by correlated and/or coherent sources can
also be mitigated via using spatial smoothing [37], and covariance matrix
differencing [38]. Moreover, bearing in mind the fact that the desired spatial
directions can be modeled via sparse representation, sparse signal
reconstruction (SSR) based methods have been developed [39, 40, 41]. Recently,
a method referred to as non-iterative subspace-based (NISB) [42] method has
been proposed where the signal subspace and noise covariance matrix are
identified by exploiting a two steps approach with the first step consisting
of using the ED of the reduced covariance matrix [36] for estimating noise
covariance matrix judiciously, followed by applying the generalized
eigendecomposition (GED) of the matrix pair of the SCM and estimated noise
covariance matrix from the first step.

In this paper, a novel algorithm that consists of three phases is devised for
addressing the problem of DOA estimation in the presence of nonuniform white
noise. The proposed algorithm remarkably provides reliable estimates in
scenarios with small snapshot size and/or relatively low SNR and/or closely
located sources. In the first phase, noise covariance matrix is estimated by
iterations via obtaining the noise subspace using the GED of two matrices,
followed by updating noise covariance matrix using LS. Only few iterations are
required to achieve an adequate estimate of the noise covariance matrix.
Furthermore, the proposed noise covariance matrix estimator is applicable to
sensor arrays with arbitrary geometry as well as any combination of correlated
and uncorrelated sources. In addition, the asymptotic variance of the
aforementioned estimator is derived. Using the noise covariance matrix
estimate, in the second phase, the nonuniform noise DOA estimation problem is
converted into the uniform noise DOA estimation problem by means of pre-
whitening. Then a rooting-based approach that takes into account the signal
subspace perturbation is used. It also utilizes the GLS for estimating the
coefficients of the desired polynomial. To deal with possibly correlated
sources, we also derive the FBA extension of the core DOA estimator. In
addition, the asymptotic variance is derived for the DOA estimator under the
high SNR assumption. In the last phase, a double number of DOA candidates is
generated, and a DOA selection strategy is developed and used to pick the
final DOAs. It is demonstrated by conducting extensive numerical simulations
for various challenging setups that the proposed DOA estimation algorithm
provides superior estimation accuracy both for the uniform and nonuniform
noise cases compared to the existing state-of-the-art techniques.

The rest of the paper is organized as follows. The signal model and problem
formulation are given in Section II. In Section III, an iterative method for
estimating diagonal nonuniform sensor noise covariance matrix is devised. The
asymptotic variance of one iteration is also given for the proposed noise
covariance estimator. The forward-only subspace-based DOA method as well as
its FBA extension are developed for the case of nonuniform noise in Section
IV. Moreover, the asymptotic performance of the proposed DOA estimator is
studied and a new DOA selection strategy is also designed to select final DOAs
in the section. Numerical simulations are provided in Section V. The paper is
concluded in Section VI, and some technical derivations are given in
Supplementary materials.

Notation: Upper-case and lower-case bold-face letters denote matrices and
vectors, respectively, while scalars are denoted by lower-case letters. The
expectation, transpose, conjugate, and Hermitian transpose are denoted by
$\mathbb{E}\\{\cdot\\}$, $\\{\cdot\\}^{T}$, $\\{\cdot\\}^{*}$, and
$\\{\cdot\\}^{H}$, respectively, while $\|\cdot\|_{2}$, $\|\cdot\|_{\rm F}$,
and $|\cdot|$ denote the Euclidean norm of a vector, the Frobenius norm of a
matrix, and the absolute value of a scalar. If the argument is a set,
$|\cdot|$ denotes the set cardinality. The Kronecker product is denoted by
$\otimes$ and $\mathrm{trace}\\{\cdot\\}$ stands for the trace of a square
matrix. The $n\times n$ identity and exchange matrices are denoted by
$\mathbf{I}_{n}$ and $\mathbf{J}_{n}$, respectively. The $n\times m$ matrix
with all elements equal zero and $n\times 1$ zero vector are denoted as
$\mathbf{0}_{m\times n}$ and $\mathbf{0}_{n}$, respectively. The $i$th entry
of the vector $\boldsymbol{\pi}$ is denoted by $[\boldsymbol{\pi}]_{i}$. The
$i$th row and $i$th column of the matrix $\boldsymbol{\Pi}$ are denoted by
$[\boldsymbol{\Pi}]_{i,:}$ and $[\boldsymbol{\Pi}]_{:,i}$, respectively, while
the entry in the interaction of the $i$th row and $j$th column is denoted as
$[\boldsymbol{\Pi}]_{i,j}$. The operator $\mathrm{vec}\\{\cdot\\}$ stacks the
columns of a matrix into a long vector. The operator
$\mathrm{diag}\\{\boldsymbol{\pi}\\}$ generates a diagonal matrix by plugging
the entries of the vector $\boldsymbol{\pi}$ into its main diagonal, while the
operator $\mathcal{D}\\{\boldsymbol{\Pi}\\}$ creates a diagonal matrix by
preserving the main diagonal of the matrix $\boldsymbol{\Pi}$ and setting all
other entries to zero. The operator $\mathrm{DFT}\\{\boldsymbol{\pi}\\}$
stands for the discrete Fourier transform (DFT) of the vector
$\boldsymbol{\pi}$, while $\mathfrak{R}\\{\cdot\\}$ returns the real part of
the bracketed argument.

##  II SIGNAL MODEL

Consider a uniform linear array (ULA) composed of $M$ omni-directional sensors
receiving $L$ ($L<M$) independent narrowband signals radiated by $L$ sources.
It is assumed that the sources are located in the far-field and have
distinguished directions, denoted as $\theta_{l}$, $l=1,\cdots,L$. Then, the
signal observed by the sensor array at the time instant $t$ is written as

$\displaystyle\mathbf{x}(t)=\mathbf{A}(\boldsymbol{\theta})\mathbf{s}(t)+\mathbf{n}(t)\
$ (1)

where $s(t)\triangleq[s_{1}(t)\cdots s_{L}(t)]^{T}\in\mathbb{C}^{L}$ denotes
the source signals, $\mathbf{n(t)}\in\mathbb{C}^{M}$ is the sensor noise
vector, the source DOAs are stacked in the vector
$\boldsymbol{\theta}\triangleq[\theta_{1}\cdots\theta_{L}]^{T}$,
$\mathbf{A}(\boldsymbol{\theta})\triangleq[\mathbf{a}(\theta_{1})\cdots\mathbf{a}(\theta_{L})]$
denotes the array manifold whose $l$th column is the steering vector
$\mathbf{a}(\theta_{l})=[1\ e^{-j2\pi\mathrm{sin}(\theta_{l})d/\lambda}\
\cdots\
e^{-j2\pi(M-1)\mathrm{sin}(\theta_{l})d/\lambda}]^{T}\in\mathbb{C}^{M}$
associated with $l$th DOA. Here $\lambda$ denotes the carrier wavelength and
$d=\lambda/2$. For notation simplicity, $\mathbf{A}$ is used instead of
$\mathbf{A}(\boldsymbol{\theta})$ hereafter, unless the argument of
$\mathbf{A}$ is different from $\boldsymbol{\theta}$.

The array covariance matrix can be written as

$\displaystyle\mathbf{R}\triangleq
E\\{\mathbf{x}(t)\mathbf{x}^{H}(t)\\}=\mathbf{A}\mathbf{P}\mathbf{A}^{H}+\mathbf{Q}\
$ (2)

where $\mathbf{P}\in\mathbb{C}^{L\times L}$ and
$\mathbf{Q}\in\mathbb{R}^{M\times M}$, respectively, stand for the signal and
noise covariance matrices defined as

$\displaystyle\mathbf{P}\triangleq E\\{\mathbf{s}(t)\mathbf{s}^{H}(t)\\},\ \
\mathbf{Q}\triangleq E\\{\mathbf{n}(t)\mathbf{n}^{H}(t)\\}.\ $ (3)

Considering the case of nonuniform sensor noise that is spatially and
temporally uncorrelated and zero-mean Gaussian, the noise covariance matrix is
the following diagonal matrix

$\displaystyle\mathbf{Q}=\mathrm{diag}\\{[\sigma_{1}^{2},\cdots,\sigma_{M}^{2}]\\}\
$ (4)

where $\sigma_{m}^{2}$, $m=1,\cdots,M$ are the noise variances, which are not
necessarily identical, i.e, $\sigma_{i}^{2}\neq\sigma_{j}^{2}$ for $i\neq j$.
When $\sigma_{1}^{2}=\sigma_{2}^{2}=\cdots=\sigma_{M}^{2}=\sigma^{2}$, the
noise covariance matrix is just a scaled identity matrix
$\mathbf{Q}=\sigma^{2}\mathbf{I}_{M}$, i.e., the sensor noise is uniform. The
latter case has been investigated to a great extent in the literature, while
the former has been drawn more attention in recent years. The critical part of
most of the methods designed for the nonuniform noise case is to estimate
$\mathbf{Q}$ efficiently, and the source DOA estimates are dependent on how
precise the estimate of $\mathbf{Q}$ is. Moreover, when the number of
available snapshots is small and/or some of the DOAs are closely located to
each other and/or SNR is relatively low, the impact of estimated $\mathbf{Q}$
is more notable.

Because $\mathbf{R}$ is unknown in practice, the SCM is considered, and it is
given by

$\displaystyle\hat{\mathbf{R}}\triangleq\frac{1}{N}\displaystyle\sum_{t=1}^{N}\mathbf{x}(t)\mathbf{x}^{H}(t)=\frac{1}{N}\mathbf{X}\mathbf{X}^{H}.\
\ $ (5)

Here, the matrix signal notation is also used

$\displaystyle\mathbf{X}=\mathbf{A}\mathbf{S}+\mathbf{N}\ \ $ (6)

with $\mathbf{X}\triangleq[\mathbf{x}(1)\cdots\mathbf{x}(N)]$,
$\mathbf{S}\triangleq~{}[\mathbf{s}(1)\cdots\mathbf{s}(N)]$,
$\mathbf{N}\triangleq[\mathbf{n}(1)\cdots\mathbf{n}(N)]$, and $N$ being the
number of snapshots.

##  III NOISE COVARIANCE MATRIX ESTIMATION

It is desirable that a method estimating $\mathbf{Q}$ would not require any
knowledge of the true DOAs, while providing an acceptable accuracy with an
affordable computational cost. Moreover, such method should be robust in
extreme scenarios. Some examples of extreme scenarios are small sample size
and presence of closely located sources. To devise such a method, we begin by
multiplying both sides of (2) by $\mathbf{U}\in\mathbb{C}^{M\times(M-L)}$
which satisfies the following condition

$\displaystyle\mathbf{A}^{H}\mathbf{U}=\boldsymbol{0}_{L\times(M-L)}\ $ (7)

where the constraint $\mathbf{U}^{H}\mathbf{U}=\mathbf{I}_{M-L}$ is also
imposed for avoiding ambiguities in defining $\mathbf{U}$. It is clear then
that the columns of $\mathbf{U}$ constitute an orthonormal basis for the noise
subspace. As $\mathbf{A}$ is unknown, finding $\mathbf{U}$ for the case of
nonuniform noise is not as simple as for the uniform noise case, when
$\mathbf{U}$ is obtained by just calculating the eigenvectors of
$\hat{\mathbf{R}}$. Multiplying (2) by $\mathbf{U}$, and exploiting (7), it
can be written that [42]

$\displaystyle\hat{\mathbf{R}}\mathbf{U}=\mathbf{Q}\mathbf{U}\ $ (8)

where $\mathbf{R}$ is replaced by the consistent estimate $\hat{\mathbf{R}}$.
Then, it can be observed from (8) that for a given estimate of $\mathbf{Q}$,
the columns of the best estimate of $\mathbf{U}$, denoted as
$\hat{\mathbf{U}}$, are given by the $M-L$ eigenvectors associated with the
$M-L$ smallest eigenvalues obtained after applying the GED to the pair of
matrices $\hat{\mathbf{R}}$ and $\mathbf{Q}$.

Furthermore, according to (8), the LS approach can be employed to estimate
$\mathbf{Q}$. In doing so, we formulate the following LS minimization problem
with respect to $\mathbf{Q}$

$\displaystyle\hat{\mathbf{Q}}=\mathrm{arg}\displaystyle\min_{\mathbf{Q}}\|(\hat{\mathbf{R}}-\mathbf{Q})\hat{\mathbf{U}}\|^{2}_{\rm
F}$ (9)

where $\mathbf{U}$ is replaced by $\hat{\mathbf{U}}$. Problem (9) needs to be
solved subject to the constraint on $\mathbf{Q}$, that is, $\mathbf{Q}$ is a
diagonal matrix. The cost function of (9) can be expressed as

$\displaystyle f(\mathbf{Q})$
$\displaystyle\triangleq\left\|(\hat{\mathbf{R}}-\mathbf{Q})\hat{\mathbf{U}}\right\|^{2}_{\rm
F}=\mathrm{trace}\left\\{\left((\hat{\mathbf{R}}-\mathbf{Q})\hat{\mathbf{U}}\right)\left((\hat{\mathbf{R}}-\mathbf{Q})\hat{\mathbf{U}}\right)^{H}\right\\}$
$\displaystyle=\mathrm{trace}\left\\{\hat{\mathbf{U}}\hat{\mathbf{U}}^{H}\hat{\mathbf{R}}^{2}\right\\}-\mathrm{trace}\left\\{\hat{\mathbf{R}}\hat{\mathbf{U}}\hat{\mathbf{U}}^{H}\mathbf{Q}\right\\}-\mathrm{trace}\left\\{\hat{\mathbf{U}}\hat{\mathbf{U}}^{H}\hat{\mathbf{R}}\mathbf{Q}\right\\}$
$\displaystyle+\mathrm{trace}\left\\{\hat{\mathbf{U}}\hat{\mathbf{U}}^{H}\mathbf{Q}^{2}\right\\}$
(10)

where the properties $\|\mathbf{X}\|^{2}_{\rm
F}=\mathrm{trace}\left\\{\mathbf{X}\mathbf{X}^{H}\right\\}$,
$\mathrm{trace}\left\\{\mathbf{X}\mathbf{Y}\right\\}=\mathrm{trace}\left\\{\mathbf{Y}\mathbf{X}\right\\}$,
$\hat{\mathbf{R}}=\hat{\mathbf{R}}^{H}$, and $\mathbf{Q}=\mathbf{Q}^{H}$ are
employed. As elaborated in Appendix A, the partial derivative of (10) with
respect to $\mathbf{Q}$ is given as

$\displaystyle\frac{\partial
f(\mathbf{Q})}{\partial\mathbf{Q}}=2\mathcal{D}\left\\{\hat{\mathbf{U}}\hat{\mathbf{U}}^{H}\right\\}\mathbf{Q}-\mathcal{D}\left\\{\hat{\mathbf{R}}\hat{\mathbf{U}}\hat{\mathbf{U}}^{H}+\hat{\mathbf{U}}\hat{\mathbf{U}}^{H}\hat{\mathbf{R}}\right\\}.$
(11)

Therefore, by equating (11) to zero, that is,

$\displaystyle
2\mathcal{D}\left\\{\hat{\mathbf{U}}\hat{\mathbf{U}}^{H}\right\\}\mathbf{Q}-\mathcal{D}\left\\{\hat{\mathbf{R}}\hat{\mathbf{U}}\hat{\mathbf{U}}^{H}+\hat{\mathbf{U}}\hat{\mathbf{U}}^{H}\hat{\mathbf{R}}\right\\}=0$
(12)

we find the optimal estimate of $\mathbf{Q}$ as

$\displaystyle\hat{\mathbf{Q}}=\frac{1}{2}\mathcal{D}\left\\{\hat{\mathbf{R}}\hat{\mathbf{U}}\hat{\mathbf{U}}^{H}+\hat{\mathbf{U}}\hat{\mathbf{U}}^{H}\hat{\mathbf{R}}\right\\}\mathcal{D}\left\\{\hat{\mathbf{U}}\hat{\mathbf{U}}^{H}\right\\}^{-1}.$
(13)

According to (8) and (13), which determine the dependencies between
$\hat{\mathbf{Q}}$ and $\hat{\mathbf{U}}$, it is natural to select an
iterative scheme for estimating $\hat{\mathbf{Q}}$ and $\hat{\mathbf{U}}$ in
alternative manner. It starts with properly initializing $\hat{\mathbf{Q}}$,
denoted as $\hat{\mathbf{Q}}^{(0)}$, followed by estimating
$\hat{\mathbf{U}}^{(0)}$ as the $M-L$ generalized eigenvectors of the pair
$\\{\hat{\mathbf{R}},\;\hat{\mathbf{Q}}^{(0)}\\}$ associated with the $M-L$
smallest eigenvalues. Then, $\hat{\mathbf{Q}}^{(1)}$ is obtained through (13)
after replacing $\hat{\mathbf{U}}$ with $\hat{\mathbf{U}}^{(0)}$. This
procedure continues until a predefined stopping criterion is satisfied.
Although any randomly generated diagonal matrix with positive diagonal entries
can be used as $\hat{\mathbf{Q}}^{(0)}$, we suggest to initialize it as
$\hat{\mathbf{Q}}^{(0)}=\mathcal{D}\\{\hat{\mathbf{R}}\\}$.

Algorithm 1: Noise Covariance Matrix Estimation
---
1: Compute $\hat{\mathbf{R}}=1/N\sum_{t=1}^{N}\mathbf{x}(t)\mathbf{x}^{H}(t)$.
2: Set $i=0$ and $\hat{\mathbf{Q}}^{(0)}=\mathcal{D}\\{\hat{\mathbf{R}}\\}$.
In addition, set the maximum number of iterations $i_{max}=5$.
while $i\leq i_{max}$
3: Carry out the GED of the pair of matrices
$\\{\hat{\mathbf{R}},\;\hat{\mathbf{Q}}^{(i)})\\}$ to obtain
$\hat{\mathbf{U}}^{(i)}$ as the $M-L$ eigenvectors corresponding to the
$M-L$ smallest eigenvalues.
4: Calculate $\hat{\mathbf{Q}}^{(i+1)}$ using (13).
5: set $i=i+1$.
end

The steps of the proposed method for noise covariance estimation are outlined
in Algorithm 1. Additionally, to demonstrate the quality of the proposed noise
covariance estimation method in (13), its asymptotic mean square error (MSE)
is found in the following proposition.

Proposition 1: The asymptotic variance of estimating each diagonal element of
$\hat{\mathbf{Q}}$ in (13), given a particular $\hat{\mathbf{U}}$, is

$\displaystyle\mathbb{E}\left\\{\left(\Delta\sigma_{m}^{2}\right)^{2}\right\\}$
$\displaystyle=\left(\frac{[\mathbf{R}]_{mm}}{2N\tau_{m}^{2}}\right)\mathfrak{R}\left\\{\mathbf{v}_{m}^{H}\mathbf{R}\left(\mathbf{v}_{m}+N\mathbf{v}_{m}^{*}\right)\right\\}$
(14)

where $\tau_{m}\triangleq[\hat{\mathbf{U}}\hat{\mathbf{U}}^{H}]_{mm}$ and
$\mathbf{v}_{m}\triangleq[\hat{\mathbf{U}}\hat{\mathbf{U}}^{H}]_{:,m}$.

Proof: See Appendix B. $\blacksquare$

Remark 1: The estimate of $\hat{\mathbf{Q}}$ in (13) is an alternative
representation of the power domain (PD) method [34], that is,

$\displaystyle\sigma_{m}^{2}=\frac{\left(\mathbf{d}_{m}^{T}\mathbf{P}_{\mathbf{A}}\hat{\mathbf{r}}_{m}+\hat{\mathbf{r}}_{m}^{H}\mathbf{P}_{\mathbf{A}}\mathbf{d}_{m}\right)}{2\mathbf{d}_{m}^{T}\mathbf{P}_{\mathbf{A}}\mathbf{d}_{m}},\
\ \ m=1,\cdots,M$ (15)

where $\mathbf{d}_{m}\in\mathbb{R}^{M}$ denotes a vector with one on the $m$th
position and zero elsewhere, $\hat{\mathbf{r}}_{m}$ is the $m$th column of
$\hat{\mathbf{R}}$, and $\mathbf{P}_{\mathbf{A}}\in\mathbb{C}^{M\times M}$ is
the orthogonal projection matrix of the signal subspace, i.e.,
$\mathbf{P}_{\mathbf{A}}=\mathbf{I}_{M}-\mathbf{A}(\mathbf{A}^{H}\mathbf{A})^{-1}\mathbf{A}^{H}$.
The main difference between (13) and (15) is the use of
$\hat{\mathbf{U}}\hat{\mathbf{U}}^{H}$ as a proper estimate of $\mathbf{P}$
instead of conducting a multidimensional search for finding an ML estimate as
in [34]. It is well-known that performing such multidimensional search is a
very computationally demanding task. Therefore, the proposed method is
advantageous compared to PD technique in the sense that it requires
significantly less computations.

Remark 2: As the actual $\mathbf{Q}$ is unknown, it is more accurate to modify
(8) for the $i$th iteration as

$\displaystyle\hat{\mathbf{R}}\hat{\mathbf{U}}^{(i)}=\hat{\mathbf{Q}}^{(i)}\hat{\mathbf{U}}^{(i)}\hat{\boldsymbol{\Lambda}}^{(i)}\
$ (16)

where $\hat{\boldsymbol{\Lambda}}^{(i)}\in\mathbb{R}^{(M-L)\times(M-L)}$ is a
diagonal matrix containing the generalized eigenvalues corresponding to the
columns of $\hat{\mathbf{U}}^{(i)}$. As a consequence, the relationship
between $\hat{\mathbf{Q}}^{(i+1)}$ and $\hat{\mathbf{Q}}^{(i)}$ can be
obtained via exploiting (13) and (16) as follows

$\displaystyle\hat{\mathbf{Q}}^{(i\\!+\\!1)}=$
$\displaystyle\frac{\mathcal{D}\\!\left\\{\hat{\mathbf{R}}\hat{\mathbf{U}}^{(i)}\hat{\mathbf{U}}^{(i)H}\\!+\\!\hat{\mathbf{U}}^{(i)}\hat{\mathbf{U}}^{(i)H}\hat{\mathbf{R}}\right\\}\\!\mathcal{D}\\!\left\\{\hat{\mathbf{U}}^{(i)}\hat{\mathbf{U}}^{(i)H}\right\\}^{-1}}{2}$
$\displaystyle=$
$\displaystyle\mathcal{D}\left\\{\hat{\mathbf{U}}^{(i)}\hat{\boldsymbol{\Lambda}}^{(i)}\hat{\mathbf{U}}^{(i)H}\right\\}\mathcal{D}\left\\{\hat{\mathbf{U}}^{(i)}\hat{\mathbf{U}}^{(i)H}\right\\}^{-1}\hat{\mathbf{Q}}^{(i)}.$
(17)

Here the properties $\hat{\mathbf{Q}}^{(i)}=(\hat{\mathbf{Q}}^{(i)})^{H}$ and
$\hat{\boldsymbol{\Lambda}}^{(i)}=(\hat{\boldsymbol{\Lambda}}^{(i)})^{H}$ have
also been used. Hence, based on (17), two consecutive estimates
$\hat{\mathbf{Q}}^{(i)}$ and $\hat{\mathbf{Q}}^{(i+1)}$ are linked to each
other through a diagonal weighting matrix $\mathbf{W}_{Q}^{(i)}$, that is,

$\displaystyle\mathbf{W}_{Q}^{(i)}\triangleq\mathcal{D}\\{\hat{\mathbf{U}}^{(i)}\hat{\boldsymbol{\Lambda}}^{(i)}\hat{\mathbf{U}}^{(i)H}\\}\mathcal{D}\\{\hat{\mathbf{U}}^{(i)}\hat{\mathbf{U}}^{(i)H}\\}^{-1}.$
(18)

Remark 3: The sufficient number of iterations for achieving a precise estimate
of the noise covariance matrix for
$\hat{\mathbf{Q}}^{(0)}=\mathcal{D}\\{\hat{\mathbf{R}}\\}$ is 3–5 as will be
shown in simulation section.

##  IV SUBSPACE-BASED DOA ESTIMATION VIA GLS

###  IV-A Forward-only Algorithm for DOA Estimation

Using the estimate $\hat{\mathbf{Q}}$ obtained based on Algorithm 1, the
received signal can be preprocessed by multiplying (6) by
$\hat{\mathbf{Q}}^{-\frac{1}{2}}$ to enforce the uniform noise. Thus, the
received signal becomes

$\displaystyle\overline{\mathbf{X}}\triangleq\hat{\mathbf{Q}}^{-\frac{1}{2}}\mathbf{X}=\hat{\mathbf{Q}}^{-\frac{1}{2}}\mathbf{A}\mathbf{S}+\hat{\mathbf{Q}}^{-\frac{1}{2}}\mathbf{N}=\hat{\mathbf{Q}}^{-\frac{1}{2}}\mathbf{A}\mathbf{S}+\overline{\mathbf{N}}\
$ (19)

where the columns of $\overline{\mathbf{N}}$ are Gaussian random vectors with
zero mean and the covariance matrix is $\mathbf{I}_{M}$.111In fact, a more
accurate $\hat{\mathbf{Q}}$ leads to the covariance matrix of
$\overline{\mathbf{N}}$ being closer to $\mathbf{I}_{M}$.

The truncated SVD of $\overline{\mathbf{X}}$ is given as

$\displaystyle\overline{\mathbf{X}}=\mathbf{U}_{\rm s}\boldsymbol{\Sigma}_{\rm
s}\mathbf{V}_{\rm s}^{H}$ (20)

where $\mathbf{U}_{\rm s}\in\mathbb{C}^{M\times L}$ and $\mathbf{V}_{\rm
s}\in\mathbb{C}^{N\times L}$ denote respectively the left and right singular
vectors associated with the $L$ principal singular values on the diagonal of
$\boldsymbol{\Sigma}_{\rm s}\in\mathbb{R}^{L\times L}$.

According to (19) and (20), the columns of $\mathbf{U}_{\rm s}$, denoted as
$\mathbf{u}_{p}$ for $p=1,\cdots,L$, and the columns of
$\hat{\mathbf{Q}}^{-\frac{1}{2}}\mathbf{A}$ span the same column space, i.e.,
$\mathrm{span}(\mathbf{U}_{\rm
s})=\mathrm{span}(\hat{\mathbf{Q}}^{-\frac{1}{2}}\mathbf{A})$. In other words,
$\mathbf{U}_{\rm s}$ and $\hat{\mathbf{Q}}^{-\frac{1}{2}}\mathbf{A}$ are
related as

$\displaystyle\mathbf{U}_{\rm
s}=\hat{\mathbf{Q}}^{-\frac{1}{2}}\mathbf{A}\mathbf{G}$ (21)

where $\mathbf{G}\in\mathbb{C}^{L\times L}$ is a non-singular matrix.
Consequently, multiplying (21) by $\hat{\mathbf{Q}}^{\frac{1}{2}}$, we have

$\displaystyle\widetilde{\mathbf{U}}_{\rm s}=\mathbf{A}\mathbf{G}$ (22)

where $\widetilde{\mathbf{U}}_{\rm
s}\triangleq\hat{\mathbf{Q}}^{\frac{1}{2}}\mathbf{U}_{\rm
s}\in\mathbb{C}^{M\times L}$.

Using (22), the $m$th entry of the $p$th column of
$\widetilde{\mathbf{U}}_{\rm s}$, denoted as
$[\widetilde{\mathbf{u}}_{p}]_{m}$, can be written as

$\displaystyle[\widetilde{\mathbf{u}}_{p}]_{m}=$
$\displaystyle[\mathbf{A}]_{m,:}\
\mathbf{g}_{p}=\displaystyle\sum_{l=1}^{L}[\mathbf{g}_{p}]_{l}\ e^{-j2\pi
d\mathrm{sin}(\theta_{l})(m-1)/\lambda},$ $\displaystyle p=1,\cdots,L,\
m=1,\cdots,M\ \ $ (23)

where $\mathbf{g}_{p}\in\mathbb{C}^{L}$ is the $p$th column of $\mathbf{G}$.

Next, the DFT can be applied to each column of $\widetilde{\mathbf{U}}_{\rm
s}$ as, for example, in [43]. Using (23) and the definition of DFT, the $k$th
bin of the DFT of $\widetilde{\mathbf{u}}_{p}$ can be expressed as

$\displaystyle[\bar{\mathbf{u}}_{p}]_{k}=\displaystyle\sum_{l=1}^{L}[\mathbf{g}_{p}]_{l}$
$\displaystyle\ \frac{1-e^{jM\beta_{l}}}{1-e^{-j\frac{2\pi
k}{M}}e^{j\beta_{l}}}=\displaystyle\sum_{l=1}^{L}\frac{\alpha_{l}}{1-\gamma_{l}W_{M}^{k}},$
$\displaystyle\ k=1,\cdots,M\ \ $ (24)

where
$\bar{\mathbf{u}}_{p}\triangleq\mathrm{DFT}\\{\widetilde{\mathbf{u}}_{p}\\}=\mathbf{W}_{D}\widetilde{\mathbf{u}}_{p},\;p=1,\cdots,L$
and

$\displaystyle\mathbf{W}_{D}$ $\displaystyle=\begin{bmatrix}1&1&\cdots&1\\\
1&e^{\frac{-j2\pi}{M}}&\cdots&e^{\frac{-j2\pi(M-1)}{M}}\\\
\vdots&\vdots&\ddots&\vdots\\\
1&e^{\frac{-j2\pi(M-1)}{M}}&\cdots&e^{\frac{-j2\pi(M-1)(M-1)}{M}}\end{bmatrix}.$

In addition,
$\alpha_{l}\triangleq[\mathbf{g}_{p}]_{l}\left(1-e^{jM\beta_{l}}\right)$,
$\gamma_{l}\triangleq e^{j\beta_{l}}$, $\beta_{l}\triangleq-2\pi d\
\mathrm{sin}(\theta_{l})/\lambda$ , and $W_{M}^{k}\triangleq e^{-j\frac{2\pi
k}{M}}$ in (24). Unifying the $L$ rational functions into one, we recast (24)
as

$\displaystyle[\bar{\mathbf{u}}_{p}]_{k}=\frac{\displaystyle\sum_{l=1}^{L}\alpha_{l}\displaystyle\prod_{\begin{subarray}{c}v=1\\\
v\neq
l\end{subarray}}^{L}(1-\gamma_{v}W_{M}^{k})}{\displaystyle\prod_{l=1}^{L}(1-\gamma_{l}W_{M}^{k})}\
\ $ (25)

where the common denominator is the product of the $L$ denominators of each
rational function. Noticing the special structure of (25), the nominator and
denominator can be expanded as two polynomials of degrees $L-1$ and $L$,
respectively, that is,

$\displaystyle\displaystyle\sum_{l=1}^{L}\alpha_{l}\displaystyle\prod_{\begin{subarray}{c}v=1\\\
v\neq
l\end{subarray}}^{L}(1-\gamma_{v}W_{M}^{k})=\displaystyle\sum_{l=1}^{L}b_{pl}(W_{M}^{k})^{l-1}=\overline{\mathbf{w}}_{k}^{\intercal}\mathbf{b}_{p}$
(26)
$\displaystyle\displaystyle\prod_{l=1}^{L}(1-\gamma_{l}W_{M}^{k})=1+\displaystyle\sum_{l=1}^{L}a_{l}(W_{M}^{k})^{l}=1+\mathbf{w}_{k}^{\intercal}\mathbf{a}$
(27)

where $\overline{\mathbf{w}}_{k}\triangleq[1\ W_{M}^{k}\ (W_{M}^{k})^{2}\
\cdots\ (W_{M}^{k})^{L-1}]^{T}$, $\mathbf{w}_{k}\triangleq[W_{M}^{k}\
(W_{M}^{k})^{2}\ (W_{M}^{k})^{3}\ \cdots\ (W_{M}^{k})^{L}]^{T}$,
$\mathbf{b}_{p}\triangleq[b_{p1}\ \cdots\ b_{pL}]^{T}$, and
$\mathbf{a}\triangleq[a_{1}\ \cdots\ a_{L}]^{T}$. Note that neither the
entries of $\mathbf{b}_{p}$’s nor the entries of $\mathbf{a}$ are dependent to
$k$.

It can be seen from (27) that estimating $\mathbf{a}$ is the key for finding
$\gamma_{l}$’s, $l=1,\cdots,L$, since $\gamma_{l}$’s, $l=1,\cdots,L$ are the
roots of the polynomial defined by the entries of $\mathbf{a}$ as

$\displaystyle\gamma^{L}+\displaystyle\sum_{l=1}^{L}[\mathbf{a}]_{l}\
\gamma^{L-l}=0.\ \ $ (28)

Finally, knowing $\gamma_{l}$’s, $\theta_{l}$’s can be extracted using the
relation $\theta_{l}=\mathrm{arcsin}(-\frac{\beta_{l}\lambda}{2\pi d})$ where
$\beta_{l}$ is the phase argument of $\gamma_{l}$ as it is defined above.

Thus, the objective now is to find an estimate of $\mathbf{a}$, denoted as
$\hat{\mathbf{a}}$. Multiplying both sides of (25) by the denominator, and
using (26) and (27), we have

$\displaystyle[\bar{\mathbf{u}}_{p}]_{k}(1+\mathbf{w}_{k}^{\intercal}\mathbf{a})=\overline{\mathbf{w}}_{k}^{\intercal}\mathbf{b}_{p}.\
\ $ (29)

Piling up all the equations that can be generated for $k=1,\cdots,M$ based on
(29), we can write

$\displaystyle\bar{\mathbf{u}}_{p}+\mathrm{diag}(\bar{\mathbf{u}}_{p})\mathbf{W}_{a}\mathbf{a}=\overline{\mathbf{W}}\mathbf{b}_{p},\
\ \ \ p=1,\cdots,L\ \ $ (30)

where
$\mathbf{W}_{a}\triangleq[\mathbf{w}_{1}\mathbf{w}_{2}\cdots\mathbf{w}_{M}]^{T}\in\mathbb{C}^{M\times
L}$ and
$\overline{\mathbf{W}}\triangleq[\overline{\mathbf{w}}_{1}\overline{\mathbf{w}}_{2}\cdots\overline{\mathbf{w}}_{M}]^{T}\in\mathbb{C}^{M\times
L}$.

We suggest to generate two sets of estimates for $\theta_{l}$’s as the most
probable candidates first, and then pick up the best candidates via using a
proper selection criteria. In doing so, we introduce the selection matrix
$\mathbf{Z}_{\mathcal{I}}\in\mathbb{R}^{|\mathcal{I}|\times M}$. Here
$\mathcal{I}$ denotes the set containing the indices of the selected
equations. The matrix $\mathbf{Z}_{\mathcal{I}}$ is used to consider different
sets of (30) for estimating $\mathbf{a}$. Since $\mathbf{Z}_{\mathcal{I}}$ is
a selection matrix, all the entries of the $i$th row of
$\mathbf{Z}_{\mathcal{I}}$ are zeros except one entry whose index is the $i$th
member of $\mathcal{I}$. The only nonzero entry of each row is set to 1. The
proposed method for selecting two sets of indices to serve as $\mathcal{I}$
will be clarified in the sequel.

Given a particular set $\mathcal{I}$, the selected subset of (30) can be
written as

$\displaystyle\mathbf{Z}_{\mathcal{I}}\bar{\mathbf{u}}_{p}+$
$\displaystyle\mathrm{diag}(\mathbf{Z}_{\mathcal{I}}\bar{\mathbf{u}}_{p})\mathbf{Z}_{\mathcal{I}}\mathbf{W}_{a}\mathbf{a}=\mathbf{Z}_{\mathcal{I}}\overline{\mathbf{W}}\mathbf{b}_{p},\quad
p=1,\cdots,L.$ (31)

To avoid estimating $\mathbf{b}_{p}$’s in (31),
$\mathbf{B}\in\mathbb{C}^{|\mathcal{I}|\times(|\mathcal{I}|-L)}$ can be
obtained using SVD such that
$\mathbf{B}^{H}\widetilde{{\mathbf{Z}}}_{\mathcal{I}}=\mathbf{0}_{(|\mathcal{I}|-L)\times
L}$ with
$\widetilde{{\mathbf{Z}}}_{\mathcal{I}}\triangleq\mathbf{Z}_{\mathcal{I}}\overline{\mathbf{W}}\in\mathbb{C}^{|\mathcal{I}|\times
L}$. Thus, multiplying both sides of (31) by $\mathbf{B}^{H}$ and using the
property
$\mathbf{B}^{H}\widetilde{{\mathbf{Z}}}_{\mathcal{I}}=\mathbf{0}_{(|\mathcal{I}|-L)\times
L}$, yields

$\displaystyle\mathbf{B}^{H}(\mathbf{Z}_{\mathcal{I}}\bar{\mathbf{u}}_{p}+$
$\displaystyle\mathrm{diag}(\mathbf{Z}_{\mathcal{I}}\bar{\mathbf{u}}_{p})\mathbf{Z}_{\mathcal{I}}\mathbf{W}_{a}\mathbf{a})=\mathbf{0}_{(|\mathcal{I}|-L)},\;p=1,\cdots,L.\
$ (32)

The impact of $\mathbf{b}_{p}$’s is eliminated in (32) independent of the
value of $p$ because of the definition of
$\widetilde{{\mathbf{Z}}}_{\mathcal{I}}$, which is independent of $p$. This
property enables us to combine the $L$ sets of linear equations generated via
(32). Rearranging the terms in (32), we get the following system of linear
equations $\mathbf{H}_{p}\mathbf{a}=\mathbf{h}_{p},\;p=1,\cdots,L$, where
$\mathbf{H}_{p}\triangleq\mathbf{B}^{H}\mathrm{diag}(\mathbf{Z}_{\mathcal{I}}\bar{\mathbf{u}}_{p})\mathbf{Z}_{\mathcal{I}}\mathbf{W}_{a}\in\mathbb{C}^{(|\mathcal{I}|-L)\times
L}$ and
$\mathbf{h}_{p}\triangleq-\mathbf{B}^{H}\mathbf{Z}_{\mathcal{I}}\bar{\mathbf{u}}_{p}\in\mathbb{C}^{(|\mathcal{I}|-L)}$.
Stacking the $L$ matrices $\mathbf{H}_{p}$ and the $L$ vectors
$\mathbf{h}_{p}$ into a larger matrix $\mathbf{H}$ and a longer vector
$\mathbf{h}$, respectively, we have

$\displaystyle\mathbf{H}\mathbf{a}=\mathbf{h}\ \ $ (33)

where
$\mathbf{H}\triangleq[\mathbf{H}_{1}^{T}\cdots\mathbf{H}_{L}^{T}]^{T}\in\mathbb{C}^{L(|\mathcal{I}|-L)\times
L}$ and
$\mathbf{h}\triangleq[\mathbf{h}_{1}^{T}\cdots\mathbf{h}_{L}^{T}]^{T}\in\mathbb{C}^{L(|\mathcal{I}|-L)}$.

However, only an estimate of $\mathbf{U}_{s}$ (which spans the columns space
of $\hat{\mathbf{Q}}^{-\frac{1}{2}}\mathbf{A}$) can be obtained via the
truncated SVD of $\overline{\mathbf{X}}$ due to the presence of noise. Hence,
a more precise in terms of notation form of (20) is

$\displaystyle\overline{\mathbf{X}}=\hat{\mathbf{U}}_{\rm
s}\hat{\boldsymbol{\Sigma}}_{\rm s}\hat{\mathbf{V}}_{\rm s}^{H}$ (34)

where $\hat{\mathbf{U}}_{\rm
s}=[\hat{\mathbf{u}}_{1}\cdots\hat{\mathbf{u}}_{L}]\in\mathbb{C}^{M\times L}$
is the matrix of $L$ left singular vectors associated with the $L$ largest
singular values on the diagonal of $\hat{\boldsymbol{\Sigma}}_{\rm
s}\in\mathbb{R}^{L\times L}$. Accordingly, $\hat{\widetilde{\mathbf{U}}}_{\rm
s}$ is defined as

$\displaystyle\hat{\widetilde{\mathbf{U}}}_{\rm s}$
$\displaystyle\triangleq\hat{\mathbf{Q}}^{\frac{1}{2}}\hat{\mathbf{U}}_{\rm
s}=\left[(\hat{\mathbf{Q}}^{\frac{1}{2}}\hat{\mathbf{u}}_{1})\cdots(\hat{\mathbf{Q}}^{\frac{1}{2}}\hat{\mathbf{u}}_{L})\right]\in\mathbb{C}^{M\times
L}.$ (35)

Replacing $\widetilde{\mathbf{u}}_{p}$ by $\hat{\widetilde{\mathbf{u}}}_{p}$
in (33) and generating $\hat{\mathbf{H}}$ and $\hat{\mathbf{h}}$, (33) turns
into an approximate equality, i.e.,

$\displaystyle\hat{\mathbf{H}}\mathbf{a}\approx\hat{\mathbf{h}}.\ $ (36)

The optimal value of $\mathbf{a}$ can now be obtained employing the GLS
technique [44], [45]. Using (36), the GLS optimization problem can be
formulated as

$\displaystyle\hat{\mathbf{a}}=\mathrm{arg}\displaystyle\min_{\mathbf{a}}(\hat{\mathbf{H}}\mathbf{a}-\hat{\mathbf{h}})^{H}\mathbf{W}(\hat{\mathbf{H}}\mathbf{a}-\hat{\mathbf{h}})$
(37)

where
$\mathbf{W}\triangleq\left(\mathbb{E}\\{\hat{\mathbf{e}}\hat{\mathbf{e}}^{H}\\}\right)^{-1}\in\mathbb{C}^{L(|\mathcal{I}|-L)\times
L(|\mathcal{I}|-L)}$ and
$\hat{\mathbf{e}}\triangleq\hat{\mathbf{H}}\mathbf{a}-\hat{\mathbf{h}}\in\mathbb{C}^{L(|\mathcal{I}|-L)}$.
The solution of (37) is given by

$\displaystyle\hat{\mathbf{a}}=(\hat{\mathbf{H}}^{H}\mathbf{W}\hat{\mathbf{H}})^{-1}\hat{\mathbf{H}}^{H}\mathbf{W}\hat{\mathbf{h}}.\
$ (38)

To use (38), an estimate of $\mathbf{W}$ is required. However, it is clear
from the definitions of $\mathbf{W}$ and $\hat{\mathbf{e}}$ that an estimate
of $\mathbf{W}$ depends on the unknown vector $\mathbf{a}$. Thus, it is
natural to utilize an iterative scheme to estimate $\hat{\mathbf{a}}$ in one
step, followed by estimating $\hat{\mathbf{W}}$ in the other step by employing
$\hat{\mathbf{a}}$ obtained in the previous step. The alternation between
these two steps is then carried on until a termination criteria is satisfied.

To figure out a principle for finding $\hat{\mathbf{W}}$, we take into
consideration the first-order subspace estimation error by expressing
$\hat{\mathbf{u}}_{p}$ as
$\hat{\mathbf{u}}_{p}\triangleq\mathbf{u}_{p}+\mathbf{\Delta}\mathbf{u}_{p}$
for $p=1,\cdots,L$. Exploiting this definition, we have

$\displaystyle\hat{\mathbf{e}}_{p}$
$\displaystyle\triangleq\hat{\mathbf{H}}_{p}\mathbf{a}-\hat{\mathbf{h}}_{p}=\mathbf{B}^{H}\left(\mathbf{Z}_{\mathcal{I}}\hat{\bar{\mathbf{u}}}_{p}+\mathrm{diag}\\{\mathbf{Z}_{\mathcal{I}}\hat{\bar{\mathbf{u}}}_{p}\\}\mathbf{Z}_{\mathcal{I}}\mathbf{W}_{a}\mathbf{a}\right)$
$\displaystyle=\mathbf{B}^{H}\left(\mathbf{Z}_{\mathcal{I}}\mathbf{W}_{D}\hat{\mathbf{Q}}^{\frac{1}{2}}\hat{\mathbf{u}}_{p}+\mathrm{diag}\\{\mathbf{Z}_{\mathcal{I}}\mathbf{W}_{D}\hat{\mathbf{Q}}^{\frac{1}{2}}\hat{\mathbf{u}}_{p}\\}\mathbf{Z}_{\mathcal{I}}\mathbf{W}_{a}\mathbf{a}\right)$
$\displaystyle=\mathbf{B}^{H}\left(\mathbf{Z}_{\mathcal{I}}\mathbf{W}_{D}\hat{\mathbf{Q}}^{\frac{1}{2}}\hat{\mathbf{u}}_{p}+\mathrm{diag}\\{\mathbf{Z}_{\mathcal{I}}\mathbf{W}_{a}\mathbf{a}\\}\mathbf{Z}_{\mathcal{I}}\mathbf{W}_{D}\hat{\mathbf{Q}}^{\frac{1}{2}}\hat{\mathbf{u}}_{p}\right)$
$\displaystyle=\mathbf{B}^{H}\left(\mathbf{I}_{|\mathcal{I}|}+\mathrm{diag}\\{\mathbf{Z}_{\mathcal{I}}\mathbf{W}_{a}\mathbf{a}\\}\right)\mathbf{Z}_{\mathcal{I}}\mathbf{W}_{D}\hat{\mathbf{Q}}^{\frac{1}{2}}\hat{\mathbf{u}}_{p}$
$\displaystyle=\mathbf{C}(\mathbf{a})\hat{\mathbf{u}}_{p}\ ,\ \ \ \
p=1,\cdots,L\ \ $ (39)

where
$\mathbf{C}(\mathbf{a})\triangleq\mathbf{B}^{H}\left(\mathbf{I}_{|\mathcal{I}|}+\mathrm{diag}\\{\mathbf{Z}_{\mathcal{I}}\mathbf{W}_{a}\mathbf{a}\\}\right)\mathbf{Z}_{\mathcal{I}}\mathbf{W}_{D}\hat{\mathbf{Q}}^{\frac{1}{2}}\in\mathbb{C}^{(|\mathcal{I}|-L)\times
M}$. In (39), we also used the property
$\mathrm{diag}\\{\mathbf{x}\\}\mathbf{y}=\mathrm{diag}\\{\mathbf{y}\\}\mathbf{x}$
as well as the equality
$\hat{\bar{\mathbf{u}}}_{p}=\mathbf{W}_{D}\hat{\widetilde{\mathbf{u}}}_{p}=\mathbf{W}_{D}\hat{\mathbf{Q}}^{\frac{1}{2}}\hat{\mathbf{u}}_{p}$.

Using the definition of $\mathbf{C}(\mathbf{a})$, (32) can be recast as

$\displaystyle\mathbf{C}(\mathbf{a})[\mathbf{u}_{1}\mathbf{u}_{2}\cdots\mathbf{u}_{L}]=\mathbf{C}(\mathbf{a})\mathbf{U}_{\rm
s}=\mathbf{0}_{(|\mathcal{I}|-L)\times L}.$ (40)

According to the definitions of $\mathbf{W}$, $\hat{\mathbf{e}}$ and
$\mathbf{C}(\mathbf{a})$, and using (39) as well as (40), $\hat{\mathbf{e}}$
can be found as

$\displaystyle\hat{\mathbf{e}}$
$\displaystyle=\mathrm{vec}\\{[\hat{\mathbf{e}}_{1}\hat{\mathbf{e}}_{2}\cdots\hat{\mathbf{e}}_{L}]\\}=\mathrm{vec}\\{\mathbf{C}(\mathbf{a})\hat{\mathbf{U}}_{\rm
s}\\}=\mathrm{vec}\\{\mathbf{C}(\mathbf{a})(\mathbf{U}_{\rm
s}+\mathbf{\Delta}\mathbf{U}_{\rm s})\\}$
$\displaystyle=\mathrm{vec}\\{\mathbf{C}(\mathbf{a})\mathbf{\Delta}\mathbf{U}_{\rm
s}\\}=(\mathbf{I}_{L}\otimes\mathbf{C}(\mathbf{a}))\mathrm{vec}\\{\mathbf{\Delta}\mathbf{U}_{\rm
s}\\}=(\mathbf{I}_{L}\otimes\mathbf{C}(\mathbf{a}))\mathbf{\Delta}\mathbf{u}_{\rm
s}$ (41)

where $\mathbf{\Delta}\mathbf{U}_{\rm
s}\triangleq[\mathbf{\Delta}\mathbf{u}_{1}\cdots\mathbf{\Delta}\mathbf{u}_{L}]\in\mathbb{C}^{M\times
L}$ and $\mathbf{\Delta}\mathbf{u}_{\rm
s}\triangleq\mathrm{vec}\\{\mathbf{\Delta}\mathbf{U}_{\rm
s}\\}\in\mathbb{C}^{ML}$. The identity
$\mathrm{vec}\\{\mathbf{X}\mathbf{Y}\mathbf{Z}\\}=(\mathbf{Z}^{T}\otimes\mathbf{X})\mathrm{vec}\\{\mathbf{Y}\\}$
has also being used in (41). Moreover, according to the definitions of
$\mathbf{W}$ and $\hat{\mathbf{e}}$, and using (41), $\mathbf{W}$ can be
expressed as

$\displaystyle\mathbf{W}$
$\displaystyle=\left[\left(\mathbf{I}_{L}\otimes\mathbf{C}(\mathbf{a})\right)\mathbb{E}\\{\mathbf{\Delta}\mathbf{u}_{\rm
s}\mathbf{\Delta}\mathbf{u}_{\rm
s}^{H}\\}\left(\mathbf{I}_{L}\otimes\mathbf{C}(\mathbf{a})\right)^{H}\right]^{-1}.$
(42)

Using the first-order perturbation expansion for SVD [24], [46], we can write
that

$\displaystyle\mathbf{\Delta}\mathbf{U}_{\rm
s}\approx(\mathbf{I}_{M}-\mathbf{U}_{\rm s}\mathbf{U}_{\rm
s}^{H})\overline{\mathbf{N}}\mathbf{V}_{\rm s}\boldsymbol{\Sigma}_{\rm
s}^{-1}.$ (43)

Applying the vectorization operator to (43), we have

$\displaystyle\mathbf{\Delta}\mathbf{u}_{\rm s}$
$\displaystyle=\mathrm{vec}\\{\mathbf{\Delta}\mathbf{U}_{\rm
s}\\}\approx\left(\boldsymbol{\Sigma}_{\rm s}^{-1}\mathbf{V}_{\rm
s}^{T}\otimes(\mathbf{I}_{M}-\mathbf{U}_{\rm s}\mathbf{U}_{\rm
s}^{H})\right)\bar{\mathbf{n}}$ (44)

where
$\bar{\mathbf{n}}\triangleq\mathrm{vec}\\{\overline{\mathbf{N}}\\}\in\mathbb{C}^{MN}$.
Inserting (44) into (42) yields

$\displaystyle\mathbf{W}$
$\displaystyle\approx[(\mathbf{I}_{L}\otimes\mathbf{C}(\mathbf{a}))(\boldsymbol{\Sigma}_{\rm
s}^{-1}\mathbf{V}_{\rm s}^{T}\otimes(\mathbf{I}_{M}-\mathbf{U}_{\rm
s}\mathbf{U}_{\rm
s}^{H}))\mathbb{E}\\{\bar{\mathbf{n}}\bar{\mathbf{n}}^{H}\\}(\mathbf{V}_{\rm
s}^{*}\boldsymbol{\Sigma}_{\rm s}^{-1}\otimes(\mathbf{I}_{M}-\mathbf{U}_{\rm
s}\mathbf{U}_{\rm
s}^{H}))(\mathbf{I}_{L}\otimes\mathbf{C}^{H}(\mathbf{a}))]^{-1}$
$\displaystyle=[(\boldsymbol{\Sigma}_{\rm s}^{-1}\mathbf{V}_{\rm
s}^{T}\otimes\mathbf{C}(\mathbf{a})(\mathbf{I}_{M}-\mathbf{U}_{\rm
s}\mathbf{U}_{\rm
s}^{H}))(\mathbf{I}_{M}\otimes\mathbf{I}_{N})(\mathbf{V}_{\rm
s}^{*}\boldsymbol{\Sigma}_{\rm s}^{-1}\otimes(\mathbf{I}_{M}-\mathbf{U}_{\rm
s}\mathbf{U}_{\rm s}^{H})\mathbf{C}^{H}(\mathbf{a}))]^{-1}$ (45)

where
$\mathbb{E}\\{\bar{\mathbf{n}}\bar{\mathbf{n}}^{H}\\}=\mathbf{I}_{MN}=(\mathbf{I}_{M}\otimes\mathbf{I}_{N})$.
The property
$(\mathbf{X}\otimes\mathbf{Y})(\mathbf{Z}\otimes\mathbf{T})=(\mathbf{X}\mathbf{Z}\otimes\mathbf{Y}\mathbf{T})$
has been used for deriving (45). Using this property again together with (40),
(45) can be further simplified as

$\displaystyle\mathbf{W}$ $\displaystyle\approx[(\boldsymbol{\Sigma}_{\rm
s}^{-1}\mathbf{V}_{\rm s}^{T}\mathbf{V}_{\rm s}^{*}\boldsymbol{\Sigma}_{\rm
s}^{-1}\otimes\mathbf{C}(\mathbf{a})\mathbf{C}^{H}(\mathbf{a}))]^{-1}$
$\displaystyle=(\boldsymbol{\Sigma}_{\rm
s}^{-2}\otimes\mathbf{C}(\mathbf{a})\mathbf{C}^{H}(\mathbf{a}))^{-1}=(\boldsymbol{\Sigma}_{\rm
s}^{2}\otimes(\mathbf{C}(\mathbf{a})\mathbf{C}^{H}(\mathbf{a}))^{-1}).$ (46)

The identities $\mathbf{V}_{\rm s}^{T}\mathbf{V}_{\rm s}^{*}=\mathbf{I}_{L}$
and
$(\mathbf{X}\otimes\mathbf{Y})^{-1}=(\mathbf{X})^{-1}\otimes(\mathbf{Y})^{-1}$
have also been used here.

However, since the matrix $\boldsymbol{\Sigma}_{\rm s}$ is unknown, we replace
it with $\hat{\boldsymbol{\Sigma}}_{\rm s}$ obtained from (34). Eventually,
$\hat{\mathbf{W}}$ is expressed as

$\displaystyle\hat{\mathbf{W}}$
$\displaystyle\approx(\hat{\boldsymbol{\Sigma}}_{\rm
s}^{2}\otimes(\mathbf{C}(\mathbf{a})\mathbf{C}^{H}(\mathbf{a}))^{-1}).$ (47)

It can now be seen from (47) that $\hat{\mathbf{W}}$ is a function of
$\mathbf{a}$. As a result, an iterative scheme should be used to estimate
$\hat{\mathbf{a}}$ and $\hat{\mathbf{W}}$ via (38) and (47) in alternative
manner. The LS solution of (36) can be used as the initial vector for
$\hat{\mathbf{a}}$, i.e.,

$\displaystyle\hat{\mathbf{a}}^{(0)}\triangleq\hat{\mathbf{a}}_{LS}=\hat{\mathbf{H}}^{\dagger}\hat{\mathbf{h}}.$
(48)

Initializing $\mathbf{C}(\mathbf{a})$ by inserting (48) into the definition of
$\mathbf{C}(\mathbf{a})$, enables us to obtain

$\displaystyle\hat{\mathbf{W}}^{(0)}=(\hat{\boldsymbol{\Sigma}}_{\rm
s}^{2}\otimes(\mathbf{C}(\hat{\mathbf{a}}^{(0)})\mathbf{C}(\hat{\mathbf{a}}^{(0)})^{H})^{-1}).$
(49)

Subsequently, a new estimate of $\hat{\mathbf{a}}$ is generated by
substituting (49) into (38). The iterations carry on until a proper
termination criteria is satisfied222Although any common termination criteria
can be adopted, we observe that performing 3 to 5 iterations are usually
sufficient to obtain a precise result. Therefore, in the numerical examples, 5
is opted as the number of iterations for implementing the proposed methods..

The remaining problem is still how to determine the members of $\mathcal{I}$,
i.e., the indices of the selected equations of (30) to be employed for
estimating DOAs. Given the cardinality $|\mathcal{I}|$, it is reasonable to
select those indices which correspond to the entries of
$\hat{\bar{\mathbf{u}}}_{1}=\mathrm{DFT}\\{\hat{\widetilde{\mathbf{u}}}_{1}\\}$
with $|\mathcal{I}|$ largest absolute values. According to (35),
$\hat{\widetilde{\mathbf{u}}}_{1}=\hat{\mathbf{Q}}^{\frac{1}{2}}\hat{\mathbf{u}}_{1}$
with $\hat{\mathbf{u}}_{1}$ denoting the left singular vector of
$\overline{\mathbf{X}}$ which corresponds to the largest singular value. The
logic of this choice is rooted in (22), where it is indicated that each
$\hat{\widetilde{\mathbf{u}}}_{p}$ can be expressed as a linear combination of
the columns of $\mathbf{A}$. Therefore, picking the indices of
$\hat{\bar{\mathbf{u}}}_{1}$ with largest absolute values is a sensible choice
because of the following three reasons. 1) The structure of the DFT basis is
completely matched with the columns of $\mathbf{A}$, and consequently
$\hat{\widetilde{\mathbf{u}}}_{p}$’s, which makes the absolute values of
$\hat{\bar{\mathbf{u}}}_{p}$’s the best option to be used for selecting the
most relevant equations. 2) Choosing the indices with largest absolute values
guarantees picking equations with the most contributions. 3) The estimation
error of finding $\hat{\mathbf{u}}_{1}$ is the smallest among
$\hat{\mathbf{u}}_{p}$’s since it associates with the largest singular values,
resulting in smaller error. Finally, the steps required for implementing the
proposed DOA estimation method are summarized in Algorithm 2.

Algorithm 2: Forward-only DOA Estimation
---
1: Compute $\hat{\mathbf{R}}=1/N\sum_{t=1}^{N}\mathbf{x}(t)\mathbf{x}^{H}(t)$
and estimate $\hat{\mathbf{Q}}$ using Algorithm 1.
2: Calculate $\overline{\mathbf{X}}=\hat{\mathbf{Q}}^{-\frac{1}{2}}\mathbf{X}$
and construct $\mathbf{W}_{a}$ as well as $\overline{{\mathbf{W}}}$.
3: For a pre-chosen $|\mathcal{I}|$, determine $\mathbf{Z}_{\mathcal{I}}$ and
obtain $\mathbf{B}$ via performing the SVD of the matrix
$\widetilde{{\mathbf{Z}}}_{\mathcal{I}}=\mathbf{Z}_{\mathcal{I}}\overline{\mathbf{W}}$
so that the
condition $\mathbf{B}^{H}\widetilde{{\mathbf{Z}}}_{\mathcal{I}}=\mathbf{0}$ is
satisfied.
4: Carry out the SVD of $\overline{\mathbf{X}}$ to obtain
$\hat{\mathbf{U}}_{\rm s}$ and $\hat{\boldsymbol{\Sigma}}_{\rm s}$, where the
former contains the $L$ left singular vectors corresponding to
the $L$ largest singular values on the diagonal of the latter.
5: Utilize (35) for computing $\hat{\widetilde{\mathbf{U}}}_{\rm s}$. Apply
DFT on the columns of $\hat{\widetilde{\mathbf{U}}}_{\rm s}$ to obtain
$\hat{\bar{\mathbf{u}}}_{p}$’s.
6: Form $\hat{\mathbf{H}}$ and $\hat{\mathbf{h}}$ using the estimates
$\hat{\bar{\mathbf{u}}}_{p}$’s instead of $\bar{\mathbf{u}}_{p}$’s.
7: Set $i=0$ and
$\hat{\mathbf{a}}^{(0)}=\hat{\mathbf{H}}^{\dagger}\hat{\mathbf{h}}$. In
addition, set the maximum number of iterations $i_{max}=5$.
while $i\leq i_{max}$
8: Compute $\mathbf{C}\left(\hat{\mathbf{a}}^{(i)}\right)$ and
$\hat{\mathbf{W}}^{(i)}$.
9: Generate a new estimate
$\hat{\mathbf{a}}^{(i+1)}=\left(\hat{\mathbf{H}}^{H}\hat{\mathbf{W}}^{(i)}\hat{\mathbf{H}}\right)^{-1}\hat{\mathbf{H}}^{H}\hat{\mathbf{W}}^{(i)}\hat{\mathbf{h}}$.
10: set $i=i+1$.
end
11: Find the $L$ roots of the polynomial
$\gamma^{L}+\sum_{l=1}^{L}[\hat{\mathbf{a}}]_{l}\ \gamma^{L-l}=0$, denoted by
$\hat{\gamma}_{l},l=1,\cdots,L$.
12: Obtain the $L$ DOA estimates as
$\hat{\theta}_{l}=\mathrm{arcsin}\left(-\frac{\beta_{l}\lambda}{2\pi
d}\right)$, where $\beta_{l}$ is the phase argument of
$\hat{\gamma}_{l},\;l=1,\cdots,L$.

To provide a theoretical measure for performance of the proposed forward-only
DOA estimation method, the asymptotic variance of the $l$th DOA estimated by
the proposed forward-only algorithm is derived under a high SNR assumption in
the following proposition.

Proposition 2: The asymptotic variance of the proposed forward-only DOA
estimation algorithm, for a particular matrix $\mathbf{Z}_{\mathcal{I}}$, is
given as

$\displaystyle\mathbb{E}\\{\Delta\theta_{l}^{2}\\}$
$\displaystyle\approx\frac{1}{2}\left(\frac{\lambda}{2\pi
d\cos(\theta_{l})}\right)^{2}\frac{\boldsymbol{\gamma}_{l}^{T}(\mathbf{H}^{H}\mathbf{W}\mathbf{H})^{-1}\boldsymbol{\gamma}_{l}^{*}}{|\phi_{l}|^{2}}$
(50)

where $\boldsymbol{\gamma}_{l}\triangleq[\gamma_{l}^{L-1}\cdots 1]^{T}$ and
$\phi_{l}\triangleq
L\gamma_{l}^{L-1}+(L-1)[\mathbf{a}]_{1}\gamma_{l}^{L-2}+\cdots+[\mathbf{a}]_{L-1}$.

Proof: See Appendix C. $\blacksquare$

We aim to generate double number of DOA candidates333In fact, the number of
DOA candidates can be arbitrary, but from diverse numerical simulations
conducted, we find that generating more candidates than $2L$ does not improve
the DOA estimation accuracy considerably. by running the proposed Algorithm 2
twice for two different values of $|\mathcal{I}|$. Then a proper DOA selection
strategy should be employed to determine the final DOA estimates. First,
consider the following example to get a more through wisdom about how
different choices of $|\mathcal{I}|$ affect the DOA estimation accuracy.

Illustrative Example 1: Consider a scenario where the signals of two
uncorrelated sources located in $\boldsymbol{\theta}=[-2^{\circ},7^{\circ}]$
are received by a ULA consisting $M=8$ sensors with half waveform adjacent
distances. The sensor noise covariance matrix is set as
$\mathbf{Q}=\mathrm{diag}\\{[10,1.2,3.5,18,2,8.5,24,6.5]\\}$, the sample size
is $N=40$, and 2000 Monte Carlo runs are conducted to calculate RMSE defined
as

$\displaystyle{\rm RMSE}=10{\rm
log}_{10}\sqrt{\frac{1}{2000L}\displaystyle\sum_{l=1}^{L}\displaystyle\sum_{i=1}^{2000}(\hat{\theta}_{l,i}-\theta_{l})^{2}}.$
(51)

We also include the deterministic CRB [32] as a benchmark. In Fig. 1, higher
accuracy in DOA estimation can be observed for larger $|\mathcal{I}|$. From
this observation, we conclude that $|I|=M-1$ and $|I|=M$ are the best choices
in the sense of providing the most precise estimates for generating double
number of DOA candidates.

Figure 1: RMSE performance of the proposed method for different
$|\mathcal{I}|$ vs. SNR for $L=2$ uncorrelated sources with
$\boldsymbol{\theta}=[-2^{\circ},7^{\circ}]$, $M=8$, and $N=40$.

###  IV-B DOA Selection Strategy

After running Algorithm 2 twice with $|\mathcal{I}|=M-1$ and
$|\mathcal{I}|=M$, $2L$ DOA candidates are generated. The natural question is
that how to select $L$ final DOA estimates. Two known conventional approaches
to DOA selection are based on conventional beamformer (CB) [16], [47, 48, 49]
and ML cost function minimization [22], [23], [50]. Recently, another method,
which has been originally proposed for joint source number detection and DOA
estimation [51], has been employed as the DOA selection scheme. It is based on
the generalized likelihood ratio (GLR), which extracts the final $L$ DOAs
sequentially [29]. The computational complexity of such DOA selection strategy
is much lower than that of based on the ML cost function minimization.
Besides, the performance provided by the GLR is comparable with that provided
by the ML-based methods. Here we aim to design a DOA selection strategy which
takes advantage of the three aforementioned approaches, i.e, it uses CB,
deterministic ML cost function, and the GLR technique.

In doing so, the following three-step selection strategy is proposed.

Step 1: Denote the vector containing all $2L$ DOA candidates by
$\boldsymbol{\theta}_{2L}$. Calculate the CB output for a proper number of
equidistant points to cover the whole interval of interest, i.e.,
$[-\frac{\pi}{2},\frac{\pi}{2}]$. Then, find the $(L+1)$th peak and define a
threshold, denoted by $\eta$, as the CB output at the $(L+1)$th peak444Note
that if the total number of peaks is smaller than $(L+1)$, select the last
peak with the lowest output.. Afterwards, calculate the CB output for entries
of $\boldsymbol{\theta}_{2L}$ and stack those entries with output larger than
$\eta$ in a new vector $\widetilde{\boldsymbol{\theta}}$. If the number of
entries in $\widetilde{\boldsymbol{\theta}}$ is less than $L$, then let
$\widetilde{\boldsymbol{\theta}}=\boldsymbol{\theta}_{2L}$.

Step 2: Determine the first DOA as that entry of
$\widetilde{\boldsymbol{\theta}}$ which maximizes the GLR, i.e.,

$\displaystyle\hat{\theta}_{1}=\mathrm{arg}\displaystyle\max_{\theta}\frac{\mathbf{a}^{H}(\theta)\hat{\mathbf{Q}}^{-1}\hat{\mathbf{R}}\hat{\mathbf{Q}}^{-1}\mathbf{a}(\theta)}{\mathbf{a}^{H}(\theta)\hat{\mathbf{Q}}^{-1}\mathbf{a}(\theta)}\
\ \ ,\ \theta\in\widetilde{\boldsymbol{\theta}}.$ (52)

The GLR in (52) is a straightforward nonuniform noise extension of the GLR in
[51] for the case of uniform noise.

Step 3: Denote the remaining entries of $\widetilde{\boldsymbol{\theta}}$ as
$\bar{\boldsymbol{\theta}}$, and the size of $\bar{\boldsymbol{\theta}}$ as
$\bar{L}$. Then, divide the $\bar{L}$ DOA candidates in
$\bar{\boldsymbol{\theta}}$ into
$\bar{G}=\frac{\bar{L}}{(L-1)!(\bar{L}-L+1)!}$ subsets containing $(L-1)$
different DOAs each. Denote these subsets as
$\Theta_{1},\cdots,\Theta_{\bar{G}}$ and associate them with
$\mathbf{A}(\Theta_{1}),\cdots,\mathbf{A}(\Theta_{\bar{G}})$. The subset that
minimizes the following deterministic ML cost function555Since the difference
of employing the deterministic ML and stochastic ML is marginal, the
deterministic ML is discussed here because it also has lower complexity.
determines the $(L-1)$ remaining DOAs

$\displaystyle\hat{\Theta}_{R}=\mathrm{arg}\displaystyle\min_{\Theta_{\mathbb{S}}}$
$\displaystyle\
\mathrm{trace}\left[\left(\mathbf{P}^{\perp}_{\widetilde{A}(\Theta_{\mathbb{S}})}-\boldsymbol{\nu}_{1}\boldsymbol{\nu}_{1}^{H}\right)\hat{\mathbf{Q}}^{-\frac{1}{2}}\hat{\mathbf{R}}\hat{\mathbf{Q}}^{-\frac{1}{2}}\right],\
\ \ \mathbb{S}\in\\{1,\cdots,\bar{G}\\}$ (53)

where
$\mathbf{P}^{\perp}_{\widetilde{A}(\Theta_{\mathbb{S}})}\triangleq\mathbf{I}_{M}-\widetilde{\mathbf{A}}(\Theta_{\mathbb{S}})\left(\widetilde{\mathbf{A}}(\Theta_{\mathbb{S}})^{H}\widetilde{\mathbf{A}}(\Theta_{\mathbb{S}})\right)^{-1}\widetilde{\mathbf{A}}(\Theta_{\mathbb{S}})^{H}$,
$\widetilde{\mathbf{A}}(\Theta_{\mathbb{S}})\triangleq\hat{\mathbf{Q}}^{-\frac{1}{2}}\mathbf{A}(\Theta_{\mathbb{S}})$,
and
$\boldsymbol{\nu}_{1}\triangleq\frac{\mathbf{P}^{\perp}_{\widetilde{A}(\Theta_{\mathbb{S}})}\hat{\mathbf{Q}}^{-\frac{1}{2}}\mathbf{a}(\hat{\theta}_{1})}{\|\mathbf{P}^{\perp}_{\widetilde{A}(\Theta_{\mathbb{S}})}\hat{\mathbf{Q}}^{-\frac{1}{2}}\mathbf{a}(\hat{\theta}_{1})\|_{2}}$.
In (53), the contribution of previously estimated $\hat{\theta}_{1}$ on the
deterministic ML cost function is isolated in $\boldsymbol{\nu}_{1}$ thanks to
the properties of the orthogonal projection matrix [52]. In addition,
$\hat{\mathbf{Q}}^{-\frac{1}{2}}\hat{\mathbf{R}}\hat{\mathbf{Q}}^{-\frac{1}{2}}$
is used instead of $\hat{\mathbf{R}}$ to consider the general case of
nonuniform noise.

Finally, $\hat{\theta}_{1}$ and $\hat{\Theta}_{R}$ obtained via (53) form
together the final $L$ DOA estimates.

###  IV-C Computational Complexity

For implementing the proposed DOA estimation method, Algorithm 1 should be
implemented first that requires $\mathcal{O}\left(M^{2}N\right)$ flops for
calculating $\hat{\mathbf{R}}$ and $\mathcal{O}\left(M^{3}\right)$ flops for
computing GED of the pair of matrices
$\left\\{\hat{\mathbf{R}},\;\hat{\mathbf{Q}}^{(i)}\right\\}$ in the $i$th
iteration. Thus, the total computational complexity of Algorithm 1 is
$\mathcal{O}\left(I_{1}(M^{3})+M^{2}N\right)$ with $I_{1}$ being the number of
iterations. The computational complexity of performing SVD of
$\overline{\mathbf{X}}$ is
$\mathcal{O}\left(\mathrm{max}(M,N)\mathrm{min}(M,N)^{2}\right)$, and the
computational complexity of computing $\mathbf{B}$ is
$\mathcal{O}\left(|\mathcal{I}|L^{2}\right)$. Applying DFT on the columns of
$\hat{\widetilde{U}}_{\rm s}$ requires
$\mathcal{O}\left(LM\mathrm{log}_{2}(M)\right)$ flops. In addition, the main
source of computational complexity for calculating $\hat{\mathbf{W}}$ is the
inversion of the matrix $\mathbf{C}(\mathbf{a})\mathbf{C}^{H}(\mathbf{a})$,
which involves $\mathcal{O}\left((|\mathcal{I}|-L)^{3}\right)$ flops, making
the computational complexity of computing $\hat{\mathbf{W}}$ be
$\mathcal{O}\left(L(|\mathcal{I}|-L)^{3}\right)$. Determining
$\hat{\mathbf{a}}$ requires
$\mathcal{O}\left(L^{3}(2(M-L)+1)+L^{2}(2(M-L)^{2}+(M-L))\right)$ flops. As
the proposed DOA estimation algorithm is run twice with $|\mathcal{I}|=M-1$
and $|\mathcal{I}|=M$, by considering $|\mathcal{I}|\approx M$, the total
complexity required to generate $2L$ DOA candidates is
$\mathcal{O}\left(M^{2}N+I_{1}(M^{3})+(\mathrm{max}(M,N)\mathrm{min}(M,N)^{2})+ML^{2}+\right.$
$\left.LM\mathrm{log}_{2}(M)+I_{2}(M^{3}L+ML^{3}-L^{4})\right)$ where $I_{2}$
denotes the number of iterations needed for finding the best
$\hat{\mathbf{a}}$. At last, the complexity of DOA selection is mainly in step
3, which is about
$\mathcal{O}\left(\bar{G}(M^{3}+3M(L-1)^{2}+(L-1)^{3})\right)$. Finally, the
total complexity of DOA estimation is reduced to
$\mathcal{O}\left(M^{2}N+\bar{G}M^{3}\right)$ in the case that
$\bar{G}\gg\mathrm{max}(I_{1},I_{2})$ and $M\gg L$. The significant point here
is that $\bar{G}$ is approximately four times smaller than the parameter $G$
in [23], resulting in the computational complexity of our method being
approximately a quarter of what is required for implementing the EPUMA.

###  IV-D FBA Extension

FBA is a natural extension/improvement of the forward-only-based DOA
estimation methods [18], [53, 30, 54, 55, 37]. The essence of FBA is to first
transform the observed signal or SCM into a new centro-Hermitian signal matrix
or centro-Hermition covariance matrix, respectively, followed by
computationally simplified and more accurate DOA estimation. Higher DOA
estimation accuracy is a consequence of decorrelating possibly correlated
source pairs and obtaining more accurate SCM estimation.

The FBA covariance matrix is given as [15]

$\displaystyle\mathbf{R}_{\rm FB}$
$\displaystyle=\frac{1}{2}\left(\mathbf{R}+\mathbf{J}_{M}\mathbf{R}^{*}\mathbf{J}_{M}\right)=\frac{1}{2}\left(\mathbf{A}\mathbf{P}\mathbf{A}^{H}+\mathbf{Q}+\mathbf{J}_{M}(\mathbf{A}^{*}\mathbf{P}^{*}\mathbf{A}^{T}+\mathbf{Q})\mathbf{J}_{M}\right).$
(54)

It can be readily verified that $\mathbf{R}_{\rm FB}$ is centro-Hermitian,
i.e., $\mathbf{R}_{\rm FB}=\mathbf{J}_{M}\mathbf{R}^{*}_{FB}\mathbf{J}_{M}$.
Rearranging (54), we have

$\displaystyle\mathbf{R}_{FB}$
$\displaystyle=\frac{1}{2}\left(\mathbf{A}\mathbf{P}\mathbf{A}^{H}+\mathbf{Q}+\mathbf{A}\mathbf{D}\mathbf{P}^{*}\mathbf{D}^{H}\mathbf{A}^{H}+\mathbf{J}_{M}\mathbf{Q}\mathbf{J}_{M}\right)=\mathbf{A}\widetilde{\mathbf{P}}\mathbf{A}^{H}+\frac{1}{2}\widetilde{\mathbf{Q}}$
(55)

where
$\widetilde{\mathbf{P}}\triangleq\frac{1}{2}(\mathbf{P}+\mathbf{D}\mathbf{P}^{*}\mathbf{D}^{H})$,
$\widetilde{\mathbf{Q}}\triangleq\mathbf{Q}+\mathbf{J}_{M}\mathbf{Q}\mathbf{J}_{M}$,
$\mathbf{D}\triangleq\mathrm{diag}\left\\{e^{-j(2\pi/\lambda)d(M-1)\sin(\theta_{1})},\cdots,e^{-j(2\pi/\lambda)d(M-1)\sin(\theta_{L})}\right\\}$.

According to (55), a proper centro-Hermitian matrix similar to
$\overline{\mathbf{X}}$ defined in (19), can be formed as

$\displaystyle\overline{\mathbf{X}}_{\rm FB}$
$\displaystyle=\left[\hat{\widetilde{\mathbf{Q}}}^{-\frac{1}{2}}\mathbf{X}\ \
\ \
\mathbf{J}_{M}\hat{\widetilde{\mathbf{Q}}}^{-\frac{1}{2}}\mathbf{X}^{*}\mathbf{J}_{N}\right]$
(56)

where
$\hat{\widetilde{\mathbf{Q}}}\triangleq\hat{\mathbf{Q}}+\mathbf{J}_{M}\hat{\mathbf{Q}}\mathbf{J}_{M}$
is used instead of $\widetilde{\mathbf{Q}}$. The matrix
$\overline{\mathbf{X}}_{\rm FB}$ can be decomposed by applying the truncated
SVD as

$\displaystyle\overline{\mathbf{X}}_{FB}=\hat{\mathbf{E}}_{\rm
s}\hat{\boldsymbol{\Pi}}_{\rm s}\hat{\mathbf{T}}_{\rm s}^{H}$ (57)

where $\hat{\mathbf{E}}_{\rm
s}=[\hat{\mathbf{e}}_{1}\cdots\hat{\mathbf{e}}_{L}]\in\mathbb{C}^{M\times L}$
is composed of $L$ left singular vectors associated with $L$ largest singular
values on the diagonal of $\hat{\boldsymbol{\Pi}}_{\rm
s}\in\mathbb{R}^{L\times L}$. The columns of $\hat{\mathbf{E}}_{\rm s}$ and
the columns of $\hat{\widetilde{\mathbf{Q}}}^{-\frac{1}{2}}\mathbf{A}$ span
the same vector space. Thus, similar to the forward-only case, the
relationship between $\hat{\mathbf{E}}_{\rm s}$ and
$\hat{\widetilde{\mathbf{Q}}}^{-\frac{1}{2}}\mathbf{A}$ can be written as
$\hat{\mathbf{E}}_{\rm
s}=\hat{\widetilde{\mathbf{Q}}}^{-\frac{1}{2}}\mathbf{A}\overline{\mathbf{G}}$.
Multiplying both sides of the latter equation by
$\hat{\widetilde{\mathbf{Q}}}^{\frac{1}{2}}$, we obtain the FBA analog of
(22), that is,

$\displaystyle\hat{\widetilde{\mathbf{E}}}_{\rm
s}\triangleq\hat{\widetilde{\mathbf{Q}}}^{\frac{1}{2}}\hat{\mathbf{E}}_{\rm
s}=\mathbf{A}\overline{\mathbf{G}}.$ (58)

Following the same steps as in Subsection IV-A, the DFT of the columns of
$\hat{\widetilde{\mathbf{E}}}_{\rm s}$ can be found as

$\displaystyle\hat{\overline{\mathbf{E}}}_{\rm s}$
$\displaystyle=[\hat{\bar{\mathbf{e}}}_{1}\cdots\hat{\bar{\mathbf{e}}}_{L}]\triangleq\mathrm{DFT}\\{\hat{\widetilde{\mathbf{E}}}_{\rm
s}\\}=\mathrm{DFT}\\{\hat{\widetilde{\mathbf{Q}}}^{\frac{1}{2}}\hat{\mathbf{E}}_{\rm
s}\\}$
$\displaystyle=\left[(\mathbf{W}_{D}\hat{\widetilde{\mathbf{Q}}}^{\frac{1}{2}}\hat{\mathbf{e}}_{1})\cdots(\mathbf{W}_{D}\hat{\widetilde{\mathbf{Q}}}^{\frac{1}{2}}\hat{\mathbf{e}}_{L})\right]\in\mathbb{C}^{M\times
L}.\ $ (59)

Considering a specific value for $|\mathcal{I}|$, a system of linear equations
similar to (36) can be formulated as

$\displaystyle\hat{\widetilde{\mathbf{H}}}$
$\displaystyle\mathbf{a}\approx\hat{\widetilde{\mathbf{h}}}$ (60)

where
$\hat{\widetilde{\mathbf{H}}}_{p}=\mathbf{B}^{H}\mathrm{diag}(\mathbf{Z}_{\mathcal{I}}\hat{\bar{\mathbf{e}}}_{p})\mathbf{Z}_{\mathcal{I}}\mathbf{W}_{a}\in\mathbb{C}^{(|\mathcal{I}|-L)\times
L}$,
$\hat{\widetilde{\mathbf{h}}}_{p}=-\mathbf{B}^{H}\mathbf{Z}_{\mathcal{I}}\hat{\bar{\mathbf{e}}}_{p}\in\mathbb{C}^{(|\mathcal{I}|-L)}$,
$\hat{\widetilde{\mathbf{H}}}=\left[\hat{\widetilde{\mathbf{H}}}_{1}^{T}\cdots\hat{\widetilde{\mathbf{H}}}_{L}^{T}\right]^{T}\in\mathbb{C}^{L(|\mathcal{I}|-L)\times
L}$,
$\hat{\widetilde{\mathbf{h}}}=\left[\hat{\widetilde{\mathbf{h}}}_{1}^{T}\cdots\hat{\widetilde{\mathbf{h}}}_{L}^{T}\right]^{T}\in\mathbb{C}^{L(|\mathcal{I}|-L)}$
with the entries of $\mathbf{a}$ being the coefficients of the polynomial
presented in (28). Finally, the GLS solution of (60) is given by

$\displaystyle\hat{\mathbf{a}}=\left(\hat{\widetilde{\mathbf{H}}}^{H}\hat{\mathbf{W}}_{FB}\hat{\widetilde{\mathbf{H}}}\right)^{-1}\hat{\widetilde{\mathbf{H}}}^{H}\hat{\mathbf{W}}_{FB}\hat{\widetilde{\mathbf{h}}}$
(61)

where $\hat{\mathbf{W}}_{\rm FB}\approx\hat{\boldsymbol{\Pi}}_{\rm
s}^{2}\otimes\left(\mathbf{C}_{\rm FB}(\mathbf{a})\mathbf{C}_{\rm
FB}(\mathbf{a})^{H}\right)^{-1}$, $\mathbf{C}_{\rm
FB}(\mathbf{a})\triangleq\mathbf{B}^{H}\left(\mathbf{I}_{|\mathcal{I}|}+\mathrm{diag}\\{\mathbf{Z}_{\mathcal{I}}\mathbf{W}_{a}\mathbf{a}\\}\right)\mathbf{Z}_{\mathcal{I}}\mathbf{W}_{D}\hat{\widetilde{\mathbf{Q}}}^{\frac{1}{2}}\in\mathbb{C}^{(|\mathcal{I}|-L)\times
M}$.

After finding $\hat{\mathbf{a}}$ using (61), the $L$ DOA estimates
$\hat{\theta}_{l}$, $l=1,\cdots,L$ are obtained as
$\hat{\theta}_{l}=\mathrm{arcsin}\left(-\frac{\beta_{l}\lambda}{2\pi
d}\right)$ with $\beta_{l}$ denoting the phase argument of $\hat{\gamma}_{l}$.
In addition, $\hat{\gamma}_{l}$ denotes the $l$th root of the polynomial
defined as $\gamma^{L}+\sum_{l=1}^{L}[\hat{\mathbf{a}}]_{l}\ \gamma^{L-l}=0$.

Remark 4: Similar to Proposition 2, it can be readily shown that the
asymptotic variance of the DOA estimation for the FBA extension of the
proposed method in the high SNR regime is given as

$\displaystyle\mathbb{E}\\{\Delta\theta_{l}^{2}\\}$
$\displaystyle\approx\frac{1}{2}\left(\frac{\lambda}{2\pi
d\cos(\theta_{l})}\right)^{2}\frac{\boldsymbol{\gamma}_{l}^{T}\left(\widetilde{\mathbf{H}}^{H}\mathbf{W}_{\rm
FB}\widetilde{\mathbf{H}}\right)^{-1}\boldsymbol{\gamma}_{l}^{*}}{|\phi_{l}|^{2}}.$
(62)

##  V SIMULATION RESULTS

The aim of this section is to evaluate the performance of the proposed method
and compare it to that of the state-of-the-art algorithms in terms of diverse
numerical simulation examples especially for challenging scenarios. Our
examples address both uniform and nonuniform sensor noise cases. For the
uniform noise case, the performance of the forward-only and FBA versions of
the proposed method is compared with that of the unitary root-MUSIC method
[15], the root-swap unitary root-MUSIC method [22], the EPUMA method [23], and
the UE GLS [24]. For the nonuniform noise case, the “NISB+MUSIC” method [42]
and the “IMLSE+MUSIC” method [35] are used for comparison. In addition, for
achieving better DOA estimations, the combinations of both the NISB and IMLSE
with the root-MUSIC framework are considered for the nonuniform noise case.
The uniform stochastic CRB [28], and the nonuniform stochastic CRB [32] are
used as the benchmarks in the corresponding examples. The number of trials
used for calculating the RMSE is 2000 in all examples. If not further
specified, a ULA with $M=10$ sensors separated by half wavelength collecting
$N=10$ snapshots is considered for the uniform noise examples, while $M=8$ for
the nonuniform noise examples. The SNR is computed as SNR
$=\frac{1}{\sigma^{2}}$ for the uniform noise case, and as SNR
$=\frac{\sigma_{\rm s}^{2}}{M}\sum_{m=1}^{M}\frac{1}{\sigma_{m}^{2}}$ for the
nonuniform noise case. Here the powers of different sources are considered to
be identical and denoted by $\sigma_{\rm s}^{2}$.

In the first example, three uncorrelated sources located at
$\boldsymbol{\theta}=[19^{\circ},34^{\circ},36^{\circ}]$ are considered. It
can be seen in Fig. 2 that the SNR threshold performance of the proposed
method is outstandingly better than that of the other methods tested. In Fig.
3, the capability of different methods to cope with correlated sources is
evaluated by letting $\boldsymbol{\theta}=[5^{\circ},8^{\circ}]$ be correlated
with the correlation coefficient $\rho=0.95$. As expected, the performance of
the forward-only version of the proposed method degrades in relatively high
SNR region, while it is as good as that of the FBA version of the proposed
method in low SNRs. In addition, the FBA version of the proposed method
possesses much better performance compared to other methods. As illustrated in
Fig. 4, the impact of the number of snapshots on the performance of the
methods tested is investigated via setting
$\boldsymbol{\theta}=[10^{\circ},34^{\circ},36^{\circ}]$, $\rho=0$ for the
fixed SNR $=5$ dB. Similarly, Fig. 5 shows the impact of the number of
snapshots for a different setup of
$\boldsymbol{\theta}=[34^{\circ},38^{\circ}]$, $\rho=0.95$, and SNR $=2$ dB.
It displays that the FBA version of the proposed method provides robust
estimates even when the number of snapshots is about one order of magnitude
smaller than that of the other methods tested. In the next setup, the
capability of the methods tested to deal with the scenario of two closely
located sources is investigated. In doing so, we regard the setup in which
$\boldsymbol{\theta}=[0^{\circ},34^{\circ},(34+\Delta\theta)^{\circ}]$,
$\rho=0$, and SNR $=10$ dB with $\Delta\theta$ varying from $0.8^{\circ}$ to
$6^{\circ}$. It can be seen in Fig. 6 that the proposed method is more
reliable in dealing with smaller angular separations compared to the other
methods tested. The results for the same setup as Fig. 6, but with $\rho=0.95$
for the last two directions, and SNR $=15$ dB, are shown Fig.7. It can be seen
that the FBA version of the proposed method has the best performance.

Figure 2: RMSE vs. SNR for $L=3$ uncorrelated sources with
$\boldsymbol{\theta}=[19^{\circ},34^{\circ},36^{\circ}]$, $M=10$, and $N=10$.
Figure 3: RMSE vs. SNR for $L=2$ correlated sources with
$\boldsymbol{\theta}=[5^{\circ},8^{\circ}]$, $\rho=0.95$, $M=10$, and $N=10$.
Figure 4: RMSE vs. the number of snapshots for $L=3$ uncorrelated sources with
$\boldsymbol{\theta}=[10^{\circ},34^{\circ},36^{\circ}]$, SNR $=5$ dB, and
$M=10$. Figure 5: RMSE vs. the number of snapshots for $L=2$ correlated
sources with $\boldsymbol{\theta}=[34^{\circ},38^{\circ}]$, $\rho=0.95$, SNR
$=2$ dB, and $M=10$. Figure 6:  RMSE vs. the angular separation for $L=3$
uncorrelated sources with
$\boldsymbol{\theta}=[0^{\circ},34^{\circ},(34+\Delta\theta)^{\circ}]$, SNR
$=10$ dB, $M=10$, and $N=10$. Figure 7:  RMSE vs. the angular separation for
$L=3$ partly correlated sources with
$\boldsymbol{\theta}=[0^{\circ},34^{\circ},(34+\Delta\theta)^{\circ}]$, SNR
$=15$ dB, $\rho=0.95$ for the last two directions, $M=10$, and $N=10$.

To study the impact of the nonuniform noise, the sensor noise covariance
matrix is set as $\mathbf{Q}=\mathrm{diag}\\{[6,2,0.5,2.5,3,1,5.5,10]\\}$ for
the follow up examples [35]. The worst noise power ratio (WNPR) used in these
examples is given as

$\displaystyle\mathrm{WNPR}=\frac{\sigma_{max}^{2}}{\sigma_{min}^{2}}=\frac{10}{0.5}=20.$

First, Fig. 8 shows the RMSE performance of the methods tested versus SNR for
the setup of $\boldsymbol{\theta}=[33^{\circ},36^{\circ}]$ and $\rho=0$. As it
can be observed, the threshold performance of the proposed methods is
substantially better than that of the other methods tested. Moreover, in Fig.
9, the same scenario as that shown in Fig. 8 is considered for the case of
correlated sources $\boldsymbol{\theta}=[33^{\circ},38^{\circ}]$ and
$\rho=0.95$. Fig. 9 demonstrates the superiority of the FBA version of the
proposed method over the other methods tested. Fig. 10 depicts how different
methods perform depending on the number of snapshots for the scenario of
$\boldsymbol{\theta}=[33^{\circ},38^{\circ}]$, $\rho=0$, and SNR $=0$ dB. Fig.
11 illustrates this dependency also for the scenario of
$\boldsymbol{\theta}=[33^{\circ},48^{\circ}]$, $\rho=0.95$, and SNR $=-4$ dB.
Based on Figs. 10 and 11, it can be concluded that the reliability of the
proposed method to the scarcity of the number of snapshots is higher for the
case of uncorrelated sources compared to the other methods tested, while the
FBA version of the proposed method copes with the correlated sources more
efficiently. Finally, Figs. 12 and 13 show the strengths of the methods tested
against the presence of closely located sources for the cases of uncorrelated
and correlated signals, respectively. The setup regarded for Fig. 12 is
$\boldsymbol{\theta}=[-10^{\circ},34^{\circ},(34+\Delta\theta)^{\circ}]$,
$\rho=0$, and SNR $=15$ dB with $\Delta\theta$ varying from $1^{\circ}$ to
$12^{\circ}$. It can be observed that the performance of the proposed method
almost achieves the CRB. Fig. 13 depicts the results obtained from conducting
the same setup as for Fig. 12 with the difference that $\rho=0.95$ for the
last two directions. The superiority of the FBA version of the proposed method
over other methods tested can be seen to be very substantial.

Figure 8:  RMSE vs. SNR for $L=2$ uncorrelated sources with
$\boldsymbol{\theta}=[33^{\circ},36^{\circ}]$, $M=8$, and $N=10$. Figure 9:
RMSE vs. SNR for $L=2$ correlated sources with
$\boldsymbol{\theta}=[33^{\circ},38^{\circ}]$, $\rho=0.95$, $M=8$, and $N=10$.
Figure 10:  RMSE vs. the number of snapshots for $L=2$ uncorrelated sources
with $\boldsymbol{\theta}=[33^{\circ},38^{\circ}]$, SNR $=0$ dB, and $M=8$.
Figure 11:  RMSE vs. the number of snapshots for $L=2$ correlated sources with
$\boldsymbol{\theta}=[33^{\circ},48^{\circ}]$, $\rho=0.95$, SNR $=-4$ dB, and
$M=8$. Figure 12:  RMSE vs. the Angular Separation for $L=3$ uncorrelated
sources with
$\boldsymbol{\theta}=[-10^{\circ},34^{\circ},(34+\Delta\theta)^{\circ}]$, SNR
$=15$ dB, $M=8$, and $N=10$. Figure 13:  RMSE vs. the Angular Separation for
$L=3$ partly correlated sources with
$\boldsymbol{\theta}=[-10^{\circ},34^{\circ},(34+\Delta\theta)^{\circ}]$, SNR
$=15$ dB, $\rho=0.95$ for the last two directions, $M=8$, and $N=10$.

##  VI Conclusion

A new algorithm for DOA estimation in the presence of nonuniform sensor noise
is introduced. The proposed algorithm works in three phases. The essence of
the first phase is to estimate the nonuniform noise covariance matrix in an
iterative manner. In each iteration, the noise subspace is estimated using GED
first and then the noise covariance estimate is updated as the solution of an
LS minimization problem. The asymptotic performance of one iteration of such
algorithm is investigated. The advantage of the proposed noise covariance
estimator is that it is applicable to any array geometry. After the noise
covariance matrix is estimated, candidate DOAs are estimated using a rooting-
based DOA estimation method based on the combination of the GLS and the first-
order signal subspace perturbation. The noise covariance estimate obtained in
the first phase is used for pre-whitening the array signal with nonuniform
noise. In addition, the forward-only DOA estimation method is extended using
FBA. Furthermore, the asymptotic performance of both the forward-only and FBA
versions of the proposed method is studied. In the third phase, the final best
$L$ DOA estimates out of $2L$ DOA candidates generated in the second phase are
selected using properly designed DOA selection strategy. Numerical simulation
examples are included to show the superiority of the proposed algorithm
compared to the state-of-the-art approaches for both cases of the uniform and
nonuniform sensor noise.

## Appendix A

Taking into account the fact that $\mathbf{Q}$ is a real-valued diagonal
matrix and also using the derivative properties [56], the partial derivatives
of the terms related to $\mathbf{Q}$ in (10) can be found to be

$\displaystyle\frac{\partial\mathrm{trace}\left\\{\hat{\mathbf{R}}\hat{\mathbf{U}}\hat{\mathbf{U}}^{H}\mathbf{Q}\right\\}}{\partial\mathbf{Q}}$
$\displaystyle=\frac{\partial\mathrm{trace}\left\\{\mathcal{D}\left\\{\hat{\mathbf{R}}\hat{\mathbf{U}}\hat{\mathbf{U}}^{H}\right\\}\mathbf{Q}\right\\}}{\partial\mathbf{Q}}=\mathcal{D}\left\\{\hat{\mathbf{R}}\hat{\mathbf{U}}\hat{\mathbf{U}}^{H}\right\\}$
(63)
$\displaystyle\frac{\partial\mathrm{trace}\left\\{\hat{\mathbf{U}}\hat{\mathbf{U}}^{H}\hat{\mathbf{R}}\mathbf{Q}\right\\}}{\partial\mathbf{Q}}$
$\displaystyle=\frac{\partial\mathrm{trace}\left\\{\mathcal{D}\left\\{\hat{\mathbf{U}}\hat{\mathbf{U}}^{H}\hat{\mathbf{R}}\right\\}\mathbf{Q}\right\\}}{\partial\mathbf{Q}}=\mathcal{D}\left\\{\hat{\mathbf{U}}\hat{\mathbf{U}}^{H}\hat{\mathbf{R}}\right\\}$
(64)
$\displaystyle\frac{\partial\mathrm{trace}\left\\{\hat{\mathbf{U}}\hat{\mathbf{U}}^{H}\mathbf{Q}^{2}\right\\}}{\partial\mathbf{Q}}$
$\displaystyle=\frac{\partial\mathrm{trace}\left\\{\mathcal{D}\left\\{\hat{\mathbf{U}}\hat{\mathbf{U}}^{H}\right\\}\mathbf{Q}^{2}\right\\}}{\partial\mathbf{Q}}=2\mathcal{D}\left\\{\hat{\mathbf{U}}\hat{\mathbf{U}}^{H}\right\\}\mathbf{Q}.$
(65)

Using (63)–(65), the partial derivative of (10) with respect to $\mathbf{Q}$
can be straightforwardly found to be (11).

## Appendix B

Using (13), the $m$th diagonal entry of $\hat{\mathbf{Q}}$, denoted by
$\hat{\sigma}_{m}^{2}$, can be written as

$\displaystyle\hat{\sigma}_{m}^{2}=\frac{\left(\mathbf{v}_{m}^{H}\hat{\mathbf{r}}_{m}+(\mathbf{v}_{m}^{H}\hat{\mathbf{r}}_{m})^{H}\right)}{2\tau_{m}}=\frac{\mathfrak{R}\\{\mathbf{v}_{m}^{H}\hat{\mathbf{r}}_{m}\\}}{\tau_{m}}$
(66)

where $\tau_{m}\triangleq[\hat{\mathbf{U}}\hat{\mathbf{U}}^{H}]_{mm}$, and
$\mathbf{v}_{m}\triangleq[\hat{\mathbf{U}}\hat{\mathbf{U}}^{H}]_{:,m}$.
Expressing $\hat{\mathbf{r}}_{m}$ as
$\hat{\mathbf{r}}_{m}=\mathbf{r}_{m}+\boldsymbol{\Delta}\mathbf{r}_{m}$, where
$\boldsymbol{\Delta}\mathbf{r}_{m}$ denotes the estimation error of the $m$th
column of the SCM, it can be written that

$\displaystyle\Delta\sigma_{m}^{2}=\frac{\mathfrak{R}\\{\mathbf{v}_{m}^{H}\boldsymbol{\Delta}\mathbf{r}_{m}\\}}{\tau_{m}}$
(67)

where $\Delta\sigma_{m}^{2}$ is the difference between the actual
$\sigma_{m}^{2}$ and the estimate $\hat{\sigma}_{m}^{2}$, i.e.,
$\Delta\sigma_{m}^{2}=\hat{\sigma}_{m}^{2}-\sigma_{m}^{2}$. As a result, the
variance of $\Delta\sigma_{m}^{2}$ can be expressed as

$\displaystyle\mathbb{E}\left\\{\left(\Delta\sigma_{m}^{2}\right)^{2}\right\\}$
$\displaystyle=\frac{1}{4\tau_{m}^{2}}\mathbb{E}\bigg{\\{}\left(\mathbf{v}_{m}^{H}\boldsymbol{\Delta}\mathbf{r}_{m}+\mathbf{v}_{m}^{T}\boldsymbol{\Delta}\mathbf{r}_{m}^{*}\right)\left(\boldsymbol{\Delta}\mathbf{r}_{m}^{H}\mathbf{v}_{m}+\boldsymbol{\Delta}\mathbf{r}_{m}^{T}\mathbf{v}_{m}^{*}\right)\bigg{\\}}$
$\displaystyle=\frac{1}{4\tau_{m}^{2}}\bigg{(}\mathbf{v}_{m}^{H}\mathbb{E}\left\\{\boldsymbol{\Delta}\mathbf{r}_{m}\boldsymbol{\Delta}\mathbf{r}_{m}^{H}\right\\}\mathbf{v}_{m}+\mathbf{v}_{m}^{H}\mathbb{E}\left\\{\boldsymbol{\Delta}\mathbf{r}_{m}\boldsymbol{\Delta}\mathbf{r}_{m}^{T}\right\\}\mathbf{v}_{m}^{*}$
$\displaystyle+\mathbf{v}_{m}^{T}\mathbb{E}\left\\{\boldsymbol{\Delta}\mathbf{r}_{m}^{*}\boldsymbol{\Delta}\mathbf{r}_{m}^{H}\right\\}\mathbf{v}_{m}+\mathbf{v}_{m}^{T}\mathbb{E}\left\\{\boldsymbol{\Delta}\mathbf{r}_{m}^{*}\boldsymbol{\Delta}\mathbf{r}_{m}^{T}\right\\}\mathbf{v}_{m}^{*}\bigg{)}.$
(68)

According to [57], the asymptotic covariance and pseudo-covariance matrices of
the vector
$\boldsymbol{\Delta}\mathbf{r}\triangleq\mathrm{vec}\\{(\hat{\mathbf{R}}-\mathbf{R})\\}\in\mathbb{C}^{M^{2}}$
are

$\displaystyle\mathbb{E}\left\\{\boldsymbol{\Delta}\mathbf{r}\boldsymbol{\Delta}\mathbf{r}^{H}\right\\}$
$\displaystyle=\frac{1}{N}(\mathbf{R}^{T}\otimes\mathbf{R})$ (69)
$\displaystyle\mathbb{E}\left\\{\boldsymbol{\Delta}\mathbf{r}\boldsymbol{\Delta}\mathbf{r}^{T}\right\\}$
$\displaystyle=\mathbf{R}^{T}\otimes\mathbf{R}.$ (70)

Using (69) and (70), it is straightforward to show that [58]

$\displaystyle\mathbb{E}\left\\{\boldsymbol{\Delta}\mathbf{r}_{m}\boldsymbol{\Delta}\mathbf{r}_{m}^{H}\right\\}$
$\displaystyle=\left(\frac{[\mathbf{R}]_{mm}}{N}\right)\mathbf{R}$ (71)
$\displaystyle\mathbb{E}\left\\{\boldsymbol{\Delta}\mathbf{r}_{m}\boldsymbol{\Delta}\mathbf{r}_{m}^{T}\right\\}$
$\displaystyle=([\mathbf{R}]_{mm})\mathbf{R}.$ (72)

Plugging (71) and (72) into (68) yields

$\displaystyle\mathbb{E}\left\\{\left(\Delta\sigma_{m}^{2}\right)^{2}\right\\}$
$\displaystyle=\left(\frac{[\mathbf{R}]_{mm}}{2N\tau_{m}^{2}}\right)\mathfrak{R}\left\\{\mathbf{v}_{m}^{H}\mathbf{R}\left(\mathbf{v}_{m}+N\mathbf{v}_{m}^{*}\right)\right\\}$
(73)

which completes the proof.

## Appendix C

The proof goes in the same steps as that in [23] (Appendix A), and is included
for the sake of completeness. As $\theta_{l}$ and $\gamma_{l}$ are related to
each other as $\gamma_{l}=e^{-j2\pi d\mathrm{sin}(\theta_{l})/\lambda}$, we
perform Taylor’s expansion and keep only the terms containing up to the first-
order perturbation terms to obtain

$\displaystyle\Delta\theta_{l}\approx-\frac{\lambda}{2\pi
d\cos(\theta_{l})}\frac{\Delta\gamma_{l}}{j\gamma_{l}}.$ (74)

To enforce $\Delta\theta_{l}$ to be a real-valued quantity, it is reasonable
to define

$\displaystyle\Delta\theta_{l}$
$\displaystyle\triangleq\frac{1}{2}(\Delta\theta_{l}+\Delta\theta_{l}^{*})=\frac{1}{2}\frac{j\lambda}{2\pi
d\cos(\theta_{l})}(\gamma_{l}^{*}\Delta\gamma_{l}-\gamma_{l}\Delta\gamma_{l}^{*}).$
(75)

Using (75), the variance of $\Delta\theta_{l}$ can be written as

$\displaystyle\mathbb{E}\\{\Delta\theta_{l}^{2}\\}$
$\displaystyle\approx\frac{1}{2}\left(\frac{\lambda}{2\pi
d\cos(\theta_{l})}\right)^{2}\left(\mathbb{E}\\{|\Delta\gamma_{l}|^{2}\\}-\mathfrak{R}\\{\mathbb{E}\\{\Delta\gamma_{l}^{2}\\}(\gamma_{l}^{*})^{2}\\}\right).$
(76)

Next we need to find an expression that connects $\Delta\gamma_{l}$ and
$\boldsymbol{\Delta}\mathbf{a}\triangleq\hat{\mathbf{a}}-\mathbf{a}$, as (76)
is dependent to the statistics of $\Delta\gamma_{l}$ which are related to the
statistics of $\boldsymbol{\Delta}\mathbf{a}$. Towards this end, using the
first-order approximation of (28) when $\gamma$ is replaced by $\gamma_{l}$,
we obtain

$\displaystyle\boldsymbol{\gamma}_{l}^{T}\boldsymbol{\Delta}\mathbf{a}+\phi_{l}\Delta\gamma_{l}\approx
0$ (77)

where
$\boldsymbol{\Delta}\mathbf{a}\triangleq[\Delta[\mathbf{a}]_{1}\cdots\Delta[\mathbf{a}]_{L}]^{T}$,
$\boldsymbol{\gamma}_{l}\triangleq[\gamma_{l}^{L-1}\cdots 1]^{T}$,
$\phi_{l}\triangleq
L\gamma_{l}^{L-1}+(L-1)[\mathbf{a}]_{1}\gamma_{l}^{L-2}+\cdots+[\mathbf{a}]_{L-1}$.
From (77), we obtain
$\Delta\gamma_{l}\approx-\frac{\boldsymbol{\gamma}_{l}^{T}\mathbf{\Delta}\mathbf{a}}{\phi_{l}}$.
Thus, it can be written that

$\displaystyle\mathbb{E}\\{|\Delta\gamma_{l}|^{2}\\}\approx\frac{\boldsymbol{\gamma}_{l}^{T}\mathbb{E}\\{\boldsymbol{\Delta}\mathbf{a}\boldsymbol{\Delta}\mathbf{a}^{H}\\}\boldsymbol{\gamma}_{l}^{*}}{|\phi_{l}|^{2}}.$
(78)

Let us define

$\displaystyle
f(\mathbf{a})=(\hat{\mathbf{H}}\mathbf{a}-\hat{\mathbf{h}})^{H}\mathbf{W}(\hat{\mathbf{H}}\mathbf{a}-\hat{\mathbf{h}}).$
(79)

Since $\hat{\mathbf{a}}$ is the vector that minimizes (79),
$f^{\prime}(\hat{\mathbf{a}})$ can be approximated under the assumption of
high SNR using Taylor’s expansion as [23], [59]

$\displaystyle 0=f^{\prime}(\hat{\mathbf{a}})\approx
f^{\prime}(\mathbf{a})+f^{\prime\prime}(\mathbf{a})\boldsymbol{\Delta}\mathbf{a}$
(80)

where $f^{\prime}(\mathbf{a})$ and $f^{\prime\prime}(\mathbf{a})$ denote
respectively the first and second derivatives of $f(\mathbf{a})$ with respect
to $\mathbf{a}$, which are given as

$\displaystyle f^{\prime}(\mathbf{a})$
$\displaystyle=2\hat{\mathbf{H}}^{H}\mathbf{W}(\hat{\mathbf{H}}\mathbf{a}-\hat{\mathbf{h}})=2\hat{\mathbf{H}}^{H}\mathbf{W}\hat{\mathbf{e}}=2\hat{\mathbf{H}}^{H}\mathbf{W}\left(\mathbf{I}_{L}\otimes\mathbf{C}(\mathbf{a})\right)\boldsymbol{\Delta}\mathbf{u}_{s}$
(81) $\displaystyle f^{\prime\prime}(\mathbf{a})$
$\displaystyle=2\hat{\mathbf{H}}^{H}\mathbf{W}\hat{\mathbf{H}}.$ (82)

Combining (80), (81) and (82), we get for high SNR that

$\displaystyle\mathbb{E}$
$\displaystyle\\{\boldsymbol{\Delta}\mathbf{a}\boldsymbol{\Delta}\mathbf{a}^{H}\\}\approx(\mathbf{H}^{H}\mathbf{W}\mathbf{H})^{-1}\mathbf{H}^{H}\mathbf{W}\left(\mathbf{I}_{L}\otimes\mathbf{C}(\mathbf{a})\right)$
$\displaystyle\times\mathbb{E}\\{\boldsymbol{\Delta}\mathbf{u}_{\rm
s}\boldsymbol{\Delta}\mathbf{u}_{\rm
s}^{H}\\}\left(\mathbf{I}_{L}\otimes\mathbf{C}^{H}(\mathbf{a})\right)\mathbf{W}\mathbf{H}(\mathbf{H}^{H}\mathbf{W}\mathbf{H})^{-1}.$
(83)

Based on (42) and (83), it can be written that

$\displaystyle\mathbb{E}\\{\boldsymbol{\Delta}\mathbf{a}\boldsymbol{\Delta}\mathbf{a}^{H}\\}\approx(\hat{\mathbf{H}}^{H}\mathbf{W}\hat{\mathbf{H}})^{-1}.$
(84)

Consequently, substituting (84) into (78), we have

$\displaystyle\mathbb{E}\\{|\Delta\gamma_{l}|^{2}\\}\approx\frac{\boldsymbol{\gamma}_{l}^{T}(\hat{\mathbf{H}}^{H}\mathbf{W}\hat{\mathbf{H}})^{-1}\boldsymbol{\gamma}_{l}^{*}}{|\phi_{l}|^{2}}.$
(85)

The final part is to compute $\mathbb{E}\\{\Delta\gamma_{l}^{2}\\}$, which has
the following form

$\displaystyle\mathbb{E}\\{\Delta\gamma_{l}^{2}\\}\approx\frac{\boldsymbol{\gamma}_{l}^{T}\mathbb{E}\\{\boldsymbol{\Delta}\mathbf{a}\boldsymbol{\Delta}\mathbf{a}^{T}\\}\boldsymbol{\gamma}_{l}}{\phi_{l}^{2}}$
(86)

where

$\displaystyle\mathbb{E}\\{\boldsymbol{\Delta}\mathbf{a}\boldsymbol{\Delta}\mathbf{a}^{T}\\}$
$\displaystyle\approx(\mathbf{H}^{H}\mathbf{W}\mathbf{H})^{-1}\mathbf{H}^{H}\mathbf{W}\left(\mathbf{I}_{L}\otimes\mathbf{C}(\mathbf{a})\right)$
$\displaystyle\times\mathbb{E}\\{\boldsymbol{\Delta}\mathbf{u}_{\rm
s}\boldsymbol{\Delta}\mathbf{u}_{\rm
s}^{T}\\}\left(\mathbf{I}_{L}\otimes\mathbf{C}(\mathbf{a})^{T}\right)\mathbf{W}^{T}\mathbf{H}^{*}(\mathbf{H}^{H}\mathbf{W}\mathbf{H})^{-T}.$
(87)

It follows from (44) that

$\displaystyle\mathbb{E}\\{\boldsymbol{\Delta}\mathbf{u}_{\rm
s}\boldsymbol{\Delta}\mathbf{u}_{\rm
s}^{T}\\}\approx\left(\boldsymbol{\Sigma}_{\rm s}^{-1}\mathbf{V}_{\rm
s}^{T}\otimes(\mathbf{I}_{M}-\mathbf{U}_{\rm s}\mathbf{U}_{\rm
s}^{H})\right)\mathbb{E}\\{\bar{\mathbf{n}}\bar{\mathbf{n}}^{T}\\}\times\left(\mathbf{V}_{\rm
s}\boldsymbol{\Sigma}_{\rm s}^{-1}\otimes(\mathbf{I}_{M}-\mathbf{U}_{\rm
s}^{*}\mathbf{U}_{\rm s}^{T})\right).$ (88)

Since
$\mathbb{E}\\{\bar{\mathbf{n}}\bar{\mathbf{n}}^{T}\\}=\mathbf{0}_{MN\times
MN}$, (88) becomes a zero matrix which gives rise to

$\displaystyle\mathbb{E}\\{\Delta\gamma_{l}^{2}\\}\approx 0.$ (89)

As a result, by combining (76), (85) and (89), we obtain (50), which completes
the proof.

## References

  * [1] H. L. Van Trees, _Optimum array processing: Part IV of detection, estimation, and modulation theory_. John Wiley & Sons, 2004.
  * [2] P.-J. Chung, M. Viberg, and J. Yu, “Doa estimation methods and algorithms,” in _Academic Press Library in Signal Processing_. Elsevier, 2014, vol. 3, pp. 599–650.
  * [3] H. Krim and M. Viberg, “Two decades of array signal processing research: the parametric approach,” _IEEE Signal Process. Mag._ , vol. 13, no. 4, pp. 67–94, Jul. 1996.
  * [4] K. Upadhya and S. A. Vorobyov, “An array processing approach to pilot decontamination for massive mimo,” in _Proc. IEEE 6th Int. Workshop Comput. Adv. Multi-Sensor Adapt. Process., Cancun, Mexico_ , Dec. 2015, pp. 453–456.
  * [5] K. Upadhya, S. A. Vorobyov, and R. W. Heath, “Low-overhead receiver-side channel tracking for mmwave mimo,” in _Proc. IEEE Int. Conf. Acoust., Speech Signal Process., Alberta, Canada_ , Apr. 2018, pp. 3859–3863.
  * [6] A. Hassanien and S. A. Vorobyov, “Transmit energy focusing for doa estimation in mimo radar with colocated antennas,” _IEEE Trans. Signal Process._ , vol. 59, no. 6, pp. 2669–2682, Jun. 2011.
  * [7] A. Khabbazibasmenj, A. Hassanien, S. A. Vorobyov, and M. W. Morency, “Efficient transmit beamspace design for search-free based doa estimation in mimo radar,” _IEEE Trans. Signal Process._ , vol. 62, no. 6, pp. 1490–1500, Mar. 2014.
  * [8] J. Zhang, D. Rakhimov, and M. Haardt, “Gridless channel estimation for hybrid mmwave mimo systems via tensor-esprit algorithms in dft beamspace,” _IEEE J. Sel. Topics Signal Process._ , vol. 15, no. 3, pp. 816–831, Apr. 2021.
  * [9] K. Ardah, S. Gherekhloo, A. L. de Almeida, and M. Haardt, “Trice: A channel estimation framework for ris-aided millimeter-wave mimo systems,” _IEEE Signal Process. Lett._ , vol. 28, pp. 513–517, Feb. 2021.
  * [10] S. Vorobyov, A. Gershman, and W. K.M., “Maximum likelihood direction of arrival estimation in unknown noise fields using sparse sensor arrays,” _IEEE Trans. Signal Process._ , vol. 53, no. 1, pp. 34–43, Jan. 2005.
  * [11] R. Schmidt, “Multiple emitter location and signal parameter estimation,” _IEEE Trans. Antennas Propag._ , vol. 34, no. 3, pp. 276–280, Mar. 1986.
  * [12] J. M. Kim, O. K. Lee, and J. C. Ye, “Compressive music: Revisiting the link between compressive sensing and array signal processing,” _IEEE Trans. Inf. Theo._ , vol. 58, no. 1, pp. 278–301, Jan. 2012.
  * [13] P. Vallet, X. Mestre, and P. Loubaton, “Performance analysis of an improved music doa estimator,” _IEEE Trans. Signal Process._ , vol. 63, no. 23, pp. 6407–6422, Dec. 2015.
  * [14] A. Barabell, “Improving the resolution performance of eigenstructure-based direction-finding algorithms,” in _Proc. IEEE Int. Conf. Acoust., Speech Signal Process._ , vol. 8, 1983, pp. 336–339.
  * [15] M. Pesavento, A. B. Gershman, and M. Haardt, “Unitary root-music with a real-valued eigendecomposition: A theoretical and experimental performance study,” _IEEE Trans. Signal Process._ , vol. 48, no. 5, pp. 1306–1314, May 2000.
  * [16] C. Qian, L. Huang, and H.-C. So, “Improved unitary root-music for doa estimation based on pseudo-noise resampling,” _IEEE Signal Process. Lett._ , vol. 21, no. 2, pp. 140–144, Feb. 2013.
  * [17] R. Roy and T. Kailath, “Esprit-estimation of signal parameters via rotational invariance techniques,” _IEEE Trans. Acoust., Speech, Signal Process._ , vol. 37, no. 7, pp. 984–995, Jul. 1989.
  * [18] M. Haardt and J. A. Nossek, “Unitary esprit: How to obtain increased estimation accuracy with a reduced computational burden,” _IEEE Trans. Signal Process._ , vol. 43, no. 5, pp. 1232–1242, May 1995.
  * [19] P. Stoica and A. Nehorai, “Music, maximum likelihood, and cramer-rao bound,” _IEEE Trans. Acoust., Speech, Signal Process._ , vol. 37, no. 5, pp. 720–741, May 1989.
  * [20] P. Stoica and K. C. Sharman, “Maximum likelihood methods for direction-of-arrival estimation,” _IEEE Trans. Acoust., Speech, Signal Process._ , vol. 38, no. 7, pp. 1132–1143, Jul. 1990.
  * [21] P. Stoica and A. Nehorai, “Performance comparison of subspace rotation and music methods for direction estimation,” _IEEE Trans. Signal Process._ , vol. 39, no. 2, pp. 446–453, Feb. 1991.
  * [22] M. Shaghaghi and S. A. Vorobyov, “Subspace leakage analysis and improved doa estimation with small sample size,” _IEEE Trans. Signal Process._ , vol. 63, no. 12, pp. 3251–3265, Jun. 2015.
  * [23] C. Qian, L. Huang, N. D. Sidiropoulos, and H. C. So, “Enhanced puma for direction-of-arrival estimation and its performance analysis,” _IEEE Trans. Signal Process._ , vol. 64, no. 16, pp. 4127–4137, Aug. 2016.
  * [24] J. Steinwandt, F. Roemer, and M. Haardt, “Generalized least squares for esprit-type direction of arrival estimation,” _IEEE Signal Process. Lett._ , vol. 24, no. 11, pp. 1681–1685, Nov. 2017.
  * [25] M. Trinh-Hoang, M. Viberg, and M. Pesavento, “Partial relaxation approach: An eigenvalue-based doa estimator framework,” _IEEE Trans. Signal Process._ , vol. 66, no. 23, pp. 6190–6203, Dec. 2018.
  * [26] M. Trinh-Hoang, M. Viberg and M. Pesavento, “Cramér-rao bound for doa estimators under the partial relaxation framework: Derivation and comparison,” _IEEE Trans. Signal Process._ , vol. 68, pp. 3194–3208, May 2020.
  * [27] M. W. Morency, S. A. Vorobyov, and G. Leus, “Joint detection and localization of an unknown number of sources using the algebraic structure of the noise subspace,” _IEEE Trans. Signal Process._ , vol. 66, no. 17, pp. 4685–4700, Sep. 2018.
  * [28] P. Stoica and A. Nehorai, “Performance study of conditional and unconditional direction-of-arrival estimation,” _IEEE Trans. Acoust., Speech, Signal Process._ , vol. 38, no. 10, pp. 1783–1795, Oct. 1990.
  * [29] M. Esfandiari and S. A. Vorobyov, “Enhanced standard esprit for overcoming imperfections in doa estimation,” in _Proc. IEEE Int. Conf. Acoust., Speech Signal Process., Toronto, Canada_ , Jun. 2021, pp. 4375–4379.
  * [30] D. A. Linebarger, R. D. DeGroat, and E. M. Dowling, “Efficient direction-finding methods employing forward/backward averaging,” _IEEE Trans. Signal Process._ , vol. 42, no. 8, pp. 2136–2145, Aug. 1994.
  * [31] S. U. Pillai and B. H. Kwon, “Forward/backward spatial smoothing techniques for coherent signal identification,” _IEEE Trans. Acoust., Speech, Signal Process._ , vol. 37, no. 1, pp. 8–15, Jan. 1989.
  * [32] M. Pesavento and A. B. Gershman, “Maximum-likelihood direction-of-arrival estimation in the presence of unknown nonuniform noise,” _IEEE Trans. Signal Process._ , vol. 49, no. 7, pp. 1310–1324, Jul. 2001.
  * [33] C. E. Chen, F. Lorenzelli, R. E. Hudson, and K. Yao, “Stochastic maximum-likelihood doa estimation in the presence of unknown nonuniform noise,” _IEEE Trans. Signal Process._ , vol. 56, no. 7, pp. 3038–3044, Jul. 2008.
  * [34] D. Madurasinghe, “A new doa estimator in nonuniform noise,” _IEEE Signal Process. Lett._ , vol. 12, no. 4, pp. 337–339, Apr. 2005.
  * [35] B. Liao, S.-C. Chan, L. Huang, and C. Guo, “Iterative methods for subspace and doa estimation in nonuniform noise,” _IEEE Trans. Signal Process._ , vol. 64, no. 12, pp. 3008–3020, Jun. 2016.
  * [36] B. Liao, L. Huang, C. Guo, and H. C. So, “New approaches to direction-of-arrival estimation with sensor arrays in unknown nonuniform noise,” _IEEE Sensors J._ , vol. 16, no. 24, pp. 8982–8989, Dec. 2016.
  * [37] J. Wen, B. Liao, and C. Guo, “Spatial smoothing based methods for direction-of-arrival estimation of coherent signals in nonuniform noise,” _Digit. Signal Process._ , vol. 67, pp. 116–122, Aug. 2017.
  * [38] C. Qi, Z. Chen, Y. Wang, and Y. Zhang, “Doa estimation for coherent sources in unknown nonuniform noise fields,” _IEEE Trans. Aerosp. Electron. Syst._ , vol. 43, no. 3, pp. 1195–1204, Jul. 2007.
  * [39] P. Stoica, P. Babu, and J. Li, “Spice: A sparse covariance-based estimation method for array processing,” _IEEE Trans. Signal Process._ , vol. 59, no. 2, pp. 629–638, Feb. 2011.
  * [40] Z.-Q. He, Z.-P. Shi, and L. Huang, “Covariance sparsity-aware doa estimation for nonuniform noise,” _Digit. Signal Process._ , vol. 28, pp. 75–81, May 2014.
  * [41] H. Wang, X. Wang, L. Wan, and M. Huang, “Robust sparse bayesian learning for off-grid doa estimation with non-uniform noise,” _IEEE access_ , vol. 6, pp. 64 688–64 697, Nov. 2018.
  * [42] M. Esfandiari, S. A. Vorobyov, S. Alibani, and M. Karimi, “Non-iterative subspace-based doa estimation in the presence of nonuniform noise,” _IEEE Signal Process. Lett._ , vol. 26, no. 6, pp. 848–852, Jun. 2019.
  * [43] S. Provencher, “Parameters estimation of complex multitone signal in the dft domain,” _IEEE Trans. Signal Process._ , vol. 59, no. 7, pp. 3001–3012, Jul. 2011.
  * [44] T. Amemiya, “Generalized least squares theory,” _Advanced Econometrics_ , Cambridge, MA, USA: Harvard Univ. Press, 1985.
  * [45] T. Kariya and H. Kurata, _Generalized Least Squares_. John Wiley & Sons, 2004.
  * [46] F. Li, H. Liu, and R. J. Vaccaro, “Performance analysis for doa estimation algorithms: unification, simplification, and observations,” _IEEE Trans. Aerosp. Electron. Syst._ , vol. 29, no. 4, pp. 1170–1184, Oct. 1993.
  * [47] A. B. Gershman, “Pseudo-randomly generated estimator banks: A new tool for improving the threshold performance of direction finding,” _IEEE Trans. Signal Process._ , vol. 46, no. 5, pp. 1351–1364, May 1998.
  * [48] V. Vasylyshyn, “Improved beamspace esprit-based doa estimation via pseudo-noise resampling,” in _Proc. IEEE 9th Europ. Radar Conf., Amsterdam, Netherlands_ , 2012, pp. 238–241.
  * [49] V. Vasylyshyn, “Removing the outliers in root-music via pseudo-noise resampling and conventional beamformer,” _Signal Process._ , vol. 93, no. 12, pp. 3423–3429, Dec. 2013.
  * [50] A. B. Gershman and P. Stoica, “New mode-based techniques for direction finding with an improved threshold performance,” _Signal Process._ , vol. 76, no. 3, pp. 221–235, Aug. 1999.
  * [51] F. Izedi, M. Karimi, and M. Derakhtian, “Joint doa estimation and source number detection for arrays with arbitrary geometry,” _Signal Process._ , vol. 140, pp. 149–160, Nov. 2017.
  * [52] S. M. Kay, _Fundamentals of statistical signal processing_. Prentice Hall PTR, 1993, Chapter 8.
  * [53] K.-C. Huarng and C.-C. Yeh, “A unitary transformation method for angle-of-arrival estimation,” _IEEE Trans. Signal Process._ , vol. 39, no. 4, pp. 975–977, Apr. 1991.
  * [54] A. B. Gershman and P. Stoica, “On unitary and forward–backward mode,” _Digit. Signal Process._ , vol. 9, no. 2, pp. 67–75, Apr. 1999.
  * [55] B. D. Rao and K. Hari, “Weighted subspace methods and spatial smoothing: analysis and comparison,” _IEEE Trans. Signal Process._ , vol. 41, no. 2, pp. 788–803, Feb. 1993.
  * [56] K. Petersen and M. Pedersen, “The matrix cookbook, version 20121115,” _Technical Univ. Denmark, Kongens Lyngby, Denmark, Tech. Rep_ , vol. 3274, 2012.
  * [57] M. Bilodeau and D. Brenner, _Theory of multivariate statistics_. Springer Science & Business Media, 2008.
  * [58] B. Ottersten, P. Stoica, and R. Roy, “Covariance matching estimation techniques for array signal processing applications,” _Digit. Signal Process._ , vol. 8, no. 3, pp. 185–210, Jul. 1998.
  * [59] H. C. So, Y. T. Chan, K. Ho, and Y. Chen, “Simple formulae for bias and mean square error computation [dsp tips and tricks],” _IEEE Signal Process. Mag._ , vol. 30, no. 4, pp. 162–165, Jul. 2013.