# Spatial Covariance Matrix Reconstruction for DOA Estimation in Hybrid
Massive MIMO Systems with Multiple Radio Frequency Chains

Yinsheng Liu, Yiwei Yan, Li You, Wenji Wang, and Hongtao Duan Yinsheng Liu is
with State Key Laboratory of Rail Traffic Control and Safety, Beijing Jiaotong
University, Beijing 100044, China, and National Mobile Communications Research
Laboratory, Southeast University, Nanjing 210096, China (e-mail:
ys.liu@bjtu.edu.cn).Yiwei Yan is with School of Electronic and Information
Engineering, Beijing Jiaotong University, Beijing 100044, China,
(e-mail:19120159@bjtu.edu.cn).Li You and Wenjin Wang are with National Mobile
Communications Research Laboratory, Southeast University, Nanjing 210096,
China, and Purple Mountain Laboratories, Nanjing 211100, China
(e-mail:liyou@seu.edu.cn,wangwj@seu.edu.cn).Hongtao Duan is with Beijing radio
monitoring station of State Radio Monitoring Center (SRMC), Beijing 100037,
China (e-mail: duanht@srrc.org.cn).Corresponding author: Yinsheng Liu.

###### Abstract

Multiple signal classification (MUSIC) has been widely applied in multiple-
input multiple-output (MIMO) receivers for direction-of-arrival (DOA)
estimation. To reduce the cost of radio frequency (RF) chains operating at
millimeter-wave bands, hybrid analog-digital structure has been adopted in
massive MIMO transceivers. In this situation, the received signals at the
antennas are unavailable to the digital receiver, and as a consequence, the
spatial covariance matrix (SCM), which is essential in MUSIC algorithm, cannot
be obtained using traditional sample average approach. Based on our previous
work, we propose a novel algorithm for SCM reconstruction in hybrid massive
MIMO systems with multiple RF chains. By switching the analog beamformers to a
group of predetermined DOAs, SCM can be reconstructed through the solutions of
a set of linear equations. In addition, based on insightful analysis on that
linear equations, a low-complexity algorithm, as well as a careful selection
of the predetermined DOAs, will be also presented in this paper. Simulation
results show that the proposed algorithms can reconstruct the SCM accurately
so that MUSIC algorithm can be well used for DOA estimation in hybrid massive
MIMO systems with multiple RF chains.

###### Index Terms:

DOA estimation, MUSIC, millimeter-wave, massive MIMO, hybrid structure.

##  I Introduction

Direction-of-arrival (DOA) estimation has been widely used in wireless
communications because it can determine the directions of unknown signal
sources [1, 2]. As an important DOA estimation approach, multiple signal
classification (MUSIC) has gained a lot of attention due to its super-
resolution property in the presence of multiple signals [3].

Massive multiple-input multiple-output (MIMO) is one of the most important
enabling technologies in 5G and Beyond 5G communication systems [4]. Due to a
large number of antennas, massive MIMO is essential to millimeter-wave bands
because the large array gain can compensate for the high path loss. With the
help of massive MIMO, the frequency resources at millimeter-wave bands can be
exploited efficiently in 5G and Beyond 5G communication systems [5, 6].

To reduce the number of radio frequency (RF) chains, hybrid structure has been
adopted for massive MIMO operating at millimeter-wave bands [7, 8, 9]. In
hybrid systems, one RF chain is connected to multiple antennas, so that the
number of RF chains can be greatly reduced. However, in hybrid massive MIMO,
the received signals are first fed to the analog phase shifters and then
combined in the analog domain before sent to the digital receiver.
Consequently, the received signals at the antennas are unavailable to the
digital receiver, and the spatial covariance matrix (SCM), which is essential
in MUSIC algorithm, cannot be obtained using the traditional sample average
approach [10]. As MUSIC algorithm is not applicable in hybrid systems, a
straightforward strategy for DOA estimation is to search for the direction
with the maximum received power [11, 12, 13], which, however, is restricted by
Rayleigh limitation [3]. Rayleigh limitation refers to the limitation on the
angle resolution in the presence of multiple signals. It is proportional to
the aperture and thus a large number of antennas are required for better
resolution.

To make use of the super-resolution property of MUSIC algorithm, we have
developed a beam sweeping algorithm for SCM reconstruction in massive MIMO
systems with single RF chain [14]. In this paper, the beam sweeping algorithm
is improved to enable SCM reconstruction in massive MIMO systems with multiple
RF chains. First, the overall SCM is divided into a number of sub-SCMs. Then,
similar to [14], by switching the directions of multiple beamformers to
predetermined DOAs in turn, each sub-SCM can be reconstructed through solving
a set of linear equations. Furthermore, it shows that there are a lot of
redundant calculation in the basic beam sweeping algorithm, because many
repeated entries exist in the SCM. Based on this observation, the algorithm
optimization is further investigated in this paper. First, a low-complexity
beam sweeping algorithm is presented, where the computation complexity can be
reduced to be linearly proportional to the number of antennas. Then, the
selection of predetermined DOAs will be also optimized in this paper. Using
the optimized selection, the number of required predetermined DOAs can be
reduced significantly, so that the SCM reconstruction procedure can be
accomplished in a shorter interval. Simulation results have also been
presented to demonstrate the proposed algorithms.

The rest of this paper is organized as follows. In Section II, signal model
for hybrid massive MIMO is introduced. In Section III, beam sweeping algorithm
is presented and the algorithm optimization is shown in Section IV. Simulation
results can be found in Section V and the conclusions are drawn in Section VI.

##  II System Model

###  II-A Signal Model

As in Fig. 1, consider a hybrid massive MIMO system composed of a uniform
linear array (ULA) with $M$ antennas and $N$ RF chains. Denote $y_{m,n}(t)$ to
be the received signal at the $m$-th antenna of the $n$-th RF chain. Since
each RF chain is connected to $\frac{M}{N}$ antennas, we have
$m=0,1,\cdots,\frac{M}{N}-1$ and $n=0,1,\cdots,N-1$. Then, the received signal
vector by the $n$-th RF chain
$\boldsymbol{y}_{n}(t)=[y_{0,n}(t),y_{1,n}(t),\cdots,y_{\frac{M}{N}-1,n}(t)]^{\mathrm{T}}$
can be represented as

$\displaystyle\boldsymbol{y}_{n}(t)=\sum_{l=0}^{L-1}\boldsymbol{a}_{n}(\theta_{l})x_{l}(t)+\boldsymbol{z}_{n}(t),$
(1)

where $x_{l}(t)$’s ($l=0,1,\cdots,L-1)$ are $L$ narrow-band signals impinging
from far field onto the array, $\theta_{l}$ is the DOA of $x_{l}(t)$,
$\boldsymbol{z}_{n}(t)$ denotes the additive Gaussian noise vector with
$\mathrm{E}\\{\boldsymbol{z}_{n}(t)\boldsymbol{z}_{n}^{\mathrm{H}}(t)\\}=N_{0}\boldsymbol{I}_{\frac{M}{N}}$
where $N_{0}$ is the noise power and $\boldsymbol{I}_{\frac{M}{N}}$ is an
$\frac{M}{N}\times\frac{M}{N}$ identity matrix, and
$\boldsymbol{a}_{n}(\theta_{l})$ is the $\frac{M}{N}\times 1$ steering vector
corresponding to the $n$-th RF chain with the $m$-th entry given by

$\displaystyle
a_{m,n}(\theta_{l})=e^{j2\pi\cdot\frac{d}{\lambda}\cdot\sin\theta_{l}\cdot(n\frac{M}{N}+m)},$
(2)

where $d=\frac{\lambda}{2}$ denotes the antenna distance and $\lambda$ is the
wave length. Take all the RF chains into account, then the overall received
signal vector, the overall steering vector, and the overall additive noise
vector can be represented as

$\displaystyle\boldsymbol{y}(t)$
$\displaystyle=[\boldsymbol{y}_{0}^{\mathrm{T}}(t),\boldsymbol{y}_{1}^{\mathrm{T}}(t),\cdots,\boldsymbol{y}_{N-1}^{\mathrm{T}}(t)]^{\mathrm{T}},$
(3) $\displaystyle\boldsymbol{a}(\theta_{l})$
$\displaystyle=[\boldsymbol{a}_{0}^{\mathrm{T}}(\theta_{l}),\boldsymbol{a}_{1}^{\mathrm{T}}(\theta_{l}),\cdots,\boldsymbol{a}_{N-1}^{\mathrm{T}}(\theta_{l})]^{\mathrm{T}},$
(4) $\displaystyle\boldsymbol{z}(t)$
$\displaystyle=[\boldsymbol{z}_{0}^{\mathrm{T}}(t),\boldsymbol{z}_{1}^{\mathrm{T}}(t),\cdots,\boldsymbol{z}_{N-1}^{\mathrm{T}}(t)]^{\mathrm{T}},$
(5)

respectively. Accordingly, we have

$\displaystyle\boldsymbol{y}(t)=\sum_{l=0}^{L-1}\boldsymbol{a}(\theta_{l})x_{l}(t)+\boldsymbol{z}(t).$
(6)

Denote
$\boldsymbol{R}=\mathrm{E}\\{\boldsymbol{y}(t)\boldsymbol{y}^{\mathrm{H}}(t)\\}$
to be the overall SCM, then using (3), the overall SCM can be divided into

$\displaystyle\boldsymbol{R}=\left[\begin{array}[]{ccc}\boldsymbol{R}_{0,0}&\cdots&\boldsymbol{R}_{0,N-1}\\\
\vdots&\ddots&\vdots\\\
\boldsymbol{R}_{N-1,0}&\cdots&\boldsymbol{R}_{N-1,N-1}\end{array}\right],$
(10)

where
$\boldsymbol{R}_{n_{1},n_{2}}=\mathrm{E}\\{\boldsymbol{y}_{n_{1}}(t)\boldsymbol{y}_{n_{2}}^{\mathrm{H}}(t)\\}$
is the $(n_{1},n_{2})$-th sub-SCM. Assuming $L$ signals are mutually
independent with zero means and the power of the $l$-th signal is
$\mathrm{E}\\{|x_{l}(t)|^{2}\\}=\sigma_{l}^{2}$, then the $(n_{1},n_{2})$-th
sub-SCM and the overall SCM will be

$\displaystyle\boldsymbol{R}_{n_{1},n_{2}}$
$\displaystyle=\sum_{l=0}^{L-1}\sigma_{l}^{2}\cdot\boldsymbol{a}_{n_{1}}(\theta_{l})\boldsymbol{a}_{n_{2}}^{\mathrm{H}}(\theta_{l})+\delta[n_{1}-n_{2}]N_{0}\boldsymbol{I}_{\frac{M}{N}},$
(11) $\displaystyle\boldsymbol{R}$
$\displaystyle=\sum_{l=0}^{L-1}\sigma_{l}^{2}\cdot\boldsymbol{a}(\theta_{l})\boldsymbol{a}^{\mathrm{H}}(\theta_{l})+N_{0}\boldsymbol{I}_{M},$
(12)

respectively, where $\delta[\cdot]$ indicates Kronecker Delta function.

Figure 1: For a hybrid massive MIMO with multiple RF chains, each RF chain is
connected to a subset of the antennas.

###  II-B Review of MUSIC Algorithm

Denote $\boldsymbol{y}[k]=\boldsymbol{y}(kT_{s})$ to be the sample of the
received signal where $T_{s}$ denotes the sampling period. In MUSIC algorithm,
$\boldsymbol{y}[k]$’s are assumed to be available in the receiver. In this
case, the overall SCM in (12) can be estimated using the sample average
approach, that is [3, 15]

$\displaystyle\boldsymbol{R}\approx\frac{1}{K}\sum_{k=0}^{K-1}\boldsymbol{y}[k]\boldsymbol{y}^{\mathrm{H}}[k],$
(13)

where $K$ denotes the number of samples. The eigenvalue decomposition of the
overall SCM can be given as
$\boldsymbol{R}=(\boldsymbol{U}_{s},\boldsymbol{U}_{n})\boldsymbol{\Lambda}_{s}(\boldsymbol{U}_{s},\boldsymbol{U}_{n})^{\mathrm{H}},$
where $\boldsymbol{U}_{s}$ and $\boldsymbol{U}_{n}$ denote the orthogonal base
vectors corresponding to the signal and the noise subspaces respectively, and
$\boldsymbol{\Lambda}_{s}$ is a diagonal matrix composed of the eigenvalues.
Then, unknown DOAs can be determined by searching the peak values of
$p(\theta)$,

$\displaystyle
p(\theta)={\|\boldsymbol{U}_{n}^{\mathrm{H}}\boldsymbol{a}(\theta)\|_{2}^{-2}},~{}~{}~{}\theta\in[-90^{\mathrm{o}},90^{\mathrm{o}}].$
(14)

For sample average in (13), $\boldsymbol{y}[k]$ is required to estimate SCM.
In this case, the received signals at all antennas should be sent via RF
chains to the digital receiver. In hybrid massive MIMO, however, Fig. 1 shows
that only the combination of the entries of $\boldsymbol{y}_{n}[k]$ can be
sent to the digital receiver because the number of RF chains is smaller than
the number of antennas. As a consequence, the sample average approach in (13)
cannot be used in hybrid massive MIMO systems. In [14], we have developed an
algorithm for SCM reconstruction in massive MIMO with single RF chain. In the
presence of multiple RF chains, however, each RF chain is connected to only a
subset of the antennas. Therefore, the overall SCM cannot be obtained directly
using the algorithm in [14].

##  III Beam Sweeping Algorithm

To reconstruct the overall SCM in the presence of multiple RF chains, we can
first reconstruct the sub-SCMs individually, and then the overall SCM can be
obtained following (10).

To reconstruct $\boldsymbol{R}_{n_{1},n_{2}}$, define
$\\{\theta^{(0)},\theta^{(1)},\cdots,\theta^{(Q-1)}\\}$ as a set of
predetermined DOA angles. The analog beamformers switch the beam directions to
the predetermined DOAs in turn. For the $q$-th sweeping beam, the
predetermined DOA is $\theta^{(q)}$, and thus the steering vector
corresponding to the $n$-th RF chain is $\boldsymbol{a}_{n}(\theta^{(q)})$.
The combination of the received signals on the $n$-th RF chain can be
represented by

$\displaystyle
c_{n}^{(q)}(t)=\boldsymbol{a}_{n}^{\mathrm{H}}(\theta^{(q)})\boldsymbol{y}_{n}(t).$
(15)

From Fig. 1, the $n$-th signal combination is sampled before sent to the
receiver, and thus the samples of the signal combination can be given by

$\displaystyle
c_{n}^{(q)}[k]=c_{n}^{(q)}(kT_{s})=\boldsymbol{a}_{n}^{\mathrm{H}}(\theta^{(q)})\boldsymbol{y}_{n}[k].$
(16)

Denote $P_{n_{1},n_{2}}^{(q)}$ to be the correlation between the outputs of
the $n_{1}$-th and the $n_{2}$-th RF chains, that is

$\displaystyle
P_{n_{1},n_{2}}^{(q)}=\frac{1}{K}\sum_{k=0}^{K-1}c_{n_{1}}^{(q)}[k]c_{n_{2}}^{(q)*}[k]$
$\displaystyle=\boldsymbol{a}_{n_{1}}^{\mathrm{H}}(\theta^{(q)})\left(\frac{1}{K}\sum_{k=0}^{K-1}\boldsymbol{y}_{n_{1}}[k]\boldsymbol{y}_{n_{2}}^{\mathrm{H}}[k]\right)\boldsymbol{a}_{n_{2}}(\theta^{(q)}).$
(17)

If the number of samples is large enough, the sample average in (III) is
equivalent to the statistical average, that is

$\displaystyle\boldsymbol{a}_{n_{1}}^{\mathrm{H}}(\theta^{(q)})\boldsymbol{R}_{n_{1},n_{2}}\boldsymbol{a}_{n_{2}}(\theta^{(q)})=P_{n_{1},n_{2}}^{(q)}.$
(18)

Using the $\mathrm{vec}(\cdot)$ operator to (18), the left-hand-side of (18)
can be given as

$\displaystyle\mathrm{vec}\\{\boldsymbol{a}_{n_{1}}^{\mathrm{H}}(\theta^{(q)})\boldsymbol{R}_{n_{1},n_{2}}\boldsymbol{a}_{n_{2}}(\theta^{(q)})\\}$
$\displaystyle=[\boldsymbol{a}_{n_{2}}^{\mathrm{T}}(\theta^{(q)})\otimes\boldsymbol{a}_{n_{1}}^{\mathrm{H}}(\theta^{(q)})]\mathrm{vec}(\boldsymbol{R}_{n_{1},n_{2}}),$
(19)

where we have used the equation (1.10.25) in [16], that is,
$\mathrm{vec}(\boldsymbol{BCD})=(\boldsymbol{D}^{\mathrm{T}}\otimes\boldsymbol{B})\mathrm{vec}(\boldsymbol{C})$
with $\otimes$ denoting the Kronecker product.

To proceed, denote
$\boldsymbol{r}_{n_{1},n_{2}}=\mathrm{vec}(\boldsymbol{R}_{n_{1},n_{2}})$ and

$\displaystyle\boldsymbol{a}_{n_{1},n_{2}}^{(q)}=\boldsymbol{a}_{n_{2}}(\theta^{(q)})\otimes\boldsymbol{a}_{n_{1}}^{*}(\theta^{(q)}),$
(20)

where $(\cdot)^{*}$ denotes the element-wise conjugation. Apparently, both
$\boldsymbol{r}_{n_{1},n_{2}}$ and $\boldsymbol{a}_{n_{1},n_{2}}^{(q)}$ are
$\frac{M^{2}}{N^{2}}\times 1$ vectors. Then, (18) can be rewritten as

$\displaystyle(\boldsymbol{a}_{n_{1},n_{2}}^{(q)})^{\mathrm{T}}\boldsymbol{r}_{n_{1},n_{2}}=P_{n_{1},n_{2}}^{(q)}.$
(21)

Considering that there are $Q$ predetermined DOA angles, then (21) can be
extended to a group of linear equations as

$\displaystyle\boldsymbol{A}_{n_{1},n_{2}}\boldsymbol{r}_{n_{1},n_{2}}=\boldsymbol{p}_{n_{1},n_{2}},$
(22)

where $\boldsymbol{A}_{n_{1},n_{2}}$ is a $Q\times\frac{M^{2}}{N^{2}}$ matrix
and $\boldsymbol{p}_{n_{1},n_{2}}$ is a $Q\times 1$ vector

$\displaystyle\boldsymbol{A}_{n_{1},n_{2}}=\left[\boldsymbol{a}_{n_{1},n_{2}}^{(0)},\boldsymbol{a}_{n_{1},n_{2}}^{(1)},\cdots,\boldsymbol{a}_{n_{1},n_{2}}^{(Q-1)}\right]^{\mathrm{T}},$
(23)
$\displaystyle\boldsymbol{p}_{n_{1},n_{2}}=\left[P_{n_{1},n_{2}}^{(0)},P_{n_{1},n_{2}}^{(1)},\cdots,P_{n_{1},n_{2}}^{(Q-1)}\right]^{\mathrm{T}}.$
(24)

Then, the $(n_{1},n_{2})$-th sub-SCM can be reconstructed by solving (22).
Similar to [14], a diagonal loading coefficient can be adopted to improve the
distribution of the eigenvalues, so that the ill-conditioned result can be
avoided. In this case, the vector form of sub-SCM can be obtained as

$\displaystyle\widehat{\boldsymbol{r}}_{n_{1},n_{2}}=(\boldsymbol{A}_{n_{1},n_{2}}^{\mathrm{H}}\boldsymbol{A}_{n_{1},n_{2}}+\sigma^{2}\boldsymbol{I}_{\frac{M^{2}}{N^{2}}})^{-1}\boldsymbol{A}_{n_{1},n_{2}}^{\mathrm{H}}\boldsymbol{p}_{n_{1},n_{2}},$
(25)

where $\sigma^{2}$ denotes the diagonal loading coefficient. Then, the
$(n_{1},n_{2})$-th sub-SCM can be reconstructed through
$\widehat{\boldsymbol{R}}_{n_{1},n_{2}}=\mathrm{unvec}(\widehat{\boldsymbol{r}}_{n_{1},n_{2}}).$
Once the sub-SCMs are obtained, the overall SCM can be reconstructed,
following (10), as

$\displaystyle\widehat{\boldsymbol{R}}=\left[\begin{array}[]{ccc}\widehat{\boldsymbol{R}}_{0,0}&\cdots&\widehat{\boldsymbol{R}}_{0,N-1}\\\
\vdots&\ddots&\vdots\\\
\widehat{\boldsymbol{R}}_{N-1,0}&\cdots&\widehat{\boldsymbol{R}}_{N-1,N-1}\end{array}\right].$
(29)

As in [14], although matrix inversion in (25) causes a huge computational
burden, the operator
$(\boldsymbol{A}_{n_{1},n_{2}}^{\mathrm{H}}\boldsymbol{A}_{n_{1},n_{2}}+\sigma^{2}\boldsymbol{I}_{\frac{M^{2}}{N^{2}}})^{-1}\boldsymbol{A}_{n_{1},n_{2}}^{\mathrm{H}}$
can be pre-calculated off-line if the predetermined DOAs are fixed. In this
case, matrix inversion can be avoided from on-line calculation, and the
computational burden is mainly caused by the matrix-vector production in (25),
which requires $QM^{2}/N^{2}$ complex multiplications.

With the reconstructed SCM in (29), unknown DOAs can obtained using (14) if
classical MUSIC algorithm is adopted. In addition to the classical MUSIC
algorithm, the variants of MUSIC algorithm, such as root-MUSIC [17], can be
also used. To use root-MUSIC and other variants, we only need to replace
sample average approach in those algorithms with the SCM reconstruction
algorithm in this paper.

##  IV Algorithm Optimization

Although the beam sweeping algorithm in Section III can reconstruct the SCM
successfully, the computation complexity is still huge due to the large
dimension of matrix-vector production in (25). In this section, a low-
complexity algorithm will be presented where the dimension of matrix-vector
production can be reduced significantly. Based on the low-complexity
algorithm, optimized selection of predetermined DOAs is further investigated.
Using the optimized predetermined DOAs, the matrix inverse operation can be
even avoided.

###  IV-A Low-complexity Implementation

Low-complexity implementation of the beam sweeping algorithm is inspired by
the fact that many repeated entries exist in the SCM. Although there are
$\frac{M^{2}}{N^{2}}$ entries in $\boldsymbol{R}_{n_{1},n_{2}}$, the number of
non-repeated unknowns is only $\frac{2M}{N}-1$. Therefore, the computational
complexity can be reduced if we only recover the non-repeated unknowns.

Denote the $(m_{1},m_{2})$-th entry of $\boldsymbol{R}_{n_{1},n_{2}}$,
$\boldsymbol{R}_{n_{1},n_{2}}[m_{1},m_{2}]$, to be
$\gamma_{n_{1},n_{2}}[m_{1}-m_{2}]$. If denote
$\boldsymbol{\gamma}_{n_{1},n_{2}}$ to be a $(\frac{2M}{N}-1)\times 1$ column
vector containing all non-repeated unknowns in $\boldsymbol{r}_{n_{1},n_{2}}$
or $\boldsymbol{R}_{n_{1},n_{2}}$, then $\boldsymbol{\gamma}_{n_{1},n_{2}}$
can be given by

$\displaystyle\boldsymbol{\gamma}_{n_{1},n_{2}}=\left(\gamma_{n_{1},n_{2}}\left[1-\frac{M}{N}\right],\cdots,\gamma_{n_{1},n_{2}}\left[\frac{M}{N}-1\right]\right)^{\mathrm{T}}.$
(30)

Accordingly, $\boldsymbol{r}_{n_{1},n_{2}}$ can be expressed by
$\boldsymbol{\gamma}_{n_{1},n_{2}}$ using

$\displaystyle\boldsymbol{r}_{n_{1},n_{2}}=\boldsymbol{E}\cdot\boldsymbol{\gamma}_{n_{1},n_{2}},$
(31)

where $\boldsymbol{E}$ is an
$\frac{M^{2}}{N^{2}}\times\left(\frac{2M}{N}-1\right)$ matrix

$\displaystyle\boldsymbol{E}=\left[\begin{array}[]{ccc}\boldsymbol{O}_{\frac{M}{N}-1}&\boldsymbol{I}_{\frac{M}{N}}&\boldsymbol{O}_{0}\\\
\boldsymbol{O}_{\frac{M}{N}-2}&\boldsymbol{I}_{\frac{M}{N}}&\boldsymbol{O}_{1}\\\
\vdots&\vdots&\vdots\\\
\boldsymbol{O}_{0}&\boldsymbol{I}_{\frac{M}{N}}&\boldsymbol{O}_{\frac{M}{N}-1}\end{array}\right]$
(36)

with $\boldsymbol{O}_{m}$ denoting an $\frac{M}{N}\times m$ all-zero matrix.
By substituting (31) into (21), (21) can be rewritten as

$\displaystyle{(\boldsymbol{b}_{n_{1},n_{2}}^{(q)})^{\mathrm{T}}}\boldsymbol{\gamma}_{n_{1},n_{2}}=P_{n_{1},n_{2}}^{(q)},$
(37)

where $(\boldsymbol{b}_{n_{1},n_{2}}^{(q)})^{\mathrm{T}}$ is a
$1\times\left(\frac{2M}{N}-1\right)$ row vector given by

$\displaystyle(\boldsymbol{b}_{n_{1},n_{2}}^{(q)})^{\mathrm{T}}=(\boldsymbol{a}_{n_{1},n_{2}}^{(q)})^{\mathrm{T}}\boldsymbol{E}.$
(38)

Then, similar to (22), if we take $Q$ predetermined DOAs into account,
equation (37) can be extended to a group of linear equations as

$\displaystyle\boldsymbol{B}_{n_{1},n_{2}}\boldsymbol{\gamma}_{n_{1},n_{2}}=\boldsymbol{p}_{n_{1},n_{2}},$
(39)

where $\boldsymbol{B}_{n_{1},n_{2}}$ is a $Q\times(\frac{2M}{N}-1)$ matrix
given by

$\displaystyle\boldsymbol{B}_{n_{1},n_{2}}=\left[\boldsymbol{b}_{n_{1},n_{2}}^{(0)},\boldsymbol{b}_{n_{1},n_{2}}^{(1)},\cdots,\boldsymbol{b}_{n_{1},n_{2}}^{(Q-1)})\right]^{\mathrm{T}}.$
(40)

Different from (25) where a diagonal loading coefficient is adopted to avoid
the ill-conditioned result, equation (39) can be solved without diagonal
loading if the predetermined DOAs are carefully selected, as will be discussed
in the next subsection. Therefore, $\boldsymbol{\gamma}_{n_{1},n_{2}}$ in (39)
can be obtained as

$\displaystyle\widehat{\boldsymbol{\gamma}}_{n_{1},n_{2}}=(\boldsymbol{B}_{n_{1},n_{2}}^{\mathrm{H}}\boldsymbol{B}_{n_{1},n_{2}})^{-1}\boldsymbol{B}_{n_{1},n_{2}}^{\mathrm{H}}\boldsymbol{p}_{n_{1},n_{2}}.$
(41)

Consequently, the sub-SCM can be reconstructed as follows

$\displaystyle\widehat{\boldsymbol{R}}_{n_{1},n_{2}}[m_{1},m_{2}]=\widehat{\gamma}_{n_{1},n_{2}}[m_{1}-m_{2}].$
(42)

Similar to (25), the operator
$(\boldsymbol{B}_{n_{1},n_{2}}^{\mathrm{H}}\boldsymbol{B}_{n_{1},n_{2}})^{-1}\boldsymbol{B}_{n_{1},n_{2}}^{\mathrm{H}}$
in (41) can be pre-calculated off-line so that the computational burden in
(41) is mainly caused by the matrix-vector production. Since
$\boldsymbol{B}_{n_{1},n_{2}}$ is a $Q\times(\frac{2M}{N}-1)$ matrix, the
number of complex multiplications required in (41) is only $Q(2M/N-1)$, which
is much lower than that required in (25). In addition to the reduction of
complexity, (39) also indicates that at least $\frac{2M}{N}-1$ predetermined
DOAs are required to achieve accurate reconstruction, that is,
$Q\geq\frac{2M}{N}-1.$ This is because there are $\frac{2M}{N}-1$ non-repeated
unknowns in $\boldsymbol{\gamma}_{n_{1},n_{2}}$, at least $\frac{2M}{N}-1$
equations are required to solve (39) with each equation corresponding to one
predetermined DOA.

###  IV-B Selection of Predetermined DOAs

In [14], the predetermined DOAs are selected as uniformly distributed from
$-90^{\circ}$ to $90^{\circ}$. Although simple, we will show in this
subsection that the selection of the predetermined DOAs can be further
optimized. Denote

$\displaystyle v^{(q)}=\frac{d}{\lambda}\sin\theta^{(q)}=0.5\sin\theta^{(q)},$
(43)

to be the space frequency corresponding to $\theta^{(q)}$. In this paper, the
predetermined DOAs are selected such that $v^{(q)}$ are uniformly distributed
from $-0.5$ to $0.5$, that is, $v^{(q)}=-0.5+q/Q$. As a result, the
predetermined DOAs are determined as

$\displaystyle\theta^{(q)}=\arcsin(-1+2q/Q).$ (44)

In following, it will show
$\boldsymbol{B}_{n_{1},n_{2}}^{\mathrm{H}}\boldsymbol{B}_{n_{1},n_{2}}$ can be
converted to a diagonal matrix using the predetermined DOAs in (44) so that
the matrix inverse operation in (41) can be avoided.

If using (20) and (36) to (38), we have

$\displaystyle(\boldsymbol{b}_{n_{1},n_{2}}^{(q)})^{\mathrm{T}}$
$\displaystyle=(\boldsymbol{a}_{n_{2}}^{\mathrm{T}}(\theta^{(q)})\otimes\boldsymbol{a}_{n_{1}}^{\mathrm{H}}(\theta^{(q)}))\cdot\boldsymbol{E}$
$\displaystyle=[a_{0,n_{2}}(\theta^{(q)})\boldsymbol{a}_{n_{1}}^{\mathrm{H}}(\theta^{(q)}),\cdots,a_{\frac{M}{N}-1,n_{2}}(\theta^{(q)})\boldsymbol{a}_{n_{1}}^{\mathrm{H}}(\theta^{(q)})]$
$\displaystyle\cdot\left[\begin{array}[]{ccc}\boldsymbol{O}_{\frac{M}{N}-1}&\boldsymbol{I}_{\frac{M}{N}}&\boldsymbol{O}_{0}\\\
\vdots&\vdots&\vdots\\\
\boldsymbol{O}_{0}&\boldsymbol{I}_{\frac{M}{N}}&\boldsymbol{O}_{\frac{M}{N}-1}\end{array}\right]$
(48) $\displaystyle=$
$\displaystyle\sum_{m=0}^{\frac{M}{N}-1}a_{m,n_{2}}(\theta^{(q)})\boldsymbol{a}_{n_{1}}^{\mathrm{H}}(\theta^{(q)})\left[\boldsymbol{O}_{\frac{M}{N}-1-m},\boldsymbol{I}_{\frac{M}{N}},\boldsymbol{O}_{m}\right]$
$\displaystyle=$
$\displaystyle\sum_{m=0}^{\frac{M}{N}-1}\left[\boldsymbol{0}_{\frac{M}{N}-1-m}^{\mathrm{T}},a_{m,n_{2}}(\theta^{(q)})\boldsymbol{a}_{n_{1}}^{\mathrm{H}}(\theta^{(q)}),\boldsymbol{0}_{m}^{\mathrm{T}}\right],$
(49)

where $\boldsymbol{0}_{m}$ is an $m\times 1$ all-zero vector. If denote
${b}_{n_{1},n_{2}}^{(q)}[m_{0}]$ as the $m_{0}$-th entry of
$\boldsymbol{b}_{n_{1},n_{2}}^{(q)}$, then we have
$m_{0}=0,1,\cdots,\frac{2M}{N}-2$ since the length of
$\boldsymbol{b}_{n_{1},n_{2}}^{(q)}$ is $\frac{2M}{N}-1$. To proceed, define
$\widetilde{a}_{m,n}$ as a sequence with infinite length where

$\displaystyle\widetilde{a}_{m,n}=\begin{cases}{a}_{m,n}(\theta^{(q)}),&m=0,1,\cdots,\frac{M}{N}-1,\\\
0,&\mathrm{others},\end{cases}.$ (50)

Then, from the last equation of (IV-B), it is easy to verify that
${b}_{n_{1},n_{2}}^{(q)}[m_{0}]$ can be obtained as

$\displaystyle{b}_{n_{1},n_{2}}^{(q)}[m_{0}]=\sum_{m=0}^{\frac{M}{N}-1}\widetilde{a}_{\frac{M}{N}-1-m_{0}+m,n_{2}}\cdot\widetilde{a}_{m,n_{1}}^{*},$
(51)

which is essentially a linear convolution between
$\widetilde{a}_{m_{0},n_{1}}^{*}$ and
$\widetilde{a}_{\frac{M}{N}-1-m_{0},n_{2}}$. Using (51), the
$(m_{1},m_{2})$-th entry of
$\boldsymbol{B}_{n_{1},n_{2}}^{\mathrm{H}}\boldsymbol{B}_{n_{1},n_{2}}$ with
$m_{1},m_{2}=0,1,\cdots,\frac{2M}{N}-2$ is given by

$\displaystyle[\boldsymbol{B}_{n_{1},n_{2}}^{\mathrm{H}}\boldsymbol{B}_{n_{1},n_{2}}]_{(m_{1},m_{2})}=\sum_{q=0}^{Q-1}[(\boldsymbol{b}_{n_{1},n_{2}}^{(q)})^{*}(\boldsymbol{b}_{n_{1},n_{2}}^{(q)})^{\mathrm{T}}]_{(m_{1},m_{2})}$
$\displaystyle=\sum_{q=0}^{Q-1}b_{n_{1},n_{2}}^{(q)*}[m_{1}]\cdot
b_{n_{1},n_{2}}^{(q)}[m_{2}]$
$\displaystyle=\sum_{i_{1}=0}^{\frac{M}{N}-1}\sum_{i_{2}=0}^{\frac{M}{N}-1}\sum_{q=0}^{Q-1}\widetilde{a}_{\frac{M}{N}-1-m_{1}+i_{1},n_{2}}^{*}\widetilde{a}_{i_{1},n_{1}}\widetilde{a}_{\frac{M}{N}-1-m_{2}+i_{2},n_{2}}\widetilde{a}_{i_{2},n_{1}}^{*}.$
(52)

Since the sequence $\widetilde{a}_{m,n}$ has non-zero values only when
$m=0,1,\cdots,\frac{M}{N}-1$, $i_{1}$ in (IV-B) should satisfy

$\displaystyle 0\leq i_{1}\leq M/N-1,$ (53) $\displaystyle 0\leq
M/N-1-m_{1}+i_{1}\leq M/N-1,$ (54)

simultaneously. Therefore, we can obtain

$\displaystyle U_{1}^{-}\leq i_{1}\leq U_{1}^{+},$ (55)

where $U_{1}^{-}=\mathrm{max}\\{0,m_{1}-(\frac{M}{N}-1)\\}$ and
$U_{1}^{+}=\mathrm{min}\\{m_{1},\frac{M}{N}-1\\}$. Similarly, we have

$\displaystyle U_{2}^{-}\leq i_{2}\leq U_{2}^{+},$ (56)

where $U_{2}^{-}=\mathrm{max}\\{0,m_{2}-(\frac{M}{N}-1)\\}$ and
$U_{2}^{+}=\mathrm{min}\\{m_{2},\frac{M}{N}-1\\}$. With the constraints of
$i_{1}$ and $i_{2}$, (IV-B) can be rewritten as

$\displaystyle[\boldsymbol{B}_{n_{1},n_{2}}^{\mathrm{H}}\boldsymbol{B}_{n_{1},n_{2}}]_{(m_{1},m_{2})}=\sum_{i_{1}=U_{1}^{-}}^{U_{1}^{+}}\sum_{i_{2}=U_{2}^{-}}^{U_{2}^{+}}\sum_{q=0}^{Q-1}e^{j2\pi
v_{q}(m_{1}-m_{2})},$ (57)

where we have used the identity

$\displaystyle
a_{\frac{M}{N}-1-m_{1}+i_{1},n_{2}}^{*}(\theta^{(q)})a_{i_{1},n_{1}}(\theta^{(q)})\cdot$
$\displaystyle
a_{\frac{M}{N}-1-m_{2}+i_{2},n_{2}}(\theta^{(q)})a_{i_{2},n_{1}}^{*}(\theta^{(q)})=e^{j2\pi
v_{q}(m_{1}-m_{2})}.$ (58)

If applying the predetermined DOAs in (44) to (57), we have

$\displaystyle\sum_{q=0}^{Q-1}e^{j2\pi(-0.5+q/Q)(m_{1}-m_{2})}=\begin{cases}Q,&m_{1}=m_{2}\\\
0,&m_{1}\neq m_{2}\end{cases}.$ (59)

As a result, we have
$[\boldsymbol{B}_{n_{1},n_{2}}^{\mathrm{H}}\boldsymbol{B}_{n_{1},n_{2}}]_{(m_{1},m_{2})}=0$
for $m_{1}\neq m_{2}$ and when $m_{1}=m_{2}=0,1,\cdots,2\frac{M}{N}-2$,

$\displaystyle[\boldsymbol{B}_{n_{1},n_{2}}^{\mathrm{H}}\boldsymbol{B}_{n_{1},n_{2}}]_{(m_{1},m_{1})}=Q(U_{1}^{+}-U_{1}^{-}+1)^{2}$
$\displaystyle=\begin{cases}Q(m_{1}+1)^{2},&m_{1}\leq\frac{M}{N}-1\\\
Q\left[2\left(\frac{M}{N}-1\right)-m_{1}+1\right]^{2},&m_{1}>\frac{M}{N}-1\end{cases}.$
(60)

Apparently,
$\boldsymbol{B}_{n_{1},n_{2}}^{\mathrm{H}}\boldsymbol{B}_{n_{1},n_{2}}$ has
been converted into a diagonal matrix by selecting predetermined DOAs as in
(44), and correspondingly, the matrix inverse operation in (41) can be avoided
completely.

##  V Simulation Results

Computer simulation is adopted in this section to investigate the proposed
algorithms. We consider a ULA with $M=64$ antennas and the distance between
antennas is $0.5\lambda$. $32$ signals are impinging onto the ULA where DOAs
of the signals are uniformly distributed from $-90^{\circ}$ to $90^{\circ}$.
The arrival signals are assumed independent with zero means and unit powers.
Without specification, the signal-to-noise ratio (SNR) is $-5$ dB, and the
predetermined DOAs are as in (44). Similar to [14], normalized-square-error
(NSE) is used in the simulation to evaluate the accuracy of reconstructed SCM,
that is
$\mathrm{NSE}=\|\widehat{\boldsymbol{R}}-\boldsymbol{R}\|_{\mathrm{F}}^{2}\cdot\|\boldsymbol{R}\|_{\mathrm{F}}^{-2}.$
To demonstrate the effectiveness of reconstructed SCM, classical MUSIC
algorithm is adopted in the simulation for DOA estimation. Accordingly, means-
quare-error (MSE) is used to evaluate the accuracy of DOA estimation, that is
$\mathrm{MSE}=\mathrm{E}\\{|\widehat{\theta}-\theta|^{2}\\}.$

Figure 2: NSE versus the number of predetermined DOAs.

For the basic reconstruction beam sweeping algorithm as in (25), the diagonal
loading coefficient is fixed as $\sigma^{2}=1$.

(a) (b) (c)

Figure 3: (a) NSE versus the sample number (b) MSE of MUSIC using
reconstructed SCM (c) number of complex multiplications versus antenna number.

Fig. 2 shows the reconstruction accuracy with different numbers of
predetermined DOAs. The number of samples is fixed as $K=5000$. For the basic
algorithm, when $Q$ is small, the reconstruction accuracy can be improved
significantly as the rising of $Q$. When $Q>\frac{2M}{N}-1$, there have been
sufficient predetermined DOAs and thus the NSE can be hardly reduced further
by increasing $Q$. For the low-complexity algorithm as in (41), when
$Q<\frac{2M}{N}-1$,
$\boldsymbol{B}_{n_{1},n_{2}}^{\mathrm{H}}\boldsymbol{B}_{n_{1},n_{2}}$ is
ill-conditioned, and thus the low-complexity algorithm is only available when
$Q\geq\frac{2M}{N}-1$. It shows that the low-complexity algorithm can achieve
minimum NSE when $Q=\frac{2M}{N}-1$. Therefore, the reconstruction procedure
can be accomplished within $(\frac{2M}{N}-1)K$ samples. For comparison, Fig. 2
also presents the existing algorithm for the single RF chain case [14].
Different from (44), the predetermined DOAs are selected as uniformly
distributed from $-90^{\circ}$ to $90^{\circ}$. In this situation, more
predetermined DOAs are required to achieve the performances of the proposed
algorithms in this paper, leading to a much longer procedure for SCM
reconstruction.

Fig. 3 (a) shows the reconstruction accuracy with different number of samples.
The number of predetermined DOAs is fixed as $Q=\frac{2M}{N}-1$. For the low-
complexity algorithm, the NSE can be reduced as the increasing of the number
of samples. For the basic algorithm ($N=4$ in particular), however, the
reconstruction accuracy can be hardly improved when $K$ is large enough. This
observation actually coincides with the result in Fig. 2. Fig. 3 (a) also
shows that better accuracy can be obtained with a smaller number of RF chains.
This is because $Q\cdot K$ samples are adopted in overall to reconstruct the
SCM. As $Q=\frac{2M}{N}-1$, more samples will be adopted when $N$ is small,
and thus the NSE will be improved. For comparison, the traditional sample
average algorithm is also included. As expected, the sample average algorithm
has the worst performance because it adopts only $K$ samples.

Fig. 3 (b) shows the performance of MUSIC algorithm based on the reconstructed
SCM. The sample-average based MUSIC algorithm [3] and the H-MUSIC algorithm
[18] are also included for comparison. For the proposed algorithm, the number
of predetermined DOAs is fixed as $Q=\frac{2M}{N}-1$. Since the NSE can be
reduced as the reduction of the number of the RF chains, the DOA estimation
accuracy can be also improved for massive MIMO systems with small number of RF
chains. It also shows that the proposed algorithm can achieve even better
performance than sample-average based MUSIC algorithm. This is because the
sample average approach is not as accurate as the proposed SCM reconstruction
algorithm, as shown in Fig. 3 (a). Although H-MUSIC algorithm can be also used
for DOA estimation in hybrid massive MIMO, the estimation accuracy is even
worse than the classic MUSIC algorithm. This observation coincides with the
result in [18].

A comparison on the computation complexity for the basic beam sweeping
algorithm and the low-complexity algorithm is shown in Fig. 3 (c). As
expected, the low-complexity algorithm in this paper can reduce the
computational burden due to the reduction of dimension of matrix-vector
production. It also shows that the complexity reduction is more significant
when the number of RF chains is smaller.

##  VI Conclusions

In this paper, a beam sweeping approach has been proposed to reconstruct the
SCM, so that MUSIC algorithm can be applied for DOA estimation in hybrid
massive MIMO systems with multiple RF chains. We have presented the basic
algorithm that can be used in the case with multiple RF chains. Then, a low-
complexity algorithm has been also introduced by removing the repeated entries
in the SCM. In addition, the selection of the predetermined DOAs has also been
optimized and the matrix inversion in the low-complexity algorithm can be
further avoided. Simulation results have shown that the proposed approach can
achieve better performance than existing baselines, and the performance of
MUSIC algorithm can be also improved accordingly.

## References

  * [1] T. E. Tuncer and B. Friedlander, _Classical and Modern Direction-of-Arrival Estimation_. Academic, Orlando, 2009.
  * [2] F. Shu, X. Wu, J. Hu, J. Li, R. Chen, and J. Wang, “Secure and precise wireless transmission for random-subcarrier-selection-based directional modulation,” _IEEE J. Sel. Areas Commun._ , vol. 36, no. 4, pp. 890–904, July 2018.
  * [3] R. O. Schmidt, “Multiple emitter location and signal parameter estimation,” _IEEE Trans. Antennas Propag._ , no. 3, pp. 276–280, Mar. 1986.
  * [4] E. G. Larsson, F. Tufvesson, O. Edfors, and T. L. Marzetta, “Massive MIMO for next generation wireless systems,” _IEEE Commun. Mag._ , vol. 52, no. 2, pp. 186–195, Feb. 2014.
  * [5] W. R. _et al_., “Millimeter-wave beamforming as an enabling technology for 5G cellular communications: Theoretical feasibility and prototype results,” _IEEE Commun. Mag._ , vol. 52, no. 2, pp. 106–113, Feb. 2014.
  * [6] L. You, X. Q. Gao, G. Y. Li, X.-G. Xia, and N. Ma, “BDMA for millimeter-wave/Terahertz massive MIMO transmission with per-beam synchronization,” _IEEE J. Sel. Areas Commun._ , vol. 35, no. 7, pp. 1550–1563, Jul. 2017.
  * [7] O. E. Ayach, S. Rajagopal, S. Abu-Surra, Z. Pi, and R. W. Heath, “Spatially sparse precoding in millimeter wave MIMO systems,” _IEEE Trans. Wireless Commun._ , vol. 13, no. 3, pp. 1499–1513, Mar. 2014.
  * [8] V. Venkateswaran and A. J. van der Veen, “Analog beamforming in MIMO communications with phase shift networks and online channel estimation,” _IEEE Trans. Signal Process._ , vol. 58, no. 8, pp. 4131–4143, Aug. 2010\.
  * [9] C. Lin and G. Y. Li, “Adaptive beamforming with resource allocation for distance-aware multi-user indoor Terahertz communications,” _IEEE Trans. Commun._ , vol. 63, no. 8, pp. 2985–2995, Aug 2015.
  * [10] S. Chuang, W. Wu, and Y. Liu, “High-resolution AoA estimation for hybrid antenna arrays,” _IEEE Trans. Antennas Propag._ , vol. 63, no. 7, pp. 2955–2968, July 2015.
  * [11] F. Shu, Y. Qin, T. Liu, L. Gui, Y. Zhang, J. Li, and Z. Han, “Low-complexity and high-resolution DOA estimation for hybrid analog and digital massive MIMO receive array,” _IEEE Trans. Commun._ , vol. 66, no. 6, pp. 2487–2501, June 2018.
  * [12] K. Aghababaiyan, V. Shah-Mansouri, and B. Maham, “High-precision OMP-based direction of arrival estimation scheme for hybrid non-uniform array,” _IEEE Commun. Lett._ , vol. 24, no. 2, Feb. 202.
  * [13] D. Hu, Y. Zhang, L. He, and J. Wu, “Low-complexity deep-learning-based DOA estimation for hybrid massive MIMO systems with uniform circular arrays,” _IEEE Wirel. Commun. Lett._ , vol. 9, no. 1, pp. 83–86, Jan. 2020.
  * [14] S. Li, Y. Liu, L. You, W. Wang, H. Duan, and X. Li, “Covariance matrix reconstruction for DOA estimation in hybrid massive MIMO systems,” _IEEE Wirel. Commun. Lett._ , vol. Early Acces.
  * [15] D. G. Manolakis, _Statistical and Adaptie Signal Processing_. ARTech House, 2005.
  * [16] X. Zhang, _Matrix Analysis and Applications (1st Edition)_. Cambridge University Press, 2017.
  * [17] B. D. Rao and K. V. S. Hari, “Performance analysis of root-MUSIC,” _IEEE Trans. Signal Process._ , vol. 37, no. 12, pp. 1939–1949, 1987.
  * [18] S. F. Chuang, W. R. Wu, and Y. T. Liu, “High-resolution AoA estimation for hybrid antenna arrays,” _IEEE Trans. Antennas Propag._ , vol. 63, no. 7, pp. 2955–2968, July 2015.