Date of publication xxxx 00, 0000, date of current version xxxx 00, 0000.
10.1109/ACCESS.2023.DOI

Corresponding author: Md Imrul Hasan<EMAIL_ADDRESS>and Mohammad
Saquib<EMAIL_ADDRESS>
# Low Complexity Iterative 2D DOA Estimation in MIMO Systems

MD IMRUL HASAN1  MOHAMMAD SAQUIB1  Department of Electrical Engineering, The
University of Texas at Dallas, Richardson, TX 75080 USA

###### Abstract

Multiple-input multiple-output (MIMO) systems play an essential role in
direction-of-arrival (DOA) estimation. A large number of antennas used in a
MIMO system imposes a huge complexity burden on the popular DOA estimation
algorithms, such as MUSIC and ESPRIT due to the implementation of eigenvalue
decomposition. This renders those algorithms impractical in applications
requiring quick DOA estimation. Consequently, we theoretically derive several
useful noise subspace vectors when the number of signal sources is less than
the number of elements in both the transmitter and receiver sides. Those noise
subspace vectors are then utilized to formulate a 2D-constrained minimization
problem, solved iteratively to obtain the DOAs of all the sources in a scene.
The convergence of the proposed iterative algorithm has been mathematically as
well as numerically demonstrated. Depending on the number of iterations, our
algorithm can provide significant complexity gain over the existing high-
resolution 2D DOA estimation algorithms in MIMO systems, such as MUSIC, while
exhibiting comparable performance for a moderate to high signal-to-noise ratio
(SNR).

###### Index Terms:

MIMO, Non-signal subspace modeling, Iterative algorithm, Low complexity, DOA
estimation.

=-15pt

##  I Introduction

Direction of arrival (DOA) is an important concept in a wide range of areas,
such as wireless communications, sonar, radar, geophysics, defense operations,
seismology, etc [1, 2, 3, 4, 5, 6]. In these applications, DOA estimation is
used, 1) to guide the beam in order to improve signal quality, suppress
interference, and reduce the power for transmission, 2) to localize a target,
which is critical for navigation, communication, surveillance, and tracking,
3) to detect the origin of seismic waves, which allows the exploration of the
Earth’s interior, predicting earthquakes, and so on.

Multiple Input Multiple Output (MIMO) technologies play an essential role in
DOA estimation by allowing diversity, a high spatial multiplexing gain, and
improved spatial resolution of array elements. Since a MIMO system utilizes
many antennas collected in one panel and permits numerous paths for the signal
to reach the receiver, it helps to reduce the impact of fading and
interference. Another advantage of a MIMO system is the formulation of the
virtual array [7, 8, 9], which describes the practical spatial locations from
which the target reflectivity gets sampled. That means if the monostatic
transmitter and receiver elements operate independently at each virtual array
element location, those will provide the exact phase measurements as the
physical configuration. The virtual array in a MIMO system enhances the
original aperture size, which as a result, increases the degrees of freedom
(DOF) [10]. These attributes offer high accuracy and robustness in the DOA
estimation with a MIMO system compared to a common phased array system.

In the literature, extensive studies have been performed, and numerous
algorithms have been developed to estimate the DOAs, i.e. multiple signal
classification (MUSIC) [11], maximum-likelihood (ML) [12], Capon [13],
estimation of signal parameters via rotational invariance techniques (ESPRIT)
[14], Min-Norm [15], etc. Among all these methods, traditional subspace-based
estimation algorithms, such as MUSIC, and ESPRIT are especially noteworthy due
to their high precision and resolution in DOA estimation[16, 17]. This noise
subspace-based estimation (NISE) algorithm [16] iteratively solves for each
source’s elevation angle, potentially yielding (depending on the number of
iterations) lower complexity than existing DoA estimation algorithms, such as
Fast Root-MUSIC (FRM) [18].

MIMO systems employ 2D antenna elements due to their 3D beamforming benefits.
Therefore, the 2D DOA estimation in a MIMO system is of great interest. For
the past few years, a number of 2D DOA estimation techniques have been
introduced based on MUSIC and ESPRIT algorithms. For example, a MUSIC-based
algorithm is proposed in [19] which estimates the DOAs in two stages, but
still requires eigenvalue decomposition and peak searching to estimate the
angles. A DOA estimation algorithm based on 2D MUSIC in MIMO systems is
proposed in [20], however, it uses spatial and eigenvalue decomposition. A 2D
unitary ESPRIT algorithm is introduced in [21], but it fails to achieve a high
precision because of the constraint concerning the array aperture. Moreover,
the eigenvalue decomposition (EVD) of the covariance matrix needs to be
computed. This leads to a high computational complexity while using large-
scale antenna arrays since the array in MIMO systems can contain hundreds of
elements[17]. Moreover, as the electromagnetic environment is getting
increasingly complex every day, the computational burden of these algorithms
in MIMO systems poses an issue of great concern, especially in a real-time
environment where a quick estimation of the DOAs is required [22, 23, 24, 25,
26]. The low complexity NISE algorithm [16] is capable of quickly estimating
1D DOA employing a uniform linear array (ULA). Unfortunately, it is not
suitable for a MIMO system.

In light of the above, we seek a low-complexity quick 2D DOA estimation
technique. For a MIMO system, closed-form expressions of several non-signal
subspace eigenvectors are derived using the algebraic expressions of the non-
signal subspace developed in [16] for a ULA. Those eigenvectors contain
information on all the DOAs in a scene and are utilized to formulate a
2D-constrained minimization problem. We solve this problem using an iterative
2D DOA estimation algorithm, namely iDEA. Our analysis shows that in a noise-
free scenario, as the number of iterations approaches infinity iDEA’s
estimates approach the true values of the DOAs. Since iDEA is an iterative
algorithm, it allows performing a trade-off between the performance and
complexity by controlling the number of iterations. An extensive numerical
study is carried out to investigate the performance of iDEA against the 2D
MUSIC algorithm. This study also includes a thorough complexity comparison
between those two.

The rest of the paper is organized as follows: Section II discusses the System
Model for MIMO systems. In Section III, eigenvectors of the non-signal
subspace for MIMO systems are formulated and later demonstrated how these
vectors could be utilized to determine the ${K}^{\underline{\text{th}}}$
source, given we have prior information of the other $K-1$ sources in the
scene. Section IV describes the iDEA algorithm for MIMO systems to estimate
all the $K$ sources sequentially while using an arbitrary number of antennas.
In Section V, the convergence proof of the algorithm is presented, and a
thorough performance analysis of iDEA w.r.t. the 2D MUSIC algorithm is
performed for a MIMO system in Section VI. Finally, Section VII contains the
concluding remarks of this paper.

Notations: We use lowercase and uppercase bold letters to denote vectors and
matrices, respectively. Lowercase letters in italics are used to represent
scalars. The notation * refers to the complex conjugate of a scalar,
$[\cdot]^{T}$ refers to transpose, and $[\cdot]^{\dagger}$ denotes the
Hermitian of a matrix. In this paper, $\otimes$ is used to portray the
Kronecker product.

Figure 1:  2D DOA estimation in MIMO system.

Properties of Kronecker Products (PKP) : In this paper, the following
properties of the Kronecker Products[27] are used to formulate and prove the
convergence of the proposed algorithm, iDEA.

  * •

PKP 1:
$\left({\mathbf{A}\otimes\mathbf{B}}\right)^{\dagger}=\mathbf{A}^{\dagger}\otimes\mathbf{B}^{\dagger}$

  * •

PKP 2: $k\mathbf{{A}}\otimes\mathbf{{B}}$ = $\mathbf{{A}}\otimes
k\mathbf{{B}}=k\left({\mathbf{{A}}\otimes\mathbf{{B}}}\right)$

  * •

PKP 3:

    * –

$\mathbf{{A}}\otimes\left({\mathbf{{B}}+\mathbf{{C}}}\right)=\mathbf{{A}}\otimes\mathbf{{B}}+\mathbf{{A}}\otimes\mathbf{{C}}$

    * –

$\left({\mathbf{{A}}+\mathbf{{B}}}\right)\otimes\mathbf{{C}}=\mathbf{{A}}\otimes\mathbf{{C}}+\mathbf{{B}}\otimes\mathbf{{C}}$

  * •

PKP 4:
$\left({\mathbf{A}\otimes\mathbf{B})(\mathbf{C}\otimes\mathbf{D}}\right)$ =
$\mathbf{(AC)}\otimes\mathbf{(BD)}$

##  II System Model

Let us consider a MIMO system, where $N_{\mathrm{tx}}$ transmitters are placed
linearly along the $\mathrm{x}-$ axis, and $N_{\mathrm{rx}}$ receivers are
positioned linearly at an angle of $\phi_{\mathrm{trx}}$ with the
transmitters; see Fig. 1. At the transmit site, $N_{\mathrm{tx}}$ different
narrowband waveforms are emitted simultaneously, which have identical
bandwidth and center frequency but are temporally orthogonal. In each
receiver, the signals are processed for all of the transmitted waveforms.
Assume the number of noncoherent far-field targets in the scene, $K$ is
previously known [28], and these signals impinge upon the MIMO with wavelength
$\lambda$, where $K<N_{\mathrm{tx}}$ and $N_{\mathrm{rx}}$. The azimuthal
angle of the $k^{\underline{\text{th}}}$ user,
$0^{\circ}\leq\phi_{k}<360^{\circ}$ is measured counterclockwise from the
$\mathrm{x}-$ axis, and its elevation angle, $0^{\circ}\leq\theta_{k}\leq
90^{\circ}$ is measured downward from the $\mathrm{z}-$ axis. Let us denote
the steering matrix for the transmitters, and the receivers as
$\mathbf{A}_{\mathrm{tx}}$, and $\mathbf{A}_{\mathrm{rx}}$, respectively.
Therefore, for $K$ number of targets, we can write

$\mathbf{A}_{\mathrm{tx}}={\left[{\mathbf{a}_{\mathrm{tx},1},\ldots,\mathbf{a}_{\mathrm{tx},K}}\right]},$

and

$\mathbf{A}_{\mathrm{rx}}={\left[{\mathbf{a}_{\mathrm{rx},1},\ldots,\mathbf{a}_{\mathrm{rx},K}}\right]}.$

Given the geometry, and the uniform array element spacing $d\leq\lambda/2$,
the ${i}^{\underline{\text{th}}}$ element of the steering vector
$\mathbf{a}_{\mathrm{tx},k}$,

$a_{\mathrm{tx},i,k}=\exp\\{j\times
i\epsilon_{\mathrm{tx},k}\\};\,i=1,\ldots,N_{\mathrm{tx}}\,\,k=1,\ldots,K,$
(1)

and the ${l}^{\underline{\text{th}}}$ element of the steering vector
$\mathbf{a}_{\mathrm{rx},k}$,

$a_{\mathrm{rx},l,k}=\exp\\{j\times
l\epsilon_{\mathrm{rx},k}\\};\,l=1,\ldots,N_{\mathrm{rx}}\,\,k=1,\ldots,K,$
(2)

where

$\epsilon_{\mathrm{tx},k}=-\left({2\pi
d/\lambda}\right)\sin\theta_{k}\cos\phi_{k},$
$\epsilon_{\mathrm{rx},k}=-\left({2\pi
d/\lambda}\right)\sin\theta_{k}\cos(\phi_{k}-\phi_{\mathrm{trx}}).$

Now, the virtual steering vector for the MIMO is given by the following
relation [7]:

$\mathbf{A}={\left[{\mathbf{A}_{1},\ldots,\mathbf{A}_{K}}\right]}\,,$ (3)

where

$\mathbf{A}_{k}=\mathbf{a}_{\mathrm{tx},k}\otimes\mathbf{a}_{\mathrm{rx},k}$

is the contribution from the ${k}^{\underline{\text{th}}}$ user. Therefore, at
the $m^{\underline{\mathrm{th}}}$ snapshot, the outputs of all the matched
filters in all the receivers can be expressed as [29]

$\mathbf{x}(m)=\mathbf{A}\mathbf{s}(m)+\mathbf{w}(m),$ (4)

where $\mathbf{x}(m)\in C^{N\times 1}$, $\mathbf{A}\in C^{N\times K}$,
$\mathbf{s}(m)\in C^{K\times 1}$, and $\mathbf{w}(m)\in C^{N\times 1}$. Here,
$\mathbf{s}(m)$ is a column vector consisting of the phases and amplitudes of
the $K$ sources, vector $\mathbf{w}$(m) contains complex additive white
Gaussian noise (AWGN) samples with average power $\sigma^{2}/2$ per dimension
and spatially and temporally independent of $\mathbf{s}(m)$, and by
definition, $N=N_{\mathrm{tx}}N_{\mathrm{rx}}$. In the next section, we will
develop the model for the eigenvectors in the non-signal subspace of the auto-
correlation matrix $\mathbf{R}$ based on the received signal vector (4).

##  III Non-Signal Subspace Analysis

The auto-correlation matrix,

$\displaystyle\mathbf{R}$
$\displaystyle=\operatorname{\mathbb{E}}\\{\mathbf{{x}}(m)\mathbf{{x}}(m)^{\dagger}\\}$
$\displaystyle=\mathbf{A}\operatorname{\mathbb{E}}\\{\mathbf{s}(m)\mathbf{s}(m)^{\dagger}\\}\mathbf{A}^{\dagger}+\operatorname{\mathbb{E}}\\{\mathbf{w}(m)\mathbf{w}(m)^{\dagger}\\}$
$\displaystyle=\mathbf{A}\mathbf{R}_{s}\mathbf{A}^{\dagger}+\mathbf{R}_{n}\,.$
(5)

In the above equation, for uncorrelated sources, the signal correlation
matrix,
$\mathbf{R}_{\mathrm{s}}=\mathrm{diag}[\sigma_{1}^{2},\ldots,\sigma_{K}^{2}]$,
where $\sigma_{k}^{2}$ denotes the ${k}^{\underline{\text{th}}}$ source power,
and the noise covariance matrix, $\mathbf{R}_{n}=\sigma^{2}\mathbf{I}$. In
practice, the auto-correlation matrix $\mathbf{R}$ will be replaced by the
sample auto-correlation matrix averaged over $M$ time samples (or snapshots).
Since $K$ sources are considered, the largest $K$ eigenvalues of this auto-
correlation matrix and their eigenvectors are related to the signal subspace.
The other $N_{\mathrm{tx}}N_{\mathrm{rx}}-K$ eigenvalues and eigenvectors
correspond to the non-signal (noise) subspace. In the development of the
iterative DOA estimation algorithm, we assume a noise-free situation (i.e.,
$\sigma^{2}=0$). Later, the performance of the developed algorithm will be
demonstrated in a noisy environment considering two different cases, where 1)
there is prior knowledge of the noise variance, and 2) the noise variance of
the system is totally unknown.

Let us consider $N_{\mathrm{tx}}=N_{\mathrm{rx}}=K+1$. By using the non-signal
subspace modeling in a uniform linear array [16], we can write, for $K$
sources, the non-signal subspace of the linear array along the axis
representing the transmitter,

$\displaystyle\mathbf{E}_{\mathrm{tx},K}$
$\displaystyle=\begin{bmatrix}\mathbf{E}_{\mathrm{tx},K-1}\\\
0\end{bmatrix}-\exp(j\epsilon_{\mathrm{tx},K})\begin{bmatrix}0\\\
\mathbf{E}_{\mathrm{tx},K-1}\end{bmatrix}$
$\displaystyle=\mathbf{\tilde{E}}_{\mathrm{tx},K-1}-\exp(j\epsilon_{\mathrm{tx},K})\mathbf{L}\mathbf{\tilde{E}}_{\mathrm{tx},K-1},$
(6)

where

$\mathbf{\tilde{E}}_{\mathrm{tx},K-1}=\begin{bmatrix}\mathbf{E}_{\mathrm{tx},K-1}\\\
0\end{bmatrix},\mathbf{L}=\begin{bmatrix}\mathbf{0}_{1\times k-1}&1\\\
\mathbf{I}_{K-1\times K-1}&\mathbf{0}^{T}_{1\times k-1}\end{bmatrix},$

and

$\mathbf{E}_{\mathrm{tx},1}=\begin{bmatrix}1\\\
-\exp(j\epsilon_{\mathrm{tx},1})\end{bmatrix}.$

The above equations show how one can recursively construct the non-signal
subspace due to the ${K}^{\underline{\text{th}}}$ source from that of the
${(K-1)}^{\underline{\text{th}}}$ source. Similarly, the non-signal subspace
of the linear array along the axis representing the receiver,

$\displaystyle\mathbf{E}_{\mathrm{rx},K}$
$\displaystyle=\begin{bmatrix}\mathbf{E}_{\mathrm{rx},K-1}\\\
0\end{bmatrix}-\exp(j\epsilon_{\mathrm{rx},K})\begin{bmatrix}0\\\
\mathbf{E}_{\mathrm{rx},K-1}\end{bmatrix}$
$\displaystyle=\mathbf{\tilde{E}}_{\mathrm{rx},K-1}-\exp(j\epsilon_{\mathrm{rx},K})\mathbf{L}\mathbf{\tilde{E}}_{\mathrm{rx},K-1},$
(7)

where

$\mathbf{\tilde{E}}_{\mathrm{rx},K-1}=\begin{bmatrix}\mathbf{E}_{\mathrm{rx},K-1}\\\
0\end{bmatrix},\quad\mathbf{E}_{\mathrm{rx},1}=\begin{bmatrix}1\\\
-\exp(j\epsilon_{\mathrm{rx},1})\end{bmatrix}.$

Next, we prove a theorem to find a vector belonging to the non-signal subspace
of the auto-correlation matrix $\mathbf{R}$ using (6) and (7).

######  Theorem 1.

When there are $K$ sources and the number of antenna elements at the
transmitter and receiver satisfies $N_{\mathrm{tx}}=N_{\mathrm{rx}}=K+1$,

$\mathbf{E}_{K}=\mathbf{E}_{\mathrm{tx},K}\otimes\mathbf{E}_{\mathrm{rx},K}$

is a non-signal subspace vector of the auto-correlation matrix
$\mathbf{R}=\operatorname{\mathbb{E}}\\{\mathbf{A}\mathbf{s}(m)\mathbf{s}^{\dagger}(m)\mathbf{A}^{\dagger}\\}$.

###### Proof:

Recall that vectors $\mathbf{E}_{\mathrm{tx},K}$ and
$\mathbf{E}_{\mathrm{rx},K}$ are the non-signal subspace of auto-correlation
matrix
$\mathbf{R}_{\mathrm{tx}}=\operatorname{\mathbb{E}}\\{\mathbf{A}_{\mathrm{tx}}\mathbf{s}(m)\mathbf{s}^{\dagger}(m)\mathbf{A}^{\dagger}_{\mathrm{tx}}\\}$,
and
$\mathbf{R}_{\mathrm{rx}}=\operatorname{\mathbb{E}}\\{\mathbf{A}_{\mathrm{rx}}\mathbf{s}(m)\mathbf{s}^{\dagger}(m)\mathbf{A}^{\dagger}_{\mathrm{rx}}\\}$,
respectively. By definition,

$\mathbf{A}_{\mathrm{tx}}^{\dagger}\mathbf{E}_{\mathrm{tx},K}={\mathbf{0}_{K\times
1}},$

and

$\mathbf{A}_{\mathrm{rx}}^{\dagger}\mathbf{E}_{\mathrm{rx},K}={\mathbf{0}_{K\times
1}}.$

The above equations yield,

$\mathbf{a}_{\mathrm{tx,k}}^{\dagger}\mathbf{E}_{\mathrm{tx},K}=\mathbf{a}_{\mathrm{rx,k}}^{\dagger}\mathbf{E}_{\mathrm{rx},K}=0,$
(8)

where $k=1,\ldots,K.$ Note that
$\mathbf{E}_{K}=\mathbf{E}_{\mathrm{tx},K}\otimes\mathbf{E}_{\mathrm{rx},K}$
will be the non-signal subspace vector of $\mathbf{R}$, if

$\mathbf{A}^{\dagger}\mathbf{E}_{K}={\mathbf{0}_{{K\times 1}}}.$

First, we use PKP 1 to get

$\displaystyle\mathbf{A}^{\dagger}\mathbf{E}_{K}$
$\displaystyle=\begin{bmatrix}\mathbf{a}^{\dagger}_{\mathrm{tx,1}}\otimes\mathbf{a}^{\dagger}_{\mathrm{rx,1}}\\\
\vdots\\\
\mathbf{a}^{\dagger}_{\mathrm{tx,K}}\otimes\mathbf{a}^{\dagger}_{\mathrm{rx,K}}\end{bmatrix}(\mathbf{E}_{\mathrm{tx},K}\otimes\mathbf{E}_{\mathrm{rx},K}),$

and then we apply PKP 4 to write

$\displaystyle\mathbf{A}^{\dagger}\mathbf{E}_{K}$
$\displaystyle=\begin{bmatrix}\mathbf{a}^{\dagger}_{\mathrm{tx,1}}\mathbf{E}_{\mathrm{tx},K}\otimes\mathbf{a}^{\dagger}_{\mathrm{rx,1}}\mathbf{E}_{\mathrm{tx},K}\\\
\vdots\\\
\mathbf{a}^{\dagger}_{\mathrm{tx,K}}\mathbf{E}_{\mathrm{tx},K}\otimes\mathbf{a}^{\dagger}_{\mathrm{rx,K}}\mathbf{\tilde{E}}_{\mathrm{rx},K}\end{bmatrix}$
$\displaystyle={\mathbf{0}_{K\times 1}},$

where the equality is due to (8) and we prove the theorem. ∎

In the next section, we will demonstrate how the non-signal subspace vector
$\mathbf{E}_{K}$ can be utilized to estimate DOAs of $K$ sources iteratively.
This demonstration is performed using the following lemma on the property of
$\mathbf{E}_{K}$ which is a straightforward application of PKP 3 with (6) and
(7).

######  Lemma 1.

Since the non-signal subspace vector

$\mathbf{E}_{K}=\mathbf{E}_{\mathrm{tx},K}\otimes\mathbf{E}_{\mathrm{rx},K}\,,$

we have

$\displaystyle\mathbf{E}_{K}$
$\displaystyle=\mathbf{u}_{0}-\mathbf{u}_{1}\exp(j\epsilon_{\mathrm{rx},K})-\mathbf{u}_{2}\exp(j\epsilon_{\mathrm{tx},K})$
$\displaystyle+\mathbf{u}_{3}\exp(j\epsilon_{\mathrm{tx},K})\exp(j\epsilon_{\mathrm{rx},K}),$
(9)

where

$\displaystyle\mathbf{u}_{0}$
$\displaystyle=\mathbf{\tilde{E}}_{\mathrm{tx},K-1}\otimes\mathbf{\tilde{E}}_{\mathrm{rx},K-1}\,,$
(10) $\displaystyle\mathbf{u}_{1}$
$\displaystyle=\mathbf{\tilde{E}}_{\mathrm{tx},K-1}\otimes(\mathbf{L}\mathbf{\tilde{E}}_{\mathrm{rx},K-1})\,,$
(11) $\displaystyle\mathbf{u}_{2}$
$\displaystyle=(\mathbf{L}\mathbf{\tilde{E}}_{\mathrm{tx},K-1})\otimes\mathbf{\tilde{E}}_{\mathrm{rx},K-1}\,,$
(12)

and

$\displaystyle\mathbf{u}_{3}$
$\displaystyle=(\mathbf{L}\mathbf{\tilde{E}}_{\mathrm{tx},K-1})\otimes(\mathbf{L}\mathbf{\tilde{E}}_{\mathrm{rx},K-1})\,.$
(13)

Next, we prove a key property of the non-signal subspace vector
$\mathbf{E}_{K}$ that will be used later to develop the proposed iterative DOA
estimation algorithm.

######  Theorem 2.

When there are $K$ sources and the number of antenna elements at the
transmitter and receiver satisfies $N_{\mathrm{tx}}=N_{\mathrm{rx}}=K+1$,

$\mathbf{u}_{1}-\exp(j\epsilon_{\mathrm{tx},K})\mathbf{u}_{3}\,,$
$\mathbf{u}_{2}-\exp(j\epsilon_{\mathrm{rx},K})\mathbf{u}_{3}\,,$

and

$\mathbf{u}_{0}-\exp(j\epsilon_{\mathrm{tx},K})\exp(j\epsilon_{\mathrm{rx},K})\mathbf{u}_{3}$

belong to the non-signal subspace.

###### Proof:

First, we will show that vector
$\\{\mathbf{u}_{1}-\exp(j\epsilon_{\mathrm{tx},K})\mathbf{u}_{3}\\}$ belongs
to the non-signal subspace. In order to claim that, we need to prove

$\displaystyle\mathbf{R}\\{\mathbf{u}_{1}-\exp(j\epsilon_{\mathrm{tx},K})\mathbf{u}_{3}\\}={\mathbf{0}_{N_{\mathrm{tx}}N_{\mathrm{rx}}\times
1}}$ (14)

Applying PKP 1 and PKP 2 to (11) and (13), we write

$\displaystyle\mathbf{u}_{1}-\exp(j\epsilon_{\mathrm{tx},K})\mathbf{u}_{3}$
$\displaystyle=\\{\mathbf{\tilde{E}}_{\mathrm{tx},K-1}-\exp(j\epsilon_{\mathrm{tx},K})(\mathbf{L}\mathbf{\tilde{E}}_{\mathrm{tx},K-1})\\}\otimes(\mathbf{L}\mathbf{\tilde{E}}_{\mathrm{rx},K-1})\,,$

and then use (6) to get

$\displaystyle\mathbf{u}_{1}-\exp(j\epsilon_{\mathrm{tx},K})\mathbf{u}_{3}=\mathbf{E}_{\mathrm{tx},K}\otimes(\mathbf{L}\mathbf{\tilde{E}}_{\mathrm{rx},K-1})\,.$

In a noiseless scenario, substituting the value of the array auto-correlation
matrix $\mathbf{R}$ from (5), we write

$\displaystyle\mathbf{R}\\{\mathbf{u}_{1}-\exp(j\epsilon_{\mathrm{tx},K})\mathbf{u}_{3}\\}=\mathbf{A}\mathbf{R}_{s}\mathbf{A}^{\dagger}\\{\mathbf{E}_{\mathrm{tx},K}\otimes(\mathbf{L}\mathbf{\tilde{E}}_{\mathrm{rx},K-1})\\}\,.$

Now, we sequentially apply PKP 1 and PKP 4 to obtain

$\displaystyle\mathbf{A}^{\dagger}\\{\mathbf{E}_{\mathrm{tx},K}\otimes(\mathbf{L}\mathbf{\tilde{E}}_{\mathrm{rx},K-1})\\}$
$\displaystyle=\begin{bmatrix}\mathbf{a}^{\dagger}_{\mathrm{tx,1}}\otimes\mathbf{a}^{\dagger}_{\mathrm{rx,1}}\\\
\vdots\\\
\mathbf{a}^{\dagger}_{\mathrm{tx,K}}\otimes\mathbf{a}^{\dagger}_{\mathrm{rx,K}}\end{bmatrix}(\mathbf{E}_{\mathrm{tx},K}\otimes\mathbf{L}\mathbf{\tilde{E}}_{\mathrm{rx},K-1})$
$\displaystyle=\begin{bmatrix}\mathbf{a}^{\dagger}_{\mathrm{tx,1}}\mathbf{E}_{\mathrm{tx},K}\otimes\mathbf{a}^{\dagger}_{\mathrm{rx,1}}\mathbf{L}\mathbf{\tilde{E}}_{\mathrm{rx},K-1}\\\
\vdots\\\
\mathbf{a}^{\dagger}_{\mathrm{tx,K}}\mathbf{E}_{\mathrm{tx},K}\otimes\mathbf{a}^{\dagger}_{\mathrm{rx,K}}\mathbf{L}\mathbf{\tilde{E}}_{\mathrm{rx},K-1}\end{bmatrix}$
$\displaystyle={\mathbf{0}_{K\times 1}}.$

Note that the last equality is due to (8) and yields the desired result in
(14). Similarly, we can also prove that

$\mathbf{R}\\{\mathbf{u}_{2}-\exp(j\epsilon_{\mathrm{rx},K})\mathbf{u}_{3}\\}={\mathbf{0}_{N_{\mathrm{tx}}N_{\mathrm{rx}}\times
1}}.$

In a noiseless scenario, by the definition of non-signal subspace, we can
write,

$\mathbf{R}\mathbf{{E}}_{K}={\mathbf{0}_{N_{\mathrm{tx}}N_{\mathrm{rx}}\times
1}}.$ (15)

Inserting $\mathbf{E}_{K}$ from (9) into (15), we find

$\mathbf{R}[\mathbf{u}_{0}-\exp(j\epsilon_{\mathrm{rx},K})\mathbf{u}_{1}-\exp(j\epsilon_{\mathrm{tx},K})\mathbf{u}_{2}\\\
+\exp(j\epsilon_{\mathrm{tx},K})\exp(j\epsilon_{\mathrm{rx},K})\mathbf{u}_{3}]={\mathbf{0}_{N_{\mathrm{tx}}N_{\mathrm{rx}}\times
1}}.$ (16)

Since $\mathbf{u}_{1}-\exp(j\epsilon_{\mathrm{tx},K})\mathbf{u}_{3}$ and
$\mathbf{u}_{2}-\exp(j\epsilon_{\mathrm{rx},K})\mathbf{u}_{3}$ belong the non-
signal subspace, the above equation yields

$\mathbf{R}{\left[{\mathbf{u}_{0}-\exp(j\epsilon_{\mathrm{tx},K})\exp(j\epsilon_{\mathrm{rx},K})\mathbf{u}_{3}}\right]}={\mathbf{0}_{N_{\mathrm{tx}}N_{\mathrm{rx}}\times
1}}\,,$ (17)

and therefore, vector
$\mathbf{u}_{0}-\exp(j\epsilon_{\mathrm{tx},K})\exp(j\epsilon_{\mathrm{rx},K})\mathbf{u}_{3}$
also belongs to the non-signal subspace. ∎

The above theorem results in the following corollary which will guide us on
how to iteratively estimate steering elements (or DOA angles) of the
${K}^{\underline{\text{th}}}$ user utilizing those of the other $K-1$ users
along with the auto-correlation matrix $\mathbf{{R}}$.

######  Corollary 1.

When there are $K$ sources in the scene and the number of antenna elements at
the transmitter and receiver satisfies $N_{\mathrm{tx}}=N_{\mathrm{rx}}=K+1$,
given the steering elements of $K-1$ sources (or
$\\{\mathbf{{u}}_{i}\\}_{i=0}^{3}$), the steering elements of the
$K^{\underline{\mathrm{th}}}$ source satisfy

$\exp(j\epsilon_{\mathrm{tx},K})=\frac{\mathbf{u}^{\dagger}_{1}\mathbf{R}\mathbf{u}_{1}}{\mathbf{u}^{\dagger}_{1}\mathbf{R}\mathbf{u}_{3}}\,,$
(18)
$\exp(j\epsilon_{\mathrm{rx},K})=\frac{\mathbf{u}^{\dagger}_{2}\mathbf{R}\mathbf{u}_{2}}{\mathbf{u}^{\dagger}_{2}\mathbf{R}\mathbf{u}_{3}}\,,$
(19)

and

$\exp(j\epsilon_{\mathrm{tx},K})\exp(j\epsilon_{\mathrm{rx},K})=\frac{\mathbf{u}^{\dagger}_{0}\mathbf{R}\mathbf{u}_{0}}{\mathbf{u}^{\dagger}_{0}\mathbf{R}\mathbf{u}_{3}}\,.$
(20)

Notice that in the above corollary (18) and (19) provide sufficient
information regarding the steering elements of the
${K}^{\underline{\text{th}}}$ user. In the other words, (20) could viewed as
redundant while (18) and (19) are available.

##  IV Iterative 2D DOA Estimation Algorithm (iDEA)

In this section, using Corollary 20, we first develop the proposed iterative
2D DOA estimation algorithm (iDEA) assuming the number of antenna elements
$N_{\mathrm{tx}}=N_{\mathrm{rx}}=K+1$ in a noise-free environment. Next, we
generalize iDEA for the number of antenna elements $N_{\mathrm{tx}}>K+1$ or
$N_{\mathrm{rx}}>K+1$.

###  IV-A When $N_{\mathrm{tx}}=N_{\mathrm{rx}}=K+1$

In one iteration, iDEA will execute the following steps.

  * •

Step 1: Initialize the ${1}^{\underline{\text{st}}}$ $(K-1)$ sources’ steering
elements and derive $\\{\mathbf{u}_{q}\\}^{q=3}_{q=0}$ assuming these $K-1$
sources form a set $\mathcal{Q}$.

  * •

Step 2: Use (18) and (19) from Corollary 20 to estimate the steering elements
of the $K^{\underline{\mathrm{th}}}$ source.

  * •

Step 3:

    * –

a. From set $\mathcal{Q}$, take out first $K-1$ users.

    * –

b. Label users 1 to $K-1$ as $2$ to $K$.

    * –

c. Label the estimated $K^{\underline{\mathrm{th}}}$ source steering elements
from Step 2 as user 1.

    * –

d. Update $\\{\mathbf{u}_{q}\\}^{q=3}_{q=0}$ from first $K-1$ users from set
$\mathcal{Q}$

    * –

e. Go TO Step 2.

When all the $K$ sources’ steering elements in the scene are updated once, and
we count that as one full iteration. We see that the complexity of iDEA will
grow with each iteration. To control its complexity, we introduce a parameter,
the maximum number of iterations $T$. The value of $T$ allows a trade-off
between the complexity of the algorithm and the accuracy of the DOA
estimation. If the iteration reaches $t=T$, iDEA immediately stops and
calculates the DOAs for all the sources in the scene using the estimated
steering elements. For the ${k}^{\underline{\text{th}}}$ source, the elevation
angle can be estimated as

$\hat{\theta}_{k}=\sin^{-1}\left({a_{d}[\\{\angle\exp(j\hat{\epsilon}_{\mathrm{tx},k})\\}^{2}+\\{\angle\exp(j\hat{\epsilon}_{\mathrm{rx},k})\\}^{2}]^{\frac{1}{2}}}\right),$
(21)

where $a_{d}=\lambda/(2\pi d)$ and $k=1,\ldots,K$. Using the above estimate of
the elevation angle, the azimuthal angle of user $k$ is estimated as

$\hat{\phi}_{k}=\cos^{-1}\left({\frac{a_{d}\angle\exp(j\hat{\epsilon}_{\mathrm{tx},k})}{\sin{\hat{\theta}_{k}}}}\right)\bigcap\\\
\sin^{-1}\left({\frac{a_{d}\angle\exp(j\epsilon_{\mathrm{rx},k})}{\sin\hat{\theta}_{k}}}\right)\,.$
(22)

Every steps of iDEA has been described with mathematical notations below:

Algorithm 1 Iterative 2D DOA Estimation Algorithm (iDEA)

1: Initialize $K-1$ sources’ DOAs, $\theta_{k}$, and $\phi_{k}$ for
$k=1,....,K-1$;

2: Calculate $\epsilon_{\mathrm{rx},k}$ and $\epsilon_{\mathrm{rx},k}$ for
$k=1,....,K-1$, and form set $\mathcal{Q}$;

3: Derive $\\{\mathbf{u}_{q}\\}^{q=3}_{q=0}$ from these $K-1$ sources’
steering elements;

4: Set the iteration index, $t=0$;

5: while $t<T$ do

6: $t=t+1$;

7: for $k=1,\ldots,K-1$ do

8: Find $K^{\underline{\mathrm{th}}}$ source steering elements,
$\exp(j\epsilon_{\mathrm{tx},K})$ and $\exp(j\epsilon_{\mathrm{rx},K})$ using
(18) and (19) from Corollary 20;

9: From set $\mathcal{Q}$, take out first $K-1$ users;

10: Label users 1 to $K-1$ as $2$ to $K$;

11: Label the estimated $K^{\underline{\mathrm{th}}}$ source steering elements
from Step 8 as user 1;

12: Update $\\{\mathbf{u}_{q}\\}^{q=3}_{q=0}$ from first $K-1$ users from set
$\mathcal{Q}$;

13: end for

14: Calculate $\exp(j\epsilon_{\mathrm{tx},K})$,
$\exp(j\epsilon_{\mathrm{rx},K})$ using Corollary 20;

15: end while

16: Using the estimated steering elements, estimate all the DOAs,
$\theta_{k}$, and $\phi_{k}$ for $k=1,....,K,$ utilizing (21)- (22);

###  IV-B When $N_{\mathrm{tx}}>K+1$ or $N_{\mathrm{rx}}>K+1$

So far, we have assumed the number of elements at the transmitter and
receiver, $N_{\mathrm{tx}}=N_{\mathrm{rx}}=K+1$. However, it is possible to
have $N_{\mathrm{tx}}>K+1$ or $N_{\mathrm{rx}}>K+1$ in the scene. In that
case, we divide the URA into different uniform sub-arrays, which have been
extensively studied to enhance the DOA estimation precision and resolution
[30, 31, 32]. These uniform sub-arrays should be of size $(K+1)\times(K+1)$,
where the number of elements in both Tx and Rx axes is $K+1$. In this way, the
total number of sub-arrays is
$N_{\mathrm{sa}}=(N_{\mathrm{tx}}-K)\times(N_{\mathrm{rx}}-K)$. Now, we
construct an averaged sample covariance matrix from all these sub-arrays,
which will be utilized to perform a covariance matrix smoothing prior to
applying the iterative algorithm. This averaged covariance matrix can be
defined as

$\hat{\mathbf{R}}_{\mathrm{sa}}=\frac{1}{N_{\mathrm{sa}}}\sum_{i=1}^{N_{\mathrm{sa}}}\hat{\mathbf{R}_{i}}.$
(23)

##  V Proof of Convergence

In this section, our goal is to prove the convergence of iDEA. Recall that
iDEA operates by updating a user’s (say user $K$’s) steering elements,
$\exp(j\epsilon_{\mathrm{tx},K,\tau})$ and
$\exp(j\epsilon_{\mathrm{rx},K,\tau})$, where $\tau$ is the update-index. At
different updates, different users are assigned which is denoted by $K$. In
order iDEA to converge, we want

$\lim_{\tau\rightarrow\infty}\epsilon_{\mathrm{m},K,\tau}=\epsilon_{\mathrm{m},K}\,,$

where notation $\mathrm{m}$ is the generalized notation that could be either
$\mathrm{tx}$ or $\mathrm{rx}$. We begin the proof of convergence by proving
the following lemmas. In order to do so, we would like to introduce following
four terms

$\displaystyle a_{\mathrm{m},k,\tau}$ $\displaystyle=$
$\displaystyle\tilde{\mathbf{{E}}}_{\mathrm{m},K-1,\tau}^{{\dagger}}\mathbf{{R}}_{\mathrm{m},k}\tilde{\mathbf{{E}}}_{\mathrm{m},K-1,\tau}\,,$
(24) $\displaystyle b_{\mathrm{m},k,\tau}$ $\displaystyle=$
$\displaystyle\tilde{\mathbf{{E}}}_{\mathrm{m},K-1,\tau}^{{\dagger}}\mathbf{{L}}^{\dagger}\mathbf{{R}}_{\mathrm{m},k}\tilde{\mathbf{{E}}}_{\mathrm{m},K-1,\tau}\,,$
(25) $\displaystyle c_{\mathrm{m},k,\tau}$ $\displaystyle=$
$\displaystyle\tilde{\mathbf{{E}}}_{\mathrm{m},K-1,\tau}^{{\dagger}}\mathbf{{R}}_{\mathrm{m},k}\mathbf{{L}}\tilde{\mathbf{{E}}}_{\mathrm{m},K-1,\tau}\,,$
(26)

and

$\displaystyle d_{\mathrm{m},k,\tau}$ $\displaystyle=$
$\displaystyle\tilde{\mathbf{{E}}}_{\mathrm{m},K-1,\tau}^{{\dagger}}\mathbf{{L}}^{\dagger}\mathbf{{R}}_{\mathrm{m},k}\mathbf{{L}}\tilde{\mathbf{{E}}}_{\mathrm{m},K-1,\tau}\,,$
(27)

where
$\mathbf{{R}}_{\mathrm{m},k}=\mathbf{{a}}_{\mathrm{m},k}\mathbf{{a}}_{\mathrm{m},k}^{\dagger}$.

######  Lemma 2.

When the sources are uncorrelated, we have

$\displaystyle\mathbf{u}^{\dagger}_{1,\tau}\mathbf{R}\mathbf{u}_{1,\tau}$
$\displaystyle=$ $\displaystyle\sum_{k=1}^{K}\sigma^{2}_{k}\times
a_{\mathrm{tx},k,\tau}\times d_{\mathrm{rx},k,\tau}\,,$ (28)
$\displaystyle\mathbf{u}^{\dagger}_{1,\tau}\mathbf{R}\mathbf{u}_{3,\tau}$
$\displaystyle=$ $\displaystyle\sum_{k=1}^{K}\sigma^{2}_{k}\times
c_{\mathrm{tx},k,\tau}\times a_{\mathrm{rx},k,\tau}\,,$ (29)
$\displaystyle\mathbf{u}^{\dagger}_{2,\tau}\mathbf{R}\mathbf{u}_{3,\tau}$
$\displaystyle=$ $\displaystyle\sum_{k=1}^{K}\sigma^{2}_{k}\times
b_{\mathrm{rx},k,\tau}\times d_{\mathrm{tx},k,\tau}\,,$ (30)
$\displaystyle\mathbf{u}^{\dagger}_{2,\tau}\mathbf{R}\mathbf{u}_{2,\tau}$
$\displaystyle=$ $\displaystyle\sum_{k=1}^{K}\sigma^{2}_{k}\times
d_{\mathrm{tx},k,\tau}\times a_{\mathrm{ax},k,\tau}\,.$ (31)

###### Proof:

Using (11), we obtain

$\displaystyle\mathbf{u}^{\dagger}_{1,\tau}\mathbf{R}\mathbf{u}_{1,\tau}$
$\displaystyle=\\{\mathbf{\tilde{E}}_{\mathrm{tx},K-1,\tau}\otimes(\mathbf{L}\mathbf{\tilde{E}}_{\mathrm{rx},K-1,\tau})\\}^{\dagger}\mathbf{R}$
$\displaystyle\times\\{\mathbf{\tilde{E}}_{\mathrm{tx},K-1,\tau}\otimes(\mathbf{L}\mathbf{\tilde{E}}_{\mathrm{rx},K-1,\tau})\\}\,.$
(32)

Since we are considering noiseless scenario with uncorrelated sources,
$\mathbf{{R}}=\sum_{k=1}^{K}\sigma^{2}_{k}\mathbf{{A}}_{k}\mathbf{{A}}_{k}^{\dagger}$.
Thus the above equation becomes

$\displaystyle\mathbf{u}^{\dagger}_{1,\tau}\mathbf{R}\mathbf{u}_{1,\tau}$
$\displaystyle=\sum_{k=1}^{K}\sigma_{k}^{2}\\{\mathbf{\tilde{E}}_{\mathrm{tx},K-1,\tau}\otimes(\mathbf{L}\mathbf{\tilde{E}}_{\mathrm{rx},K-1,\tau})\\}^{\dagger}\mathbf{A}_{k}$
$\displaystyle\times\mathbf{A}_{k}^{\dagger}\\{\mathbf{\tilde{E}}_{\mathrm{tx},K-1,\tau}\otimes(\mathbf{L}\mathbf{\tilde{E}}_{\mathrm{rx},K-1},\tau)\\}\,.$
(33)

Recall that
$\mathbf{A}_{k}=\mathbf{a}_{\mathrm{tx},k}\otimes\mathbf{a}_{\mathrm{rx},k}$.
Thus, PKP-1 and PKP-4 jointly allow us to write

$\\{\mathbf{\tilde{E}}_{\mathrm{tx},K-1,\tau}\otimes(\mathbf{L}\mathbf{\tilde{E}}_{\mathrm{rx},K-1,\tau})\\}^{\dagger}\mathbf{A}_{k}$

equals to

$\mathbf{\tilde{E}}_{\mathrm{tx},K-1,\tau}^{\dagger}\mathbf{a}_{\mathrm{tx},k}\otimes\mathbf{\tilde{E}}_{\mathrm{rx},K-1,\tau}^{\dagger}\mathbf{L}^{\dagger}\mathbf{a}_{\mathrm{rx},k}$

and

$\mathbf{A}_{k}^{\dagger}\\{\mathbf{\tilde{E}}_{\mathrm{tx},K-1,\tau}\otimes(\mathbf{L}\mathbf{\tilde{E}}_{\mathrm{rx},K-1,\tau})\\}$

equals to

$\mathbf{a}_{\mathrm{tx},k}^{\dagger}\mathbf{\tilde{E}}_{\mathrm{tx},K-1,\tau}\otimes\mathbf{a}_{\mathrm{rx},k}^{\dagger}\mathbf{L}\mathbf{\tilde{E}}_{\mathrm{rx},K-1,\tau}\,.$

Therefore, after applying PKP-4 we rewrite

$\displaystyle\mathbf{u}^{\dagger}_{1,\tau}\mathbf{R}\mathbf{u}_{1,\tau}$ (34)
$\displaystyle=\sum_{k=1}^{K}\sigma_{k}^{2}\mathbf{\tilde{E}}_{\mathrm{tx},K-1,\tau}^{\dagger}\mathbf{R}_{\mathrm{tx},k}\mathbf{\tilde{E}}_{\mathrm{tx},K-1,\tau}$
$\displaystyle\qquad\qquad\otimes\mathbf{\tilde{E}}_{\mathrm{rx},K-1,\tau}^{\dagger}\mathbf{L}^{\dagger}\mathbf{R}_{\mathrm{rx},k}\mathbf{L}\mathbf{\tilde{E}}_{\mathrm{rx},K-1,\tau}$
$\displaystyle=\sum_{k=1}^{K}\sigma_{k}^{2}\times
a_{\mathrm{tx},k,\tau}\otimes d_{\mathrm{rx},k,\tau}$
$\displaystyle=\sum_{k=1}^{K}\sigma_{k}^{2}\times a_{\mathrm{tx},k,\tau}\times
d_{\mathrm{rx},k,\tau}\,.$ (35)

The last equality is due to the fact that $\otimes$ operation of two scalars
is equivalent to their multiplication. The equality before the last one is due
to (24) and (27). Similarly, we can derive (29)-(31). ∎

Following the proof of Theorem 1 of [16], we get

$\tilde{\mathbf{{E}}}_{\mathrm{m},K-1,\tau}^{{\dagger}}\mathbf{{a}}_{\mathrm{m},k}=\Pi_{\iota=1}^{K-1}\left({1-\exp\left\\{j\left({\epsilon_{\mathrm{m},\iota,\tau}-\epsilon_{\mathrm{m},k}}\right)\right\\}}\right)\,,$

and

$\displaystyle\tilde{\mathbf{{E}}}_{\mathrm{m},K-1,\tau}^{{\dagger}}\mathbf{{L}}^{\dagger}\mathbf{{a}}_{\mathrm{m},k}$
$\displaystyle=\exp\left({j\epsilon_{\mathrm{m},k}}\right)\times$
$\displaystyle\Pi_{\iota=1}^{K-1}\left({1-\exp\left\\{j\left({\epsilon_{\mathrm{m},\iota,\tau}-\epsilon_{\mathrm{m},k}}\right)\right\\}}\right)\,,$

which we apply to (24) - (27) for obtaining the following lemma.

######  Lemma 3.

When the sources are uncorrelated, (24), (25), and (26) can be expressed in
terms of the following trigonometric functions:

$\displaystyle a_{\mathrm{m},k,\tau}$ $\displaystyle=$ $\displaystyle
4^{K-1}\Pi_{\iota=1}^{{K-1}}\sin^{2}\left(\frac{\epsilon_{\mathrm{m},k}-\hat{\epsilon}_{\mathrm{m},\iota,\tau}}{2}\right)\,,$
(36) $\displaystyle b_{\mathrm{m},k,\tau}$ $\displaystyle=$ $\displaystyle
4^{K-1}e^{j{\epsilon}_{\mathrm{m},k}}\Pi_{\iota=1}^{K-1}\sin^{2}\left(\frac{\epsilon_{\mathrm{m},k}-\hat{\epsilon}_{\mathrm{m},\iota,\tau}}{2}\right)$
(37) $\displaystyle=$ $\displaystyle
e^{j{\epsilon}_{\mathrm{m},k}}a_{\mathrm{m},k,\tau}\,,$ $\displaystyle
c_{\mathrm{m},k,\tau}$ $\displaystyle=$ $\displaystyle
4^{K-1}e^{-j{\epsilon}_{\mathrm{m},k}}\Pi_{\iota=1}^{K-1}\sin^{2}\left(\frac{\epsilon_{\mathrm{m},k}-\hat{\epsilon}_{\mathrm{m},\iota,\tau}}{2}\right)$
(38) $\displaystyle=$ $\displaystyle
e^{-j{\epsilon}_{\mathrm{m},k}}a_{\mathrm{m},k,\tau}\,,$

and $d_{\mathrm{m},k,\tau}=a_{\mathrm{m},k,\tau}$.

Now we are almost ready to prove the convergence of iDEA. Before doing so, we
need the following two cost functions for updating
$\exp(j\epsilon_{\mathrm{tx},K,\tau})$

$\displaystyle Q_{\mathrm{tx},K,\tau}$ $\displaystyle=$
$\displaystyle\\{\mathbf{u}_{1,\tau}-\exp(j\epsilon_{\mathrm{tx},K,\tau})\mathbf{u}_{3,\tau}\\}^{\dagger}\mathbf{R}$
(39) $\displaystyle\times$
$\displaystyle\\{\mathbf{u}_{1,\tau}-\exp(j\epsilon_{\mathrm{tx},K,\tau})\mathbf{u}_{3,\tau}\\}\,,$

and

$\displaystyle{Q}_{\mathrm{tx},K,\tau+1}$ $\displaystyle=$
$\displaystyle\\{\mathbf{u}_{1,\tau}-\exp(j\epsilon_{\mathrm{tx},K,\tau+1})\mathbf{u}_{3,\tau}\\}^{\dagger}\mathbf{R}$
(40) $\displaystyle\times$
$\displaystyle\\{\mathbf{u}_{1,\tau}-\exp(j\epsilon_{\mathrm{tx},K,\tau+1})\mathbf{u}_{3,\tau}\\}\,.$

Similarly, for updating $\exp(j\epsilon_{\mathrm{rx},K,\tau})$ we define

$\displaystyle Q_{\mathrm{rx},K,\tau}$ $\displaystyle=$
$\displaystyle\\{\mathbf{u}_{2,\tau}-\exp(j\epsilon_{\mathrm{rx},K,\tau})\mathbf{u}_{3,\tau}\\}^{\dagger}\mathbf{R}$
(41) $\displaystyle\times$
$\displaystyle\\{\mathbf{u}_{2,\tau}-\exp(j\epsilon_{\mathrm{rx},K,\tau})\mathbf{u}_{3,\tau}\\}\,,$

and

$\displaystyle{Q}_{\mathrm{rx},K,\tau+1}$ $\displaystyle=$
$\displaystyle\\{\mathbf{u}_{2,\tau}-\exp(j\epsilon_{\mathrm{rx},K,\tau+1})\mathbf{u}_{3,\tau}\\}^{\dagger}\mathbf{R}$
(42) $\displaystyle\times$
$\displaystyle\\{\mathbf{u}_{2,\tau}-\exp(j\epsilon_{\mathrm{rx},K,\tau+1})\mathbf{u}_{3,\tau}\\}\,,$

######  Theorem 3.

Since

$Q_{\mathrm{m},K,\tau+1}<{Q}_{\mathrm{m},K,\tau}\longrightarrow\lim_{\tau\rightarrow\infty}Q_{\mathrm{m},K,\tau}=0\,,$

for both m = tx and m = rx. Therefore,

$\lim_{\tau\rightarrow\infty}\epsilon_{\mathrm{m},K,\tau}=\epsilon_{\mathrm{m},K}\,,$

which ensures the convergence of iDEA.

###### Proof:

First, we demonstrate the desired result for m = tx by proving that

${Q}_{\mathrm{tx},K,\tau}-Q_{\mathrm{tx},K,\tau+1}>0\,.$

After some simplifications, (39) and (40) jointly yield

${Q}_{\mathrm{tx},K,\tau}-Q_{\mathrm{tx},K,\tau+1}\\\
=-\exp(-\epsilon_{\mathrm{tx},K,\tau+1})\mathbf{u^{\dagger}}_{3,\tau}\mathbf{R}\mathbf{u}_{1,\tau}-\exp(\epsilon_{\mathrm{tx},K,\tau+1})\mathbf{u}_{1,\tau}^{\dagger}\mathbf{R}\mathbf{u}_{3,\tau}\\\
+\exp(-\epsilon_{\mathrm{tx},K,\tau})\mathbf{u}_{3,\tau}^{\dagger}\mathbf{R}\mathbf{u}_{1,\tau}+\exp(\epsilon_{\mathrm{x},K,\tau})\mathbf{u}_{1,\tau}^{\dagger}\mathbf{R}\mathbf{u}_{3,\tau}.$
(43)

According to Corollary 20,

$\exp(\epsilon_{\mathrm{tx},K,\tau+1})=\frac{\mathbf{u}^{\dagger}_{1,\tau}\mathbf{R}\mathbf{u}_{1,\tau}}{\mathbf{u}^{\dagger}_{1,\tau}\mathbf{R}\mathbf{u}_{3,\tau}},$

which yields the first line of the right-hand side of (43) as

$\displaystyle-\exp(-\epsilon_{\mathrm{tx},K,\tau+1})\mathbf{u^{\dagger}}_{3,\tau}\mathbf{R}\mathbf{u}_{1,\tau+1}-\exp(\epsilon_{\mathrm{tx},K,\tau})\mathbf{u}_{1,\tau}^{\dagger}\mathbf{R}\mathbf{u}_{3,\tau}$
$\displaystyle=-2\mathbf{u}^{\dagger}_{1,\tau}\mathbf{R}\mathbf{u}_{1,\tau}=2\sum_{k=1}^{K}\sigma^{2}_{k}a_{\mathrm{tx},k,\tau}a_{\mathrm{rx},k,\tau}$
(44)

Last equality is due to Lemmas 2 and 3. These two lemmas also yield the second
line of the right-hand side of (43) as

$\displaystyle-\exp(-\epsilon_{\mathrm{tx},K,\tau})\mathbf{u}_{3,\tau}^{\dagger}\mathbf{R}\mathbf{u}_{1,\tau}-\exp(\epsilon_{\mathrm{tx},K,\tau})\mathbf{u}_{1,\tau}^{\dagger}\mathbf{R}\mathbf{u}_{3,\tau}$
(45)
$\displaystyle=-\sum_{k=1}^{K}\sigma^{2}_{k}a_{\mathrm{tx},k,\tau}a_{\mathrm{rx},k,\tau}(\exp(\epsilon_{\mathrm{tx},k}-\epsilon_{\mathrm{tx},K,\tau})$
$\displaystyle\qquad\qquad+\exp(-\epsilon_{\mathrm{tx},k}+\epsilon_{\mathrm{tx},K,\tau}))\,,$
(46)

which can be further simplified to write

$\displaystyle-\exp(-\epsilon_{\mathrm{tx},K,\tau})\mathbf{u}_{3,\tau}^{\dagger}\mathbf{R}\mathbf{u}_{1,\tau}-\exp(\epsilon_{\mathrm{tx},K,\tau})\mathbf{u}_{1,\tau}^{\dagger}\mathbf{R}\mathbf{u}_{3,\tau}$
$\displaystyle=-2\sum_{k=1}^{K}\sigma^{2}_{k}a_{\mathrm{tx},k,\tau}a_{\mathrm{rx},k,\tau}\cos\left({\epsilon_{\mathrm{tx},k}-\epsilon_{\mathrm{tx},K,\tau}}\right)$
(47)

Using (44) and (47) with (43), we get

${Q}_{\mathrm{tx},K,\tau}-Q_{\mathrm{tx},K,\tau+1}\\\
=2\sum_{k=1}^{K}\sigma^{2}_{k}a_{\mathrm{tx},k,\tau}a_{\mathrm{rx},k,\tau}\left({1-\cos\left({\epsilon_{\mathrm{tx},k}-\epsilon_{\mathrm{tx},K,\tau}}\right)}\right)\\\
=4^{2K-1}\sum_{k=1}^{K}\sigma^{2}_{k}\sin^{2}\left(\frac{\epsilon_{\mathrm{tx},k}-\epsilon_{\mathrm{tx},K,\tau}}{2}\right)\times\\\
\Pi_{m=1}^{{K-1}}\sin^{2}\left(\frac{\epsilon_{\mathrm{tx},k}-\epsilon_{\mathrm{tx},m,\tau}}{2}\right)\times\\\
\Pi_{n=1}^{{K-1}}\sin^{2}\left(\frac{\epsilon_{\mathrm{rx},k}-\epsilon_{\mathrm{rx},n,\tau}}{2}\right)$
(48)

Every term from $k=1$ to $k=K-1$ in then summation of (48) is non-negative,
thus

${Q}_{\mathrm{tx},K,\tau}-Q_{\mathrm{tx},K,\tau+1}\\\
\geq\sigma^{2}_{K}\sin^{2}\left(\frac{\epsilon_{\mathrm{tx},K}-\epsilon_{\mathrm{tx},K,\tau}}{2}\right)\times\\\
\Pi_{m=1}^{{K-1}}\sin^{2}\left(\frac{\epsilon_{\mathrm{tx},K}-\epsilon_{\mathrm{tx},m,\tau}}{2}\right)\times\\\
\Pi_{n=1}^{{K-1}}\sin^{2}\left(\frac{\epsilon_{\mathrm{rx},K}-\epsilon_{\mathrm{rx},n,\tau}}{2}\right)$
(49)

In a wireless environment with random DOA, the lower bound in (49) is positive
definite (i.e. strictly $>0$) unless
$\epsilon_{\mathrm{tx},K}=\epsilon_{\mathrm{tx},K,\tau}$. Since the cost
function is bounded from below by zero and has been found to be strictly
monotonically decreasing, it is deduced that the algorithm will converge [33],
i.e.

$\lim_{\tau\rightarrow\infty}\epsilon_{\mathrm{tx},K,\tau}=\epsilon_{\mathrm{tx},K}\,.$

Similarly, using the following positive definite upper-bound

${Q}_{\mathrm{rx},K,\tau}-Q_{\mathrm{rx},K,\tau+1}\\\
\geq\sigma^{2}_{K}\sin^{2}\left(\frac{\epsilon_{\mathrm{rx},K}-\epsilon_{\mathrm{rx},K,\tau}}{2}\right)\times\\\
\Pi_{m=1}^{{K-1}}\sin^{2}\left(\frac{\epsilon_{\mathrm{tx},K}-\epsilon_{\mathrm{tx},m,\tau}}{2}\right)\times\\\
\Pi_{n=1}^{{K-1}}\sin^{2}\left(\frac{\epsilon_{\mathrm{rx},K}-\epsilon_{\mathrm{rx},n,\tau}}{2}\right)>0\,,$
(50)

we claim

$\lim_{\tau\rightarrow\infty}\epsilon_{\mathrm{rx},K,\tau}=\epsilon_{\mathrm{rx},K}\,.$

∎

##  VI Numerical Results

In this section, our objectives are, 1) to demonstrate the convergence of the
proposed algorithm, iDEA, 2) to provide numerical simulation results to
exhibit the performance of iDEA, and 3) to perform a complexity analysis
between iDEA, and the 2D MUSIC algorithm in a MIMO system. In all the
numerical examples, iDEA is implemented without any prior knowledge of the
noise power, $\sigma^{2}$. However, its performance can be further improved by
estimating $\sigma^{2}$ [34, 35, 36], and subtracting it from the auto-
correlation matrix. Of course, this estimation will increase algorithm
complexity. In addition, the implementation of the baseline algorithm, 2D
MUSIC in a MIMO system also does not hinge on the knowledge of $\sigma^{2}$;
thus, the above implementation of iDEA will make the comparison fair.

For the simulations, we consider $K=2$ non-coherent far-field targets in a
MIMO system with $N_{\mathrm{tx}}=3$ transmitters, and $N_{\mathrm{rx}}=3$
receivers, which are uniformly placed (spacing, $d=\lambda/2$) along the
$\mathrm{x}-$ , and the $\mathrm{y}-$ axis, respectively. Therefore, the angle
between Tx and Rx axes, $\phi_{\mathrm{trx}}=90^{\circ}.$ Targets’ DOA angels
$(\theta_{1},\phi_{1})=(30^{\circ},25^{\circ})$, and
$(\theta_{2},\phi_{2})=(70^{\circ},80^{\circ})$. Lastly, in order to run iDEA,
we consider the initial values of the DOAs,
$(\theta_{1,\mathrm{in}},\phi_{1,\mathrm{in}})=(10^{\circ},10^{\circ})$. For
the 2D MUSIC algorithm in the MIMO system, we have used the search with
$0.033^{\circ}$ angle precision for the DOAs. The above parameters are kept
the same throughout the entire numerical study unless otherwise specified.

###  VI-A Convergence as a Function of the Iterations

iDEA is an iterative algorithm, thus the demonstration of its convergence is
important. In Figure 2, and 3, we plot the cost function $Q_{\mathrm{tx},K}$
in (39), and $Q_{\mathrm{rx},K}$ in (41) as a function of the number of
iterations ($T$) for two different number of samples, $M=10$, and $M=500$ at
SNR = 30 dB. As expected the convergence of iDEA depends on the number of
samples, and a higher number of samples yields faster convergence. For
example, when $M=10$, iDEA takes 4 iterations to get closer to the convergence
value, and for $M=500$, this number is 2. Note that, after reaching the
convergence, any increase in the number of iterations, $T$ may not contribute
to the performance enhancement of the algorithm.

Figure 2: Convergence of iDEA for $M=10$. Figure 3: Convergence of iDEA for
$M=500$.

###  VI-B Complexity Analysis

Before we go any further evaluating the performance of our algorithm, we
compare the computational complexity of iDEA against 2D MUSIC for estimating
DOAs in MIMO environments. We consider the complexity incurred after
performing the matched filtering, and it is measured by the number of complex
multiplications associated with the major operations. Note that one complex
multiplication is equivalent to four real multiplications.

For iDEA, the cost of estimating the auto-correlation matrix is given by
$(K+1)^{2}M+M^{2}+(K+1)^{2}$ when $N_{\mathrm{tx}}=N_{\mathrm{rx}}=K+1$ [16].
if $N>K+1$, then $N_{\mathrm{sa}}$ sub arrays can be formed using (23). Hence,
the total calculation associated with the auto-correlation matrix is
$\\{(K+1)^{2}M+M^{2}+(K+1)^{2}\\}N_{\mathrm{sa}}$. Recall that iDEA operates
by estimating the ${K}^{\underline{\text{th}}}$ source steering elements using
(18) and (19) from Corollary 20. The cost of estimating these elements is
$2KT\\{(K+1)\times(K+2)+(K+1)+1\\}$, where $K$ is the number of sources and
$T$ is the number of iterations. Finally, after estimating all the steering
elements, the cost associated with the estimation of the DOAs using (21), and
(22) is $(7+4+4+1)/4=4$ complex multiplications [37, 38]. Therefore, the total
cost associated with iDEA is
$\\{(K+1)^{2}M+M^{2}+(K+1)^{2}\\}N_{\mathrm{sa}}+2KT\\{(K+1)\times(K+2)+(K+1)+1\\}+4$.
Notice that the complexity of the iDEA algorithm is primarily dictated by the
number of samples ($M$), the number of iterations ($T$) and the number of
sources ($K$).

Now, we focus on calculating the complexity of 2D MUSIC which performs
eigenvalue decomposition (EVD) of the estimated auto-correlation matrix. EVD
often is obtained from a singular value decomposition (SVD). Here, the cost
associated with the estimation of the auto-correlation matrix is given by
$N^{2}M+M^{2}+N^{2}$, and the complexity associated with an SVD is $12N^{3}$,
where $N=N_{\mathrm{tx}}N_{\mathrm{rx}}$[16]. Then MUSIC estimates the DOAs by
2D angle search using the null space of the auto-correlation matrix. This cost
is calculated as $K\times N_{\theta}N_{\phi}\\{N(N-K)+1\\}$, where
$N_{\theta}$ and $N_{\phi}$ represent the searching point number on the
azimuthal and elevation plane. Therefore, the total cost associated with the
2D MUSIC algorithm is $N^{2}M+M^{2}+N^{2}+12N^{3}+K\times
N_{\theta}N_{\phi}\\{N(N-K)+1\\}$. Notice that the virtual array size ($N$),
the number of samples ($M$), the number of sources ($K$), and the product of
the searching points ($N_{\theta}$ and $N_{\phi}$) dominate the complexity of
the 2D MUSIC algorithm. The costs associated with implementing iDEA and 2D
MUSIC are tabulated in Table I for comparison.

TABLE I: Complexity analysis Algorithm | Complexity (Number of complex multiplications)
---|---
iDEA | $\\{(K+1)^{2}(M+1)+M^{2}\\}N_{\mathrm{sa}}$
$+2KT\\{(K+1)(K+3)+1\\}+4$
2D MUSIC | $N^{2}(M+1)+M^{2}+12N^{3}$
$+KN_{\theta}N_{\phi}\\{N(N-K)+1\\}$
Figure 4: Complexity of 2D MUSIC w.r.t. iDEA algorithm, where
$N_{\mathrm{tx}}=N_{\mathrm{rx}}=K+1$, $M=50$.

Now, we use Table I to demonstrate the complexities of the DOA estimating
algorithms by using a numerical example; see Figure 4, where the complexity of
2D MUSIC (w.r.t the iDEA algorithm) is plotted as a function of the total
number of sources $K$, considering the number of Tx and Rx antennas,
$N_{\mathrm{tx}}=N_{\mathrm{rx}}=K+1$, number of samples, $M=50$, and the
number of iterations in iDEA for three different cases, $T=K$, $T=2K$, and
$T=4K$. Here, the search in 2D MUSIC is conducted with $1^{\circ}$ precision
of the DOA angles in both the elevation and the azimuthal planes. The results
in Figure 4 suggest that depending on the system parameters our proposed iDEA
is capable of offering significant complexity gain over 2D MUSIC. For
instance, the offered gain is over $36$ dB for $K=3$, and $T=4K$.

###  VI-C Performance as a function of SNR/the number of samples

We now shift our focus to the performance iDEA. First, we demonstrate the root
mean square error (RMSE) of iDEA and 2D MUSIC as a function of SNR for the
number of iteration $T=6$ and the number of samples $M=50$; see Figure 5 and 6
for the elevation and azimuthal planes, respectively. Here, it can be noticed
that, at low SNRs, the 2D MUSIC algorithm performs slightly better than iDEA.
However, as SNR increases, the performance difference between the algorithms
starts to diminish. As expected, similar observations are made after comparing
the performance of the above two algorithms by varying the number of samples
for $M$ for given SNR = 10 dB; see Figure 9, and 10. Note that although iDEA
does not have any prior knowledge of noise power in the environment, it
performs well even in the presence of noise.

Figure 5: RMSE vs SNR for the elevation angle estimation. Figure 6: RMSE vs
SNR for the azimuthal angle estimation. Figure 7: RMSE vs SNR for the
elevation angle in two different virtual array scenarios. Figure 8: RMSE vs
SNR for the azimuthal angle in two different virtual array scenarios.

###  VI-D Performance as a function of the number of elements

Finally, we investigate the performance of iDEA by varying the number of array
elements. Here the same two source scenario ($K=2$) are considered for two
different numbers of elements, $\\{N_{\mathrm{tx}}=3,N_{\mathrm{tx}}=3\\}$ ,
and $\\{N_{\mathrm{tx}}=5,N_{\mathrm{tx}}=4\\}$. Recall that when
$N_{\mathrm{tx}}>K+1$ or $N_{\mathrm{rx}}>K+1$, we form $N_{\mathrm{sa}}$ sub-
arrays and average all $N_{\mathrm{sa}}$ auto-correlation matrices according
to (23) for obtaining the auto-correlation matrix for iDEA to operate.
Therefore, for the first case, we obtain, the number of sub-arrays
$N_{\mathrm{sa}}=1$, and for the second case $N_{\mathrm{sa}}=6$. Figure 7,
and 8 depict RMSE for those 2 cases on both the elevation and the azimuthal
planes while varying the SNR from 5 to 30 dB. Here, we have used the number of
iterations, $T=6$, and the total number of samples, $M=50$. It can be seen
that the use of sub-arrays can significantly reduce the RMSE of the proposed
algorithm iDEA. For example, at SNR = 5 dB, these reductions in RMSE are
$55.33\%,44.70\%,29.65\%,37.09\%$ while estimating
$\theta_{1},\,\theta_{2},\,\phi_{1},\,\phi_{2}$, respectively. Note that the
maximum number of sources iDEA can handle is
$\min{\\{N_{\mathrm{tx}},N_{\mathrm{rx}}\\}}-1$. Thus, when
$K<\min{\\{N_{\mathrm{tx}},N_{\mathrm{rx}}\\}}-1$, the use of sub-arrays can
be used to improve the accuracy of the DOA estimation. However, this
performance gain will be achieved at the cost of hardware and software
complexity.

Figure 9: RMSE vs $N$ for the elevation angle estimation. Figure 10: RMSE vs
$N$ for the azimuthal angle estimation.

##  VII Conclusion

In this paper, a low-complexity iterative 2D DOA estimation algorithm, namely
iDEA, has been developed for a MIMO system. This algorithm exploits the
properties of theoretically derived eigenvectors associated with the non-
signal subspace in a MIMO system when the number of signal sources is less
than the number of elements in both transmitters and receivers. The
convergence of the algorithm has been established both mathematically and
numerically. Numerical results showed that the performance of the proposed
algorithm was comparable to that of the 2D MUSIC algorithm in a MIMO system
for a moderate to high SNR while exhibiting a significantly low complexity.
For instance, this complexity gain of the algorithm was over 36 dB for a
number of sources, $K=3$, and a number of iterations, $T=4K$. In all the
numerical examples, iDEA was implemented without the knowledge of the noise
power. Its performance can be improved if the prior knowledge of the noise
power is used or the sub-arrays are formed by using additional antennas at the
transmitter or receiver or both. Using only two extra elements in the
transmitter, and one element in the receiver, in a two-user system the RMSE of
iDEA was reduced by $55.33\%,44.70\%,29.65\%,37.09\%$ while estimating
$\theta_{1},\,\theta_{2},\,\phi_{1},\,\phi_{2}$, respectively at SNR = 5 dB.

## References

  * [1] Z. Chen, G. Gokeda, and Y. Yu, _Introduction to Direction-of-arrival Estimation_. Artech House, 2010.
  * [2] N. Ruan, H. Wang, F. Wen, and J. Shi, “Doa estimation in b5g/6g: Trends and challenges,” _Sensors_ , vol. 22, no. 14, p. 5125, 2022.
  * [3] M. I. Hasan and M. Saquib, “Robust 2-d doa estimation in a polarization sensitive single user environment,” _IEEE Access_ , vol. 10, pp. 115 914–115 929, 2022.
  * [4] ——, “Low complexity single source 2-d doa estimation based on reduced dimension svr,” in _2022 IEEE 22nd Annual Wireless and Microwave Technology Conference (WAMICON)_. IEEE, 2022, pp. 1–4.
  * [5] M. I. Hasan, M. Dautta, and M. A. Haque, “A study of ferromagnetic object detection in ocean environment using empirical mode decomposition,” in _TENCON 2017-2017 IEEE Region 10 Conference_. IEEE, 2017, pp. 2132–2135.
  * [6] N. Zimerman, J. D. Rosenblatt, and T. Routtenberg, “Colored noise in doa estimation from seismic data: an empirical study,” in _2020 54th Asilomar Conference on Signals, Systems, and Computers_. IEEE, 2020, pp. 1240–1244.
  * [7] A. Di Serio, P. Hügler, F. Roos, and C. Waldschmidt, “2-d mimo radar: A method for array performance assessment and design of a planar antenna array,” _IEEE Transactions on Antennas and Propagation_ , vol. 68, no. 6, pp. 4604–4616, 2020.
  * [8] H. Lee and J. Chun, “Virtual array response vector for angle estimation of mimo radar with a wide-band interleaved ofdm signal,” _IEEE Communications Letters_ , vol. 25, no. 5, pp. 1539–1543, 2021.
  * [9] H. Chen, Z. Pan, L. Tian, J. Shi, G. Yang, and M. Suzuki, “A novel awsf algorithm for doa estimation in virtual mimo systems,” _IEEE journal on selected areas in communications_ , vol. 31, no. 10, pp. 1994–2003, 2013.
  * [10] D. Bliss and K. Forsythe, “Multiple-input multiple-output (mimo) radar and imaging: degrees of freedom and resolution,” in _The Thrity-Seventh Asilomar Conference on Signals, Systems & Computers, 2003_, vol. 1. IEEE, 2003, pp. 54–59.
  * [11] R. Schmidt, “Multiple emitter location and signal parameter estimation,” _IEEE transactions on antennas and propagation_ , vol. 34, no. 3, pp. 276–280, 1986.
  * [12] A. G. Jaffer, “Maximum likelihood direction finding of stochastic sources: A separable solution,” in _ICASSP-88., International Conference on Acoustics, Speech, and Signal Processing_. IEEE Computer Society, 1988, pp. 2893–2894.
  * [13] J. Capon, “High-resolution frequency-wavenumber spectrum analysis,” _Proceedings of the IEEE_ , vol. 57, no. 8, pp. 1408–1418, 1969.
  * [14] A. Paulraj, R. Roy, and T. Kailath, “A subspace rotation approach to signal parameter estimation,” _Proceedings of the IEEE_ , vol. 74, no. 7, pp. 1044–1046, 1986.
  * [15] R. Kumaresan and D. W. Tufts, “Estimating the angles of arrival of multiple plane waves,” _IEEE Transactions on Aerospace and electronic Systems_ , no. 1, pp. 134–139, 1983.
  * [16] E. A. Santiago and M. Saquib, “Noise subspace-based iterative technique for direction finding,” _IEEE Transactions on Aerospace and Electronic Systems_ , vol. 49, no. 4, pp. 2281–2295, 2013.
  * [17] C. Stoeckle, J. Munir, A. Mezghani, and J. A. Nossek, “Doa estimation performance and computational complexity of subspace-and compressed sensing-based methods,” in _WSA 2015; 19th International ITG Workshop on Smart Antennas_. VDE, 2015, pp. 1–6.
  * [18] Q. Ren and A. Willis, “Fast root MUSIC algorithm,” _Electronics Letters_ , vol. 33, no. 6, pp. 450 –451, Mar. 1997.
  * [19] K.-Y. Yang, J.-Y. Wu, and W.-H. Li, “A low-complexity direction-of-arrival estimation algorithm for full-dimension massive mimo systems,” in _2014 IEEE International Conference on Communication Systems_. IEEE, 2014, pp. 472–476.
  * [20] H. Meng, Z. Zheng, Y. Yang, K. Liu, and Y. Ge, “A low-complexity 2-d doa estimation algorithm for massive mimo systems,” in _2016 IEEE/CIC international conference on communications in China (ICCC)_. IEEE, 2016, pp. 1–5.
  * [21] C. Jian, S. Wang, and L. Lin, “Two-dimensional doa estimation of coherent signals based on 2d unitary esprit method,” in _2006 8th international Conference on Signal Processing_ , vol. 1. IEEE, 2006.
  * [22] G. Y. Suk, Y.-G. Lim, H. B. Yilmaz, J.-N. Shim, D. K. Kim, and C.-B. Chae, “Low complexity doa estimation in millimeter wave mimo with rf lens,” in _2018 IEEE Wireless Communications and Networking Conference (WCNC)_. IEEE, 2018, pp. 1–6.
  * [23] Y. Xiao, L. Licheng, and Y. Wang, “A new now complexity doa estimation algorithm for massive mimo systems,” in _IEEE International Conference on Consumer Electronics China_ , 2016.
  * [24] A. Wang, L. Liu, and J. Zhang, “Low complexity direction of arrival (doa) estimation for 2d massive mimo systems,” in _2012 IEEE Globecom Workshops_. IEEE, 2012, pp. 703–707.
  * [25] Y. Liu, J. Chai, Y. Zhang, Z. Liu, M. Jin, and T. Qiu, “Low-complexity neural network based doa estimation for wideband signals in massive mimo systems,” _AEU-International Journal of Electronics and Communications_ , vol. 138, p. 153853, 2021.
  * [26] R. Feng, F. Uysal, and A. Yarovoy, “Target localization using mimo-monopulse: Application on 79 ghz fmcw automotive radar,” in _2018 15th European Radar Conference (EuRAD)_. IEEE, 2018, pp. 59–62.
  * [27] K. Schacke, “On the kronecker product,” _Master’s thesis, University of Waterloo_ , 2004.
  * [28] J. Zhen and X. Si, “A method for determining number of coherent signals with arbitrary plane array,” in _The 2010 IEEE International Conference on Information and Automation_. IEEE, 2010, pp. 240–244.
  * [29] C. Duofang, C. Baixiao, and Q. Guodong, “Angle estimation using esprit in mimo radar,” _Electronics Letters_ , vol. 44, no. 12, p. 1, 2008.
  * [30] T. Wu, Y. Li, X. Zhang, Y. Huang, Q. Gu, and C. Fu, “Doa estimation of two-dimensional coherently distributed sources based on spatial smoothing of uniform rectangular arrays,” _International Journal of Antennas and Propagation_ , vol. 2019, 2019.
  * [31] J. Shi, G. Hu, X. Zhang, F. Sun, and Y. Xiao, “Computationally efficient 2d doa estimation with uniform rectangular array in low-grazing angle,” _Sensors_ , vol. 17, no. 3, p. 470, 2017.
  * [32] B. Cantrell, W. Gordon, and G. Trunk, “Maximum likelihood elevation angle estimates of radar targets using subapertures,” _IEEE Transactions on Aerospace and Electronic Systems_ , no. 2, pp. 213–221, 1981.
  * [33] E. M. Stein and R. Shakarchi, _Real analysis: measure theory, integration, and Hilbert spaces_. Princeton University Press, 2009.
  * [34] M.-W. Wu, Y. Jin, Y. Li, T. Song, and P.-Y. Kam, “Maximum-likelihood, magnitude-based, amplitude and noise variance estimation,” _IEEE Signal Processing Letters_ , vol. 28, pp. 414–418, 2021.
  * [35] J. Fei, R. Guangliang, and Z. Zhe, “A new noise variance and post detection snr estimation method for mimo ofdm systems,” in _2008 11th IEEE International Conference on Communication Technology_. IEEE, 2008, pp. 179–182.
  * [36] M. Villano, “Snr and noise variance estimation in polarimetric sar data,” _IEEE Geoscience and Remote Sensing Letters_ , vol. 11, no. 1, pp. 278–282, 2013.
  * [37] HandWiki, ’Computational complexity of mathematical operations’ [Online]. Available: https://handwiki.org/wiki/Computational_complexity_of_mathematical_operations [Accessed: 08- October- 2021].
  * [38] R. P. Brent and P. Zimmermann, _Modern computer arithmetic_. Cambridge University Press, 2010, vol. 18.

|  Md Imrul Hasan (Student Member, IEEE) received a B.Sc. degree from
Bangladesh University of Engineering and Technology (BUET), in March 2016 with
a major in Electrical and Electronic Engineering. He completed his master’s
degree in May 2022 and now pursuing his doctoral degree with a major in
Electrical Engineering from The University of Texas at Dallas. MD IMRUL HASAN
joined LitePoint Corporation, a Teradyne company, in October 2022 as a
Wireless Signal Processing Engineer, where he is currently an integral part of
the engineering development team. Here, he develops and implements advanced
signal-processing algorithms for wireless technologies.
---|---
|  Mohammad Saquib (Senior Member, IEEE) received the B.Sc. degree from the
Bangladesh University of Engineering and Technology, Bangladesh, in 1991, and
the M.S. and Ph.D. degrees from Rutgers University, New Brunswick, NJ, USA, in
1995 and 1998, respectively, all in electrical engineering. He worked as a
member of the technical staff with the Massachusetts Institute of Technology
Lincoln Laboratory, and as an Assistant Professor with Louisiana State
University, Baton Rouge, LA, USA. He is a Professor with Electrical
Engineering Department, The University of Texas at Dallas, Richardson, TX,
USA. His current research interests include the various aspects of wireless
data transmission, radio resource management, and signal processing techniques
for low-cost radar applications.
---|---