8954 IEEE TRANSACTIONS ON VEHICULAR TECHNOLOGY, VOL. 67, NO.

9, SEPTEMBER 2018

Super-Resolution Channel Estimation for MmWave codebooks. However, the channel estimation resolution is limited by
Massive MIMO With Hybrid Precoding the codebook size. [7] was able to achieve better angle estimation
by performing an amplitude comparison with respect to the auxiliary
Chen Hu , Linglong Dai , Senior Member, IEEE, Talha Mir, beam pair. On the other hand, by exploiting the angular channel spar-
Zhen Gao , Member, IEEE, and Jun Fang, Senior Member, IEEE sity, 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)
Abstract—Channel estimation is challenging for millimeter-wave mas-
sive MIMO with hybrid precoding, since the number of radio frequency
lie in discrete points in the angle domain (i.e., “on-grid” AoAs/AoDs),
chains is much smaller than that of antennas. Conventional compressive while the actual AoAs/AoDs are continuously distributed (i.e., “off-
sensing based channel estimation schemes suffer from severe resolution grid” AoAs/AoDs) in practice. The assumption of on-grid AoAs/AoDs
loss due to the channel angle quantization. To improve the channel esti- results in the power leakage problem, which severely degrades the
mation accuracy, we propose an iterative reweight-based superresolution channel estimation accuracy. To solve this resolution limitation caused
channel estimation scheme in this paper. By optimizing an objective func-
tion through the gradient descent method, the proposed scheme can itera- by the on-grid angle estimation, we propose an iterative reweight (IR)-
tively move the estimated angle of arrivals/departures towards the optimal based super-resolution channel estimation scheme to estimate the off-
solutions, and finally realize the superresolution channel estimation. In grid AoAs/AoDs.1
the optimization, a weight parameter is used to control the tradeoff be- Specifically, we iteratively optimize the estimates of AoAs/AoDs,
tween the sparsity and the data fitting error. In addition, a singular value
decomposition-based preconditioning is developed to reduce the compu-
to decrease the weighted summation of the sparsity and the data fit-
tational complexity of the proposed scheme. Simulation results verify the ting error. The weight controlling the tradeoff between the sparsity
better performance of the proposed scheme than conventional solutions. 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
Index Terms—Millimeter-wave (mmWave), massive MIMO, hybrid pre-
coding, angle of arrival (AoA), angle of departure (AoD), superresolution initial angle-domain grids towards the actual off-grid AoAs/AoDs, the
channel estimation. proposed scheme is able to achieve the super-resolution channel esti-
mation. In addition, we propose a singular value decomposition (SVD)-
based preconditioning method to reduce the computational complexity
I. INTRODUCTION of the proposed scheme, which is realized by reducing the number of
Millimeter-Wave (mmWave) massive MIMO has been recognized initial candidates of AoAs/AoDs in the IR procedure. Simulation results
as a promising technology for future 5G wireless communications [1]. show that the proposed IR-based super-resolution channel estimation
To reduce the hardware cost and power consumption, hybrid precoding can achieve better performance than conventional solutions.
has been proposed for practical mmWave massive MIMO systems, The contributions of this paper are the follows. We propose a novel
where hundreds of antennas are driven by a much smaller number of IR-based super-resolution channel estimation scheme for mmWave
radio frequency (RF) chains [2], [3]. The analog and digital co-design in massive MIMO with hybrid precoding. Comparing with the state-of-art
hybrid precoding requires accurate channel state information. However, schemes schemes such as those in [6]–[8], we can achieve super-
the digital baseband cannot directly access all antennas due to the small resolution channel estimation, which means substantially improved es-
number of RF chains, so it is difficult to accurately estimate the high- timation accuracy. Moreover, the proposed SVD based preconditioning
dimensional MIMO channel [4], [5]. significantly reduces the computational complexity of the IR procedure,
Several novel channel estimation schemes have been recently pro- and makes the method practical in mmWave channel estimation.
posed for mmWave massive MIMO with hybrid precoding [5]–[9]. Notation: In this paper, the boldface lower and upper-case symbols
Specifically, [5], [6] proposed the adaptive codebook-based channel denote vectors and matrices. (·)T , (·)H and (·)−1 denote the trans-
sounding scheme, where the transmitter and receiver search for the pose, the conjugate transpose, and the inverse of a matrix, respectively.
best beam pair by adjusting the predefined precoding and combining diag(x) is the diagonal matrix with the vector x on its diagonal. The
0 -norm, 2 -norm, and Frobenius norm are given by  · 0 ,  · 2 , and
Manuscript received October 28, 2017; revised March 18, 2018; accepted May
 · F , respectively.
27, 2018. Date of publication June 1, 2018; date of current version September 17,
2018. This work was supported in part by the National Natural Science Founda- II. SYSTEM MODEL
tion of China for Outstanding Young Scholars under Grant 61722109, in part by
the National Natural Science Foundation of China under Grant 61571270, and We consider a hybrid-precoding mmWave massive MIMO with ar-
in part by the Royal Academy of Engineering through the U.K.–China Industry bitrary array geometry. Let NT , NR , NTR F , and NRR F be the number of
Academia Partnership Programme Scheme under Grant UK-CIAPP\49. The
review of this paper was coordinated by Prof. S. Coleri Ergen. (Corresponding
transmit antennas, receive antennas, transmitter RF chains, and receiver
author: Linglong Dai.) RF chains, respectively. For practical mmWave massive MIMO with
C. Hu, L. Dai, and T. Mir are with the Department of Electronic Engi- hybrid precoding, the number of RF chains is much smaller than that
neering, Tsinghua University, Beijing 100084, China (e-mail:, huc16@mails. of antennas, i.e., NTR F < NT , NRR F < NR [1]–[3]. The system model
tsinghua.edu.cn; [email protected]; [email protected]). can be given by
Z. Gao is with the Advanced Research Institute for Multidisciplinary
Science, Beijing Institute of Technology, Beijing 100081, China (e-mail:, r = QH HPs + n, (1)
[email protected]).
where r ∈ C N R ×1 is the received signal, Q ∈ C N R ×N R is the hy-
RF RF

J. Fang is with the National Key Laboratory of Science and Technology on

Communications, University of Electronic Science and Technology of China, brid combining matrix, H ∈ C N R ×N T is the channel matrix, P ∈
Chengdu 611731, China (e-mail:,[email protected]).
Color versions of one or more of the figures in this paper are available online
at http://ieeexplore.ieee.org. 1 Simulation codes are provided to reproduce the results presented in this
Digital Object Identifier 10.1109/TVT.2018.2842724 paper: http://oa.ee.tsinghua.edu.cn/dailinglong/publications/publications.html

0018-9545 © 2018 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission.
See http://www.ieee.org/publications standards/publications/rights/index.html for more information.
IEEE TRANSACTIONS ON VEHICULAR TECHNOLOGY, VOL. 67, NO. 9, SEPTEMBER 2018 8955

RF
C N T ×N T is the hybrid precoding matrix, s ∈ C N T ×1 is the trans-
RF
of the channel matrix H, the sparse channel estimation problem can be
RF
mitted signal, and n ∈ C N R ×1 is the received noise after combining. formulated as

The channel model
min ẑ0 , s.t. Y − W H ĤX ≤ ε, (8)
ẑ , θ̂ R , θ̂ T F

L
   
H= zl aR φaRz,il , φeRle, l aH
T φaTz,il , φeTle, l (2) where ẑ0 is the number of non-zero elements of ẑ, which means the
l= 1 estimated number of paths L̂, Ĥ is the estimated channel matrix, and
ε is the error tolerance parameter [9].
is widely adopted in mmWave massive MIMO systems, and it is re-
garded almost unchanged within the channel coherence time for chan-
nel estimation [6]–[8], where L is the number of propagation paths, III. PROPOSED IR-BASED SUPER-RESOLUTION
L  min (NR , NT ), zl , φaRz,il (φeRle, l ) and φaTz,il (φeTle, l ) are the complex CHANNEL ESTIMATION
path gain, the azimuth (elevation) AoA and AoD of the l-th path, re-
A. Proposed Optimization Formulation
spectively. aR(φaRz,il , φeRle, l ) and aT(φaTz,il , φeTle, l ) are the steering vector
at the receiver and the steering vector at the transmitter, respectively. The main difficulty in solving (8) lies in the fact that the l0 -norm
These steering vectors depend on the array geometry. Ignoring the is not computationally efficient for finding the optimal solution. By
 for the typical N1 × N2 uniform
subscripts without loss of replacing the l0 -norm with a log-sum function [10], we have
 generality,
planar arrays (UPAs), a φal z i , φel le is given by [2]

Δ

L

log |zl |2 + δ , s.t. Y −WH ĤX ≤ ε, (9)
 i e le   j 2π d sin φa z isin φe l e/λ T min F (z) =
z,θR ,θT F
, · · ·, ej 2π (N1−1)d sin φ sin φ /λ
azi e le
a φa z,φ = 1, e l= 0

 T where δ > 0 ensures that the logarithmic function is well-defined [10],

⊗ 1, ej 2π d c o s φ /λ , · · ·, ej 2π (N2 −1)d c o s φ /λ ,
e le e le
(3) Ĥ is determined by the parameters z, θ R and θ T defined in (5). By
adding a regularization parameter λ > 0, we can further formulate the
where d is the antenna spacing, λ is the wavelength, ⊗ denotes the problem (9) as a unconstrained optimization problem:
Kronecker product. For uniform linear arrays (ULAs), the steering 
L

Δ
vector is only determined by one angle [6] min G (z, θ R , θ T ) = log |zl |2 +δ
z,θR ,θT
l= 1
 j 2π d sin φ / λ

j 2π (N −1)d sin φ / λ T 2
a (φ) = 1, e , · · ·, e . (4)
+ λ Y − WH ĤX . (10)
F
azi Δ Moreover, by using an iterative surrogate function instead of the log-
By defining the normalized spacial angles by θ =
Δ
d sin φa z i sin φe le /λ and θ e le = d cos φe le /λ, the mmWave channel sum function, the minimization of G (z, θ R , θ T ) is equivalent to the
matrix H in (2) can be also written as minimization of the surrogate function [10]:
2

min S (i )(z, θ R , θ T ) = λ−1 zH D(i ) z + Y − W H ĤX , (11)
Δ

H = AR (θ R ) diag (z) AH
T (θ T ) , (5) z,θR ,θT F

where D(i ) is defined as

where z = [z1 , z2 , . . . , zL ]T , θ R = [θR, azi e le azi e le
1 , θR, 1 , θR, 2 , θR, 2 , . . ., ⎛ ⎞
 T  
e le T
azi
θR, e le
L , θR, L  , θ T = θT, azi e le azi e le azi
1 , θT, 1 , θT, 2 , θT, 2 , . . ., θT, L , θT, L , Δ ⎜ 1 1 1 ⎟
D(i ) = diag ⎝     ···   ⎠, (12)
AR (θ R ) = aR (θRa z,i1 , θRe le, 1 ) aR (θRa z,i2 , θRe le, 2 ) · · · aR (θRa z,iL , θRe le, L ) ,  (i ) 2  (i ) 2  (i ) 2
  ẑ1  +δ ẑ2  +δ ẑL  +δ
AT (θ T ) = aT (θTa z, i1 , θTe le, 1 ) aT (θTa z, i2 , θTe le, 2 ) · · · aT (θTa z, iL , θTe le, L ) .
N T ×1
Denote x = Ps ∈ C , where the i-th element of x is the trans- where ẑ(i ) is the estimate of z at the i-th iteration.
mitted signal at the i-th transmit antenna. Suppose that the transmit- Then, as proved in Appendix A, we can optimize (11) with regard
ter sends NX (NX < NT ) different pilot sequences x1 , x2 , . . . , xN X . to the path gains z, to find the optimal point of ẑ and the corresponding
Since the number of RF chains is smaller than the required dimen- optimal value of S (i ) as follows:
sion of received pilot sequence, for each transmit pilot sequence
Δ
xp (1 ≤ p ≤ NX ), we use M time slots to obtain an NY -dimensional (i )
zo p t (θ R , θ T ) = arg min S (i ) (z, θ R , θ T )
z
received pilot sequence yp , where NY = M NRR F . Thus, the training
 −1 N 
overhead is T = M NX . In the m-th time slot, we use the combining 
NX
 X
−1
matrix W m to obtain an NRR F -dimensional received pilot sequence = λ D (i )
+ KH
p Kp KH
p yp , (13)
p=1 p=1
yp , m = W m
H
Hxp + np , m . (6) Δ
(i )
So p t (θ R , θ T ) = min S (i ) (z, θ R , θ T )
z
By collecting the received pilots in the M time slots, we
N   −1
have yp = W H Hxp + np , where yp = [ypT, 1 , ypT, 2 , . . . , ypT, M ]T ∈  X 
NX
−1 (i )
C N Y ×1 , W = [W 1 , W 2 , . . . , W M ] ∈ C N R ×N Y , np ∈ C N Y ×1 is the =− yp K p · λ D +
H
Kp Kp
H

noise. By defining Y = [y1 , y2 , . . . , yN X ], X = [x1 , x2 , . . . , xN X ], p=1 p=1

N = [n1 , n2 , ldots, nN X ], we have N 

 X 
NX
· KH
p yp + ypH yp , (14)
Y = W H HX + N. (7) p=1 p=1
 
The estimation of the channel matrix H in (7) is equivalent to the where Kp = W H AR diag AH T xp . After that, we only need to opti-
estimation of the number of paths, the normalized spacial angles (θ T , mize the normalized spatial angles θ R and θ T in (14), which will be
θ R ), and path gains z for all L paths. Due to the angle-domain sparsity discussed in the next section.
8956 IEEE TRANSACTIONS ON VEHICULAR TECHNOLOGY, VOL. 67, NO. 9, SEPTEMBER 2018

Algorithm 1: IR-based Super-Resolution Channel Algorithm 2: SVD-based preconditioning.

Estimation. Input: Noisy received signals Y, transmit pilot signals X,
Input: Noisy received signals Y, transmit pilot signals X, combining matrix W, and Nin it , the number of
combining matrix W, initial on-grid AoAs and paths to detect.
(0) (0)
AoDs θ̂ R , θ̂ T , pruning threshold zth and Output: Coarse AoAs/AoDs estimates of the Nin it paths.
termination threshold εth . 1: [U, Σ, V] = SVD(Y).
Output: Estimated AoAs/AoDs andpath gains of all paths. 2: Take the first Nin it columns, {u1 , u2 , · · · , uN i n i t } from

(0)
1: Initialize ẑ(0) = zo p t θ̂ R , θ̂ T according to (13).
(0) U, and {v1 , v2 , · · · , vN i n i t } from V, which are
correspondent to the Nin it largest singular values.
2: repeat
i = 1, 2, · · · , Nin it do
3: for

3: Update λ by (15). a z i(0) e le (0)

(i ) 4: θ̂R , i , θ̂R , i = arg max uH i W aR
H
4: Construct the function So p t (θ R , θ T ) by (14).
(θ Ra z i , θ Re l e )∈Ω R
(i + 1) (i + 1)  a z i e le 
5: Search for new angle estimates θ̂ R , θ̂ T by (17). θ , θR .

R 
6: Estimate the path gains ẑ(i + 1) according to (13). a z i(0) e le (0)
(i + 1)
5: θ̂T , i , θ̂T , i = arg max viH XH aT
7: Prune path l if ẑl <z . (θ Ta z i , θ Te l e )∈Ω T
th  a z i e le 
8: until L(i ) = L(i + 1) and z(i + 1) − z(i ) 2 < εth . θT , θT .
(la st) (la st) 6: end for
9: θ̂ R = θ̂ R , θ̂ T = θ̂ T , ẑ = ẑ(la st) .

B. IR-Based Super-Resolution Channel Estimation estimates are almost the same as the previous ones. With our proposed
IR-based super-resolution channel estimation scheme, the estimates of
In the previous section, we have already simplified the constrained
(θ R , θ T ) can be moved from the initial on-grid coarse estimates to its
optimization problem (8) to an unconstrained angle optimization prob-
actual off-grid positions, thus the super-resolution channel estimation
lem (14). To solve this reformulated problem, now we propose an
can be realized.
IR-based super-resolution channel estimation scheme as described in
It is worthy to point out that the sparsity level L is unknown in
Algorithm 1.
practice. In the proposed scheme, the sparsity level can be initialized
The objective function S (i ) (z, θ R , θ T ) is the weighted sum of
to be larger than the real channel sparsity. During the iteration process,
two parts: zH Dz controlling the sparsity of the estimation result and
the paths with too small path gains will be regarded as noise instead
Y − W H ĤXF denoting the residue. In addition, λ is the regular-
of real paths. Then, our algorithm prune these paths to make the result
ization parameter that controls the tradeoff between the sparsity and
sparser. By iteratively pruning these paths, the estimated sparsity level
the data fitting error.
will decrease to the real number of paths.
In the iterative reweighted method [10], λ is not fixed but updated in
The computational complexity in each iteration lies in calculating
each iteration. To be specific, if the previous iteration is poorly-fitted,
the gradient in  Step 5. The computational
 complexity to calculate the
we will choose a smaller λ to make the estimate sparser. On the other
gradient is O NX NY (NR +NT )L2 . As a result, the number of initial
hand, if the previous iteration returns a well-fitted estimate and leads
candidates L(0) is critical, and it should be as small as possible to make
to a small residue, our method will choose a larger λ to accelerate the
the computation affordable. The problem how to effectively select the
searching for the best-fitting estimate. In the proposed algorithm, λ is (0) (0)
updated by initial θ̂ R and θ̂ T before the iteration will be discussed in the next
λ = min(d/r (i ) , λm a x ), (15) section.

where d is a constant scaling factor, and λm a x is selected to make the

problem well-conditioned, r (i ) is the squared residue in the previous C. SVD-Based Preconditioning
step, i.e.,
In this section, we propose a singular value decomposition (SVD)-
(i ) (i ) based preconditioning as shown in Algorithm 2, to reduce the com-
r (i )
= Y − W H
AR (θ̂ R ) diag (ẑ(i ) )AH
T (θ̂ T )X2F . (16)
putational complexity of the IR procedure in the proposed IR-based
The updation of λ was discussed in [10] with more details. super-resolution channel estimation scheme. The proposed scheme can
The proposed algorithm starts iteration at the angle domain find the angle-domain grids nearest to the real AoAs/AoDs. Compar-
grids. In the i-th iteration, our task is to search for new ing to using all NR NT angle-domain grids as initial candidates, the
(i + 1) (i + 1) preconditioning can significantly reduce the computational complexity
estimates θ̂ R and θ̂ T in the neighborhood of the previous esti-
(i ) (i ) of the IR-based super-resolution channel estimation scheme.
mates θ̂ R and θ̂ T
to make the objective function S (i ) become smaller. to the matrix Y, we
This searching can be accomplished via gradient descent method:
Specifically, by applying SVD
  have Y =
UΣV H , where Σ = diag σ1 , σ2 , . . . , σm in (N X , N Y ) ∈ RN Y ×N X
(i + 1) (i )

(i ) (i )  whose diagonal entries σ1 ≥ σ2 ≥ · · · ≥ σm in (N X , N Y ) ≥ 0 are the sin-
(i )
θ̂ R = θ̂ R − η · ∇θ R So p t θ̂ R , θ̂ T , gular values of Y, and UH U = IN Y ×N Y , V H V = IN X ×N X . From
(i + 1) (i )

(i ) (i )  (17) (5) and (7), we have
(i )
θ̂ T = θ̂ T − η · ∇θ T So p t θ̂ R , θ̂ T ,
   H
where the gradients can be calculated according to Y = W H AR (θ R ) diag (z) XH AT (θ T ) + N. (18)
Appendix B, and η is the chosen step-length to make sure
(i ) (i + 1) (i + 1) (i ) (i ) (i )
So p t(θ̂ R , θ̂ T ) ≤ So p t(θ̂ R , θ̂ T ). The estimates become more As the noise is small, the largest L singular values and their cor-
and more accurate during the iterative searching, until the new responding singular vectors are approximately determined by the L
IEEE TRANSACTIONS ON VEHICULAR TECHNOLOGY, VOL. 67, NO. 9, SEPTEMBER 2018 8957

Fig. 1. NMSE performance comparison of different channel estimation Fig. 2. NMSE performance comparison of different channel estimation
schemes under NLoS channel. schemes under LoS channel.

paths, i.e., for i = 1, 2, · · · , L, we have

   
σi ≈ |zl i | WH aR θRa z,il i ,θRe le, l i 2 XH aT θTa z, il i ,θTe le, l i 2 ,
   
ui ≈ WH aR θRa z,il i , θRe le, l i / WH aR θTa z, il i , θTe le, l i 2 , (19)
   
vi ≈ XH aT θRa z,il i , θRe le, l i / XH aT θTa z, il i , θTe le, l i 2 ,

where ui and vi are the i-th column of U and V, respectively,

{l1 , l2 , · · · , lL } is a permutation of {1, 2, · · · , L}.
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 ΩR (ΩT ). Take UPA as an example. For an N1 ×
N2 receiver array, the set of grids can be defined by ΩR =
{(i/N1 , j/N2 )|i = 0, 1, . . . , N1 − 1; j = 0, 1, . . . , N2 − 1}. We can
similarly define ΩT for the transmitter.
In [10], the initial candidates of Algorithm 1 are set to be Fig. 3. Average spectral efficiency when different channel estimation schemes
allthe grids, i.e., L(0) = NRNT . The computational complexity is are used.
O NX NY (NR +NT )NR2 NT2 , which is unaffordable when NR and
NT are very large. Fortunately, with the proposed SVD-based pre-
conditioning shown in Algorithm 2, the coarse estimates will be dB in the LoS scenario. In both cases, the proposed scheme achieves
used as the initial candidates of Algorithm 1, i.e., L(0) = Nin it ≈ L. much better NMSE performance when SNR becomes large. Moreover,
Thus,
 the computational complexity after SVD preconditioning is we show the performance of the proposed scheme when UPA is consid-
O NX NY (NR +NT )L2 , which is much lower than directly applying ered. Both the transmitter and the receiver adopt 64-antenna UPA with
the scheme in [10]. 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
IV. SIMULATION RESULTS AND DISCUSSIONS achieve super-resolution channel estimation when UPA is used. Since
the estimation errors of both azimuth and elevation angles contribute
In this section, simulation results are provided to investigate the to the NMSE, under the same number of antennas and number of pilot
performance. We consider the mmWave massive MIMO system with overhead, the NMSE performance of UPA is higher than that of ULA.
hybrid precoding, where L = 3, d = λ/2, NR = NT = 64, NRR F = Fig. 3 compares the average spectral efficiency when different chan-
NTR F = 4 and NX = NY = 32. The path gains are assumed Gaussian, nel estimation schemes are used. The spectral efficiency is evaluated in
i.e., αl ∼ CN(0, σα ).j ωEach element of the transmitted pilots X satis-
2
the hybrid precoding system [2]. The case with ideal CSI was adopted
fies xi , j = ρ/NT e i , j , where ρ is the transmitted power, ωi , j is as the upper bound for performance comparison. It can be observed that
the random phase uniformly distributed in [0, 2π). The signal-to-noise the proposed super-resolution channel estimation is able to approach
2
ratio (SNR) is defined by SNR = ρσσ2α , where σn2 is the noise variance. this upper bound. This is because the angle resolution of the proposed
n
We consider the ULA geometry, so that the adaptive codebook-based scheme does not suffer from limited codebook size or angle quantiza-
channel estimation [6], the auxiliary beam pair based channel estima- tion. Thus, we can conclude that the proposed scheme can achieve the
tion [7], and the OMP-based channel estimation [8] can be adopted for super-resolution channel estimation.
performance comparison. There is a tradeoff between the channel estimation accuracy and the
Figs. 1 and 2 compares the normalized mean square error (NMSE) computation complexity. The proposed scheme with SVD precondi-
performance against SNR, under none-line-of-sight (NLoS) and tioning is able to achieve much higher channel estimation accuracy, but
line-of-sight (LoS) channels, respectively. The Rician K-factor is 20 has a higher computational complexity. The computational complexity
8958 IEEE TRANSACTIONS ON VEHICULAR TECHNOLOGY, VOL. 67, NO. 9, SEPTEMBER 2018

of the proposed super-resolution

 channel estimation
 scheme with θ T can be obtained as
SVD precondistioning is O NX NY (NR +NT )L2 . In comparison,  −1 N 
the computational complexity of the OMP-based channel estimation 
NX
 X

zo p t (θ R , θ T ) = λ−1 D+ KH
p Kp KH
p yp , (22)
[8] is O (NX NY (NR +NT )L). In order to achieve higher estimation
p=1 p=1
accuracy, the increase in computational complexity is acceptable since N H  −1
 X NX

L is usually small for mmWave channels. −1
So p t (θ R , θ T ) = − KH
p yp λ D+ KH
p Kp
p=1 p=1
N 
V. CONCLUSION  X 
NX
· KH
p yp + ypH yp . (23)
In this correspondence, we have proposed an IR-based super-
p=1 p=1
resolution channel estimation scheme for mmWave massive MIMO
with hybrid precoding. Specifically, we have transformed the channel
APPENDIX B
estimation problem to the optimization problem of a new objective
GRADIENT OF Sopt (θ R , θ T )
function, which is the weighted summation of the sparsity and the data
 X H −1
N X H
fitting error. The proposed scheme starts from the on-grid points in Denote v = N p = 1 K p yp , A = λ D + p = 1 K p K p , we have
the angle domain, and iteratively moves them to the neighboring off- −1
 NX
So p t = −v A v + p = 1 yp yp . Take partial derivative with respect
H H
grid actual positions via gradient descent method. In addition, we have to θR , l , we have
proposed an SVD-based preconditioning to reduce the computational
complexity. Simulation results have confirmed that the proposed super- ∂So p t ∂v H −1 ∂A−1 ∂v
=− A v−v H v−v H A−1
resolution channel estimation scheme can advance the state-of-art by ∂θR , l ∂θR , l ∂θR , l ∂θR , l
estimating the off-grid AoAs/AoDs with much increased accuracy. An- (24)
gle estimation is the key of channel estimation for mmWave massive ∂v H −1 ∂A −1 ∂v
=− A v+v H A−1 A v−v H A−1 ,
MIMO. Estimating the AoAs/AoDs with higher resolution is a practical ∂θR , l ∂θR , l ∂θR , l
way to realize higher spectral efficiency. For future work, it would be
where
interesting to study other super-resolution channel estimation schemes  
with reduced complexity. In addition, super-resolution channel estima- ∂v NX
∂Kp H
∂A  X N
∂KHp ∂Kp
tion under high mobility is an important yet challenging topic to be = yp , = Kp +KH
p ,
∂θR , l p = 1 ∂θR , l ∂θR , l p = 1 ∂θR , l ∂θR , l
investigated. (25)
∂Kp  
∂ a R (θ R , l ) H
= 0 · · · 0 WH ∂ θR ,l
aT (θT , l ) xp 0 ··· 0 .
APPENDIX A ∂θR , l
OPTIMIZATION OF S IN (11) WITH REGARD TO z
For notational conciseness, we ignore the superscript (i) of S (i )
REFERENCES
and D(i ) in (11), and use AR, AT for  AR (θ R ), AT (θ T ) respec-
tively. Let Kp = WH Ar diag AH t xp . In order to find the optimal
[1] S. Mumtaz, J. Rodriquez, and L. Dai, MmWave Massive MIMO: A
S (z, θ R , θ T ) with regard to z, we can expand the objective function Paradigm for 5G. New York, NY, USA: Academic, 2016.
[2] O. E. Ayach, S. Rajagopal, S. Abu-Surra, Z. Pi, and R. W. Heath, “Spa-
S as tially sparse precoding in millimeter wave MIMO systems,” IEEE Trans.
Wireless Commun., vol. 13, no. 3, pp. 1499–1513, Mar. 2014.
[3] X. Gao, L. Dai, S. Han, C.-L. I, and R. Heath, “Energy-efficient hy-
S (z, θ R , θ T ) = λ−1 zH Dz brid analog and digital precoding for mmWave MIMO systems with
large antenna arrays,” IEEE J. Sel. Areas Commun., vol. 34, no. 4,

NX
2
+ yp − W H AR diag (z) AH
T xp 2
pp. 998–1009, Apr. 2016.
[4] L. Dai, Z. Wang, and Z. Yang, “Spectrally efficient time-frequency train-
p=1 ing OFDM for mobile large-scale MIMO systems,” IEEE J. Sel. Areas
Commun., vol. 31, no. 2, pp. 251–263, Feb. 2013.

NX
[5] S. Hur, T. Kim, D. J. Love, J. V. Krogmeier, T. A. Thomas, and A. Ghosh,
= λ−1 zH Dz + (yp − Kp z)H (yp − Kp z) “Millimeter wave beamforming for wireless backhaul and access in small
p=1 cell networks,” IEEE Trans. Commun., vol. 61, no. 10, pp. 4391–4403,
  N  N  Oct. 2013.

NX
 X  X
[6] A. Alkhateeb, O. E. Ayach, G. Leus, and R. W. Heath, “Channel estimation
−1
=z H
λ D+ KH
p Kp z−z H
KH
p yp − yp K p z
H
and hybrid precoding for millimeter wave cellular systems,” IEEE J. Sel.
p=1 p=1 p=1 Topics Signal Process., vol. 8, no. 5, pp. 831–846, Oct. 2014.
[7] D. Zhu, J. Choi, and R. W. Heath, “Auxiliary beam pair enabled AoD
 NX
and AoA estimation in closed-loop large-scale millimeter-wave MIMO
+ yH
p yp , (20) systems,” IEEE Trans. Wireless Commun., vol. 16, no. 7, pp. 4770–4785,
p=1 May 2017.
[8] J. Lee, G. T. Gil, and Y. H. Lee, “Channel estimation via orthogonal
matching pursuit for hybrid MIMO systems in millimeter wave commu-
Then, we can obtain the partial derivative by nications,” IEEE Trans. Commun., vol. 64, no. 6, pp. 2370–2386, Jun.
2016.
  N  [9] Z. Marzi, D. Ramasamy, and U. Madhow, “Compressive channel estima-
∂S (z, θ R , θ T ) 
NX
 X
−1 tion and tracking for large arrays in mm-wave picocells,” IEEE J. Sel.
= z λ D+
H
Kp Kp −
H
yp K p .
H
∂z Topics Signal Process., vol. 10, no. 3, pp. 514–527, Apr. 2016.
p=1 p=1
[10] J. Fang, F. Wang, Y. Shen, H. Li, and R. S. Blum, “Super-resolution
(21) compressed sensing for line spectral estimation: An iterative reweighted
By setting the derivative to zero, the minimum point z and the corre- approach,” IEEE Trans. Signal Process., vol. 64, no. 18, pp. 4649–4662,
sponding minimum value of S (z, θ R , θ T ) as the function of θ R and Sep. 2016.