Beamspace Channel Estimation for Millimeter Wave Massive MIMO System with Hybrid Precoding and Combining

TL;DR

This paper proposes a beamspace channel estimation framework for mmWave massive MIMO systems with hybrid precoding, introducing three new schemes that outperform existing methods in efficiency and accuracy.

Contribution

The paper introduces three novel beamspace channel estimation schemes and a comprehensive framework for hybrid precoding and combining in mmWave massive MIMO systems.

Findings

Proposed schemes outperform existing methods in estimation accuracy.

Total training time slots are significantly reduced.

Simulation results approach ideal performance levels.

Abstract

In this paper, a framework of beamspace channel estimation in millimeter wave (mmWave) massive MIMO system is proposed. The framework includes the design of hybrid precoding and combining matrix as well as the search method for the largest entry of over-sampled beamspace receiving matrix. Then based on the framework, three channel estimation schemes including identity matrix approximation (IA)-based scheme, scattered zero off-diagonal (SZO)-based scheme and concentrated zero off-diagonal (CZO)-based scheme are proposed. These schemes together with the existing channel estimation schemes are compared in terms of computational complexity, estimation error and total time slots for channel training. Simulation results show that the proposed schemes outperform the existing schemes and can approach the performance of the ideal case. In particular, total time slots for channel training can be substantially reduced.

Tables

Table 1. TABLE I: Comparisons of different schemes.

	Computational Complexity	Estimation Error without Noise	Estimation Error with Noise	Total Time Slots
IA	low	$>2/K$	small	flexible
SZO	low	$2/K$	small	fixed
CZO	low	$2/K$	large	flexible
DCS	high	$>2/K$	large	flexible
OCS	high	$>2/K$	large	flexible
ECS	high	$>2/K$	large	flexible
HMC	low	depend on $T$	small	flexible
JOINT	low	depend on $T$	small	flexible

Equations

$$\mathbf{y}={\mathbf{W}}_B{\mathbf{W}}_R\sum _{u=1}^U{\mathbf{H}}_u{\mathbf{F}}_{R,u}{\mathbf{F}}_{B,u}{\mathbf{s}}_u+{\mathbf{W}}_B{\mathbf{W}}_R\mathbf{n}$$

$${\mathbf{H}}_u=\sqrt{\frac{N_A{M}_A}{L_u}}\sum _{i=1}^{L_u}{g}_{u,i}\mathbf{\alpha }(N_A,{\theta }_{u,i}){\mathbf{\alpha }}^H(M_A,{\phi }_{u,i})$$

$$\mathbf{\alpha }(N,\theta )=\frac{1}{\sqrt{N}}{\left[1,e^{-j\pi \theta },...,e^{-j\pi \theta (N-1)}\right]}^T.$$

$${\mathbf{Y}}^{t_1,t_2}={\mathbf{W}}_B^{t_2}{\mathbf{W}}_R^{t_2}\sum _{u=1}^U{\mathbf{H}}_u{\mathbf{F}}_{R,u}^{t_1}{\mathbf{F}}_{B,u}^{t_1}{\mathbf{P}}_u+{\mathbf{W}}_B^{t_2}{\mathbf{W}}_R^{t_2}{\mathbf{N}}^{t_1,t_2}$$

$${\mathbf{r}}_u^{t_1,t_2}={\mathbf{Y}}^{t_1,t_2}{\mathbf{p}}_u={\mathbf{W}}^{t_2}{\mathbf{H}}_u{\mathbf{f}}_u^{t_1}+{\tilde{\mathbf{n}}}^{t_1,t_2}$$

$$\begin{array}{rl}{\mathbf{W}}^{t_2}\triangleq {\mathbf{W}}_B^{t_2}{\mathbf{W}}_R^{t_2},& {\mathbf{f}}_u^{t_1}\triangleq {\mathbf{F}}_{R,u}^{t_1}{\mathbf{F}}_{B,u}^{t_1}{\text{1}}_{M_R},\\ {\tilde{\mathbf{n}}}^{t_1,t_2}& \triangleq {\tilde{\mathbf{N}}}^{t_1,t_2}{\mathbf{p}}_u.\end{array}$$

$${\mathbf{r}}_u^{t_1}=\mathbf{W}{\mathbf{H}}_u{\mathbf{f}}_u^{t_1}+{\tilde{\mathbf{n}}}^{t_1}$$

$$\begin{array}{rl}{\mathbf{r}}_u^{t_1}& \triangleq [({\mathbf{r}}_u^{t_1,1}{)}^T,({\mathbf{r}}_u^{t_1,2}{)}^T,\dots ,({\mathbf{r}}_u^{t_1,T_2}{)}^T{]}^T\in {\mathbb{C}}^{T_3},\\ \mathbf{W}& \triangleq [({\mathbf{W}}^1{)}^T,({\mathbf{W}}^2{)}^T,\dots ,({\mathbf{W}}^{T_2}{)}^T{]}^T\in {\mathbb{C}}^{T_3\times N_A},\\ {\tilde{\mathbf{n}}}^{t_1}& \triangleq [({\tilde{\mathbf{n}}}^{t_1,1}{)}^T,({\tilde{\mathbf{n}}}^{t_1,2}{)}^T,\dots ,({\tilde{\mathbf{n}}}^{t_1,T_2}{)}^T{]}^T\in {\mathbb{C}}^{T_3}.\end{array}$$

$$\begin{array}{rl}{\mathbf{R}}_u& \triangleq [{\mathbf{r}}_u^1,{\mathbf{r}}_u^2,\dots ,{\mathbf{r}}_u^{T_1}]\in {\mathbb{C}}^{T_3\times T_1},\\ {\mathbf{F}}_u& \triangleq [{\mathbf{f}}_u^1,{\mathbf{f}}_u^2,\dots ,{\mathbf{f}}_u^{T_1}]\in {\mathbb{C}}^{M_A\times T_1},\\ \tilde{\mathbf{n}}& \triangleq [{\tilde{\mathbf{n}}}^1,{\tilde{\mathbf{n}}}^2,\dots ,{\tilde{\mathbf{n}}}^{T_1}]\in {\mathbb{C}}^{T_3\times T_1}.\end{array}$$

$${\mathbf{R}}_u=\mathbf{W}{\mathbf{H}}_u{\mathbf{F}}_u+\tilde{\mathbf{n}}.$$

$$\begin{split}&\|\boldsymbol{D}(N_{A},K)^{H}\boldsymbol{W}^{H}\boldsymbol{W}-\gamma_{N}\boldsymbol{D}(N_{A},K)^{H}\|_{F}^{2}=\|\boldsymbol{D}(N_{A},K)^{H}(\boldsymbol{W}^{H}\boldsymbol{W}-\gamma_{N}\boldsymbol{I}_{N_{A}})\|_{F}^{2}\\ =&\textrm{Tr}\big{(}(\boldsymbol{W}^{H}\boldsymbol{W}-\gamma_{N}\boldsymbol{I}_{N_{A}})^{H}\boldsymbol{D}(N_{A},K)\boldsymbol{D}(N_{A},K)^{H}(\boldsymbol{W}^{H}\boldsymbol{W}-\gamma_{N}\boldsymbol{I}_{N_{A}})\big{)}\\ =&K\textrm{Tr}\big{(}(\boldsymbol{W}^{H}\boldsymbol{W}-\gamma_{N}\boldsymbol{I}_{N_{A}})^{H}(\boldsymbol{W}^{H}\boldsymbol{W}-\gamma_{N}\boldsymbol{I}_{N_{A}})\big{)}/N_{A}=K\|\boldsymbol{W}^{H}\boldsymbol{W}-\gamma_{N}\boldsymbol{I}_{N_{A}}\|_{F}^{2}/N_{A}\end{split}$$

$${\bar{\mathbf{H}}}_u^v=\mathbf{D}(N_A,N_A{)}^H{\mathbf{H}}_u\mathbf{D}(M_A,M_A)$$

$$\mathbf{D}($$

$$\mathbf{\alpha }(N,-1+4/K),\dots ,\mathbf{\alpha }(N,-1+2(K-1)/K)].$$

$${\mathbf{H}}_u^v=\mathbf{D}(N_A,K{)}^H{\mathbf{H}}_u\mathbf{D}(M_A,K)$$

$${\mathbf{R}}_u^v=\mathbf{D}(N_A,K{)}^H{\mathbf{W}}^H\mathbf{W}{\mathbf{H}}_u{\mathbf{F}}_u{\mathbf{F}}_u^H\mathbf{D}(M_A,K)+{\tilde{\mathbf{n}}}^v$$

$$\begin{array}{rl}\mathbf{D}(N_A,K{)}^H{\mathbf{W}}^H\mathbf{W}& ={\gamma }_N\mathbf{D}(N_A,K{)}^H,\\ {\mathbf{F}}_u{\mathbf{F}}_u^H\mathbf{D}(M_A,K)& ={\gamma }_M\mathbf{D}(M_A,K).\end{array}$$

$$\begin{array}{rl}\underset{\mathbf{W}}{\min }& \Vert \mathbf{D}(N_A,K{)}^H{\mathbf{W}}^H\mathbf{W}-{\gamma }_N\mathbf{D}(N_A,K{)}^H{\Vert }_F\\ s.t.& {\mathbf{W}}_R^{t_2}\in {\mathcal{W}}_{\mathcal{R}},t_2=1,2,\dots ,T_2,\\ & \Vert {\mathbf{W}}_B^{t_2}{\mathbf{W}}_R^{t_2}{\Vert }_F^2=P_W,t_2=1,2,\dots ,T_2,\end{array}$$

$$\begin{array}{rl}\underset{\mathbf{W}}{\min }& \Vert {\mathbf{W}}^H\mathbf{W}-{\gamma }_N{\mathbf{I}}_{N_A}{\Vert }_F\\ s.t.& {\mathbf{W}}_R^{t_2}\in {\mathcal{W}}_{\mathcal{R}},t_2=1,2,\dots ,T_2,\\ & \Vert {\mathbf{W}}_B^{t_2}{\mathbf{W}}_R^{t_2}{\Vert }_F^2=P_W,t_2=1,2,\dots ,T_2.\end{array}$$

$$\begin{array}{rl}\underset{{\mathbf{W}}_B^{t_2},{\mathbf{W}}_R^{t_2}}{\min }& \Vert {\mathbf{W}}_B^{t_2}{\mathbf{W}}_R^{t_2}-{\tilde{\mathbf{W}}}^{t_2}{\Vert }_F\\ s.t.& {\mathbf{W}}_R^{t_2}\in {\mathcal{W}}_{\mathcal{R}},\Vert {\mathbf{W}}_B^{t_2}{\mathbf{W}}_R^{t_2}{\Vert }_F^2=P_W.\end{array}$$

$${\mathbf{W}}_B^{t_2}\leftarrow \sqrt{P_W}\frac{{\mathbf{W}}_B^{t_2}}{\Vert {\mathbf{W}}_B^{t_2}{\mathbf{W}}_R^{t_2}{\Vert }_F},t_2=1,2,\dots ,T_2$$

$$\Vert {\mathbf{W}}_B^{t_2}{\mathbf{W}}_R^{t_2}-{\tilde{\mathbf{W}}}^{t_2}{\Vert }_F^2$$

$$=\Vert {\mathbf{W}}_B^{t_2}({\mathbf{W}}_R^{t_2}-{\beta }^{-1}{{\mathbf{W}}_B^{t_2}}^H{\tilde{\mathbf{W}}}^{t_2}){\Vert }_F^2$$

$$\displaystyle=\textrm{Tr}\big{(}(\boldsymbol{W}_{R}^{t_{2}}-{\boldsymbol{W}_{D}^{t_{2}}}^{H}\widetilde{\boldsymbol{W}}^{t_{2}})^{H}{\boldsymbol{W}_{B}^{t_{2}}}^{H}\boldsymbol{W}_{B}^{t_{2}}(\boldsymbol{W}_{R}^{t_{2}}-{\boldsymbol{W}_{D}^{t_{2}}}^{H}\widetilde{\boldsymbol{W}}^{t_{2}})\big{)}$$

$$\displaystyle=\textrm{Tr}\big{(}(\boldsymbol{W}_{R}^{t_{2}}-{\boldsymbol{W}_{D}^{t_{2}}}^{H}\widetilde{\boldsymbol{W}}^{t_{2}})^{H}\beta(\boldsymbol{W}_{R}^{t_{2}}-{\boldsymbol{W}_{D}^{t_{2}}}^{H}\widetilde{\boldsymbol{W}}^{t_{2}})\big{)}$$

$$\displaystyle=\beta\textrm{Tr}\big{(}(\boldsymbol{W}_{R}^{t_{2}}-{\boldsymbol{W}_{D}^{t_{2}}}^{H}\widetilde{\boldsymbol{W}}^{t_{2}})^{H}(\boldsymbol{W}_{R}^{t_{2}}-{\boldsymbol{W}_{D}^{t_{2}}}^{H}\widetilde{\boldsymbol{W}}^{t_{2}})\big{)}$$

$$=\beta \Vert {\mathbf{W}}_R^{t_2}-{{\mathbf{W}}_D^{t_2}}^H{\tilde{\mathbf{W}}}^{t_2}{\Vert }_F^2.$$

$$\begin{array}{rl}\underset{{\mathbf{W}}_B^{t_2},{\mathbf{W}}_R^{t_2}}{\min }& \Vert {\mathbf{W}}_R^{t_2}-{{\mathbf{W}}_D^{t_2}}^H{\tilde{\mathbf{W}}}^{t_2}{\Vert }_F^2\\ s.t.& {\mathbf{W}}_R^{t_2}\in {\mathcal{W}}_{\mathcal{R}}.\end{array}$$

$${\mathbf{W}}_R^{t_2}=\text{arg}({{\mathbf{W}}_D^{t_2}}^H{\tilde{\mathbf{W}}}^{t_2},{\mathcal{W}}_{\mathcal{R}})$$

$$\begin{array}{rl}\underset{{\mathbf{W}}_D^{t_2}}{\min }& \Vert {\mathbf{W}}_R^{t_2}-{{\mathbf{W}}_D^{t_2}}^H{\tilde{\mathbf{W}}}^{t_2}{\Vert }_F^2\\ s.t.& {{\mathbf{W}}_D^{t_2}}^H{\mathbf{W}}_D^{t_2}={\beta }^{-1}{\mathbf{I}}_{N_R}.\end{array}$$

$${\mathbf{W}}_D^{t_2}={\beta }^{-1/2}{\mathbf{VU}}^H,$$

$$ϵ\triangleq \frac{\Vert {\mathbf{W}}_R^{t_2,i}-{\mathbf{W}}_R^{t_2,i-1}{\Vert }_F^2+\Vert {\mathbf{W}}_B^{t_2,i}-{\mathbf{W}}_B^{t_2,i-1}{\Vert }_F^2}{\Vert {\mathbf{W}}_R^{t_2,i-1}{\Vert }_F^2+\Vert {\mathbf{W}}_B^{t_2,i-1}{\Vert }_F^2},$$

$$\begin{array}{rl}\underset{{\mathbf{F}}_u}{\min }& \Vert {\mathbf{F}}_u{\mathbf{F}}_u^H\mathbf{D}(M_A,K)-{\gamma }_M\mathbf{D}(M_A,K){\Vert }_F\\ s.t.& {\mathbf{F}}_{R,u}^{t_1}\in {\mathcal{F}}_{\mathcal{R}},t_1=1,2,\dots ,T_1,\\ & \Vert {\mathbf{F}}_{R,u}^{t_1}{\mathbf{F}}_{B,u}^{t_1}{\Vert }_F^2=P_F,t_1=1,2,\dots ,T_1,\end{array}$$

Full text

I Introduction

Millimeter wave (mmWave) communication is a promising technology for next generation wireless communication owing to its abundant frequency spectrum resource[1, 2]. However, realizing mmWave massive MIMO in practice is not a trivial task, which faces the problem of high propagation loss due to the high carrier frequency [3]. To compensate for the propagation loss, antenna arrays are usually used to form directional beamforming. Fortunately, thanks to the short wavelength of the mmWave frequency, large antenna arrays are possible to be packed into small form factors.

In order to exploit the spatial degree of freedom, the hybrid analog and digital precoding is usually employed [4, 5, 6]. A small number of RF chains are tied to a large antenna array. This structure enables parallel transmission, and thus provides the potential to approach the capacity bound that can be achieved by digital precoding. On the other hand, the large antenna arrays challenge the low-complexity design of hybrid precoding and channel estimation [7]. In particular, the hybrid precoding may require matrix operations with a scale of antenna size, which is generally large in mmWave communication [4]. Moreover, the channel estimation is also rather time consuming due to the large number of antennas at both transmitting and receiving sides [8].

To reduce the complexity of channel estimation in mmWave massive MIMO system, some advanced schemes based on the beamspace channel have been proposed very recently [8, 9, 10, 11, 12, 13, 14]. A hybrid and multiresolutional codebook (HMC)-based channel estimation method and a JOINT-based channel estimation method are proposed in [8] and [9], respectively. The common idea of [8] and [9] is to use the hierarchical codebook, where the precision of channel estimation relies on the number of the layers in the hierarchical codebook. The key ideas of [11, 12, 13, 14] are to efficiently explore the sparsity of beamspace channel by sparse signal processing techniques. An enhanced compressive sensing (ECS)-based channel estimation scheme is proposed in [11], which explores the channel sparsity in beamspace. With higher resolution of phase shifters, improved channel estimation accuracy can be achieved since higher freedom is available for the design of measurement matrix for the sparse recovery[4, 5].

However, considering the limited beamspace resolution, the sparsity of beamspace channel may be impaired by power leakage [15], indicating that the beamspace channel is not ideally sparse and there are many small nonzero entries [16]. Therefore, it brings extra challenge for the sparse recovery [17, 18]. To solve this problem, an adaptive support detection (ASD)-based channel estimation scheme is proposed to iteratively detect and adjust the channel support to find one with the largest channel power [14]. Note that the beamspace resolution is proportional to the reciprocal of the antenna number. Therefore, ASD-based scheme cannot obtain high precision channel estimation regarding the channel power leakage [8]. A discrete compressive sensing (DCS)-based channel estimation scheme is proposed in [15], where the beamspace resolution can be set arbitrarily and can be higher than the reciprocal of the antenna number, leading to better channel estimation performance than ASD-based scheme. An over-sampling compressive sensing (OCS)-based channel estimation scheme is proposed in [17], where the measurement matrix is consisted of a large number of over-sampling steering vectors and is capable of estimating the angle of arrival (AoA) and angle of departure (AoD) more accurate than conventional measurement matrix. However, both DCS-based and OCS-based schemes use random measurement matrices, which can not always achieve the optimal performance. In particular, the estimated AoA and AoD can not be always within the scale of quantization error even without any noise. Therefore, it is better to use deterministic measurement matrix.

In this paper, we first propose a framework of beamspace channel estimation, which is divided into two subproblems, including the design of the hybrid precoding and combining matrix, and the search method for the largest entry of over-sampled beamspace receiving matrix. Note that the design of the hybrid precoding and combining matrix as well as the search method for the largest entry is original. Then based on the framework, we propose three channel estimation schemes.

1. We propose an identity matrix approximation (IA)-based channel estimation scheme. We formulate the design of hybrid combining and hybrid precoding as two optimization problems with the constraint of total power and constraint of the constant envelope required by phase shifters. Due to the non-convexity of the problems, we decouple the design of analog combining and digital combining, and then obtain closed-form solutions. We propose an algorithm for the design of hybrid combining matrix. The algorithm repeatedly fixes the analog combining matrix to obtain digital combining matrix, and then fixes the digital combining matrix to obtain the analog combining matrix in turn, until the stop condition is satisfied. Detailed steps summarized in an algorithm table are provided. Since the design of the hybrid precoding matrix is similar, we briefly describe the design of the hybrid precoding matrix. Note that the hybrid combining matrix and the hybrid precoding matrix can be designed off-line before the channel training. After that, we propose an algorithm to search the largest entry of the over-sampled beamspace receiving matrix. The algorithm is based on trichotomy search and includes two stages, where we find the main lobe in the first stage, and we find the largest entry corresponding to the peak within the main lobe in the second stage.

2. We design hybrid precoding matrix and hybrid combining matrix so that the coordinates of the largest entry of over-sampled beamspace receiving matrix are the AoA and AoD within the quantization error. Then we convert the discrete problem into a continuous problem to obtain the derivative. By setting the derivative zero, we get two solutions where the first solution is found to be meaningless. Based on the second solution, we propose two zero off-diagonal (ZO) beamspace channel estimation schemes, namely, scattered zero off-diagonal (SZO)-based scheme and concentrated zero off-diagonal (CZO)-based scheme.

• •

In the SZO-based channel estimation scheme, the nonzero diagonal entries are uniformly distributed with the same interval. Since the main lobe and the side lobes have the same envelope, we first make beam training based on codebook to find the main lobe and then use complementary channel estimation to further estimate the channel AoA and AoD within the main lobe. An integration-based codebook design method which results in closed-form expression of codewords is proposed and compared with the existing sparse-based codebook design method. Note that the proposed integration-based codebook design method is original.

• •

In the CZO-based channel estimation scheme, the nonzero diagonal entries are concentrated on the upper left corner of the matrix. Since the envelope of the main lobe and the side lobe is different, we can directly employ the algorithm proposed in the IA-based scheme to search the largest entry of over-sampled beamspace receiving matrix corresponding to the channel AoA and AoD within the main lobe.

Additionally, we also compare the above three schemes together with the existing channel estimation schemes in terms of computational complexity, estimation error and total time slots for channel training.

The rest of the paper is organized as follows. In Section II, the system model and problem formulation of beamspace channel estimation with hybrid precoding and combining are provided. In Section III, we propose three beamspace channel estimation schemes. The simulation results are provided in Section IV. Finally, Section V concludes the paper.

The notations are defined as follows. Symbols for matrices (upper case) and vectors (lower case) are in boldface. $(\cdot)^{T}$, $(\cdot)^{H}$, $(\cdot)^{-1}$, $\boldsymbol{I}_{L}$, $\textbf{1}_{L}$, $\mathbb{C}^{M\times{N}}$, $\otimes$, $\textrm{vec}(\cdot)$, $\operatorname{E}\{\cdot\}$, $\mathcal{O}(\cdot)$, $\boldsymbol{0}^{M}$, $\boldsymbol{0}_{M\times N}$, $\|\cdot\|_{0}$, $\|\cdot\|_{2}$, $\|\cdot\|_{F}$, $\boldsymbol{A}[p,q]$, $\langle\cdot\rangle$, $\lfloor\cdot\rfloor$, $\textrm{Tr}(\cdot)$, $\mathbb{Z}$, $\mathbb{Re}\{\cdot\}$ and $\mathcal{CN}$, denote the transpose, conjugate transpose (Hermitian), inverse, identity matrix of size $L$, vector of size $L$ with all entries being 1, the set of $M\times{N}$ complex-valued matrices, kronecker product, vectorization, expectation, order of complexity, zero vector of size $M$, $M\times{N}$ zero matrix, $l_{0}$-norm, $l_{2}$-norm, Frobenius norm, entry of $\boldsymbol{A}$ at the $p$th row and $q$th column, round function, floor function, trace, set of integer, real part and complex Gaussian distribution, respectively.

II System Model and Problem Formulation

II-A System Model

We consider a time division duplexing (TDD) multi-user mmWave massive MIMO system comprising a base station (BS) and $U$ users. We focus on the uplink transmission. Both the BS and users are equipped with an uniform linear array (ULA) [19]. Let $N_{A}$, $M_{A}$, $N_{R}$ and $M_{R}$ denote the number of antennas at the BS, number of antennas at each user, number of RF chains at the BS and number of RF chains at each user. In practical mmWave massive MIMO with hybrid precoding and combining, the number of RF chains is much smaller than that of antennas, i.e., $N_{R}\ll N_{A}$ and $M_{R}\ll M_{A}$. This is because we use large antenna array to form directional beamforming, which can compensate for the propagation loss caused by the high carrier frequency [20].

For uplink transmission, each user performs analog precoding in RF and digital precoding in the baseband, while the BS performs analog combining in RF and digital combining in the baseband. The received signal vector at the BS can be represented as

[TABLE]

where $\boldsymbol{F}_{B,u}\in{\mathbb{C}^{M_{R}\times{M_{R}}}}$, $\boldsymbol{F}_{R,u}\in{\mathbb{C}^{M_{A}\times{M_{R}}}}$, $\boldsymbol{W}_{B}\in{\mathbb{C}^{N_{R}\times{N_{R}}}}$, and $\boldsymbol{W}_{R}\in{\mathbb{C}^{N_{R}\times{N_{A}}}}$ are the digital precoding matrix, analog precoding matrix, digital combining matrix, and analog combining matrix for the $u(u=1,2,\ldots,U)$th user, respectively. $\boldsymbol{s}_{u}\in{\mathbb{C}^{M_{R}}}$ denotes the signal vector satisfying $\operatorname{E}\{\boldsymbol{s}_{u}\boldsymbol{s}_{u}^{H}\}=\boldsymbol{I}_{M_{R}}$. $\boldsymbol{n}\in{\mathbb{C}^{N_{A}}}$ denotes additive white Gaussian noise (AWGN) vector satisfying $\boldsymbol{n}\sim\mathcal{CN}(0,\sigma^{2}\boldsymbol{I}_{N_{A}})$. $\boldsymbol{H}_{u}\in{\mathbb{C}^{N_{A}\times{M_{A}}}}$ denotes the channel matrix between the BS and the $u$th user and can be expressed according to the widely used Saleh-Valenzuela channel model [1] as

[TABLE]

where $L_{u}$ and $g_{u,i}$ denote the total number of resolvable paths and the channel fading coefficient of the $i$th path for the $u$th user, respectively. Usually, there is a strong line-of-sight (LOS) path ($i=1$) and $L_{u}-1$ much weaker non-line-of-sight (NLOS) paths ($2\leq i\leq L_{u}$). The mmWave transmission essentially relies on the strong LOS path. The steering vector $\boldsymbol{\alpha}(N,\theta)$ is defined as

[TABLE]

Define the AoA and AoD of the $i$th path of the $u$th user as $\Theta_{u,i}$ and $\Phi_{u,i}$, respectively. Further define $\theta_{u,i}\triangleq\frac{2d_{BS}}{\lambda}\sin{\Theta_{u,i}}$ and $\varphi_{u,i}\triangleq\frac{2d_{UE}}{\lambda}\sin{\Phi_{u,i}}$, where $d_{BS}$ and $d_{UE}$ denote the antenna interval of the BS and users, respectively. We usually set $d_{BS}=d_{UE}=\lambda/2$, where $\lambda$ is the wavelength of mmWave signal. In practice, both $\Theta_{u,i}$ and $\Phi_{u,i}$ obey the uniform distribution $[-\pi,\pi]$ [8, 21].

II-B Problem Formulation

Note that $\boldsymbol{y}$ in (1) is a combination of signal from different users. We use $T_{1}$ different digital precoding matrices and analog precoding matrices, denoted as $\boldsymbol{F}_{B,u}^{t_{1}}\in{\mathbb{C}^{M_{R}\times{M_{R}}}}$ and $\boldsymbol{F}_{R,u}^{t_{1}}\in{\mathbb{C}^{M_{A}\times{M_{R}}}}$, respectively, $t_{1}=1,2,\ldots,T_{1}$, at the $u(u=1,2,\ldots,U)$th user. We use $T_{2}$ different digital combining matrices and analog combining matrices, denoted as $\boldsymbol{W}_{B}^{t_{2}}\in{\mathbb{C}^{N_{R}\times{N_{R}}}}$ and $\boldsymbol{W}_{R}^{t_{2}}\in{\mathbb{C}^{N_{R}\times{N_{A}}}}$, respectively, $t_{2}=1,2,\ldots,T_{2}$, at the BS. To distinguish different user signal at the BS, each user repeatedly transmits an orthogonal pilot sequence $\boldsymbol{p}_{u}\in{\mathbb{C}^{U}}$ for $T_{1}T_{2}$ times. For simplicity, we suppose each user transmit the same pilot sequence for all $M_{R}$ RF chains, where the pilot matrix for the $u$th user can be defined as $\boldsymbol{P}_{u}\triangleq[\boldsymbol{p}_{u},\boldsymbol{p}_{u},\ldots,\boldsymbol{p}_{u}]^{H}=\textbf{1}_{M_{R}}\boldsymbol{p}_{u}^{H}\in{\mathbb{C}^{M_{R}\times{U}}}$. The channel keeps constant during $T\triangleq T_{1}T_{2}U$ time slots [16]. During the $T_{1}$ repetitive transmission of pilot sequence from the $((t_{2}-1)T_{1}+1)$th transmission to $(t_{2}T_{1})$th transmission, we use $T_{1}$ different $\boldsymbol{F}_{B,u}^{t_{1}}$ and $\boldsymbol{F}_{R,u}^{t_{1}}$ for hybrid precoding while using the same $\boldsymbol{W}_{B}^{t_{2}}$ and $\boldsymbol{W}_{R}^{t_{2}}$ for hybrid combining, where the received pilot matrix $\boldsymbol{Y}^{t_{1},t_{2}}\in{\mathbb{C}^{N_{R}\times{U}}}$ can be denoted as

[TABLE]

with $\boldsymbol{N}^{t_{1},t_{2}}\in{\mathbb{C}^{N_{A}\times{U}}}$ representing the AWGN matrix. Each entry of $\boldsymbol{N}^{t_{1},t_{2}}$ independently obeys complex Gaussian distribution with zero mean and variance of $\sigma^{2}$. To ease the notation, we define $\widetilde{\boldsymbol{N}}^{t_{1},t_{2}}\triangleq\boldsymbol{W}_{B}^{t_{1}}\boldsymbol{W}_{R}^{t_{1}}\boldsymbol{N}^{t_{1},t_{2}}$. Due to the orthogonality of $\boldsymbol{p}_{u}$, i.e., \boldsymbol{p}_{u}^{H}$$\boldsymbol{p}_{u}=1 and \boldsymbol{p}_{u}^{H}$$\boldsymbol{p}_{i}=0, $\forall u,i\in\{1,2,\ldots,U\}$, $i\neq{u}$ [14], we can obtain the measurement vector $\boldsymbol{r}_{u}^{t_{1},t_{2}}\in{\mathbb{C}^{N_{R}}}$ for the $u$th user by multiplying $\boldsymbol{Y}^{t_{1},t_{2}}$ with $\boldsymbol{p}_{u}$ as

[TABLE]

where

[TABLE]

Note that we can distinguish different user signal utilizing the orthogonality of pilot sequences. In this way, the spatial distribution of users has no impact on the performance of our proposed schemes. Define $T_{3}\triangleq T_{2}N_{R}$. We stack the $T_{2}$ received pilot sequences together and have

[TABLE]

where

[TABLE]

Further define

[TABLE]

We have

[TABLE]

Considering the LOS channel path concentrates most channel power in mmWave massive MIMO systems, we usually use LOS channel path to transmit data for the $u$th user [15]. Therefore, it is important to design $\boldsymbol{W}$ and $\boldsymbol{F}_{u}$ to estimate the AoA and the AoD of the LOS path of $\boldsymbol{H}_{u}$ in (10) for the $u$th user, which will be discussed in the following sections.

III Beamspace Channel Estimation

In this section, we first propose a framework of beamspace channel estimation. Then based on this framework, three channel estimation schemes are proposed. Finally, the comparisons of these three schemes together with the existing HMC-based [8], JOINT-based [9], ECS-based [11], DCS-based [15] and OCS-based [17] channel estimation schemes are also presented.

III-A Framework of Beamspace Channel Estimation

The beamspace channel matrix for the $u$th user $\bar{\boldsymbol{H}}_{u}^{v}\in{\mathbb{C}^{N_{A}\times{M_{A}}}}$ can be represented as [14]

[TABLE]

where $\boldsymbol{D}(N,K)\in{\mathbb{C}^{N\times{K}}}$ is the sampling matrix, which is defined as

[TABLE]

In fact, $\boldsymbol{D}(N,K)$ samples the beamspace [-1,1] in an interval of $2/K$ by $K$ steering vectors. For $\bar{\boldsymbol{H}}_{u}^{v}$, the AoA and AoD is sampled in an interval of $2/N_{A}$ and $2/M_{A}$, respectively. Therefore the quantization error for the estimated AoA and AoD is $2/N_{A}$ and $2/M_{A}$, respectively. In order to decrease the quantization error, we introduce the over-sampled beamspace channel matrix for the $u$th user $\boldsymbol{H}_{u}^{v}\in{\mathbb{C}^{K\times{K}}}$ as

[TABLE]

where $K$ is the number of steering vectors with $K>N_{A}$ and $K>M_{A}$. Then the coordinates of the largest entry of $\boldsymbol{H}_{u}^{v}$ are the AoA and AoD of the LOS path with the quantization error of $2/K$. To reduce the quantization error and improve channel estimation, we can use a large $K$ by finding the largest entry of $\boldsymbol{H}_{u}^{v}$.

However, we cannot directly obtain $\boldsymbol{H}_{u}^{v}$ based on $\boldsymbol{R}_{u}$ in (10), due to the hybrid precoding and combining operations of $\boldsymbol{F}_{u}$ and $\boldsymbol{W}$, respectively. Note that the dimension of $\boldsymbol{H}_{u}$ is $N_{A}\times M_{A}$, while the dimension of $\boldsymbol{W}\boldsymbol{H}_{u}\boldsymbol{F}_{u}$ is $T_{3}\times T_{1}$. In order to obtain the over-sampled beamspace channel matrix as described in (14), we multiply $\boldsymbol{R}_{u}$ with $\boldsymbol{W}^{H}$ on the left and $\boldsymbol{F}_{u}^{H}$ on the right, which can make the dimension of $\boldsymbol{W}^{H}\boldsymbol{W}\boldsymbol{H}_{u}\boldsymbol{F}_{u}\boldsymbol{F}_{u}^{H}$ the same as that of $\boldsymbol{H}_{u}$. Now we can obtain an over-sampled beamspace receiving matrix $\boldsymbol{R}_{u}^{v}\in{\mathbb{C}^{K\times{K}}}$ as

[TABLE]

where $\widetilde{\boldsymbol{n}}^{v}\triangleq\boldsymbol{D}(N_{A},K)^{H}\boldsymbol{W}^{H}\widetilde{\boldsymbol{n}}\boldsymbol{F}_{u}^{H}\boldsymbol{D}(M_{A},K)$. It is expected that

[TABLE]

where $\gamma_{N}\triangleq\|\boldsymbol{W}^{H}\boldsymbol{W}\|_{F}/\sqrt{N_{A}}$ and $\gamma_{M}\triangleq\|\boldsymbol{F}_{u}\boldsymbol{F}_{u}^{H}\|_{F}/\sqrt{M_{A}}$. However, in this case, it requires that $T_{3}\geq N_{A}$ and $T_{1}\geq M_{A}$, leading to huge pilot overhead. For example, in an mmWave massive MIMO system with $N_{A}=64$, $M_{A}=16$, $N_{R}=4$ and $M_{R}=1$, the pilot sequence should be repetitively transmitted for $T_{1}T_{2}=256$ times. Therefore, it is required to reduce the pilot overhead in practice, where none of $\boldsymbol{D}(N_{A},K)^{H}\boldsymbol{W}^{H}\boldsymbol{W}=\gamma_{N}\boldsymbol{D}(N_{A},K)^{H}$ and $\boldsymbol{F}_{u}\boldsymbol{F}_{u}^{H}\boldsymbol{D}(M_{A},K)=\gamma_{M}\boldsymbol{D}(M_{A},K)$ can be satisfied. Now we have the following two subproblems.

Subproblem 1 (Hybrid Precoding and Combining Matrix Design): $\boldsymbol{W}$ and $\boldsymbol{F}_{u}$ should be well designed so that the coordinates of the largest entry of $\boldsymbol{R}_{u}^{v}$ are the AoA and AoD of the LOS path with the quantization error of $2/K$. Therefore, the AoA and AoD of the LOS path can be estimated by finding the largest entry of $\boldsymbol{R}_{u}^{v}$.

Subproblem 2 (Search the Largest Entry): Due to the large dimension of $\boldsymbol{R}_{u}^{v}$, finding the largest entry of $\boldsymbol{R}_{u}^{v}$ is computationally expensive. Therefore, it is better to design a low-complexity search algorithm fully regarding the structure of $\boldsymbol{R}_{u}^{v}$.

III-B IA-based Beamspace Channel Estimation Scheme

III-B1 Hybrid Precoding and Combining Matrix Design

It is observed that $\boldsymbol{R}_{u}^{v}$ in (15) is the over-sampled beamspace channel matrix with noise if (16) is satisfied. However, due to the fact that $T_{3}<N_{A}$, $T_{1}<M_{A}$, the rank of $\boldsymbol{W}^{H}\boldsymbol{W}$ and $\boldsymbol{F}_{u}\boldsymbol{F}_{u}^{H}$ is less than $N_{A}$ and $M_{A}$, respectively. So we cannot find a proper $\boldsymbol{W}$ and $\boldsymbol{F}_{u}$ satisfying (16).

The optimization problem for hybrid combining matrix design is

[TABLE]

where $\mathcal{W_{R}}$ is the set of all feasible analog combining matrix and $P_{W}=1$ to normalize the hybrid combining matrix. Note that the design of hybrid precoder and hybrid combiner is independent. We can well design the hybrid precoding and combining matrices before the transmission of pilot sequences. Note that $\boldsymbol{D}(N_{A},K)\boldsymbol{D}(N_{A},K)^{H}=K\boldsymbol{I}_{N_{A}}/N_{A}$. We have (11). Then (17) can be further rewritten as

[TABLE]

It is seen that $\boldsymbol{W}$ defined in (8) is a flat matrix where the columns are more than the rows. Therefore, it is infeasible that $\boldsymbol{W}^{H}\boldsymbol{W}$ equals $\gamma_{N}\boldsymbol{I}_{N_{A}}$. Instead, it is important to design $\boldsymbol{W}$ so that $\boldsymbol{W}^{H}\boldsymbol{W}/\gamma_{N}$ approximates the identity matrix, a.k.a, identity matrix approximation (IA). To minimize $\|\boldsymbol{W}^{H}\boldsymbol{W}-\gamma_{N}\boldsymbol{I}_{N_{A}}\|_{F}$, $\boldsymbol{W}$ can be a submatrix of $\sqrt{\gamma_{N}}\boldsymbol{U}$ by selecting the first $T_{3}$ rows of $\sqrt{\gamma_{N}}\boldsymbol{U}$, where $\boldsymbol{U}$ is any $N_{A}\times N_{A}$ unitary matrix [22]. For example, we obtain $\boldsymbol{U}$ by singular value decomposition (SVD) of a $N_{A}\times N_{A}$ random matrix $\boldsymbol{A}$, i.e., $\boldsymbol{A}=\boldsymbol{U\Sigma V}^{H}$, where each entry of $\boldsymbol{A}$ obeys the uniform distribution [0,1]. In this way, we obtain $\widetilde{\boldsymbol{W}}$. According to (8), we can obtain $\widetilde{\boldsymbol{W}}^{t_{2}},~{}t_{2}=1,2,\ldots,T_{2}$, which is essentially dividing $\widetilde{\boldsymbol{W}}$ into $T_{2}$ submatrices. Then (18) is converted into $T_{2}$ subproblems, where each subproblem can be expressed as

[TABLE]

We cannot directly obtain solutions for (19) due to the non-convexity of the constraints. Note that [23] has proved that temporarily neglecting the second power constraint of (19) during the optimization of $\boldsymbol{W}_{B}^{t_{2}}$ and $\boldsymbol{W}_{R}^{t_{2}}$ will have little impact on the optimality of the hybrid precoding problem. After $\boldsymbol{W}_{B}^{t_{2}}$ and $\boldsymbol{W}_{R}^{t_{2}}$ are obtained, we may set $\boldsymbol{W}_{B}^{t_{2}}$ as

[TABLE]

to satisfy the second constraint of (19).

To mitigate the interference among different data streams, we impose a common constraint that the columns of the digital combining matrix are mutually orthogonal, i.e., ${\boldsymbol{W}_{B}^{t_{2}}}^{H}\boldsymbol{W}_{B}^{t_{2}}=\beta\boldsymbol{I}_{N_{R}}$, where $\beta\triangleq P_{W}/(N_{A}N_{R})$. Define $\boldsymbol{W}_{D}^{t_{2}}\triangleq\beta^{-1}\boldsymbol{W}_{B}^{t_{2}}$. Then we have

[TABLE]

Therefore (19) can be expressed as

[TABLE]

It shows that $\boldsymbol{W}_{R}^{t_{2}}$ and $\boldsymbol{W}_{D}^{t_{2}}$ are decoupled. Given $\boldsymbol{W}_{D}^{t_{2}}$, the solution of $\boldsymbol{W}_{R}^{t_{2}}$ can be expressed as

[TABLE]

where $\textrm{arg}(\boldsymbol{A},\mathcal{R})$ first normalizes each entry of $\boldsymbol{A}$ and then quantizes the normalized matrix in terms of $\mathcal{R}$. Note that the quantization is required since the resolution of phase shifters is limited in practice, e.g., the resolution is 2/64 if the phase shifter is 6 bits and the range of the angle is [-1,1].

Similarly, given $\boldsymbol{W}_{R}^{t_{2}}$, the optimization of $\boldsymbol{W}_{D}^{t_{2}}$ based on (22) can be expressed as

[TABLE]

It is seen that (24) is similar to the orthogonal Procrustes problem [24]. Then the solution to (24) can be obtained as

[TABLE]

where $\boldsymbol{W}_{R}^{t_{2}}(\widetilde{\boldsymbol{W}}^{t_{2}})^{H}=\boldsymbol{U\Sigma V}^{H}$ represents the SVD of $\boldsymbol{W}_{R}^{t_{2}}(\widetilde{\boldsymbol{W}}^{t_{2}})^{H}$. Then we obtain $\boldsymbol{W}_{B}^{t_{2}}=\beta\boldsymbol{W}_{D}^{t_{2}}$. Note that both (23) and (25) are optimal closed-form solutions to the problems expressed in (22) and (24), respectively.

As shown in $\mathbf{Algorithm~{}\ref{alg1}}$, we propose an algorithm of hybrid combining matrix design for the IA-based channel estimation. We repeatedly fix $\boldsymbol{W}_{R}^{t_{2}}$ to obtain $\boldsymbol{W}_{B}^{t_{2}}$ via (25), and then fix $\boldsymbol{W}_{R}^{t_{2}}$ to obtain $\boldsymbol{W}_{B}^{t_{2}}$ via (23) in turn. Define the normalized iteration error $\epsilon$ as

[TABLE]

the stop condition is that the iterative update of both $\boldsymbol{W}_{R}^{t_{2}}$ and $\boldsymbol{W}_{B}^{t_{2}}$ is stable, i.e., $\epsilon<\delta$, where $\delta$ is the threshold.

The optimization problem for hybrid precoding matrix design is

[TABLE]

where $\mathcal{F_{R}}$ is the set of all feasible analog precoding matrix and $P_{F}$ is the given power for the hybrid precoding. Similar to (18), (27) can be rewritten as

[TABLE]

We can obtain a solution as $\widetilde{\boldsymbol{F}}_{u}$ by selecting the first $T_{1}$ columns of $\sqrt{\gamma_{M}}\boldsymbol{U}$, where $\boldsymbol{U}$ is any $M_{A}\times M_{A}$ unitary matrix. According to (9), we can obtain $\widetilde{\boldsymbol{f}}^{t_{1}}_{u}$ as the $t_{1}$-th column of $\widetilde{\boldsymbol{F}}_{u}$, $t_{1}=1,2,\ldots,T_{1}$. Define

[TABLE]

Then we have

[TABLE]

Similar to (19), (28) can be converted to $T_{1}$ subproblems, where each subproblem is expressed as

[TABLE]

Similar to (19), the second constraint of (31) can also be temporarily neglected. Therefore, we may replace $\boldsymbol{W}_{B}^{t_{2}}$, $\boldsymbol{W}_{R}^{t_{2}}$, $\widetilde{\boldsymbol{W}}^{t_{2}}$, $\mathcal{W_{R}}$ and $P_{W}$ in (19) with $(\boldsymbol{f}_{B,u}^{t_{1}})^{H}$, $(\boldsymbol{F}_{R,u}^{t_{1}})^{H}$, $(\widetilde{\boldsymbol{f}}^{t_{1}}_{u})^{H}$, $\mathcal{F_{R}}$ and $P_{F}$, respectively.

In order to run $\mathbf{Algorithm~{}\ref{alg1}}$ to obtain $\boldsymbol{F}_{u}$, we have to further replace $N_{A}$, $T_{2}$, $\widetilde{\boldsymbol{W}}$ and $\boldsymbol{W}$ with $M_{A}$, $T_{1}$, $\widetilde{\boldsymbol{F}}_{u}$ and $\boldsymbol{F}_{u}$. The routine of the algorithm is exactly the same, except that an additional operation to obtain $\boldsymbol{F}_{B,u}^{t_{1}}$ as

[TABLE]

is required after finishing step 9; and $\boldsymbol{W}_{B}^{t_{2}}$ should be replaced by $(\boldsymbol{F}_{B,u}^{t_{1}})^{H}$ at step 10.

In summary, in the first half of the IA-based channel estimation, we design of hybrid combining matrix and hybrid precoding matrix; while in the other half to be discussed, we will search the largest entry of the over-sampled beamspace channel matrix.

III-B2 Search the Largest Entry

Instead of exhaustively search the largest entry from $\boldsymbol{R}_{u}^{v}$ defined in (15), we can improve the efficiency of the search algorithm by analyzing the structure of $\boldsymbol{R}_{u}^{v}$. Neglecting the term from the additive noise and assuming there is single path, we rewrite $\boldsymbol{R}_{u}^{v}$ as

[TABLE]

where

[TABLE]

[TABLE]

Since $\boldsymbol{r}_{N}$ is a column vector and $\boldsymbol{r}_{M}^{H}$ is a row vector, the largest entry of $\boldsymbol{R}_{u}^{v}$ essentially depends on the largest entry of $\boldsymbol{r}_{N}$ and $\boldsymbol{r}_{M}$. Therefore we will analyze the structure of $\boldsymbol{r}_{N}$ and $\boldsymbol{r}_{M}$.

In the ideal case, i.e., $\boldsymbol{F}_{u}\boldsymbol{F}_{u}^{H}\boldsymbol{D}(M_{A},K)=\gamma_{M}\boldsymbol{D}(M_{A},K)$, $\boldsymbol{r}_{M}$ is an over-sampled transmit steering vector of $\boldsymbol{\alpha}(M_{A},\varphi_{u,1})$ with an interval of $2/K$. As shown in Fig. 1, we illustrate the envelope of $\boldsymbol{r}_{M}$, where the position corresponding to the peak of $\boldsymbol{r}_{M}$ is the channel AoD with quantization error of $2/K$. In order to apply fast algorithms such as trichotomy search to find the peak of the curve, we have to first find the main lobe with the width of $4/M_{A}$; otherwise these algorithms may stop the search at the peak of side lobes.

Now we propose $\mathbf{Algorithm~{}\ref{alg2}}$ to fast search the largest entry of $\boldsymbol{R}_{u}^{v}$, considering the structure of the steering vectors. $\mathbf{Algorithm~{}\ref{alg2}}$ includes two stages. In the first stage, we find the main lobe of $\boldsymbol{r}_{M}$ and $\boldsymbol{r}_{N}$ without oversampling. In the second stage with oversampling, we apply the trichotomy search to find the peak of the main lobe.

In the first stage from step 3 to step 6, we search the main lobe, which is formed by two adjacent columns and two adjacent rows of $\bar{\boldsymbol{H}}_{u}^{v}$ defined in (12) [25]. We obtain the beamspace receiving matrix $\bar{\boldsymbol{R}}_{u}^{v}\in{\mathbb{C}^{N_{A}\times{M_{A}}}}$ at step 3 as

[TABLE]

where $\bar{\boldsymbol{n}}^{v}\triangleq\boldsymbol{D}(N_{A},N_{A})^{H}\boldsymbol{W}^{H}\widetilde{\boldsymbol{n}}\boldsymbol{F}_{u}^{H}\boldsymbol{D}(M_{A},M_{A})$ is a noise term. We first find two adjacent columns indexed by {$s_{q}$, $s_{q}+1$} with the largest channel power from $\bar{\boldsymbol{H}}_{u}^{v}$ at step 4 via

[TABLE]

where $\bar{\boldsymbol{R}}_{u,q}^{v}\in{\mathbb{C}^{N_{A}\times{2}}}$ represents a submatrix consisted of two consecutive columns of $\bar{\boldsymbol{R}}_{u}^{v}$, with column indices denoted as $q$ and $q+1$. Similarly, we find two adjacent rows indexed by {$s_{p}$,$s_{p}+1$} with the largest channel power from $\bar{\boldsymbol{R}}_{u}^{v}$ at step 4 via

[TABLE]

where $\bar{\boldsymbol{R}}_{u,p}^{v}\in{\mathbb{C}^{2\times{M_{A}}}}$ represents a submatrix consisted of two consecutive rows of $\bar{\boldsymbol{R}}_{u}^{v}$, with row indices denoted as $p$ and $p+1$. By finding out the largest two adjacent columns {$s_{q}$, $s_{q}+1$} and rows {$s_{p}$, $s_{p}+1$} of beamspace channel matrix $\bar{\boldsymbol{H}}_{u}^{v}$, the search of AoA and AoD can be limited to the range of

[TABLE]

respectively, where

[TABLE]

In this way, we narrow down the search space of the AoA and AoD from $[-1,1]$ to $\boldsymbol{\Gamma}$ and $\boldsymbol{\Upsilon}$, respectively. As shown in Fig. 2, the two adjacent rows and two adjacent columns found in the first stage are illustrated in a grey square area.

In the second stage from step 8 to step 14, we find the coordinates of the largest entry of $\boldsymbol{R}_{u}^{v}$ corresponding to the AoA in $\boldsymbol{\Gamma}$ and AoD in $\boldsymbol{\Upsilon}$. Note that the quantization error of the AoD and AoA is reduced from $2/M_{A}$ and $2/N_{A}$ to both $2/K$ by oversampling. We apply the trichotomy search. The two points that divide $\boldsymbol{\Gamma}$ into three equal parts are $2\Gamma_{1}/3+\Gamma_{2}/3$ and $\Gamma_{1}/3+2\Gamma_{2}/3$. Similarly, the two points that divide $\boldsymbol{\Upsilon}$ into three equal parts are $2\Upsilon_{1}/3+\Upsilon_{2}/3$ and $\Upsilon_{1}/3+2\Upsilon_{2}/3$. These four points marked in light blue in Fig. 2 can divide the area of $\boldsymbol{R}_{u}^{v}$ corresponding to $\boldsymbol{\Gamma}$ and $\boldsymbol{\Upsilon}$ into nine smaller areas. The entries corresponding to these four points are

[TABLE]

where ${quan}()$ is the quantization function to quantize the consecutive $\theta$ into $K$ discrete points, which is defined as

[TABLE]

Then we compare the amplitude of $Q_{1}$, $Q_{2}$, $Q_{3}$ and $Q_{4}$ to find the smallest one, which is expressed as

[TABLE]

We delete the parts that include $Q_{min}$ and update $\boldsymbol{\Gamma}$ and $\boldsymbol{\Upsilon}$. As shown in Fig. 2, $Q_{min}$ is marked in dark blue, where the entries within the area of the updated $\boldsymbol{\Gamma}$ and $\boldsymbol{\Upsilon}$ are marked in dark grey. For example, if $Q_{min}=Q_{1}$, the updated $\Gamma_{1}$, $\Gamma_{2}$, $\Upsilon_{1}$ and $\Upsilon_{2}$, denoted as $\bar{\Gamma}_{1}$, $\bar{\Gamma}_{2}$, $\bar{\Upsilon}_{1}$ and $\bar{\Upsilon}_{2}$, respectively, can be represented as

[TABLE]

We repeat the procedures until $\Gamma_{2}-\Gamma_{1}\leq 2/K$ and $\Upsilon_{2}-\Upsilon_{1}\leq 2/K$, which means the resolution $2/K$ of over-sampling search is reached for both the AoA and AoD. Finally we output the estimated AoA and AoD of LOS path at step 15.

III-C ZO-based Beamspace Channel Estimation Schemes

In the IA-based beamspace channel estimation, we first solve two optimization problems (17) and (27) to design the hybrid combining matrix and hybrid precoding matrix. However, considering that $\boldsymbol{W}^{H}\boldsymbol{W}$ does not exactly equal $\gamma_{N}\boldsymbol{I}_{N_{A}}$ and $\boldsymbol{F}_{u}\boldsymbol{F}_{u}^{H}$ does not exactly equal $\gamma_{M}\boldsymbol{I}_{M_{A}}$, there are possible errors for the search of the largest entry even without noise. Also note that the off-diagonal entries of $\boldsymbol{W}^{H}\boldsymbol{W}$ and $\boldsymbol{F}_{u}\boldsymbol{F}_{u}^{H}$ are not restricted to be zero, which may introduce some interference for beamspace channel estimation.

In this subsection, we design hybrid precoding matrix $\boldsymbol{F}_{u}$ and hybrid combining matrix $\boldsymbol{W}$ so that the coordinates of the largest entry of $\boldsymbol{R}_{u}^{v}$ are the AoA and AoD of the LOS path with the quantization error of $2/K$ for single path. Since the largest entry of $\boldsymbol{R}_{u}^{v}$ essentially depends on the largest entry of $\boldsymbol{r}_{N}$ in (34) and $\boldsymbol{r}_{M}$ in (35), respectively, we will first analyze $\boldsymbol{r}_{N}$ and $\boldsymbol{r}_{M}$.

The $k$th entry of $\boldsymbol{r}_{M}$ can be represented as $\boldsymbol{r}_{M}[k]=\boldsymbol{\alpha}^{H}(M_{A},\varphi_{u,1})\boldsymbol{F}_{u}\boldsymbol{F}_{u}^{H}\boldsymbol{\alpha}(M_{A},-1+2(k-1)/K)$. Therefore, $\boldsymbol{F}_{u}$ should be designed so that the largest entry of $\boldsymbol{r}_{M}$ is $\boldsymbol{r}_{M}[{quan}(\varphi_{u,1})]$, which indicates the quantization error is $2/K$. However, it is difficult to solve this problem for discrete $K$ points. Now we convert this discrete problem into continuous problem to obtain the derivative. To ease the notation, we define $\boldsymbol{\alpha}_{M}(\varphi)\triangleq\boldsymbol{\alpha}(M_{A},\varphi)$ and $\varphi_{g}\triangleq\varphi_{u,1}$. We further define a function of $\varphi$ as \mathcal{R}(\varphi_{g},\varphi)\triangleq\big{|}\boldsymbol{\alpha}_{M}^{H}(\varphi_{g})\boldsymbol{F}_{u}\boldsymbol{F}_{u}^{H}\boldsymbol{\alpha}_{M}(\varphi)\big{|}^{2}, which can be regarded as the continuous version of $\boldsymbol{r}_{M}$. The problem can be formulated as

[TABLE]

where $\varphi_{g}$ is the genuine channel parameter defined in (2) but is unknown to the receiver. It is expected that one of the solutions to (45) is $\varphi_{g}$.

Define $\boldsymbol{F}^{+}_{u}\triangleq\boldsymbol{F}_{u}\boldsymbol{F}_{u}^{H}$. Note that $\boldsymbol{F}^{+}_{u}$ is a Hermitian matrix, i.e., ${\boldsymbol{F}^{+}_{u}}^{H}=\boldsymbol{F}^{+}_{u}$. The partial derivative of $\mathcal{R}(\varphi_{g},\varphi)$ over $\varphi$ is

[TABLE]

Define a diagonal matrix $\boldsymbol{B}\in{\mathbb{C}^{M_{A}\times{M_{A}}}}$ as

[TABLE]

which is a diagonal matrix with $M_{A}-1$ nonzero entries. Then (III-C) can be further written as

[TABLE]

Since $\varphi_{g}$ is one of the solutions to (45), we have

[TABLE]

Then we have

[TABLE]

Since $\boldsymbol{F}^{+}_{u}$ is a Hermitian matrix, $\boldsymbol{\alpha}^{H}_{M}(\varphi_{g})\boldsymbol{F}^{+}_{u}\boldsymbol{\alpha}_{M}(\varphi_{g})$ is a real number. Therefore we have

[TABLE]

There are two solutions of $\boldsymbol{F}^{+}_{u}$ to (51), denoted as Solution 1 and Solution 2. Solution 1 of $\boldsymbol{F}^{+}_{u}$ satisfies

[TABLE]

while Solution 2 of $\boldsymbol{F}^{+}_{u}$ satisfies

[TABLE]

In terms of Solution 1, if (52) is satisfied, we set $\varphi_{g}=-1,-1+2/M_{A},\ldots,-1+2(M_{A}-1)/M_{A}$ and obtain $M_{A}$ equations. We combine these $M_{A}$ equations together, having

[TABLE]

where $\boldsymbol{A}_{1}\in{\mathbb{C}^{M_{A}\times{M_{A}}}}$ is a matrix with zero diagonal entries, leading to

[TABLE]

But the trace of the expression on the left side of (54) is

[TABLE]

where $(a)$ is true because $\boldsymbol{F}^{+}_{u}$ is positive semi-definite, and the equality of $(a)$ holds only when $\boldsymbol{F}^{+}_{u}=\boldsymbol{0}_{M_{A}}$. Simultaneously satisfying (55) and (III-C) leads to

[TABLE]

However, in practice $\boldsymbol{F}^{+}_{u}$ can not be zero matrix due to the power constraint in (27). Therefore Solution 1 is meaningless.

In terms of Solution 2, we first rewrite the expression on the left side of (53) as

[TABLE]

Therefore (53) can be further written as

[TABLE]

where $\xi$ is a constant. Similar to (54), we have

[TABLE]

where $\boldsymbol{A}_{2}\in{\mathbb{C}^{M_{A}\times{M_{A}}}}$ is a matrix whose diagonal entries all equal $\xi$, leading to

[TABLE]

The trace of the expression on the left of (60) is \textrm{Tr}\big{(}\boldsymbol{F}^{+}_{u}\big{)}/M_{A} according to (III-C). Combining (61) and (III-C), we have

[TABLE]

Now we solve (59) to obtain the solution of $\boldsymbol{F}^{+}$. The expression on the left side of (59) is essentially in a quadratic form, which can be written in detail as

[TABLE]

Based on (59), (62) and (63), we have

[TABLE]

Since in practice $\varphi_{g}$ can be any value in [-1,1], $\boldsymbol{F}^{+}_{u}$ should satisfy

[TABLE]

which means the summation of off-diagonal entries on the same line parallel to the main diagonal line of $\boldsymbol{F}^{+}_{u}$ is zero. To simplify the problem, we set all off-diagonal entries of $\boldsymbol{F}^{+}_{u}$ to be zero, a.k.a., zero off-diagonal (ZO). In this way, we have to shut down some antenna ports for ZO-based schemes, while the number of active RF chains keeps the same. In fact, according to the existing literature [26, 27, 28], some antenna ports are shut down in mmWave massive MIMO systems. In [26] and [27], some antenna ports are shut down to achieve the antenna selection aiming at the sum-rate maximization. In [28], some antenna ports are shut down to achieve wide mainlobe of beams for mmWave beam training based on hierarchical codebooks. Since $\boldsymbol{F}_{u}\in{\mathbb{C}^{M_{A}\times T_{1}}}$ is a tall matrix with $M_{A}>T_{1}$ and $\boldsymbol{F}^{+}_{u}=\boldsymbol{F}_{u}\boldsymbol{F}_{u}^{H}$, the rank of $\boldsymbol{F}^{+}_{u}$ is no more than $T_{1}$, indicating that there are only $T_{1}$ nonzero entries of the diagonal matrix $\boldsymbol{F}^{+}_{u}$. We set the $T_{1}$ nonzero entries to be the same $\gamma_{M}\sqrt{M_{A}/T_{1}}$ and the left $M_{A}-T_{1}$ diagonal entries to be zero.

Given $\boldsymbol{F}^{+}_{u}$, we can obtain $\boldsymbol{F}_{u}$ by making SVD of $\boldsymbol{F}^{+}_{u}$, i.e., $\boldsymbol{F}^{+}_{u}=\boldsymbol{U}_{M}\boldsymbol{\Sigma}_{M}\boldsymbol{U}_{M}^{H}$, where $\boldsymbol{U}_{M}\in{\mathbb{C}^{M_{A}\times{T_{1}}}}$ is an unitary matrix and $\boldsymbol{\Sigma}_{M}\in{\mathbb{C}^{T_{1}\times{T_{1}}}}$ is a real diagonal matrix. Then we can obtain $\boldsymbol{F}_{u}$ by

[TABLE]

Similar to $\boldsymbol{F}^{+}$, we define $\boldsymbol{W}^{+}\triangleq\boldsymbol{W}^{H}\boldsymbol{W}$, where $\boldsymbol{W}^{+}$ should satisfy

[TABLE]

We set $\boldsymbol{W}^{+}$ to be a diagonal matrix, where $T_{3}$ nonzero diagonal entries are set to be the same $\gamma_{N}\sqrt{N_{A}/T_{3}}$ and the left $N_{A}-T_{3}$ diagonal entries are zero. Given $\boldsymbol{W}^{+}$, we can obtain $\boldsymbol{W}$ by making SVD of $\boldsymbol{W}^{+}$, i.e., $\boldsymbol{W}^{+}=\boldsymbol{U}_{N}\boldsymbol{\Sigma}_{N}\boldsymbol{U}_{N}^{H}$, where $\boldsymbol{U}_{N}\in{\mathbb{C}^{N_{A}\times{T_{3}}}}$ is an unitary matrix and $\boldsymbol{\Sigma}_{N}\in{\mathbb{C}^{T_{3}\times{T_{3}}}}$ is a diagonal matrix. Then we can obtain $\boldsymbol{W}$ by

[TABLE]

According to different layout of nonzero diagonal entries, now we propose a scattered zero off-diagonal (SZO) and a concentrated zero off-diagonal (CZO) based beamspace channel estimation scheme. In the SZO scheme, the nonzero diagonal entries of $\boldsymbol{F}_{u}^{+}$ and $\boldsymbol{W}^{+}$ are uniformly distributed with the same interval. In the CZO scheme, the nonzero diagonal entries of $\boldsymbol{F}_{u}^{+}$ and $\boldsymbol{W}^{+}$ are concentrated on the upper left corner of the matrix.

III-C1 SZO-based Beamspace Channel Estimation Scheme

We design $\boldsymbol{F}^{+}_{u}$ as

[TABLE]

where $\boldsymbol{Z}_{N}$ denotes an $N\times{N}$ matrix where only the entry on the upper-left corner is one and all the other entries are zero. We design $\boldsymbol{W}^{+}$ as

[TABLE]

Then we analyze the largest entry of $\boldsymbol{r}_{M}$. The amplitude of $\boldsymbol{r}_{M}[k]$ can be derived by (III-C1).

It is verified by (III-C1) that the off-diagonal entries of $\boldsymbol{F}_{u}^{+}$ are all zero. We observe that \big{|}\boldsymbol{r}_{M}[k]\big{|} is essentially the inner product between a steering vector $\boldsymbol{\alpha}(T_{1},M_{A}(-1+2(k-1)/K)/T_{1})$ and a steering vector $\boldsymbol{\alpha}(T_{1},\varphi_{g}M_{A}/T_{1})$, indicating that \big{|}\boldsymbol{r}_{M}[k]\big{|} is maximized when the angles of these two steering vectors are the closest, i.e.,

[TABLE]

Note that the angles of these two steering vectors may be out of $[-1,1]$ because $M_{A}/T_{1}>1$. Therefore we need to add the term of $2l$ in (73) to guarantee that $(-1+2(k-1)/K-\varphi_{g})M_{A}/T_{1}+2l$ is within $[-1,1]$. From (73), we obtain

[TABLE]

It is seen that the number of solutions to (74) is $\lfloor M_{A}/T_{1}\rfloor$. Define $\widehat{\varphi}_{g}\triangleq-1+2(\widehat{k}-1)/K$ as an estimation of $\varphi_{g}$. We have

[TABLE]

Suppose $T_{1}K/M_{A}$ to be an integer. Then (75) can be rewritten as

[TABLE]

where $2\langle(\varphi_{g}+1)K/2\rangle/K-1$ is essentially the quantization of $\varphi_{g}$ with resolution of $2/K$. It seen from (76) that $\widehat{\varphi}_{g}$ is periodic with $l$, indicating that the main lobe and the side lobes of $|\boldsymbol{r}_{M}|$ have the same envelope, which is illustrated in the third sub-figure of Fig. 1.

Similarly, we define $\theta_{g}\triangleq\theta_{u,1}$ to ease the notation, and further define $\widehat{\theta}_{g}$ as an estimation of $\theta_{g}$, we have

[TABLE]

It seen from (77) that $\widehat{\theta}_{g}$ is periodic with $i$.

In order to eliminate the uncertainty of $l$ in (76) and $i$ in (77), we resort to beam training based on codebook [29] to find the main lobes of $|\boldsymbol{r}_{M}|$ and $|\boldsymbol{r}_{N}|$, and then design $\boldsymbol{F}_{u}$ and $\boldsymbol{W}$ to estimate the AoA and AoD of the LOS path within the main lobe.

(i) Codebook Design and Beam Training: As shown in Fig. 3, a typical hierarchical codebook has $S$ layers, which satisfy $N=2^{S}$, where $N$ is the number of antennas. In the $s(s=1,2,\ldots,S)$th layer, there are $2^{s}$ codewords with the same beam width but different steering angles. The union of beam angle of all the codewords in each layer is $[-1,1]$. Denote $\boldsymbol{c}(s,n)$ as the $n(n=1,2,\ldots,2^{s})$th codeword in the $s$th layer, covering the beam angle

[TABLE]

where $\omega_{1}(s,n)\triangleq-1+(n-1)/2^{s-1}$ and $\omega_{2}(s,n)\triangleq-1+n/2^{s-1}$. Now we propose an integration-based codebook design method. The codeword $\boldsymbol{c}(s,n)$ is the integration of steering vectors $\boldsymbol{\alpha}(N,\theta)$ with $\theta$ from $\omega_{1}(s,n)$ to $\omega_{2}(s,n)$, which is expressed in (III-C1). In Fig. 4, we compare the integration-based codebook design method with the existing sparse-based codebook design method [8], with respect to the codewords in the first and second layer of the codebook. It shows the codewords designed by the integration-based method have the same beam pattern as those designed by the sparse-based method. However, our method has a closed-form expression, which is much easier for beam generation than the sparse-based method.

After transmitting the beams formed by the codewords of the $S$th layer, the AoA and AoD can be coarsely estimated within $[-1+2(n_{W}-3/2)/N_{A},-1+2(n_{W}-1/2)/N_{A}$ and $[-1+2(n_{F}-3/2)/M_{A},-1+2(n_{F}-1/2)/M_{A}]$, respectively, where $n_{F}\in\{1,2,\ldots,M_{A}\}$ and $n_{W}\in\{1,2,\ldots,N_{A}\}$ are the codeword indices of the $S$th layer. After finishing the search of the codebook, the resolution of the AoA and AoD is $2/N_{A}$ and $2/M_{A}$, respectively.

(ii) Complementary Channel Estimation: After the codebook training in (i), the main lobes of $|\boldsymbol{r}_{M}|$ and $|\boldsymbol{r}_{N}|$ are identified. However, the estimation precision of the AoA and AoD is not satisfied. Therefore, we make complementary channel estimation after the codebook training. During the complementary channel estimation, limited pilot transmission with small pilot overhead is needed.

As illustrated in the third sub-figure of Fig. 1, the width of the main lobe of $|\boldsymbol{r}_{M}|$ is $4/M_{A}$, which requires the period of $2T_{1}/M_{A}$ no smaller than $4/M_{A}$, resulting in $T_{1}\geq 2$. Similarly, we have $T_{3}=T_{2}N_{R}\geq 2$, resulting in $T_{2}\geq 1$. We set $T_{1}=2$ and $T_{2}=1$, meaning that we only need two different hybrid precoding matrices $\boldsymbol{F}_{u}$ and one hybrid combining matrix $\boldsymbol{W}$, where the overhead of pilot training is substantially reduced compared to the IA-based scheme.

Finally, we run the second stage of $\mathbf{Algorithm~{}\ref{alg2}}$ with $\Gamma_{1}\triangleq-1+2(n_{W}-3/2)/N_{A}$, $\Gamma_{2}\triangleq-1+2(n_{W}-1/2)/N_{A}$, $\Upsilon_{1}\triangleq-1+2(n_{F}-3/2)/M_{A}$ and $\Upsilon_{2}\triangleq-1+2(n_{F}-1/2)/M_{A}$. Note that we cannot run the first stage of $\mathbf{Algorithm~{}\ref{alg2}}$ to find the main lobe of $|\boldsymbol{r}_{M}|$ and $|\boldsymbol{r}_{N}|$, since the main lobe and the side lobes have the same envelope, as shown in the third sub-figure of Fig. 1. We can only use the beam training based on codebook to find the main lobe.

III-C2 CZO-based Beamspace Channel Estimation Scheme

In this scheme, the nonzero diagonal entries of $\boldsymbol{F}_{u}^{+}$ and $\boldsymbol{W}^{+}$ are concentrated on the upper left corner of the matrix. We design $\boldsymbol{F}^{+}_{u}$ and $\boldsymbol{W}^{+}$ as

[TABLE]

[TABLE]

Then we analyze the largest entry of $\boldsymbol{r}_{M}$. Similar to (III-C1), the amplitude of $\boldsymbol{r}_{M}[k]$ can be derived by (III-C1). We also observe that \big{|}\boldsymbol{r}_{M}[k]\big{|} is the inner product between a steering vector $\boldsymbol{\alpha}(T_{1},-1+2(k-1)/K))$ and a steering vector $\boldsymbol{\alpha}(T_{1},\varphi_{g})$, indicating that \big{|}\boldsymbol{r}_{M}[k]\big{|} is maximized when the angles of these two steering vectors are the closest, i.e.,

[TABLE]

Note that different with (73), there is no term of $2l$ in (82), because both $-1+2(k-1)/K$ and $\varphi_{g}$ are within $[-1,1]$. From (82), we obtain

[TABLE]

Define $\widehat{\varphi}_{g}\triangleq-1+2(\widehat{k}-1)/K$ as an estimation of $\varphi_{g}$. We have

[TABLE]

where $2\langle(\varphi_{g}+1)K/2\rangle/K-1$ is essentially the quantization of $\varphi_{g}$ with resolution of $2/K$. Unlike (75), $\widehat{\varphi}_{g}$ in (84) is not periodic. Therefore, the envelope of the main lobe and the side lobe is different, as shown in the last sub-figure of Fig. 1, where the beam training based on codebook to identify the main lobe is not necessary. We can directly employ $\mathbf{Algorithm~{}\ref{alg2}}$ to find the largest entry, which corresponds to the AoD and AoA of the LOS path.

III-D Comparisons

Now we compare the proposed three schemes together with the HMC-based, JOINT-based, ECS-based, DCS-based and OCS-based channel estimation scheme in terms of computational complexity, estimation error and total time slots for channel training, which are summarized in Table I.

III-D1 Computational Complexity

As shown in (65), (67) and $\mathbf{Algorithm~{}\ref{alg1}}$, the proposed hybrid precoding and combining schemes are independent of the channel matrix, so we can design the hybrid precoding and combining matrix well before the transmission of pilot sequences. In this way, we do not need to consider the channel coherence time when designing the hybrid precoding and combing matrix. In other words, the computational complexity mainly comes from the search of the largest entry in $\mathbf{Algorithm~{}\ref{alg2}}$. For both the IA-based scheme and CZO-based scheme, in the first stage, we need to compute Frobenius norm of the $N_{A}\times 2$ matrix $\bar{\boldsymbol{R}}_{u,q}^{v}$ in (37) for $M_{A}-1$ times and Frobenius norm of the $2\times M_{A}$ matrix $\bar{\boldsymbol{R}}_{u,p}^{v}$ in (38) for $N_{A}-1$ times, resulting in the complexity to be $\mathcal{O}(4N_{A}(M_{A}-1)+4M_{A}(N_{A}-1))$. In the second stage, we use the trichotomy search to find the largest entry among $2K/N_{A}\times 2K/M_{A}$ entries. Since the length of current $\boldsymbol{\Gamma}$ and $\boldsymbol{\Upsilon}$ is $2/3$ of that of the previous $\boldsymbol{\Gamma}$ and $\boldsymbol{\Upsilon}$ in (44), the number of iterations is $\log_{3/2}(2K/M_{A})$. In the first $\log_{3/2}(2K/N_{A})$ iterations, we compute four trichotomy points in each iteration. The entry on each point is the multiplication of three parts, including a $1\times N_{A}$ row vector of $\boldsymbol{D}(N_{A},K)^{H}$, the $N_{A}\times M_{A}$ matrix $\boldsymbol{W}^{H}\boldsymbol{W}\boldsymbol{H}_{u}\boldsymbol{F}_{u}\boldsymbol{F}_{u}^{H}$, and an $M_{A}\times 1$ column vector of $\boldsymbol{D}(M_{A},K)$, leading to the complexity to be $\mathcal{O}(8(\log_{3/2}(2K/N_{A}))(N_{A}+1)M_{A})$. In the following $\log_{3/2}(2K/M_{A})-\log_{3/2}(2K/N_{A})$ iterations, since the AoA has been estimated, we only need to compute two trichotomy points in each iteration. The entry on each point is the multiplication of two parts, including a $1\times M_{A}$ row vector of $\boldsymbol{D}(N_{A},K)^{H}\boldsymbol{W}^{H}\boldsymbol{W}\boldsymbol{H}_{u}\boldsymbol{F}_{u}\boldsymbol{F}_{u}^{H}$ and an $M_{A}\times 1$ column vector of $\boldsymbol{D}(M_{A},K)$, leading to the complexity to be $\mathcal{O}(4(\log_{3/2}(2K/M_{A})-\log_{3/2}(2K/N_{A}))M_{A})$. Therefore the total computational complexity for both the IA-based scheme and CZO-based scheme is

[TABLE]

For the SZO-based scheme, where only the second stage of $\mathbf{Algorithm~{}\ref{alg2}}$ is needed, we use the trichotomy search to find the largest entry among $K/N_{A}\times K/M_{A}$ entries, resulting in the total computational complexity to be

[TABLE]

For the DCS-based channel estimation scheme [15], where the main lobe is first searched and then the largest entry within the main lobe is further searched by the exhaustive search method, the computational complexity is

[TABLE]

For the ECS-based [11] and OCS-based [17] channel estimation schemes, the largest entry is directly searched by the exhaustive search method, the computational complexity is

[TABLE]

Since $K>N_{A}$ and $K>M_{A}$, the computational complexity of the IA-based scheme, SZO-based scheme and CZO-based scheme is much lower than that of the ECS-based, OCS-based and DCS-based channel estimation scheme. Moreover, for the HMC-based and JOINT-based schemes, since they use hierarchical codebook for beam training instead of computing multiplication in (15) to obtain $\boldsymbol{R}_{u}^{v}$, the computational complexity is almost zero.

III-D2 Estimation Error

For four schemes including the IA-based scheme, ECS-based scheme, DCS-based scheme and OCS-based scheme, the estimation error of the AoA and AoD is larger than $2/K$ without noise. However, as shown in (76) and (84), the estimation error of the AoA and AoD is $2/K$ without noise for both SZO-based and CZO-based scheme. For the HMC-based scheme and JOINT-based scheme, the precision of channel estimation relies on the number of the layers in the hierarchical codebook. In fact, the number of time slots for training is proportional to the layers in the codebook. However, for the IA-based scheme, ECS-based scheme, DCS-based scheme, OCS-based scheme, SZO-based scheme and CZO-based scheme, $K$ can be set independent of the number of time slots for training.

As shown in the last sub-figure of Fig. 1 for the CZO-based scheme, the width of main lobe of $\boldsymbol{r}_{M}$ and $\boldsymbol{r}_{N}$ is $4/T_{1}$ and $4/T_{3}$, respectively, which are larger than that of the IA-based scheme and SZO-based scheme. Therefore, $\boldsymbol{r}_{N}$ and $\boldsymbol{r}_{M}$ in the CZO-based scheme are flatter than those of the IA-based scheme and SZO-based scheme, meaning that the CZO-based scheme has the largest estimation error among the three schemes with the same noise power.

III-D3 Total Time Slots for Channel Training

For the SZO-based scheme where the beam training based on codebook is used, the number of time slots consumed by the beam training is $U(5\log_{2}M_{A}+2\log_{2}(N_{A}/M_{A}))$. Additionally, we need $2U$ time slots to transmit pilot sequences. Therefore, the number of total times slots for channel training is $U(5\log_{2}M_{A}+2\log_{2}(N_{A}/M_{A})+2)$ for the SZO-based scheme. However, for the IA-based scheme, CZO-based scheme, ECS-based scheme, DCS-based scheme, OCS-based scheme, HMC-based scheme and JOINT-based scheme, the total time slots are flexible, depending on $T_{1}$ and $T_{2}$.

IV Simulation Results

Now we evaluate the performance of three proposed schemes. We consider uplink transmission of a multi-user mmWave massive MIMO system. The BS serving $U=4$ users has $N_{A}=64$ antennas and $N_{R}=4$ RF chains, while each user has $M_{A}=16$ antennas and $M_{R}=1$ RF chain. Suppose we use 6 bit and 4 bit digital phase shifters at the BS and users, respectively. The number of resolvable paths in mmWave channel is set to be $L_{u}=3$, while $g_{u,1}\thicksim\mathcal{CN}(0,1)$ and $g_{u,i}\thicksim\mathcal{CN}(0,0.01)$ for $i=2,3$. We use $K=1024$ over-sampling steering vectors. For SZO-based scheme, the number of total time slots for channel training is fixed to be $U(5\log_{2}M_{A}+2\log_{2}(N_{A}/M_{A})+2)=144$. For the IA-based scheme, SZO-based scheme, ECS-based scheme, HMC-based scheme and JOINT-based scheme, the total time slots for channel training are flexible.

As shown in Fig. 5, we compare the channel estimation performance in terms of normalized mean square error (NMSE) for the proposed three schemes, the DCS-based scheme [15] and the OCS-based scheme [17] with different SNR. The NMSE is defined as $\frac{1}{U}\operatorname{E}\left\{\sum_{u=1}^{U}\sqrt{(\hat{\theta}_{u,1}-\theta_{u,1})^{2}+(\hat{\varphi}_{u,1}-\varphi_{u,1})^{2}}\right\}$, which reflects the estimation accuracy of AoA and AoD. We set $T_{1}=4$ and $T_{2}=4$. Then $T_{3}=T_{2}N_{R}=16$. In order to make fair comparisons, we set the total time slots of the DCS-based scheme, OCS-based scheme, ECS-based scheme, HMC-based scheme and JOINT-based scheme the same as the proposed schemes. The number of total time slots for pilot training is fairly set to be $UT_{1}T_{2}=64$. Given $T_{1}T_{2}$, the number of codewords in the last layer of hierarchical codebook is $2^{\langle T_{1}T_{2}/5\rangle}=8$. We also include the ideal case where $\boldsymbol{W}^{H}\boldsymbol{W}=\gamma_{N}\boldsymbol{I}_{N_{A}}$ and $\boldsymbol{F}_{u}\boldsymbol{F}_{u}^{H}=\gamma_{M}\boldsymbol{I}_{M_{A}}$ for comparisons.

It is observed from Fig. 5 that both the IA-based scheme and SZO-based scheme outperform the DCS-based, OCS-based, ECS-based, HMC-based and JOINT-based schemes. At SNR of 15dB, the IA-based scheme has 92.8%, 66.9%, 94.2%, 96.5% and 96.5% performance improvement compared with the DCS-based, OCS-based, ECS-based, HMC-based and JOINT-based schemes, respectively, while the SZO-based scheme has 96.8%, 85.2%, 97.4%, 98.4% and 98.4% improvement compared with the DCS-based, OCS-based, ECS-based, HMC-based and JOINT-based schemes, respectively. The reason for the performance of the DCS-based and OCS-based schemes is that both the DCS-based and OCS-based schemes employ random precoding and random combining matrix, which are not optimal. Although both the DCS-based and OCS-based schemes outperform the CZO-based scheme in low SNR region, at high SNR region such as SNR of 20dB, the CZO-based scheme has 86.8% and 34.5% improvement compared with the DCS-based and OCS-based schemes, respectively. The reason is that the width of the main lobe in CZO-based scheme is larger than DCS-based and OCS-based scheme. Besides, both the SZO-based scheme and CZO-based scheme can approach the ideal performance as SNR increases, which shows the superiority of SZO-based scheme and CZO-based scheme in their small estimation error. The reason for the unsatisfactory performance of the ECS-based scheme is that the ECS-based scheme does not consider the power leakage due to the limited beamspace resolution. Therefore the ideal sparse property of beamspace channel is impaired. The reason for the unsatisfactory performance of the HMC-based and JOINT-based schemes is the number of codewords in the last layer of hierarchical codebook is only 8 and much smaller than $K=1024$, where the total time slots for training are set to be the same with other schemes for fair comparisons. In addition, compared with the CZO-based scheme, the IA-based scheme has smaller estimation error in the low SNR region. The IA-based scheme can achieve better NMSE performance than the CZO-based scheme when SNR $<$ 20 dB, which means we can use the IA-based scheme for better NMSE performance when SNR $<$ 20dB.

As shown in Fig. 6, we compare the sum-rate for the proposed three schemes, the DCS-based, OCS-based, ECS-based, HMC-based and JOINT-based schemes. It is seen that the proposed three schemes achieve better performance than the other schemes, especially in low SNR region. At SNR of 10dB, the performance gap between the proposed three schemes and the ideal case is around 1dB.

As shown in Fig. 7, we compare the channel estimation performance in terms of NMSE for different schemes with different number of total time slots for channel training, which is $T=UT_{1}T_{2}$. Since the number of total time slots for channel training is fixed to be 144 for the SZO-based scheme, SZO-based scheme is not included for comparison. We fix SNR to be 15dB. Given $T$, we set $T_{1}=T_{2}=\sqrt{T/U}$. It is seen that when $T$ is small, the CZO-based scheme performs the best; when $T$ is large, the IA-based scheme performs the best. When $T=16$, the CZO-based scheme has 78.2%, 77.8%, 83.5%, 76.0% and 76.1% performance improvement compared with the DCS-based, OCS-based, ECS-based, HMC-based and JOINT-based schemes, respectively. When $T=100$, the IA-based scheme has 76.3%, 11.2%, 81.4%, 91.4% and 91.5% improvement compared with the DCS-based, OCS-based, ECS-based, HMC-based and JOINT-based schemes, respectively. In the ideal case that both $\boldsymbol{W}^{H}\boldsymbol{W}=\gamma_{N}\boldsymbol{I}_{N_{A}}$ and $\boldsymbol{F}_{u}\boldsymbol{F}_{u}^{H}=\gamma_{M}\boldsymbol{I}_{M_{A}}$ can be achieved, the number of total time slots for channel training is $UM_{A}N_{A}=4096$. It is seen from Fig. 7 that with substantially reduced training overhead, i.e., $T\ll 4096$, satisfactory performance can almost be achieved, e.g., 0.0027 of NMSE with $T=100$ for the IA-based scheme. The number of total time slots for channel training is flexible and fixed for the IA-based scheme and the SZO-based scheme, respectively, which indicates the flexibility of the IA-based scheme. To achieve the NMSE of 0.0027, $T=100$ is enough for the IA-based scheme, while $T=144$ is fixed for the SZO-based scheme.

As shown in Fig. 8, we compare the sum-rate for different number of total time slots for channel training. We fix SNR to be 15dB. We observe that the CZO-based scheme can achieve better performance than the DCS-based, OCS-based, ECS-based, HMC-based and JOINT-based schemes when $T$ is small. When $T$ is close to 100, the sum-rate of the IA-based scheme and CZO-based scheme is almost invariant, indicating that $T=100$ is enough to achieve the maximal sum-rate.

As shown in Fig. 9, we verify the convergence of $\mathbf{Algorithm~{}\ref{alg1}}$. Suppose we use 5, 6 and 7 bit digital phase shifters at the BS, respectively. It is seen that $\epsilon$ in (26) decreases rapidly as the number of iterations grows, and $\epsilon$ decreases faster with higher resolution of digital phase shifters. Using 6 bit digital phase shifters at the BS, $\epsilon$ is smaller than $10^{-3}$ when the number of iterations is larger than 25, which means 25 iterations is enough for $\delta=10^{-3}$.

V Conclusions

This paper has investigated beamspace channel estimation for multi-user mmWave massive MIMO system. A framework of beamspace channel estimation has been proposed. Then based on the this framework, three channel estimation schemes have been proposed. These schemes together with the existing channel estimation schemes have been compared in terms of computational complexity, estimation error and total time slots for channel training. Simulation results have shown that the proposed schemes outperform the existing schemes and can approach the performance of the ideal case with substantially reduced training overhead. Future work will focus on the design of hybrid precoding and hybrid combining for multi-user data transmission regarding the energy efficiency and the theoretical proof of the convergence of iterative hybrid precoding and combining design.