Super-resolution channel estimation for mmWave massive MIMO with hybrid precoding

TL;DR

This paper introduces an iterative reweighting method for super-resolution channel estimation in mmWave massive MIMO systems with hybrid precoding, overcoming resolution loss of traditional compressive sensing approaches.

Contribution

It proposes a novel IR-based super-resolution channel estimation scheme with gradient descent optimization and SVD-based preconditioning, improving accuracy over conventional methods.

Findings

Enhanced channel estimation accuracy demonstrated in simulations

Outperforms traditional compressive sensing schemes

Reduced computational complexity through SVD preconditioning

Abstract

Channel estimation is challenging for millimeter-wave (mmWave) massive MIMO with hybrid precoding, since the number of radio frequency (RF) chains is much smaller than that of antennas. Conventional compressive sensing based channel estimation schemes suffer from severe resolution loss due to the channel angle quantization. To improve the channel estimation accuracy, we propose an iterative reweight (IR)-based super-resolution channel estimation scheme in this paper. By optimizing an objective function through the gradient descent method, the proposed scheme can iteratively move the estimated angle of arrivals/departures (AoAs/AoDs) towards the optimal solutions, and finally realize the super-resolution channel estimation. In the optimization, a weight parameter is used to control the tradeoff between the sparsity and the data fitting error. In addition, a singular value decomposition (SVD)-based preconditioning is developed to reduce the computational complexity of the proposed scheme. Simulation results verify the better performance of the proposed scheme than conventional solutions.

Equations

$$\text{vspace}\ast -1mm\mathbf{r}={\mathbf{Q}}^H\mathbf{HPs}+\mathbf{n},$$

$$\mathbf{H}=\sum _{l=1}^L{z}_l{\mathbf{a}}_{\mathrm{R}}\left({\varphi }_{\mathrm{R},l}^{\mathrm{azi}},{\varphi }_{\mathrm{R},l}^{\mathrm{ele}}\right){\mathbf{a}}_{\mathrm{T}}^H\left({\varphi }_{\mathrm{T},l}^{\mathrm{azi}},{\varphi }_{\mathrm{T},l}^{\mathrm{ele}}\right)\text{vspace}\ast -1.5mm$$

$$\mathbf{a}\left({\varphi }^{\mathrm{azi}},{\varphi }^{\mathrm{ele}}\right)={\left[1,e^{j2\pi d\sin {\varphi }^{\mathrm{azi}}\sin {\varphi }^{\mathrm{ele}}/\lambda },\cdots ,e^{j2\pi \left(N_1-1\right)d\sin {\varphi }^{\mathrm{azi}}\sin {\varphi }^{\mathrm{ele}}/\lambda }\right]}^T$$

$$\otimes {\left[1,e^{j2\pi d\cos {\varphi }^{\mathrm{ele}}/\lambda },\cdots ,e^{j2\pi \left(N_2-1\right)d\cos {\varphi }^{\mathrm{ele}}/\lambda }\right]}^T,$$

$$\text{vspace}\ast -1mm\mathbf{a}\left(\varphi \right)={\left[1,e^{j2\pi d\sin \varphi /\lambda },\cdots ,e^{j2\pi \left(N-1\right)d\sin \varphi /\lambda }\right]}^T.\text{vspace}\ast -0.5mm$$

$$\mathbf{H}={\mathbf{A}}_{\mathrm{R}}\left({\mathbf{\theta }}_{\mathrm{R}}\right)\mathrm{diag}\left(\mathbf{z}\right){\mathbf{A}}_{\mathrm{T}}^H\left({\mathbf{\theta }}_{\mathrm{T}}\right),\text{vspace}\ast -1mm$$

$${\mathbf{y}}_{p,m}={\mathbf{W}}_m^H\mathbf{H}{\mathbf{x}}_p+{\mathbf{n}}_{p,m}.\text{vspace}\ast -0.5mm$$

$$\mathbf{Y}={\mathbf{W}}^H\mathbf{HX}+\mathbf{N}.\text{vspace}\ast -0.5mm$$

$$\underset{\hat{\mathbf{z}},{\hat{\mathbf{\theta }}}_{\mathrm{R}},{\hat{\mathbf{\theta }}}_{\mathrm{T}}}{\min }{\Vert \hat{\mathbf{z}}\Vert }_0,\mathrm{s}.t.{\Vert \mathbf{Y}-{\mathbf{W}}^H\hat{\mathbf{H}}\mathbf{X}\Vert }_{\mathrm{F}}\le \epsilon ,$$

$$\frac{\underset{\mathbf{z},{\mathbf{\theta }}_{\mathrm{R}},{\mathbf{\theta }}_{\mathrm{T}}}{\min }F\left(\mathbf{z}\right)\text{buildrel}\Delta }{=\sum _{l=0}^L\log \left({\left|z_l\right|}^2+\delta \right),\mathrm{s}.t.{\Vert \mathbf{Y}-{\mathbf{W}}^H\hat{\mathbf{H}}\mathbf{X}\Vert }_{\mathrm{F}}\le \epsilon ,}$$

$$\frac{\underset{\mathbf{z},{\mathbf{\theta }}_{\mathrm{R}},{\mathbf{\theta }}_{\mathrm{T}}}{\min }G\left(\mathbf{z},{\mathbf{\theta }}_{\mathrm{R}},{\mathbf{\theta }}_{\mathrm{T}}\right)\text{buildrel}\Delta }{=\sum _{l=1}^L\log \left({\left|z_l\right|}^2+\delta \right)+\lambda {\Vert \mathbf{Y}-{\mathbf{W}}^H\hat{\mathbf{H}}\mathbf{X}\Vert }_{\mathrm{F}}^2.}$$

$$\frac{\underset{\mathbf{z},{\mathbf{\theta }}_{\mathrm{R}},{\mathbf{\theta }}_{\mathrm{T}}}{\min }{S}^{(i)}\left(\mathbf{z},{\mathbf{\theta }}_{\mathrm{R}},{\mathbf{\theta }}_{\mathrm{T}}\right)\text{buildrel}\Delta }{={\lambda }^{-1}{\mathbf{z}}^{\mathrm{H}}{\mathbf{D}}^{(i)}\mathbf{z}+{\Vert \mathbf{Y}-{\mathbf{W}}^H\hat{\mathbf{H}}\mathbf{X}\Vert }_{\mathrm{F}}^2,}$$

$$\frac{{\mathbf{D}}^{(i)}\text{buildrel}\Delta }{=\mathrm{diag}\left(\frac{1}{{\left|{\hat{z}}_1^{(i)}\right|}^2+\delta }\frac{1}{{\left|{\hat{z}}_2^{(i)}\right|}^2+\delta }\cdots \frac{1}{{\left|{\hat{z}}_L^{(i)}\right|}^2+\delta }\right),}$$

$$\frac{{\mathbf{z}}_{\mathrm{opt}}^{(i)}\left({\mathbf{\theta }}_{\mathrm{R}},{\mathbf{\theta }}_{\mathrm{T}}\right)\text{buildrel}\Delta }{=}$$

$$=$$

$$\frac{S_{\mathrm{opt}}^{(i)}\left({\mathbf{\theta }}_{\mathrm{R}},{\mathbf{\theta }}_{\mathrm{T}}\right)\text{buildrel}\Delta }{=\underset{\mathbf{z}}{\min }{S}^{(i)}\left(\mathbf{z},{\mathbf{\theta }}_{\mathrm{R}},{\mathbf{\theta }}_{\mathrm{T}}\right)}$$

$$=$$

$$\text{vspace}-1mm\lambda =\min \left(d/r^{(i)},{\lambda }_{\max }\right),$$

$$\text{vspace}-1mm{r}^{(i)}={\Vert \mathbf{Y}-{\mathbf{W}}^H{\mathbf{A}}_{\mathrm{R}}\left({\hat{\mathbf{\theta }}}_{\mathrm{R}}^{(i)}\right)\mathrm{diag}\left({\hat{\mathbf{z}}}^{(i)}\right){\mathbf{A}}_{\mathrm{T}}^H\left({\hat{\mathbf{\theta }}}_{\mathrm{T}}^{(i)}\right)\mathbf{X}\Vert }_{\mathrm{F}}^2.$$

$${\hat{\mathbf{\theta }}}_{\mathrm{R}}^{(i+1)}={\hat{\mathbf{\theta }}}_{\mathrm{R}}^{(i)}-\eta \cdot {\nabla }_{{\mathbf{\theta }}_{\mathrm{R}}}{S}_{\mathrm{opt}}^{(i)}\left({\hat{\mathbf{\theta }}}_{\mathrm{R}}^{(i)},{\hat{\mathbf{\theta }}}_{\mathrm{T}}^{(i)}\right),$$

$${\hat{\mathbf{\theta }}}_{\mathrm{T}}^{(i+1)}={\hat{\mathbf{\theta }}}_{\mathrm{T}}^{(i)}-\eta \cdot {\nabla }_{{\mathbf{\theta }}_{\mathrm{T}}}{S}_{\mathrm{opt}}^{(i)}\left({\hat{\mathbf{\theta }}}_{\mathrm{R}}^{(i)},{\hat{\mathbf{\theta }}}_{\mathrm{T}}^{(i)}\right),$$

$$\mathbf{Y}=\left({\mathbf{W}}^H{\mathbf{A}}_{\mathrm{R}}\left({\mathbf{\theta }}_{\mathrm{R}}\right)\right)\mathrm{diag}\left(\mathbf{z}\right){\left({\mathbf{X}}^H{\mathbf{A}}_{\mathrm{T}}\left({\mathbf{\theta }}_{\mathrm{T}}\right)\right)}^H+\mathbf{N}.\text{vspace}-1mm$$

$${\sigma }_i\approx \left|z_{l_i}\right|{\Vert {\mathbf{W}}^H{\mathbf{a}}_{\mathrm{R}}\left({\theta }_{\mathrm{R},l_i}^{\mathrm{azi}},{\theta }_{\mathrm{R},l_i}^{\mathrm{ele}}\right)\Vert }_2{\Vert {\mathbf{X}}^H{\mathbf{a}}_{\mathrm{T}}\left({\theta }_{\mathrm{T},l_i}^{\mathrm{azi}},{\theta }_{\mathrm{T},l_i}^{\mathrm{ele}}\right)\Vert }_2,$$

$${\mathbf{u}}_i\approx {\mathbf{W}}^H{\mathbf{a}}_{\mathrm{R}}\left({\theta }_{\mathrm{R},l_i}^{\mathrm{azi}},{\theta }_{\mathrm{R},l_i}^{\mathrm{ele}}\right)/{\Vert {\mathbf{W}}^H{\mathbf{a}}_{\mathrm{R}}\left({\theta }_{\mathrm{T},l_i}^{\mathrm{azi}},{\theta }_{\mathrm{T},l_i}^{\mathrm{ele}}\right)\Vert }_2,$$

$${\mathbf{v}}_i\approx {\mathbf{X}}^H{\mathbf{a}}_{\mathrm{T}}\left({\theta }_{\mathrm{R},l_i}^{\mathrm{azi}},{\theta }_{\mathrm{R},l_i}^{\mathrm{ele}}\right)/{\Vert {\mathbf{X}}^H{\mathbf{a}}_{\mathrm{T}}\left({\theta }_{\mathrm{T},l_i}^{\mathrm{azi}},{\theta }_{\mathrm{T},l_i}^{\mathrm{ele}}\right)\Vert }_2,$$

$$S\left(\mathbf{z},{\mathbf{\theta }}_{\mathrm{R}},{\mathbf{\theta }}_{\mathrm{T}}\right)={\lambda }^{-1}{\mathbf{z}}^H\mathbf{Dz}+\sum _{p=1}^{N_{\mathrm{X}}}{\Vert {\mathbf{y}}_p-{\mathbf{W}}^H{\mathbf{A}}_{\mathrm{R}}\mathrm{diag}\left(\mathbf{z}\right){\mathbf{A}}_{\mathrm{T}}^H{\mathbf{x}}_p\Vert }_2^2$$

$$={\lambda }^{-1}{\mathbf{z}}^H\mathbf{Dz}+\sum _{p=1}^{N_{\mathrm{X}}}{\left({\mathbf{y}}_p-{\mathbf{K}}_p\mathbf{z}\right)}^H\left({\mathbf{y}}_p-{\mathbf{K}}_p\mathbf{z}\right)$$

$$={\mathbf{z}}^H\left({\lambda }^{-1}\mathbf{D}+\sum _{p=1}^{N_{\mathrm{X}}}{\mathbf{K}}_p^H{\mathbf{K}}_p\right)\mathbf{z}-{\mathbf{z}}^H\left(\sum _{p=1}^{N_{\mathrm{X}}}{\mathbf{K}}_p^H{\mathbf{y}}_p\right)-\left(\sum _{p=1}^{N_{\mathrm{X}}}{\mathbf{y}}_p^H{\mathbf{K}}_p\right)\mathbf{z}+\sum _{p=1}^{N_{\mathrm{X}}}{\mathbf{y}}_p^H{\mathbf{y}}_p,$$

$$\frac{\partial S\left(\mathbf{z},{\mathbf{\theta }}_{\mathrm{R}},{\mathbf{\theta }}_{\mathrm{T}}\right)}{\partial \mathbf{z}}=$$

$${\mathbf{z}}_{\mathrm{opt}}\left({\mathbf{\theta }}_{\mathrm{R}},{\mathbf{\theta }}_{\mathrm{T}}\right)={\left({\lambda }^{-1}\mathbf{D}+\sum _{p=1}^{N_{\mathrm{X}}}{\mathbf{K}}_p^H{\mathbf{K}}_p\right)}^{-1}\left(\sum _{p=1}^{N_{\mathrm{X}}}{\mathbf{K}}_p^H{\mathbf{y}}_p\right),$$

$$S_{\mathrm{opt}}\left({\mathbf{\theta }}_{\mathrm{R}},{\mathbf{\theta }}_{\mathrm{T}}\right)=-{\left(\sum _{p=1}^{N_{\mathrm{X}}}{\mathbf{K}}_p^H{\mathbf{y}}_p\right)}^H{\left({\lambda }^{-1}\mathbf{D}+\sum _{p=1}^{N_{\mathrm{X}}}{\mathbf{K}}_p^H{\mathbf{K}}_p\right)}^{-1}$$

$$\cdot \left(\sum _{p=1}^{N_{\mathrm{X}}}{\mathbf{K}}_p^H{\mathbf{y}}_p\right)+\sum _{p=1}^{N_{\mathrm{X}}}{\mathbf{y}}_p^H{\mathbf{y}}_p.$$

$$\frac{\partial S_{\mathrm{opt}}}{\partial {\theta }_{\mathrm{R},l}}$$

$$=-\frac{\partial {\mathbf{v}}^H}{\partial {\theta }_{\mathrm{R},l}}{\mathbf{A}}^{-1}\mathbf{v}+{\mathbf{v}}^H{\mathbf{A}}^{-1}\frac{\partial \mathbf{A}}{\partial {\theta }_{\mathrm{R},l}}{\mathbf{A}}^{-1}\mathbf{v}-{\mathbf{v}}^H{\mathbf{A}}^{-1}\frac{\partial \mathbf{v}}{\partial {\theta }_{\mathrm{R},l}},$$

$$\frac{\partial \mathbf{v}}{\partial {\theta }_{\mathrm{R},l}}$$

Full text

I Introduction

Millimeter-wave (mmWave) massive MIMO has been recognized as a promising technology for future 5G wireless communications[1]. To reduce the hardware cost and power consumption, hybrid precoding has been proposed for practical mmWave massive MIMO systems, where hundreds of antennas are driven by a much smaller number of radio frequency (RF) chains[3, 2]. The analog and digital co-design in hybrid precoding requires accurate channel state information. However, the digital baseband cannot directly access all antennas due to the small number of RF chains, so it is difficult to accurately estimate the high-dimensional MIMO channel [4, 5].

Several novel channel estimation schemes have been recently proposed for mmWave massive MIMO with hybrid precoding[5, 6, 7, 8, 9]. Specifically, [5, 6] proposed the adaptive codebook-based channel sounding scheme, where the transmitter and receiver search for the best beam pair by adjusting the predefined precoding and combining codebooks. However, the channel estimation resolution is limited by the codebook size. [7] was able to achieve better angle estimation by performing an amplitude comparison with respect to the auxiliary beam pair. On the other hand, by exploiting the angular channel sparsity, the on-grid compressive sensing based methods [8, 9] could estimate the channel with reduced training overhead. However, such solutions assumed that the angle of arrivals/departures (AoAs/AoDs) lie in discrete points in the angle domain (i.e., “on-grid” AoAs/AoDs), while the actual AoAs/AoDs are continuously distributed (i.e., “off-grid” AoAs/AoDs) in practice. The assumption of on-grid AoAs/AoDs results in the power leakage problem, which severely degrades the channel estimation accuracy. To solve this resolution limitation caused by the on-grid angle estimation, we propose an iterative reweight (IR)-based super-resolution channel estimation scheme to estimate the off-grid AoAs/AoDs.111 Simulation codes are provided to reproduce the results presented in this paper: http://oa.ee.tsinghua.edu.cn/dailinglong/publications/publications.html.

Specifically, we iteratively optimize the estimates of AoAs/AoDs, to decrease the weighted summation of the sparsity and the data fitting error. The weight controlling the tradeoff between the sparsity and the data fitting error, is iteratively updated to avoid over-fitting or under-fitting. Since the estimated AoAs/AoDs can be moved from the initial angle-domain grids towards the actual off-grid AoAs/AoDs, the proposed scheme is able to achieve the super-resolution channel estimation. In addition, we propose a singular value decomposition (SVD)-based preconditioning method to reduce the computational complexity of the proposed scheme, which is realized by reducing the number of initial candidates of AoAs/AoDs in the IR procedure. Simulation results show that the proposed IR-based super-resolution channel estimation can achieve better performance than conventional solutions.

The contributions of this paper are the follows. We propose a novel IR-based super-resolution channel estimation scheme for mmWave massive MIMO with hybrid precoding. Comparing with the state-of-art schemes schemes such as those in [6, 7, 8], we can achieve super-resolution channel estimation, which means substantially improved estimation accuracy. Moreover, the proposed SVD based preconditioning significantly reduces the computational complexity of the IR procedure, and makes the method practical in mmWave channel estimation.

*Notation: *In this paper, the boldface lower and upper-case symbols denote vectors and matrices. $(\cdot)^{T}$, $(\cdot)^{H}$ and $(\cdot)^{-1}$ denote the transpose, the conjugate transpose, and the inverse of a matrix, respectively. ${\rm{diag}}({\bf{x}})$ is the diagonal matrix with the vector ${\bf{x}}$ on its diagonal. The $\ell_{0}$-norm, $\ell_{2}$-norm, and Frobenius norm are given by $\|\cdot\|_{0}$, $\|\cdot\|_{2}$, and $\|\cdot\|_{\rm F}$, respectively.

II System Model

We consider a hybrid-precoding mmWave massive MIMO with arbitrary array geometry. Let $N_{\rm T}$, $N_{\rm R}$, $N_{\rm T}^{\rm RF}$, and $N_{\rm R}^{\rm RF}$ be the number of transmit antennas, receive antennas, transmitter RF chains, and receiver RF chains, respectively. For practical mmWave massive MIMO with hybrid precoding, the number of RF chains is much smaller than that of antennas, i.e., $N_{\rm T}^{\rm RF}<N_{\rm T}$, $N_{\rm R}^{\rm RF}<N_{\rm R}$ [1, 3, 2]. The system model can be given by

[TABLE]

where ${\bf{r}}\in{\mathbb{C}^{N_{\rm R}^{\rm RF}\times 1}}$ is the received signal, ${\bf{Q}}\in{\mathbb{C}^{N_{\rm R}\times N_{\rm R}^{\rm RF}}}$ is the hybrid combining matrix, ${\bf{H}}\in{\mathbb{C}^{N_{\rm R}\times N_{\rm T}}}$ is the channel matrix, ${\bf{P}}\in{\mathbb{C}^{N_{\rm T}\times N_{\rm T}^{\rm RF}}}$ is the hybrid precoding matrix, ${\bf{s}}\in{\mathbb{C}^{N_{\rm T}^{\rm RF}\times 1}}$ is the transmitted signal, and ${\bf{n}}\in{\mathbb{C}^{N_{\rm R}^{\rm RF}\times 1}}$ is the received noise after combining.

The channel model

[TABLE]

is widely adopted in mmWave massive MIMO systems, and it is regarded almost unchanged within the channel coherence time for channel estimation [6, 7, 8], where $L$ is the number of propagation paths, $L\ll{\rm min}\left(N_{\rm R},N_{\rm T}\right)$, $z_{l}$, $\phi^{\rm azi}_{{\rm R},l}~{}(\phi^{\rm ele}_{{\rm R},l})$ and $\phi^{\rm azi}_{{\rm T},l}~{}(\phi^{\rm ele}_{{\rm T},l})$ are the complex path gain, the azimuth (elevation) AoA and AoD of the $l$-th path, respectively. ${{\bf{a}}_{\rm R}}\!\left(\phi^{\rm azi}_{{\rm R},l},\phi^{\rm ele}_{{\rm R},l}\right)$ and ${{\bf{a}}_{\rm T}}\!\left(\phi^{\rm azi}_{{\rm T},l},\phi^{\rm ele}_{{\rm T},l}\right)$ are the steering vector at the receiver and the steering vector at the transmitter, respectively. These steering vectors depend on the array geometry. Ignoring the subscripts without loss of generality, for the typical $N_{1}\times N_{2}$ uniform planar arrays (UPAs), ${{\bf{a}}}\!\left(\phi^{\rm azi}_{l},\phi^{\rm ele}_{l}\right)$ is given by [2]

[TABLE]

where $d$ is the antenna spacing, $\lambda$ is the wavelength, $\otimes$ denotes the Kronecker product. For uniform linear arrays (ULAs), the steering vector is only determined by one angle [6]

[TABLE]

By defining the normalized spacial angles by ${\theta^{\rm azi}}\buildrel\Delta\over{=}d\sin{\phi^{\rm azi}}\sin{\phi^{\rm ele}}/\lambda$ and ${\theta^{\rm ele}}\buildrel\Delta\over{=}d\cos{\phi^{\rm ele}}/\lambda$, the mmWave channel matrix $\bf H$ in (2) can be also written as

[TABLE]

where ${\bf z}\!=\!\left[z_{1},z_{2},\cdots,z_{L}\right]^{T}$, ${{\bm{\theta}_{\rm R}}}\!=\!\left[{\theta^{\rm azi}_{{\rm R}\!,1}},{\theta^{\rm ele}_{{\rm R}\!,1}},{\theta^{\rm azi}_{{\rm R}\!,2}},{\theta^{\rm ele}_{{\rm R}\!,2}},{\cdots},\right.$

$\left.{\theta^{\rm azi}_{{\rm R}\!,L}},{\theta^{\rm ele}_{{\rm R}\!,L}}\right]^{T}$, ${{\bm{\theta}_{\rm T}}}\!=\!\left[{\theta^{\rm azi}_{{\rm T}\!,1}},{\theta^{\rm ele}_{{\rm T}\!,1}},{\theta^{\rm azi}_{{\rm T}\!,2}},{\theta^{\rm ele}_{{\rm T}\!,2}},{\cdots},{\theta^{\rm azi}_{{\rm T}\!,L}},{\theta^{\rm ele}_{{\rm T}\!,L}}\right]^{T}$, ${{\bf{A}}_{\!\rm R}}\!\left({{\bm{\theta}_{\rm R}}}\right)\!=\!\left[{{{\bf{a}}_{\rm R}}({\theta^{\rm azi}_{{\rm R},1}},{\theta^{\rm ele}_{{\rm R},1}})}~{}{{{\bf{a}}_{\rm R}}({\theta^{\rm azi}_{{\rm R},2}},{\theta^{\rm ele}_{{\rm R},2}})}~{}{\cdots}~{}{{{\bf{a}}_{\rm R}}({\theta^{\rm azi}_{{\rm R},L}},{\theta^{\rm ele}_{{\rm R},L}})}\right]$, ${{\bf{A}}_{\!\rm T}}\!\left({{\bm{\theta}_{\rm T}}}\right)\!=\!\left[{{{\bf{a}}_{\rm T}}({\theta^{\rm azi}_{{\rm T},1}},{\theta^{\rm ele}_{{\rm T},1}})}~{}{{{\bf{a}}_{\rm T}}({\theta^{\rm azi}_{{\rm T},2}},{\theta^{\rm ele}_{{\rm T},2}})}~{}{\cdots}~{}{{{\bf{a}}_{\rm T}}({\theta^{\rm azi}_{{\rm T},L}},{\theta^{\rm ele}_{{\rm T},L}})}\right]$.

Denote ${\bf{x}}={\bf{Ps}}\in{\mathbb{C}^{N_{\rm T}\times 1}}$, where the $i$-th element of $\bf{x}$ is the transmitted signal at the $i$-th transmit antenna. Suppose that the transmitter sends $N_{\rm X}~{}(N_{\rm X}<N_{\rm T})$ different pilot sequences ${\bf{x}}_{1},{\bf{x}}_{2},\cdots,{\bf{x}}_{N_{\rm X}}$. Since the number of RF chains is smaller than the required dimension of received pilot sequence, for each transmit pilot sequence ${\bf{x}}_{p}~{}(1\!\leq\!p\!\leq\!N_{\rm X})$, we use $M$ time slots to obtain an $N_{\rm Y}$-dimensional received pilot sequence ${\bf{y}}_{p}$, where $N_{\rm Y}=MN_{\rm R}^{\rm RF}$. Thus, the training overhead is $T=MN_{\rm X}$. In the $m$-th time slot, we use the combining matrix ${\bf W}_{m}$ to obtain an $N_{\rm R}^{\rm RF}$-dimensional received pilot sequence

[TABLE]

By collecting the received pilots in the $M$ time slots, we have ${{\bf y}_{p}}={\bf W}^{H}{\bf Hx}_{p}+{\bf n}_{p}$, where ${{\bf y}_{p}}=[{{\bf y}_{p,1}^{T}},~{}{{\bf y}_{p,2}^{T}},~{}\cdots,~{}{{\bf y}_{p,M}^{T}}]^{T}\in{\mathbb{C}^{N_{\rm Y}\times 1}}$, ${{\bf W}}=[{{\bf W}_{1}},{{\bf W}_{2}},\cdots,{{\bf W}_{M}}]\in{\mathbb{C}^{N_{\rm R}\times N_{\rm Y}}}$, ${\bf n}_{p}\in{\mathbb{C}^{N_{\rm Y}\times 1}}$ is the noise. By defining ${\bf{Y}}=[{\bf{y}}_{1},{\bf{y}}_{2},\cdots,{\bf{y}}_{N_{\rm X}}]$, ${{\bf{X}}=\left[{\bf{x}}_{1},{\bf{x}}_{2},\cdots,{\bf{x}}_{N_{\rm X}}\right]}$, ${{\bf{N}}=\left[{\bf{n}}_{1},{\bf{n}}_{2},\cdots,{\bf{n}}_{N_{\rm X}}\right]}$, we have

[TABLE]

The estimation of the channel matrix $\bf H$ in (7) is equivalent to the estimation of the number of paths, the normalized spacial angles $({{\bm{\theta}}_{{\rm T}}}$, ${{\bm{\theta}}_{{\rm R}}})$, and path gains ${\bf z}$ for all $L$ paths. Due to the angle-domain sparsity of the channel matrix $\bf H$, the sparse channel estimation problem can be formulated as

[TABLE]

where ${\left\|{\bf{\hat{z}}}\right\|_{0}}$ is the number of non-zero elements of $\bf\hat{z}$, which means the estimated number of paths $\hat{L}$, ${\bf{\hat{H}}}$ is the estimated channel matrix, and $\varepsilon$ is the error tolerance parameter [9].

III Proposed IR-Based Super-Resolution Channel Estimation

III-A Proposed Optimization Formulation

The main difficulty in solving (8) lies in the fact that the $l_{0}$-norm is not computationally efficient for finding the optimal solution. By replacing the $l_{0}$-norm with a log-sum function [10], we have

[TABLE]

where $\delta>0$ ensures that the logarithmic function is well-defined [10], ${\bf\hat{H}}$ is determined by the parameters ${\bf{z}},{{\bm{\theta}}_{\rm R}}$ and ${{\bm{\theta}}_{\rm T}}$ defined in (5). By adding a regularization parameter $\lambda>0$, we can further formulate the problem (9) as a unconstrained optimization problem:

[TABLE]

Moreover, by using an iterative surrogate function instead of the log-sum function, the minimization of $G\left({{\bf{z}},{{\bm{\theta}}_{\rm R}},{{\bm{\theta}}_{\rm T}}}\right)$ is equivalent to the minimization of the surrogate function[10]:

[TABLE]

where ${{\bf{D}}^{(i)}}$ is defined as

[TABLE]

where ${\hat{\bf{z}}^{(i)}}$ is the estimate of ${\bf{z}}$ at the $i$-th iteration.

Then, as proved in Appendix A, we can optimize (11) with regard to the path gains $\bf{z}$, to find the optimal point of $\bf\hat{z}$ and the corresponding optimal value of $S^{(i)}$ as follows:

[TABLE]

[TABLE]

where ${{\bf{K}}_{p}}={\bf{W}}^{H}{{\bf{A}}_{\rm R}}{\rm{diag}}\left({{\bf{A}}_{\rm T}^{{H}}{{\bf{x}}_{p}}}\right)$. After that, we only need to optimize the normalized spatial angles ${{\bm{\theta}}_{\rm R}}$ and ${{\bm{\theta}}_{\rm T}}$ in (14), which will be discussed in the next subsection.

III-B IR-Based Super-Resolution Channel Estimation

In the previous subsection, we have already simplified the constrained optimization problem (8) to an unconstrained angle optimization problem (14). To solve this reformulated problem, now we propose an IR-based super-resolution channel estimation scheme as described in Algorithm 1.

The objective function ${{{S}}^{(i)}}\left({{\bf z}},{{{\bm{\theta}}_{\rm R}},{{\bm{\theta}}_{\rm T}}}\right)$ is the weighted sum of two parts: ${{\bf{z}}^{\rm{H}}}{\bf{Dz}}$ controlling the sparsity of the estimation result and ${\left\|{{{\bf{Y}}-{\bf{W}}^{H}{\bf{{\hat{H}}X}}}}\right\|_{\rm F}}$ denoting the residue. In addition, $\lambda$ is the regularization parameter that controls the tradeoff between the sparsity and the data fitting error.

In the iterative reweighted method [10], $\lambda$ is not fixed but updated in each iteration. To be specific, if the previous iteration is poorly-fitted, we will choose a smaller $\lambda$ to make the estimate sparser. On the other hand, if the previous iteration returns a well-fitted estimate and leads to a small residue, our method will choose a larger $\lambda$ to accelerate the searching for the best-fitting estimate. In the proposed algorithm, $\lambda$ is updated by

[TABLE]

where $d$ is a constant scaling factor, and $\lambda_{\rm max}$ is selected to make the problem well-conditioned, $r^{(i)}$ is the squared residue in the previous step, i.e.,

[TABLE]

The updation of $\lambda$ was discussed in [10] with more details.

The proposed algorithm starts iteration at the angle domain grids. In the $i$-th iteration, our task is to search for new estimates ${\bm{\hat{\theta}}}_{\rm R}^{(i+1)}$ and ${\bm{\hat{\theta}}}_{\rm T}^{(i+1)}$ in the neighborhood of the previous estimates ${\bm{\hat{\theta}}}_{\rm R}^{(i)}$ and ${\bm{\hat{\theta}}}_{\rm T}^{(i)}$ to make the objective function ${S}^{(i)}$ become smaller. This searching can be accomplished via gradient descent method:

[TABLE]

where the gradients can be calculated according to Appendix B, and $\eta$ is the chosen step-length to make sure $S_{\rm opt}^{(i)}\!\!\left({\bm{\hat{\theta}}}_{\rm R}^{(i+1)},{\bm{\hat{\theta}}}_{\rm T}^{(i+1)}\right)\!\leq\!S_{\rm opt}^{(i)}\!\!\left({\bm{\hat{\theta}}}_{\rm R}^{(i)},{\bm{\hat{\theta}}}_{\rm T}^{(i)}\right)$. The estimates become more and more accurate during the iterative searching, until the new estimates are almost the same as the previous ones. With our proposed IR-based super-resolution channel estimation scheme, the estimates of $({\bm{\theta}}_{\rm R},{\bm{\theta}}_{\rm T})$ can be moved from the initial on-grid coarse estimates to its actual off-grid positions, thus the super-resolution channel estimation can be realized.

It is worthy to point out that the sparsity level $L$ is unknown in practice. In the proposed scheme, the sparsity level can be initialized to be larger than the real channel sparsity. During the iteration process, the paths with too small path gains will be regarded as noise instead of real paths. Then, our algorithm prune these paths to make the result sparser. By iteratively pruning these paths, the estimated sparsity level will decrease to the real number of paths.

The computational complexity in each iteration lies in calculating the gradient in Step 5. The computational complexity to calculate the gradient is $\mathcal{O}\left(N_{\rm X}N_{\rm Y}(N_{\rm R}\!+\!N_{\rm T})L^{2}\right)$. As a result, the number of initial candidates $L^{(0)}$ is critical, and it should be as small as possible to make the computation affordable. The problem how to effectively select the initial ${\bm{\hat{\theta}}}_{\rm R}^{(0)}$ and ${\bm{\hat{\theta}}}_{\rm T}^{(0)}$ before the iteration will be discussed in the next section.

III-C SVD-based Preconditioning

In this section, we propose a singular value decomposition (SVD)-based preconditioning as shown in Algorithm 2, to reduce the computational complexity of the IR procedure in the proposed IR-based super-resolution channel estimation scheme. The proposed scheme can find the angle-domain grids nearest to the real AoAs/AoDs. Comparing to using all $N_{\rm R}N_{\rm T}$ angle-domain grids as initial candidates, the preconditioning can significantly reduce the computational complexity of the IR-based super-resolution channel estimation scheme.

Specifically, by applying SVD to the matrix $\bf Y$, we have ${{\bf Y}}={\bf U\Sigma V}^{H}$, where ${\bf\Sigma}={\rm diag}\left(\sigma_{1},\sigma_{2},\cdots,\sigma_{{\rm min}(N_{\rm X},N_{\rm Y})}\right)\in{\mathbb{R}^{N_{\rm Y}\times N_{\rm X}}}$ whose diagonal entries $\sigma_{1}\geq\sigma_{2}\geq\cdots\geq\sigma_{{\rm min}(N_{\rm X},N_{\rm Y})}\geq 0$ are the singular values of $\bf Y$, and ${\bf U}^{H}{\bf U}={\bf I}_{N_{\rm Y}\times N_{\rm Y}}$, ${\bf V}^{H}{\bf V}={\bf I}_{N_{\rm X}\times N_{\rm X}}$. From (5) and (7), we have

[TABLE]

As the noise is small, the largest $L$ singular values and their corresponding singular vectors are approximately determined by the $L$ paths, i.e., for $i=1,2,\cdots,L$, we have

[TABLE]

where ${\bf u}_{i}$ and ${\bf v}_{i}$ are the $i$-th column of $\bf U$ and $\bf V$, respectively, $\left\{l_{1},l_{2},\cdots,l_{L}\right\}$ is a permutation of $\left\{1,2,\cdots,L\right\}$.

Then, in steps 4-5 of Algorithm 2, we search for the coarse estimate of normalized AoAs (AoDs) in the finite set of angle-domain grids $\Omega_{\rm R}$ ($\Omega_{\rm T}$). Take UPA as an example. For an $N_{1}\times N_{2}$ receiver array, the set of grids can be defined by $\Omega_{\rm R}=\left\{(i/N_{1},j/N_{2})|i=0,1,\cdots,N_{1}-1;j=0,1,\cdots,N_{2}-1\right\}$. We can similarly define $\Omega_{\rm T}$ for the transmitter.

In [10], the initial candidates of Algorithm 1 are set to be all the grids, i.e., $L^{(0)}=N_{\rm R}N_{\rm T}$. The computational complexity is $\mathcal{O}\left(N_{\rm X}N_{\rm Y}(N_{\rm R}\!+\!N_{\rm T})N_{\rm R}^{2}N_{\rm T}^{2}\right)$, which is unaffordable when $N_{\rm R}$ and $N_{\rm T}$ are very large. Fortunately, with the proposed SVD-based preconditioning shown in Algorithm 2, the coarse estimates will be used as the initial candidates of Algorithm 1, i.e., $L^{(0)}=N_{\rm init}\approx L$. Thus, the computational complexity after SVD preconditioning is $\mathcal{O}\left(N_{\rm X}N_{\rm Y}(N_{\rm R}\!+\!N_{\rm T})L^{2}\right)$, which is much lower than directly applying the scheme in [10].

IV Simulation Results and Discussions

In this section, simulation results are provided to investigate the performance. We consider the mmWave massive MIMO system with hybrid precoding, where $L=3$, $d={\lambda}/{2}$, $N_{\rm R}=N_{\rm T}=64$, $N_{\rm R}^{\rm RF}=N_{\rm T}^{\rm RF}=4$ and $N_{\rm X}=N_{\rm Y}=32$. The path gains are assumed Gaussian, i.e., $\alpha_{l}\sim{\cal CN}(0,{\sigma_{\alpha}^{2}})$. Each element of the transmitted pilots $\bf X$ satisfies $x_{i,j}={\sqrt{{\rho}/{N_{\rm T}}}}e^{j\omega_{i,j}}$, where $\rho$ is the transmitted power, $\omega_{i,j}$ is the random phase uniformly distributed in $[0,2\pi)$. The signal-to-noise ratio (SNR) is defined by ${\rm SNR}=\frac{\rho\sigma_{\alpha}^{2}}{\sigma_{n}^{2}}$, where $\sigma_{n}^{2}$ is the noise variance. We consider the ULA geometry, so that the adaptive codebook-based channel estimation [6], the auxiliary beam pair based channel estimation [7], and the OMP-based channel estimation [8] can be adopted for performance comparison.

Fig. 1 and Fig. 2 compares the normalized mean square error (NMSE) performance against SNR, under none-line-of-sight (NLoS) and line-of-sight (LoS) channels, respectively. The Rician K-factor is $20$ dB in the LoS scenario. In both cases, the proposed scheme achieves much better NMSE performance when SNR becomes large. Moreover, we show the performance of the proposed scheme when UPA is considered. Both the transmitter and the receiver adopt 64-antenna UPA with 8 rows and 8 columns. We estimate the azimuth and elevation angles at both sides. We can observe that the proposed scheme is also able to achieve super-resolution channel estimation when UPA is used. Since the estimation errors of both azimuth and elevation angles contribute to the NMSE, under the same number of antennas and number of pilot overhead, the NMSE performance of UPA is higher than that of ULA.

Fig. 3 compares the average spectral efficiency when different channel estimation schemes are used. The spectral efficiency is evaluated in the hybrid precoding system [2]. The case with ideal CSI was adopted as the upper bound for performance comparison. It can be observed that the proposed super-resolution channel estimation is able to approach this upper bound. This is because the angle resolution of the proposed scheme does not suffer from limited codebook size or angle quantization. Thus, we can conclude that the proposed scheme can achieve the super-resolution channel estimation.

There is a tradeoff between the channel estimation accuracy and the computation complexity. The proposed scheme with SVD preconditioning is able to achieve much higher channel estimation accuracy, but has a higher computational complexity. The computational complexity of the proposed super-resolution channel estimation scheme with SVD precondistioning is $\mathcal{O}\left(N_{\rm X}N_{\rm Y}(N_{\rm R}\!+\!N_{\rm T})L^{2}\right)$. In comparison, the computational complexity of the OMP-based channel estimation [8] is $\mathcal{O}\left(N_{\rm X}N_{\rm Y}(N_{\rm R}\!+\!N_{\rm T})L\right)$. In order to achieve higher estimation accuracy, the increase in computational complexity is acceptable since $L$ is usually small for mmWave channels.

V Conclusions

In this paper, we have proposed an IR-based super-resolution channel estimation scheme for mmWave massive MIMO with hybrid precoding. Specifically, we have transformed the channel estimation problem to the optimization problem of a new objective function, which is the weighted summation of the sparsity and the data fitting error. The proposed scheme starts from the on-grid points in the angle domain, and iteratively moves them to the neighboring off-grid actual positions via gradient descent method. In addition, we have proposed an SVD-based preconditioning to reduce the computational complexity. Simulation results have confirmed that the proposed super-resolution channel estimation scheme can advance the state-of-art by estimating the off-grid AoAs/AoDs with much increased accuracy. Angle estimation is the key of channel estimation for mmWave massive MIMO. Estimating the AoAs/AoDs with higher resolution is a practical way to realize higher spectral efficiency. For future work, it would be interesting to study other super-resolution channel estimation schemes with reduced complexity. In addition, super-resolution channel estimation under high mobility is an important yet challenging topic to be investigated.

Appendix A Optimization of $S$ in (11) with regard to $\bf{z}$

For notational conciseness, we ignore the superscript $(i)$ of $S^{(i)}$ and ${\bf D}^{(i)}$ in (11), and use ${\bf A}_{\rm R}$, ${\bf A}_{\rm T}$ for ${\bf A}_{\rm R}\left({\bm{\theta}}_{\rm R}\right)$, ${\bf A}_{\rm T}\left({\bm{\theta}}_{\rm T}\right)$ respectively. Let ${{\bf{K}}_{p}}={\bf{W}}^{\!H\!}{{\bf{A}}_{r}}{\rm{diag}}\left({{\bf{A}}_{t}^{\rm{H}}{{\bf{x}}_{p}}}\right)$. In order to find the optimal $S\left({{\bf{z}},{{\bm{\theta}}_{\rm R}},{{\bm{\theta}}_{\rm T}}}\right)$ with regard to ${\bf{z}}$, we can expand the objective function $S$ as

[TABLE]

Then, we can obtain the partial derivative by

[TABLE]

By setting the derivative to zero, the minimum point ${\bf{z}}$ and the corresponding minimum value of $S\left({{\bf{z}},{{\bm{\theta}}_{\rm R}},{{\bm{\theta}}_{\rm T}}}\right)$ as the function of ${{\bm{\theta}}_{\rm R}}$ and ${{\bm{\theta}}_{\rm T}}$ can be obtained as

[TABLE]

[TABLE]

Appendix B Gradient of $S_{\rm opt}\left({{{\bm{\theta}}_{\rm R}},{{\bm{\theta}}_{\rm T}}}\right)$

Denote ${\bf v}\!=\!{\sum_{p=1}^{N_{\rm X}}{{\bf{K}}_{p}^{H}{{\bf{y}}_{p}}}}$, ${\bf A}\!=\!{\lambda^{\!-\!1}}{\bf{D}}+{\sum_{p=1}^{N_{\rm X}}{{\bf{K}}_{p}^{H}{{\bf{K}}_{p}}}}$, we have $S_{\rm opt}\!=\!-{\bf v}^{H}{\bf A}^{\!-1}{\bf v}+{\sum_{p=1}^{N_{\rm X}}{{\bf{y}}_{p}^{H}{{\bf{y}}_{p}}}}$. Take partial derivative with respect to ${\theta_{{\rm R},l}}$, we have

[TABLE]

where

[TABLE]