See discussions, stats, and author profiles for this publication at: [Link]

net/publication/316538670

Frequency-domain Compressive Channel Estimation for Frequency-Selective

Hybrid mmWave MIMO Systems

Article  in  IEEE Transactions on Wireless Communications · April 2017

DOI: 10.1109/TWC.2018.2804943

CITATIONS READS
90 137

4 authors:

Javier Rodríguez Nuria González Prelcic

University of Vigo University of Vigo
23 PUBLICATIONS   222 CITATIONS    129 PUBLICATIONS   4,523 CITATIONS   

SEE PROFILE SEE PROFILE

Kiran Venugopal Robert W. Heath

Qualcomm Flarion Technologies, Inc. United States, Bridgewater North Carolina State University
23 PUBLICATIONS   748 CITATIONS    922 PUBLICATIONS   51,437 CITATIONS   

SEE PROFILE SEE PROFILE

Some of the authors of this publication are also working on these related projects:

Limited Feedback in Mobile Communications View project

EIC IEEE Signal Processing Magazine View project

All content following this page was uploaded by Nuria González Prelcic on 15 May 2017.

The user has requested enhancement of the downloaded file.

 1

Frequency-domain Compressive Channel

Estimation for Frequency-selective Hybrid
mmWave MIMO Systems
Javier Rodrı́guez-Fernández† , Nuria González-Prelcic† , Kiran Venugopal‡ , and
arXiv:1704.08572v1 [[Link]] 27 Apr 2017

Robert W. Heath Jr.‡

†
Universidade de Vigo, Email: {jrodriguez,nuria}@[Link]
‡
The University of Texas at Austin, Email: {kiranv,rheath}@[Link]

Abstract

Channel estimation is useful in millimeter wave (mmWave) MIMO communication systems. Chan-
nel state information allows optimized designs of precoders and combiners under different metrics such
as mutual information or signal-to-interference-noise (SINR) ratio. At mmWave, MIMO precoders and
combiners are usually hybrid, since this architecture provides a means to trade-off power consumption
and achievable rate. Channel estimation is challenging when using these architectures, however, since
there is no direct access to the outputs of the different antenna elements in the array. The MIMO channel
can only be observed through the analog combining network, which acts as a compression stage of the
received signal. Most of prior work on channel estimation for hybrid architectures assumes a frequency-
flat mmWave channel model. In this paper, we consider a frequency-selective mmWave channel and
propose compressed-sensing-based strategies to estimate the channel in the frequency domain. We
evaluate different algorithms and compute their complexity to expose trade-offs in complexity-overhead-
performance as compared to those of previous approaches.

This work was partially funded by the Agencia Estatal de Investigacin (Spain) and the European Regional Development
Fund (ERDF) under project MYRADA (TEC2016-75103-C2-2-R), the U.S. Department of Transportation through the Data-
Supported Transportation Operations and Planning (D-STOP) Tier 1 University Transportation Center, by the Texas Department
of Transportation under Project 0-6877 entitled Communications and Radar-Supported Transportation Operations and Planning
(CAR-STOP) and by the National Science Foundation under Grant NSF-CCF-1319556 and NSF-CCF-1527079.
 2

Index Terms

Wideband channel estimation; millimeter wave MIMO; hybrid architecture.

I. I NTRODUCTION

MIMO architectures with large arrays are a key ingredient of mmWave communication sys-
tems, providing gigabit-per-second data rates [1]. Hybrid MIMO structures have been proposed
to operate at mmWave because the cost and power consumption of an all-digital achitecture
is prohibitive at these frequencies [2]. Optimally configuring the digital and analog precoders
and combiners requires channel knowledge when the design goal is maximizing metrics such as
the achievable rate or the SINR. Acquiring the mmWave channel is challenging with a hybrid
architecture, however, because the channel is seen through the analog combining network, the
SNR is low before beamforming, and the size of the channel matrices is large [3].
A significant number of papers have proposed solutions to the problem of channel estimation
at mmWave with hybrid architectures, but under a frequency-flat channel model [4]–[12]. These
strategies exploit the spatially sparse structure in the mmWave MIMO channel, formulating the
estimation of the channel as a sparse recovery problem. The support of the estimated sparse
vector identifies the pairs of direction-of-arrival/direction-of-departure (DoA/DoD) for each path
in the mmWave channel, while the amplitudes of the non-zero coefficients provide the channel
gains for each path. Compressive estimations lead to a reduction in the channel training length
when compared to conventional approaches such as those based on least squares (LS) estimation
[12].
Recently, some approaches for channel estimation in frequency-selective mmWave channels
have been proposed. In a recent paper [13], we designed a time-domain approach to estimate
the wideband mmWave channel assuming a hybrid MIMO architecture. This algorithm exploits
the sparsity of the wideband millimeter wave channel in both the angular and delay domains.
The sparse formulation of the problem in [13] includes the effect of non-integer sampling of the
transmit pulse shaping filter, with the subsequent leakage effect and increase of sparsity level
in the channel matrix. The main limitation of [13] is that the algorithm can be applied only to
single carrier systems. A frequency-domain strategy to estimate frequency-selective mmWave
channels was also proposed in [14]. A sparse reconstruction problem was formulated there to
 3

estimate the channel independently for every subcarrier, without exploiting spatial congruence
between subbands. Another approach in the frequency domain was designed in [15], but only
exploting the information from a reduced number of subcarriers. A different algorithm operating
in the frequency domain for a MIMO-OFDM system was proposed in [16]. Exploiting the fact
that spatial propagation characteristics do not change significantly within the system bandwidth,
[16] assumed spatially common sparsity between the channels corresponding to the different
subcarriers. A structured sparse recovery algorithm was then considered to reconstruct the
channels in the frequency domain. Thus, [16] is an interesting initial solution to the problem,
but has several limitations when applied to a mmWave communication system:
1) The effect of sampling the pulse shaping filter delayed by a non integer factor was not
considered in the channel model for a given delay tap. As shown in [18], not accounting
for this effect leads to virtual MIMO matrices with an artificially enhanced sparsity.
2) The algorithm was evaluated only for medium and high SNR regimes (larger than 10 dBs),
very unlikely at mmWave, where the expected SNR is below 0 dB.
3) The reconstruction algorithm provides accurate results when Gaussian measurement matri-
ces are employed; to generate the Gaussian matrices, unquantized phases were considered
in the training precoders, which is unrealistic in a practical implementation of a mmWave
system based on a hybrid architecture.
4) The algorithm is based on the strong assumption that the SNR is perfectly known at the
receiver before explicit channel estimation.
Another algorithm exploiting common sparsity in the frequency domain at mmWave was
proposed in [17]. Unlike the algorithm proposed in [16], the algorithm in [17] is proposed to
estimate mmWave wideband MU-MIMO channels. Besides the limitations 1)-3) described above,
which also hold in this case, this algorithm exhibits another problem that makes it less feasible
to be applied in a real mmWave communication system. A Line-of-Sight (LoS) Rician channel
model with Kfactor = 20 dB is considered, which is only applicable when there is a strong LoS
path. Owing to this artifact of the channel model, the algorithm in [17] estimates a single path
for each user, such that the task of channel estimation can not be successfully accomplished.
This paper proposes two novel frequency-domain approaches for the estimation of frequency-
selective mmWave MIMO channels. These approaches overcome the limitations of prior work and
 4

provide different trade-offs between complexity and achievable rate for a fixed training length.
As in recent work on hybrid architectures for frequency-selective mmWave channels [16], [19],
we also consider a MIMO-OFDM communications system. Similar to [14], we introduce zero-
padding (ZP) as a means to avoid loss and/or distortion of training data during reconfiguration
of RF circuitry. A geometric channel is considered to model the different scattering clusters as
in [13], [17], [18], including the band-limiting property in the overall channel response.
Our two proposed approaches exploit the spatially common sparsity within the system band-
width. The first algorithm aims at exploiting the information on the support coming from every
subcarrier in the MIMO-OFDM system and provides the best performance. In contrast, the second
algorithm uses less information to estimate the different frequency-domain subchannels, thereby
managing to significantly reduce the computational complexity. We show that both strategies are
asymptotically efficient, since they both attain the Cramér-Rao Lower Bound (CRLB). Further,
we show that asymptotic efficiency can be achieved without using frequency-selective baseband
precoders and combiners during the training stage, thereby reducing computational complexity.
Simulation results in the low SNR regime show that the two proposed algorithms significantly
outperform the approach in the frequency domain developed in [14]. Comparisons with the
algorithms proposed in [16] and [17] are also provided to show their performance in terms
of estimation error in the SNR regime that mmWave systems are expected to work. The two
proposed channel estimation algorithms provide a good trade-off between performance and
overhead. Results show that using a reasonably small training length, approximately in the
range of 60 − 100 frames, leads to low estimation errors. The computational complexity of the
proposed algorithms and previous strategies is also analyzed to compare the trade-offs between
performance-complexity provided by the different algorithms. Finally, we also show that it is
not necessary to exploit the information on the support coming from every OFDM subcarrier to
estimate the different mmWave subchannels. Yet, a reduced number of subcarriers is enough to
asymptotically attain the CRLB.
The paper is organized as follows: Section II introduces the system model that is used
throughout this work. Section III proposes two frequency-domain compressive channel estimation
approaches, including the computation of the CRLB and suitable estimators derived from it in
order to obtain the different frequency-domain channel matrices for each subcarrier. Thereafter,
Section IV provides the main simulation results for the two proposed algorithms, the OMP-
 5

.   .  
.   .  
.   .  

RF   .  
.  
.  
.  
.  
.  
RF  
chain   .. .. chain  
. .
Baseband   H Baseband  
.. Precoder   .. .  .   .  .   .. Combiner   ..
Ns . . Lt .   Nt Nr .   . Lr . Ns
FBB[k] WBB[k]
RF   RF  
chain   .. .. chain  
. .
.   .  
.   .  
.   .  

RF   RF  
Precoder   Combiner  
FRF WRF

Fig. 1. Illustration of the structure of a hybrid MIMO architecture, which include analog and digital precoders and combiners.

based compressive approach proposed in [13] and the SSAMP and DGMP algorithms proposed
in [16] and [17], respectively. Finally, Section V collects the main conclusions derived from the
simulation results and also describes future work to be conducted.
Notation: We use the following notation throughout this paper: bold uppercase A is used to
denote matrices, bold lowercase a denotes a column vector and non-bold lowercase a denotes a
scalar value. We use A to denote a set. Further, kAkF is the Frobenius norm and A∗ , A, AT and
A† denote the conjugate transpose, conjugate, transpose and Moore-Penrose pseudoinverse of a
matrix A, respectively. The (i, j)-th entry of a matrix A is denoted using {A}i,j . Similarly, the
i-th entry of a column vector a is denoted as {a}i . The identity matrix of order N is denoted as
IN . If A and B are two matrices, A◦B is the Khatri-Rao product of A and B and A⊗B is their
Kronecker product. We use N (µ, C) to denote a circularly symmetric complex Gaussian random
vector with mean µ and covariance matrix C. We use E to denote expectation. Discrete-time
signals are represented as x[n], while frequency-domain signals are denoted using x[k].

II. S YSTEM MODEL

We consider an OFDM based mmWave MIMO link employing K subcarriers to send Ns data
streams using a transmitter with Nt antennas and a receiver with Nr antennas. The system is
based on a hybrid MIMO architecture as shown in Fig. 1, with Lt and Lr RF chains at the
transmitter and receiver sides. For a general exposition, a frequency-selective hybrid precoder is
used, with F[k] = FRF FBB [k] ∈ CNt ×Ns , k = 0, . . . , K−1, where FRF is the analog precoder and
 6

FBB [k] the digital one. Note that the analog precoder is frequency-flat, while the digital precoder
is different for every subcarrier. The RF precoder and combiner are implemented using a fully
connected network of phase shifters, as described in [12]. The symbol blocks are transformed
into the time domain using Lr parallel K-point IFFTs. As in [14], [20], we consider Zero-
Padding (ZP) to both suppress Inter Symbol Interference (ISI) and account for the RF circuitry
reconfiguration time. The discrete-time complex baseband signal at subcarrier k can be written
as
x[k] = FRF FBB [k]s[k], (1)

where the transmitted symbol sequence at subcarrier k of size Ns × 1 is denoted as s[k].

The MIMO channel between the transmitter and the receiver is assumed to be frequency-
selective, having a delay tap length Nc in the time domain. The d-th delay tap of the channel
is represented by a Nr × Nt matrix denoted as Hd , d = 0, 1, ..., Nc − 1, which, assuming a
geometric channel model [18], can be written as
s
L
Nt Nr X
Hd = α` prc (dTs − τ` )aR (φ` )a∗T (θ` ), (2)
LρL `=1
where ρL denotes the path loss between the transmitter and the receiver, L denotes the number
of paths, prc (τ ) is a filter that includes the effects of pulse-shaping and other lowpass filtering
evaluated at τ , α` ∈ C is the complex gain of the `th path, τ` ∈ R is the delay of the `th
path, φ` ∈ [0, 2π) and θ` ∈ [0, 2π) are the angles-of-arrival and departure (AoA/AoD), of the
`th path, and aR (φ` ) ∈ CNr ×1 and aT (θ` ) ∈ CNt ×1 are the array steering vectors for the receive
and transmit antennas. Each one of these matrices can be written in a more compact way as

Hd = AR ∆d A∗T , (3)

where ∆d ∈ CL×L is diagonal with non-zero complex entries, and AR ∈ CNr ×L and AT ∈ CNt ×L
contain the receive and transmit array steering vectors aR (φ` ) and aT (θ` ), respectively. The
channel Hd can be approximated using the extended virtual channel model defined in [2] as

Hd ≈ ÃR ∆vd Ã∗T , (4)

where ∆vd ∈ CGr ×Gt is a sparse matrix which contains the path gains of the quantized spatial
frequencies in the non-zero elements. The dictionary matrices ÃT and ÃR contain the transmitter
and receiver array response vectors evaluated on a grid of size Gr for the AoA and a grid of size
 7

Gt for the AoD. Due to the few scattering clusters in mmWave channels, the sparse assumption
for ∆vd ∈ CGr ×Gt is commonly accepted.
Finally, the channel at subcarrier k can be written in terms of the different delay taps as
Nc −1
2πk
X
H[k] = Hd e−j K
d
. (5)
d=0

It is also useful to write this matrix in terms of the sparse matrices ∆vd and the dictionaries
 NXc −1 
v −j 2πk
H[k] ≈ ÃR ∆d e K d
Ã∗T = ÃR ∆[k]Ã∗T (6)
d=0

to help expose the sparse structure later.

Assuming that the receiver applies a hybrid combiner W[k] = WRF WBB [k] ∈ CNr ×Ns , the
received signal at subcarrier k can be written as
∗ ∗
y[k] = WBB [k]WRF H[k]FRF FBB [k]s[k]
(7)
∗ ∗
+WBB [k]WRF n[k],
where n[k] ∼ N (0, σ 2 I) is the circularly symmetric complex Gaussian distributed additive noise
vector. The receive signal model in (7) corresponds to the data transmission phase. As we will
see in Section III, during the channel acquisition phase, we will consider analog-only training
precoders and combiners to reduce complexity.

III. C OMPRESSIVE C HANNEL E STIMATION IN THE F REQUENCY D OMAIN

In this section, we formulate a compressed sensing problem to estimate the vectorized sparse
channel vector in the frequency domain. We also propose two algorithms to solve this problem
that leverage the common support between the channel matrices for every subcarrier, providing
different trade-offs complexity-performance. The first algorithm leverages the common support
between the K different subchannels providing a very good performance, while the second one
only exploits information from a reduced number of subcarriers, thereby keeping computational
complexity at a lower level.

A. Problem formulation

We assume that Lt RF chains are used at the transmitter. During the training phase, for the
(m) (m)
m-th frame we use a training precoder Ftr and a training combining matrix Wtr . This means
 8

that during the training phase, only analog precoders and combiners are considered to keep the
complexity of the sparse recovery algorithm low. We assume that the transmitted symbols satisfy
P
E{s(m) [k]s(m)∗ [k]} = I ,
Ns Ns
with P the total transmitted power and Ns = Lt . Furthermore, each
(m) (m)
entry in Ftr , Wtr is normalized to have squared-modulus Nt−1 and Nr−1 , respectively. Then,
the received samples in the frequency domain for the m-th training frame can be written as
(m) ∗ (m)
y(m) [k] = Wtr H[k]Ftr s(m) [k] + n(m)
c [k], (8)

where H[k] ∈ CNr ×Nt is the frequency-domain MIMO channel response at the k-th subcarrier
(m) (m) (m) ∗ (m)
and nc [k] ∈ CLr ×1 , nc [k] = Wtr n [k], is the frequency-domain combined noise vector
received at the k-th subcarrier. The average received SNR per transmitted symbol is given by
P
SNR = ρ` σ 2
. We assume that the channel coherence time is larger than the frame duration and
that the same channel can be considered for several consecutive frames.
Using the result vec{AXC} = (CT ⊗ A) vec{X}, the vectorized received signal is
(m) T (m) ∗
vec{y(m) [k]} = (s(m)T [k]Ftr ⊗ Wtr ) vec{H[k]} + n(m)
c [k]. (9)

Taking into account the expression in (6), the vectorized channel matrix can be written as
¯ ⊗ Ã ) vec{∆[k]}. Therefore, if we define the measurement matrix Φ(m) [k] ∈
vec{H[k]} = (Ã T R

CLr ×Nt Nr as
(m) T (m) ∗
Φ(m) [k] = (s(m)T [k]Ftr ⊗ Wtr ) (10)

and the dictionary Ψ ∈ CNt Nr ×Gt Gr as

¯ ⊗ Ã ),
Ψ = (Ã (11)
T R

(9) can be rewritten as

vec{y(m) [k]} = Φ(m) [k]Ψhv [k] + n(m)

c [k], (12)

where hv [k] = vec{∆[k]} ∈ CGr Gt ×1 is the sparse vector containing the complex channel gains.
To have enough measurements and accurately reconstruct the sparse vector hv [k], it is necessary
to use several training frames, especially in the very-low SNR regime. If the transmitter and
receiver communicate during M training steps using different pseudorandomly built precoders
 9

and combiners, (12) can be extended to

     
(1)
y(1) [k] Φ(1) [k] nc [k]
     
 (2)  (2)
 y [k]   Φ(2) [k]  nc [k]
   
 
=  Ψhv [k] +  . (13)
.. .. ..

     

 .  
  . 


 . 

(M )
y(M ) [k] Φ(M ) [k] nc [k]
| {z } | {z } | {z }
y[k] Φ[k] nc [k]

Finally, the vector hv [k] can be found by solving the sparse reconstruction problem

min khv [k]k1 subject to ky[k] − Φ[k]Ψhv [k]k22 < , (14)

where  is a tunable parameter defining the maximum error between the measurement and the
received signal assuming the reconstructed channel between the transmitter and the receiver.
Since the sparsity (number of channel paths) is usually unknown, this parameter can be set to
the noise variance [13].
There are a great variety of algorithms to solve (14). We could use, for example, Orthogonal
Matching Pursuit (OMP) to find the sparsest approximation of the vectors containing the channel
gains. Since the vector hv [k] depends on the frequency bin, it would be necessary to run the
algorithm as many times as the number of subcarriers at which the MIMO channel response
is to be estimated. In the next subsections we consider an additional assumption to solve this
problem, which avoids the need to run K OMP algorithms in parallel as proposed in [14].
The formulation in (8)-(14) is provided as a general framework to leverage the sparse structure
of the frequency-selective mmWave channel. This formulation, however, leads to subcarrier-
dependent measurement matrices, resulting in high computational complexity for channel esti-
mation. Using more than a single RF chain at the transmitter offers more degrees of freedom
to design the measurement matrices and, consequently, obtain better estimation performance.
Nonetheless, exploiting these degrees of freedom results in additional computational complexity
with regard to the case of using a single RF chain at the transmitter. Furthermore, the difference
in performance between using Lt = 1 and Lt > 1 was observed to be less than 0.5 dB when
solving the problem with the algorithm proposed in this section, which comes from the fact that a
reasonable number of RF chains must be used at a given transceiver to keep power consumption
low. For this reason, our focus is on designing algorithms exhibiting reasonable performance
and low computational complexity.
 10

For this reason, we consider the special case of Lt = 1 to present our algorithms with reduced
complexity. With this choice, the frequency-selective transmitted symbols are scalar and its effect
can be easily inverted at the receiver (i.e., by multiplying by (s(m) [k])−1 or multiplying by s(m) [k]
if the transmitted symbols belong to a QPSK constellation with energy-normalized symbols).
As shown in Section IV, the use of a single frequency-flat measurement matrix is enough to
asymptotically attain the CRLB with on-grid channel parameters. Therefore, a single RF chain
and analog-only precoders are used during training, thereby managing a reasonable trade-off
between performance and computational complexity. Nonetheless, the algorithms we propose
in the next subsections can be easily extended to work with subcarrier-dependent measurement
matrices.
By compensating for s(m) [k], for a given training step m, the measurement matrix Φ(m) [k]
(m) T (m)∗
reduces to Φ(m) = ftr ⊗ Wtr . At this point, it is important to highlight that: 1) Analog-only
training precoders and combiners are used to keep computational complexity low, and 2) A single
data stream is transmitted to guarantee that the same measurement matrix is enough to estimate
the different frequency-domain subchannels. In fact, these choices enable both online and offline
complexity reduction when the noise statistics are used to estimate the MIMO channel at the
different subcarriers.
The matrices ∆[k] exhibit an interesting property that can be exploited when solving the
compressed channel estimation problems defined in (14). Let us define the Gt Gr × 1 vectorized
virtual channel matrix for a given delay tap as

hvd , vec{∆vd }. (15)

Let us denote the supports of the virtual channel matrices ∆vd as T0 , T1 , . . . , TNc −1 , d = 0, . . . , Nc −
P c −1 v −j 2πk d
1. Then, since hv [k] = vec{∆[k]}, with ∆[k] = N d=0 ∆d e
N , k = 0, . . . , K − 1, it is clear
that
c −1
N[
v
supp{h [k]} = supp{hvd } k = 0, . . . , K − 1, (16)
d=0
where the union of the supports of the time-domain virtual channel matrices comes from the
additive nature of the Fourier transform. The sparse assumption on the vectorized channel matrix
for a given delay tap hvd is commonly accepted, since in mmWave channels L << Gr Gt . The
vectorized channel matrix hv [k] will have, in the worst case, Nc L non-zero coefficients. Typical
values for Nc in mmWave channels are usually lower than 64 symbols (for example IEEE
 11

802.11ad has been designed to work robustly for a maximum of 64 delay taps in the channel),
while the number of measured paths usually satisfies L < 30 for outdoor and indoor scenarios
[21]. From these values, using dictionaries of size Gr ≥ 64 Gt ≥ 64, allows us to assume a
spatially sparse structure for hv [k] as well. Furthermore, this sparse structure is the same for all
k, since from (6), ÃR and ÃT do not depend on k, which means that the AoAs and AoDs do
not change with frequency in the transmission bandwidth.

B. Simultaneous weighted estimation exploiting common support between subcarriers (SW-OMP)

To develop a channel estimation algorithm that leverages the sparse nature and the common
support property for all hv [k], we propose to modify the S-OMP algorithm proposed in [22].
For a given iteration, this algorithm aims at finding a new index of the support such that its
likelihood is larger than if the signals for the different subbands were individually processed. In
this way, the K different signals contribute to decide which is the most likely index belonging
to the support in a cooperative fashion. The S-OMP algorithm in [22] computes the gains of
the sparse vector using a LS approach once the support is obtained, assuming that the noise
covariance matrix is the identity matrix IM Lr . In this paper, we propose two modifications of
the S-OMP for the estimation of the support and the channel gains. The first one accounts for
the correlated nature of the noise at the output of the RF combiner, and leads to a more precise
estimation of the support in this particular application. The second one consists of a different
approach to compute the channel gains. In particular, we develop the minimum variance unbiased
estimator for the channel gains, which is shown to be a weighted LS estimator. We show that
this estimator attains the CRLB.
1) Support computation with correlated noise: Before explicit estimation of the channel gains,
it is necessary to compute the atom, i.e., vector in the measurement matrix, which yields the
largest sum-projection onto the received signals, provided that the support of the different sparse
vectors is the same. The S-OMP algorithm is based on the assumption that the perturbation (noise)
covariance matrix is diagonal, such that no correlation between the different noise components
is present. The projection vector c[k] ∈ CM Lr is defined as

c[k] = Υ∗ y[k], (17)

in which Υ ∈ CM Lr ×Gt Gr , Υ = ΦΨ is the measurement matrix and y[k] ∈ CM Lr ×1 is the received

signal for a given k, k = 0, . . . , K − 1. If there is correlation between noise components, the
 12

atom estimated from the projection in (17) is not likely to be the actual one. To introduce the
appropriate correction in the projection, the specific form the noise covariance matrix takes needs
to be used. It is important to highlight that depending on the design of the combiners for each
of the M training steps, the noise statistics will be different. Specifically, the noise covariance
matrix contains the (complex) inner product between the different hybrid combiners, for each
(m) (i)
training step m, m = 1, . . . , M . Let us consider two arbitrary (hybrid) combiners Wtr ,
(m) (j)
Wtr ∈ CNr ×Lr for two arbitrary training steps i,j and a given subcarrier k. Therefore, if the
(i) (m) (i)∗
combined noise at a given training step i and subcarrier k is denoted as nc [k] = Wtr n(i) [k],
with n(i) [k] ∼ N (0, σ 2 ILr ), the noise cross-covariance matrix is given by
(m) (i)∗ (i) (m) (j)
E{n(i) (j)∗
c [k]nc [k]} = E{Wtr n [k]n(j)∗ [k]Wtr }
(m) (i)∗ (m) (j)
= Wtr E{n(i) [k]n(j)∗ [k]}Wtr
(m) (i)∗ 2 (m) (j)
= Wtr σ δ[i − j]Wtr , (18)

such that it is the same for all subcarriers during training. It can be written as a block diagonal
matrix C ∈ CM Lr ×M Lr whose i-th diagonal block is the Gram matrix of the (hybrid) combiner
W(i) , that is, C = σ 2 Cw , where
(m) (1)∗ (m) (1) (m) (2)∗ (m) (2) (m) (M )∗ (m) (M )
Cw = blkdiag{Wtr Wtr , Wtr Wtr , . . . , Wtr Wtr }. (19)

Now, we can resort to the Cholesky decomposition of Cw

Cw = D∗w Dw , (20)

with Dw ∈ CM Lr ×M Lr an upper triangular matrix. The subscript in Cw and Dw indicates that

(m)
these matrices only depend on the combiners Wtr used during M consecutive training frames.
Therefore, the projection step is performed as

c[k] = Υ∗w y[k], (21)

where Υw ∈ CM Lr ×Gt Gr is the whitened measurement matrix given by Υw = D−1

w Υ. In this way,

the resulting projection simultaneously whitens the spatial noise components and estimates the
projection of the received signal onto the subspace spanned by the columns in the measurement
matrix.
 13

2) Computation of the channel gains: If the vectorized channel matrix for a given subcarrier
k, k = 0, 1, . . . , K − 1, is given by vec{H[k]} = (ĀT ◦ ĀR )ξ[k], with ξ[k] the L × 1 vector
of channel gains, then the received signal y[k] is distributed according to N (ΦΨξ[k], C). Once
the different AoAs/AoDs are found, the signal model for the k-th subcarrier can be written as

y[k] = Υ{:,T̂ } ξ̃[k] + ñc [k], (22)

where ñc [k] ∈ CM Lr ×1 is the residual noise in our linear model after estimating the channel
support. If the estimation of the support is accurate, ñc [k] will be close to the post-combining
noise vector nc [k]. The L̂ × 1 vector ξ̃[k] is the vector of channel gains to be estimated after
sparse recovery, where L̂ = |T̂ | is the estimate of the sparsity level. The matrix Υ{:,T̂ } ∈ CM Lr ×L̂
is defined as
 
Υ{:,T̂ } = ΦΨ {:,T̂ } . (23)

It is important to remark that the support estimated by the proposed algorithm may be different
from the actual channel support. Provided that, in general, T̂ can be different from the actual
support, the vector ξ[k] ∈ CL̂×1 containing the complex channel gains can be also different from
the vector containing the true channel gains. Since the model in (22) is linear on the parameter
vector ξ[k], there is a Minimum Variance Unbiased (MVU) Estimator that happens to be the
Best Linear Unbiased Estimator (BLUE) as well [23].
The mathematical equation in (22) is usually referred to as the General Linear Model (GLM),
for which the solution of ξ̃[k] for real parameters is provided in [23]. The extension for a complex
vector of parameters is straightforward and given by
ˆ = (Υ∗ C−1 Υ
ξ̃[k] −1 ∗ −1
{:,T̂ } {:,T̂ } ) Υ{:,T̂ } C y[k]. (24)
ˆ is the MVU estimator for our parameter vector ξ̃[k], k = 0, . . . , K − 1. Hence,
Therefore, ξ̃[k]
it is unbiased and attains the CRLB is the support is estimated correctly. It is interesting to
note that this corresponds to a Weighted Least-Squares (WLS) estimator, with the corresponding
weights given by the inverse covariance matrix of the noise. We are interested in finding the
CRLB for the estimation of the channel matrices at each subcarrier assuming perfect sparse
reconstruction. To that end, taking into account only the non-zero entries in hv [k], the Fisher
Information Matrix (FIM) is derived from the GLM in (13) as

I(ξ[k]) = Υ∗{:,T } C−1 Υ{:,T } . (25)

 14

Note that the previous expression gives the FIM for the vector ξ[k], which contains the actual
channel gains. Accordingly, the overall variance of the estimator for the vector of channel gains
is given by the sum of the individual CRLBs for each of the complex gains. Therefore, if we
denote the overall variance of the estimator for ξ[k] as γ(ξ[k]), k = 0, 1, . . . , K − 1, then

CRLB(ξ[k]) = I−1 (ξ[k]) (26)

and
γ(ξ[k]) = trace{CRLB(ξ[k])}. (27)

The last equation only takes into account the channel gains for each path in the extended virtual
channel representation. Our interest, however, is to compare the channel estimation performance
at the antenna level. To that end, the frequency-domain channel matrix can be vectorized as
(subcarrier-wise)
vec{H[k]} = (AT ◦ AR )ξ[k], (28)

where AT ∈ CNt ×L , AR ∈ CNr ×L are the array response matrices evaluated on the AoAs/AoD.
Notice that the decomposition in (28) is expressed with equality, which follows from the assump-
tion that the estimation of the support is perfect and only the channel gains are to be estimated.
The assumption on the perfect estimation of the support comes from the fact that the channel
vectors hv [k] are assumed to be sparse and the AoA/AoD are estimated perfectly. Since the
Fisher Information requires that the mean of the received signal is continuous and differentiable
on the parameters to estimate, it cannot be applied to the estimation of the angular parameters.
Owing to our lacking space, the solution to the problem of finding the CRLB in a more general
case accounting for continuous AoA and AoD is left for future work.
The overall minimum variance for an unbiased estimator of the Nt Nr entries in H[k] is
given by the sum of the variances for each element {H[k]}{i,j} , i = 1, . . . , Nr , j = 1, . . . , Nt .
Mathematically, following the same notation as with γ(ξ[k]), we will denote the overall variance
of the estimator for H[k] as γ(vec{H[k]}), for a given subcarrier k. Then, γ(vec{H[k]}) is
derived as
∂ vec{H∗ [k]}
 
∂ vec{H[k]}
γ(vec{H[k]}) = trace CRLB(ξ[k])
∂ξ[k] ∂ξ[k]
 
∗
= trace (AT ◦ AR )CRLB(ξ[k])(AT ◦ AR ) . (29)
 15

An important feature of this estimator is that the difference in performance given by the LS
and the WLS estimators is more accentuated as the number of RF chains grows (if and only if
the hybrid combiner is not built from orthonormal vectors). Furthermore, this estimator works
even though the noise variance is unknown. This is because the noise covariance matrix C
can be written as C = σ 2 Cw , where Cw models the coupling between the different combining
vectors. Since the combining vectors are known at the receiver, this matrix Cw is known as well.
ˆ ∈ CL̂×1
Therefore, we can write our estimator for the vector containing the channel gains, ξ̃[k]
as
ˆ = (Υ∗ C−1 Υ
ξ̃[k] −1 ∗ −1
{:,T̂ } w {:,T̂ } ) Υ{:,T̂ } Cw y[k]. (30)

When the noise variance is unknown, it follows from the Slepian-Bangs formula that the CRLB
for the vector ξ[k] does not change, since the noise mean does not depend on the variance σ 2 .
Thus, the elements in the FIM modeling the coupling between the noise variance and the channel
gains are zero. Therefore, the estimator of ξ[k] is efficient. In fact, it is also the ML estimator for
this complex parameter vector, since it can be seen to maximize the Log-Likelihood Function
(LLF) of the received signal. Nonetheless, there is no efficient estimator for both ξ[k] and σ 2 ,
since the MVUE of σ 2 depends on the true value for ξ[k]. Despite this, the MLE for σ 2 can
still be found by setting the partial derivative of the LLF to zero. If we use the identities
 
∂ ln det{C} −1 ∂C
= trace C (31)
∂σ 2 ∂σ 2
and
ˆ ∗ C−1 (y[k] − Υ
∂(y[k] − Υ{:,T̂ } ξ̃[k]) ˆ
{:,T̂ } ξ̃[k]) ˆ ∗ C−1 ×
= −(y[k] − Υ{:,T̂ } ξ̃[k])
∂σ 2 (32)
∂C −1 ˆ
C × (y[k] − Υ{:,T̂ } ξ̃[k]),
∂σ 2
it is possible to find that
1 ˆ ∗ C−1 (y[k] − Υξ̃[k]).
ˆ
σˆ2 = (y[k] − Υξ̃[k]) w (33)
M Lr
Finally, we can consider the information coming from all the subcarriers to enhance the
estimation of the noise variance. For the different subcarriers, we have the model

σˆ2 [k] = σ 2 + ν[k], k = 0, . . . , K − 1, (34)

where ν[k] models the estimation error for the variance in the k-th subcarrier. We known that,
in high SNR regime, the ML estimator is Gaussian distributed. Therefore, we have a set of K
 16

linear models for the estimation of the variance. Since the model is linear on the parameter σ 2 ,
we find that the ML of the variance is also the BLUE estimator, given by
1
σˆ2 ML = 1T σ, (35)
K
where σ ∈ RK×1 is defined as σ = [σˆ2 [0], σˆ2 [1], . . . , σˆ2 [K − 1]]T and 1 is the K × 1 column
vector containing a 1 in each position.

Algorithm 1 Simultaneous Weighted Orthogonal Matching Pursuit (SW-OMP)

1: procedure SW-OMP(z[k],Φ,Ψ,)

2: Compute the equivalent observation matrix

3: Υ = ΦΨ
4: Initialize the residual vectors to the input signal vectors and support estimate
5: r[k] = y[k], k = 0, . . . , K − 1, T̂ = {∅}
6:
7: while MSE >  do
8: Distributed Projection
9: c[k] = Υ∗ D−∗
w r[k], k = 0, . . . , K − 1
10: Find the maximum projection along the different spaces
p? = arg max K−1
P
11: k=0 |{c[k]}p |
p
12: Update the current guess of the common support
13:
14: T̂ = T̂ ∪ p?
15: Project the input signal onto the subspace given by the support using WLS
16: xT̂ [k] = (Υ∗{:,T̂ } Cw−1 Υ{:,T̂ } )−1 Υ∗{:,T̂ } C−1
w y[k], k = 0, . . . , K − 1

17:
18: T̂ ∪ T̂ = {1, 2, . . . , Gt Gr } , T̂ ∩ T̂ = {∅}
19: Update residual
20:
ˆ , where ξ̃[k]
r[k] = y[k] − Υ{:,T̂ } ξ̃[k] ˆ = x [k], k = 0, . . . , K − 1
21: T̂

22: Compute the current MSE

MSE = KM1 Lr K−1 ∗ −1
P
23:
24: k=0 r [k]Cw r[k]

25: end while

26: end procedure

The derived estimator for the channel gains ξ̂[k] can be used in the S-OMP algorithm. In
fact, the larger the number of subcarriers, the smaller the estimation variance the ML estimator
 17

can achieve. Thereby, if the number of averaging subcarriers K is large enough, the lack of
knowledge of the sparsity level is not so critical because of two reasons: 1) the computation
of the support is more precise due to noise averaging during the projection process, and 2)
if the support is estimated correctly, a particular estimate of σ 2 will be very close to the true
noise variance, such that the chosen halting criterion is optimal from the Maximum Likelihood
perspective. It should be clear that the higher the covariance between adjacent noise components,
the larger the performance gap between the S-OMP and the SW-OMP algorithm will be, which
actually depends on the ratio between Nr (Nt ) and Lr (Lt ).
The modification of the S-OMP algorithm to include the MVU estimator for the channel
gains, as well as the whitening matrix to estimate the support is provided in Algorithm 1. Notice
that the proposed algorithm performs both noise whitening and channel estimation, such that
interferences due to noise correlation are mitigated.

C. Subcarrier Selection Simultaneous Weighted-OMP + Thresholding

Despite the use of a single, subcarrier-independent measurement matrix Υ to estimate the

frequency-domain channels, the algorithm presented in the previous section exhibits high com-
putational complexity. The SW-OMP algorithm considers the distributed projection coming
from every subcarrier; however, a trade-off between performance and computational complexity
can be achieved if a small number of subcarriers Kp << K is used, instead. The problem
amounts as to how to choose those subcarriers, since no quality measure is available beforehand.
The ideal situation would require knowledge of the Signal-to-Noise Ratio (SNR), which is
unknown so far. Nonetheless, we have access to different frequency-domain received vectors
y[k], k = 0, 1, . . . , K − 1. Therefore, the quality measure to be used could be the l2 -norm of the
different vectors. Thereby, the Kp subcarriers having largest l2 -norm can be used to derive an
estimate of the support of the already defined sparse channel vectors hv [k], k = 0, . . . , K − 1.
The main problem concerning Matching Pursuit (MP) algorithms comes from the lack of
knowledge of the channel sparsity L. For this reason, there is usually an iteration at which L
paths have been detected but the estimate of the average residual energy is a little larger than
the noise variance itself. This makes the algorithm find additional paths which are not actually
contained in the MIMO channel. These paths usually have low power, and a pruning procedure
is needed to filter out these undesired components.
 18

Algorithm 2 Subcarrier Selection Simultaneous Weighted Orthogonal Matching Pursuit +

Thresholding (SS-SW-OMP+Th)
1: procedure SS-SW-OMP(y[k],Φ,Ψ,KP ,β,)

2: Initialize counter, set of subcarriers and residual vectors

3: i = 0, K = {∅}, r[k] = y[k], k = 0, . . . , K − 1
4: Find the Kp strongest subcarriers
5: while i ≤ Kp do
6: K = K ∪ arg max ky[k]k22
k6∈K
7: i=i+1
8: end while
9: Compute the equivalent observation matrix
10: Υ = ΦΨ
11: while MSE >  do
12: Distributed Projection
13: c[k] = Υ∗ D−∗
w r[k], k∈K
14: Find the maximum projection along the different spaces
p? = arg max k∈K |{c[k]}p |
P
15:
p
16: Update the current guess of the common support
17: T̂ = T̂ ∪ p?
18: Project the input signal onto the subspace given by the support using WLS
19: xT̂ [k] = (Υ∗{:,T̂ } Cw−1 Υ{:,T̂ } )−1 Υ∗{:,T̂ } C−1
w y[k], k = 0, . . . , K − 1

20: xT̂ = 0, where T̂ ∪ T̂ = {1, 2, . . . , Gt Gr } , T̂ ∩ T̂ = {∅}

21: Update residual
22:
ˆ , where ξ̃[k]
r[k] = y[k] − Υ{:,T̂ } ξ̃[k] ˆ = x [k], k = 0, . . . , K − 1
T̂

23: Compute the current MSE

MSE = KM1 Lr K−1 ∗ −1
P
24: k=0 r [k]Cw r[k]

25: end while

26: Thresholding based on maximum average power
P ? = max K1 K−1
P 2 1
PK−1 2
27: k=0 |{ ξ̂[k]}` | , p̂ av,i = K k=0 |{ξ̂[k]}i | , i = 1, . . . , L̂
`
T̂Th = i / p̂av,i ≥ βP ? , i ∈ T̂
S
28:

29: ξ̂[k] = {ξ̂[k]}T̂Th , k = 0, . . . , K − 1

30: end procedure
 19

A reasonable way to prune the undesired paths is to remove those components whose power
falls below a given threshold, which can be related to the average power of the component
in the estimated sparse vectors having maximum average power. Let us denote this power by
P ? . Then, the threshold can be defined as η = βP ? , β ∈ (0, 1). The value P ? is taken as
P ? = max K1 K−1 2
P
k=0 |{ξ̂[k]}` | . To keep the common sparsity property we must ensure that
`
the channel support after thresholding remains invariant. For this purpose, we define a signal
p̂av ∈ CL̂×1 whose i-th component is given by p̂av,i = K1 K−1 2
P
k=0 |{ξ̂[k]}i | , i = 1, . . . , L̂, such
that p̂av measures the average power along the different subbands of each spatial component in the
quantized angle grid. The final support after thresholding T̂Th is defined as T̂Th = L̂i=1 i / p̂av,i ≥
S

βP ? . Therefore, the components in ξ̂[k] indexed by T̂Th are the final channel gains estimates for
each subcarrier. The modification of the proposed SW-OMP algorithm to reduce computational
complexity and implement this pruning procedure is provided in Algorithm 2.

IV. R ESULTS

This section includes the main empirical results obtained with the two proposed algorithms,
SW-OMP and SS-SW-OMP + Thresholding, and comparisons with other the frequency-domain
channel estimation algorithms including SSAMP [16] and DGMP [17]. To obtain these results,
we perform Monte Carlo simulations averaged over many trials to evaluate the normalized mean
squared error (NMSE) and the ergodic rate as a function of SNR and number of training frames
M . We also provide calculations of the computational complexity for the proposed algorithms
in Table I and prior work in Table II.
The typical parameters for our system configuration are summarized as follows. Both the
transmitter and the receiver are assumed to use Uniform Linear Arrays (ULAs) with half-
wavelength separation. Such a ULA has steering vectors obeying the expressions {aT (θ` )}n =
q q
1 jnπ cos (θ` ) 1 jmπ cos (φ` )
Nt
e , n = 0, . . . , Nt − 1 and {aR (φ` )}m = Nr
e , m = 0, . . . , Nr − 1.
We take Nt = Nr = 32 and Gt = Gr = 64. The phase-shifters used in both the transmitter and
the receiver are assumed to have NQ quantization bits, so that the entries of the training vectors
2π(2NQ −1)
n o
(m) (m)
ftr , wtr , m = 1, 2, . . . , M are drawn from a set A = 0, 22π
NQ , . . . , 2NQ
. The number of
RF chains is set to Lt = 1 at the transmitter and Lr = 4 at the receiver. The number of OFDM
subcarriers is set to K = 16.
We generate channels according to (2) with the following parameters:
 20

• The L = 4 channel paths are assumed to be independent and identically distributed, with
1
delay τ` chosen uniformly from [0, (Nc − 1)Ts ], with Ts = 1760
s, as in the IEEE 802.11ad
wireless standard.
• The AoAs/AoDs are assumed to be uniformly distributed in [0, π).
• The gains of each path are zero-mean complex Gaussian distributed such that Ek {kH[k]k2F } =
Nr Nt
ρL
.
• The band-limiting filter prc (t) is assumed to be a raised-cosine filter with roll-off factor of
0.8.
• The number of delay taps of the channel is set to Nc = 4 symbols.
The simulations we perform consider channel realizations in which the AoA/AoD are off-grid,
i.e. do not correspond to the angles used to build the dictionary, and also the on-grid case, to
analyze the loss due to the model mismatch.

A. NMSE Comparison

One performance metric is the Normalized Mean Squared Error (NMSE) of a channel estimate
Ĥ[k] for a given realization, defined as
PK−1
k=0 kĤ[k] − H[k]k2F
NMSE = PK−1 . (36)
k=0 kH[k]k2F
The NMSE will be our baseline metric to compute the performance of the different estimators,
and will be averaged over many channel realizations. The normalized CRLB (NCRLB) is also
provided to compare the average performance of each algorithm with the lowest achievable
NMSE, and will also be averaged over many channel realizations
PK−1
CRLB(vec{H[k]})
NCRLB = k=0PK−1 2
. (37)
k=0 kH[k]kF

We compare the average NMSE versus SNR obtained for the different channel estimation
algorithms in Fig. 2 for a practical SNR range of −15 dB to 10 dB, on-grid AoA/AoDs, and two
different values for the number of training frames. SW-OMP performs the best, achieving NMSE
values very close to the NCRLB. SS-SW-OMP+Th performs similarly to SW-OMP, although
there is some performance loss due to the fact that SS-SW-OMP+Th does not employ every
subcarrier to estimate the common support of the sparse channel vectors. In SS-SW-OMP+Th,
the number of selected subcarriers for the estimation of the support was set to Kp = 4 for
 21

illustration and the parameter β was chosen as β = 0.025, which is a reasonably small value to
filter out undesired components in the sparse channel estimate. From the curves shown in Fig. 2,
OMP performs poorly for all the SNR range due to the fact that it is not designed to process
several vectors which are sparse in a common vector basis. Exploiting common sparsity provides
an NMSE reduction of approximately 7 dB, although there are slight variations depending on
the SNR. This improvement comes at the cost of a higher offline computational complexity in
the proposed algorithms in comparison with OMP, as we show in Section IV-C.
We also observe that the DGMP algorithm from [17] performs the worst, due to the fact that
it has been designed to estimate near-LoS mmWave MIMO channels. Since it only estimates
a single path, the estimation error for NLOS channels is large. The algorithm SSAMP is also
shown for comparison. We can see that, at low SNR regime, the information on the common
support is enough to outperform the OMP algorithm, but not at the high SNR regime. This comes
from the fact that the algorithm is estimating more than a single path per iteration. Provided that
the dictionary matrices are not square in this setting, the redundancy between columns in the
transmit and receive array matrices makes it difficult to properly estimate more than a support
index per iteration.

0
SW-OMP
0
SW-OMP SW-OMP+Th
-5
SW-OMP+Th OMP
-5 OMP CRLB
NMSE (dB)

CRLB -10 SSAMP

NMSE (dB)

SSAMP DGMP
-10
DGMP
-15
-15
-20
-20
-25
-25

-15 -10 -5 0 5 10 -15 -10 -5 0 5 10

SNR (dB) SNR (dB)

(a) (b)
Fig. 2. Evolution of the NMSE versus SNR for the different frequency-domain algorithms when the AoA/AoDs are assumed
to lie on the dictionary grid. The number of training frames is set to M = 80 (a) and M = 120 (b).

The previous simulations showed the performance of the algorithms when the channel fits
 22

the on-grid model, but it is also important to analyze the performance in a practical scenario,
when the AoA/AoD do not fall within the quantized spatial grid. Fig. 3 shows the performance
of the different algorithms for off-grid angular parameters. It can be observed that when using
Gt = Gr = 64 there is a considerable loss in performance due to the use of a fixed dictionary
to estimate the AoA/AoD. However, the estimation error is below −10 dB for values of SNR
in the order of 0 and beyond. On the other hand, provided that the SNR expected in mmWave
communication systems is in the order of −20 dB up to 0 dB, the large gap between the proposed
algorithms and the NCRLB needs to be reduced. Increasing the size of the dictionary is one of the
possible solutions, as shown by the curves corresponding to Gt = Gr = 128 and Gt = Gr = 256.

SW-OMP SW-OMP
0 SW-OMP+Th 0 SW-OMP+Th
OMP OMP
CRLB
-5 CRLB -5
NMSE (dB)

NMSE (dB)

-10
-10

-15
-15
Gt = Gr = 64
-20 Gr = Gt = 64
-20 Gt = Gr = 128
Gr = Gt = 128
-25
-25 Gt = Gr = 256 Gr = Gt = 256

-15 -10 -5 0 5 10 -15 -10 -5 0 5 10

SNR (dB) SNR (dB)

(a) (b)
Fig. 3. Evolution of the NMSE versus SNR for the different frequency-domain algorithms. The number of training frames is
set to M = 80 (a) and M = 120 (b). The AoD/AoA are generated from a continuous uniform distribution.

We show in Fig. 4 the performance of the different frequency-domain algorithms when

increasing the number of subcarriers. The parameters for the simulation scenario are the same
as in Fig. 2, however, the number of subcarriers is set to K = 64 in this case. Kp is set
to 32 subcarriers and β = 0.025σ 2 . Interestingly, both SW-OMP and SS-SW-OMP+Th are
asymptotically efficient since they are both unbiased and attain the NCRLB. A magnified plot
for SNR = −5 dB is also shown to clearly see the performance gap between the different
algorithms and also the NCRLB.
Fig. 5 shows the average NMSE vs number of training frames M . The number of training
frames M is increased from 20 to 100. The remaining parameters in the simulation scenario are
 23

SW-OMP
0 SW-OMP+Th
OMP
-5 CRLB

NMSE (dB)
-10

-15 -10

-12
-20
-14
-25
-10 -5 0

-15 -10 -5 0 5 10
SNR (dB)

Fig. 4. Comparison of evolution of the NMSE versus SNR for the different frequency-domain algorithms. The number of
training frames is set to M = 80. The number of subcarriers is set to K = 64.

the same as in Fig. 2. Results are shown for channel realizations in which the angular parameters
fall within the quantized angle grid and when they do not. The average performance of OMP
algorithm is poor for all the considered cases, which comes from its inability to exploit the
common support property shared by the different subchannels. SW-OMP can be observed to
provide the best performance for all the values of M and SNR. Furthermore, its performance
lies closer to the NCRLB at higher SNR. The larger the number of training frames and the
higher the SNR, the estimation of the support is more robust and gets closer to the actual one.
The NCRLB provided assumes perfect sparse recovery, which is the reason why the estimation
error does not lie within the bound. Nonetheless, in the on-grid case, if the number of training
frames is large enough and the SNR is not low, the performance gap between SW-OMP and
the NCRLB is lower than 1 dB. The difference in performance between SS-SW-OMP+Th and
SW-OMP reduces when either M or SNR is increased. As in Fig. 2, there is a big difference
in performance between OMP and SW-OMP, depending on both the SNR and the number of
frames. We also observe a high loss in performance when the angular parameters do not fall
within the grid.

B. Spectral efficiency comparison

Another performance metric is the spectral efficiency, computed assuming fully-digital precod-
ing and combining using the Ns dominant left and right singular vectors of the channel estimate.
 24

-4
4 0
SW-OMP SW-OMP SW-OMP
SW-OMP+Th SW-OMP+Th -6 SW-OMP+Th
2 -2
OMP OMP OMP
CRLB CRLB -8 CRLB
0 -4
NMSE (dB)

NMSE (dB)

NMSE (dB)
-10
-2 -6
-12
-8
-4 -14
-10
-6 -16
-12
-8 -18

20 40 60 80 100 20 40 60 80 100 20 40 60 80 100

Training frames (M) Training frames (M) Training frames (M)

(a) (b) (c)

SW-OMP 0 SW-OMP
6 SW-OMP 2 SW-OMP+Th SW-OMP+Th
SW-OMP+Th
OMP OMP
4 OMP 0 CRLB CRLB
CRLB
-2 -5
NMSE (dB)

NMSE (dB)

NMSE (dB)
2

0 -4

-6 -10
-2
-8
-4
-10 -15
-6
20 40 60 80 100 20 40 60 80 100 20 40 60 80 100
Training frames (M) Training frames (M) Training frames (M)

(d) (e) (f)

Fig. 5. Comparison of evolution of the NMSE versus number of training frames M at different SNR for the different channel
estimation algorithms. The SNR is set to −10 dB (a) and (d), −5 dB (b) and (e) and 0 dB (c) and (f). The curves in the first
row consider on-grid angular parameters on the channel realizations, while the ones in the second row consider the off-grid
case.

This gives Ns parallel channels each with gain λn (Heff [k]), leading to average spectral efficiency
K−1 Ns  
1 XX ρ 2
R= log 1 + λn (Heff [k]) . (38)
K k=0 n=1 KNs σ 2
In Fig. 6(a), we show the achievable spectral efficiency as a function of the SNR for the
different channel estimation algorithms when off-grid AoA/AoDs are considered. The simulation
parameters are the same as in Fig. 2. The difference in performance between OMP and either
of our two proposed algorithms is approximately constant as the SNR increases, which comes
from the fact that OMP does not force the channel estimates to share the same support. The two
proposed algorithms perform very similarly for all the range of SNR, which is an indicator
that Kp < K subcarriers are enough to obtain a reliable channel estimate. Therefore, SS-
SW-OMP+Th can be claimed to be a good trade-off between performance and computational
complexity.
 25

5.5

5
12
Spectral Efficiency (bps/Hz)
SNR = 0 dB
4.5

Spectral Efficiency (bps/Hz)

10 4 SW-OMP
SW-OMP+Th
SW-OMP 3.5 OMP
8 Perfect CSI
SW-OMP+Th 3
OMP
6 Perfect CSI 2.5
SNR = -5 dB
2
4
1.5 SNR = -10 dB

2 1

0.5

-15 -10 -5 0 5 10 20 30 40 50 60 70 80 90 100

SNR (dB) Training frames (M)

(a) (b)
Fig. 6. (a) Evolution of the spectral efficiency versus SNR for the different frequency-domain algorithms. The number of training
frames is set to M = 80. (b) Evolution of the spectral efficiency versus number of training frames M at different SNR for the
different frequency-domain algorithms.

Finally, Fig. 6(b) shows the spectral efficiency as a function of the number of training
frames M for the different channel estimation algorithms. Comparisons are provided for SNR =
{−10, −5, 0} dB. We observe that SW-OMP is the algorithm providing best performance, fol-
lowed closely by SS-SW-OMP+Th, whilst OMP performs the worst. For low values of SNR,
there is a noticeable gap between SW-OMP and the perfect CSI case, which becomes smaller
as M increases. For small values of M , the performance of the sparse recovery algorithms is
very poor. Conversely, using M ≥ 60 frames does not result in a significant improvement in
performance. Simulations also show that near-optimal achievable rates can be achieved by using
a reasonable number of frames, i.e., 60 ≤ M ≤ 100.

C. Computational complexity

The computational complexity for each step in the different algorithms is also provided in j
in Table I. Since some steps can be performed before running the channel estimation algorithms,
we will distinguish between on-line and off-line operations. Values are provided for a single
iteration. In the case of the SSAMP algorithm in [16] and the DGMP algorithm in [17], we take
the notation used in the corresponding papers for frequency-domain vectors and measurement
matrices.
 26

TABLE I
O NLINE C OMPUTATIONAL COMPLEXITY OF PROPOSED ALGORITHMS

SW-OMP

Operation Complexity
∗
(K) × c[k] = Υ D−∗
w r[k] O(K(Gr Gt − (j − 1))Lr M )
P −1
Maximum of N k=0 |{c[k]}p | O(K(Gr Gt − (j − 1)))
(K) × xT̂ [k] = (Υ∗{:,T̂ } C−1
w Υ{:,T̂ } )
−1 ∗
Υ{:,T̂ } C−1
w z[k] O(2j 2 Lr M + j 3 + KjLr M )
(K) × r[k] = z[k] − Υ{:,T̂ } ξ̂[k] O(KLr M )
1
PK−1 ∗ −1
MSE = KM Lr k=0 r [k]Cw [k]r[k] O(2KLr M )
Overall O(Lr M (K(Gr Gt − (j − 1)) + 2j 2 + Kj))

SS-SW-OMP+Thresholding

Operation Complexity
Find the Kp strongest subcarriers O(KLr M )
∗
(Kp ) × c[k] = Υ D−∗
w r[k] O(Kp (Gr Gt − (j − 1))Lr M )
P
Maximum of k∈K |{c[k]}p | O(Kp (Gr Gt − (j − 1)))
(K) × xT̂ [k] = (Υ∗{:,T̂ } C−1
w Υ{:,T̂ } )
−1 ∗
Υ{:,T̂ } C−1
w z[k] O(2j 2 Lr M + j 3 + KjLr M )
(K) × r[k] = z[k] − Υ{:,T̂ } ξ̂[k] O(KLr M )
1
PK−1 ∗ −1
MSE = KM Lr k=0 r [k]Cw [k]r[k] O(2KLr M )
Thresholding O(K L̂)
Overall O(Lr M (Kp (Gr Gt − (j − 1)) + 2j 2 + Kj))

The computational complexity of SS-SW-OMP+Th is lower than in SW-OMP, since only

a smaller number of projections are computed to estimate the channel support. It must be
noticed that the complexity of SS-SW-OMP+Th is lower than that of OMP, since the matrix
product Υ∗ C−1
w can be computed before explicit channel estimation. The online computational

complexity of SW-OMP is the same as that of OMP, since OMP performs the same projections
as SW-OMP. Conversely, the offline computational complexities of both SW-OMP and SS-SW-
OMP+Th are higher than those for the other algorithms, since the matrix product Υ∗ C−1
w must be

computed before explicit channel estimation. During estimation, only those rows corresponding
to the estimate of the support must be selected to compute the channel gains. Nonetheless,
OMP works in an isolated manner, independently solving each sparse recovery problem for
each subcarrier, whereas SW-OMP exploits information that is common to every subcarrier.
The complexity of the SSAMP algorithm proposed in [16] is also shown for comparison. It is
 27

TABLE II
O NLINE C OMPUTATIONAL COMPLEXITY OF PREVIOUSLY PROPOSED ALGORITHMS

OMP from [13]

Operation Complexity
∗
(K) × c[k] = Υ [k]r[k] O(K(Gr Gt − (j − 1))Lr M )
(K)× Maximum of |{c[k]}p | O(K(Gr Gt − (j − 1)))
(K) × xT̂ [k] = Υ†{:,T̂ } [k]z[k] O(K(2j 2 Lr M + j 3 ))
(K) × r[k] = z[k] − Υ{:,T̂ } [k]ξ̂[k] O(KLr M )
(K) × MSE = 1
M Lr
r∗ [k]r[k] O(KLr M )
Overall O(KLr M (Gr Gt − (j − 1)) + 2j 2 )

SSAMP from [16]

Operation Complexity
(K) × ap = Φ∗p bi−1
p O(KGr Gt Lr M )
PN −1 2
Maximum of p=0 kap k2 O(N Gt Gr )
(K) × {tp }Ωi−1 ∪Γ = ((Φp )†Ωi−1 ∪Γ rp O(K(4j 2 M Lr + (2j)3 ))
Prune support Ω = arg max { P 2
P
p=1 k{tp }Ω̃ k2 }, |Ω̃| =T O(2Kj)
Ω̃

(K) × {cp }Ω = ({Φp }Ω )† rp O(K(2j 2 M Lr + j 3 ))

(K) × bp = rp − Φp cp O(KLr M Gt Gr )
PK−1 2
Computation of total error p=0 kbp k2 O(KM Lr )
Overall O(2KLr M (Gr Gt + 6j 2 ))

DGMP from [17]

Operation Complexity
(K) × ap = Υ∗p rp O(KGr Gt Lr M )
PN −1 2
Maximum of ρ = arg max p=0 kap k2 O(KGt Gr )
ρ̃

(K) × {αp }ρ = ({Φp }†ρ rp O(K(2j 2 M Lr + j 3 ))

Overall O(KLr M (Gr Gt + 2j 2 ))

important to remark that this algorithm is based on the assumption that the frequency-selective
hybrid precoders and combiners used during training are semiunitary. Therefore, this algorithm
does not take into account the noise covariance matrix Cw , neither for the computation of the
channel gains nor for the estimation of the support of the sparse channel matrices.
We observe that the computational complexity of SW-OMP is lower than its SSAMP coun-
terpart. This is because SSAMP exhibits an increase in complexity of at most O(4j 2 ) owing
 28

to the estimation of j paths at the j-th iteration of SSAMP. This is because this algorithm
uses an iteration index i to estimate the sparsity level L, and a stage index j to estimate the
j channel paths found at the current iteration. Afterwards, the support of the channel estimate
is pruned to select the j most likely channel paths. Therefore, at a given iteration i and stage
j, at most 2j = |Ω̃i−1 ∪ Γ| paths are estimated and then pruned, such that only j paths are
selected among the 2j candidates. The union of the sets Ω and Γ comes from the possibility
of finding new potential paths at the i-th iteration and the j-th stage. This is done by jointly
considering the paths found at (i − 1)-th iteration and the ones found in the j-th stage within
the i-th iteration. Therefore, whereas both SW-OMP, SS-SW-OMP+Th and OMP estimate a
single path at a given iteration j, SSAMP estimates at most 2j different paths by using LS.
When computing the pseudoinverse during LS estimation, this results in an additional increase
in complexity of O(4j 2 ), as shown in Table II. By contrast, as shown in Table I, our proposed
ML estimator for the channel gains exhibits computational complexity in the order of O(2j 2 ),
thereby slightly reducing the number of operations. Further, both SW-OMP, SS-SW-OMP+Th
and OMP do not compute the projections between every column of the measurement matrix
and the received frequency-domain signals. Instead, a path that was already found in a previous
iteration is not estimated again, thereby keeping computational complexity lower than that of
SSAMP.
Last, SSAMP considers the use of K different measurement matrices for channel estimation.
As discussed in Section III, this work considers analog-only training precoders and combiners
to estimate the channel, which are constrained to be frequency flat. Consequently, the noise
covariance matrix is constant across the different subcarriers. Besides, a frequency-flat data
stream s(m) ∈ CNs ×1 is transmitted at every subcarrier for a given training step m, 1 ≤ m ≤ M .
Thereby, a single frequency-flat measurement matrix is used to jointly estimate the K different
subchannels. This results in a complexity reduction by a factor of K at every iteration, both
during the estimation of the support and during estimation of the channel gains. Thereby, the
entire process of estimating the channel is greatly simplified.
While OMP does not require any off-line operation, both SW-OMP and SS-SW-OMP+Th need
to compute the whitened measurement matrices Υw = D−1
w Υ for the different subcarriers. The

offline computation of D−1 M 3 −1

w has complexity of O( 3 Lr ), since Cw is a block diagonal matrix

containing M hermitian matrices. Therefore, the Cholesky decomposition of C−1

w is nothing but
 29

the block diagonal matrix containing the Cholesky decompositions of the covariance matrices for
3
the different training frames. The cost of such decomposition for a matrix A ∈ Ck×k is O( k3 ).
The overall cost is calculated taking this individual cost into account. It is important to remark
that this cost comes from the use of frequency-flat precoders/combiners and training symbols.
This entails a reduction in computational complexity with respect to the case in which frequency-
selective baseband combiners were used during the channel estimation stage. Nonetheless, in our
simulations, only analog precoders and combiners are used during the training stage, therefore
reducing the computational complexity by a factor of K.

V. C ONCLUSIONS

In this paper, we proposed two compressive channel estimation approaches suitable for OFDM-
based communication systems. These two strategies are based on jointly-sparse recovery to
exploit information on the common basis that is shared for every subcarrier. Our compressive
approaches enable MIMO operation in mmWave systems since the different subchannels are
simultaneously estimated during the training phase. Further, if there is no grid quantization error
and the estimation of the support is correct, we showed that our algorithms are asymptotically
efficient since they asymptotically attain the CRLB. In simulations, we found that only a small
number of subcarriers provide a high probability of correct support detection, thus our estimators
approach the CRLB with reduced computational complexity. The approaches were also found
to work well even with off-grid parameters, and to outperform competitive frequency-domain
channel estimation approaches. For future work, it would be interesting to analytically calculate
the minimum required number of subcarriers to guarantee a high probability of correctly recov-
ering the support of the sparse channel vectors. It would also be interesting to study the effects
of other impairments including array miscalibration, beam squint, and synchronization errors.

R EFERENCES
[1] T. Rappapport et al, Millimeter Wave Wireless Communications. Pearson Education, Inc., 2014.
[2] R.W. Heath, N. González-Prelcic, S. Rangan, W. Roh, and A. M. Sayeed, ”An overview of signal processing techniques
for millimeter wave MIMO systems”, IEEE J. Sel. Areas Commun., vol. 10, no. 3, pp. 436-453, April 2016.
[3] A. Alkhateeb, J. Mo, N. González-Prelcic, and R. W. Heath Jr., ”MIMO precoding and combining solutions for millimeter-
wave systems”, IEEE Commun. Mag., vol. 52, no. 12, pp. 122-131, Dec 2014
[4] M. L. Malloy and R. D. Nowak, ”Near-optimal adaptive compressed sensing”, in Proc. Asil. Conf. Signals, Syst. Comp.
(ASILOMAR), Pacific Grove, CA, 2012, pp. 1935-1939.
 30

[5] -, ”Near-optimal compressive binary search”, arXiv preprint arXiv:1306.6239, 2012.

[6] M. Iwen and A. Tewfik, ”Adaptive strategies for target detection and localization in noisy environments”, IEEE Trans.
Signal Process., vol. 60, no. 5, pp. 2344-2353, 2012.
[7] A. Alkhateeb, O. E. Ayach, G. Leus, and R. W. Heath Jr. ”Channel estimation and hybrid precoding for millimeter wave
cellular systems”, IEEE J. Sel. Topics Signal Process., vol. 8, no. 5, pp. 831-846, Oct. 2014.
[8] D. Ramasamy, S. Venkateswaran, and U. Madhow, ”Compressive adaptation of large steerable arrays”, in Information Theory
and Applications Workshop (ITA), 2012, Feb. 2012, pp. 234-239.
[9] -, ”Compressive tracking with 1000-element arrays: A framework for multi-Gbps mm wave cellular downlinks”, in Proc.
Annual Allerton Conference on Communication, Control, and Computing, Oct. 2012, pp. 690-697.
[10] D. E. Berraki, S. M. D. Armour, and A. R. Nix, ”Application of compressive sensing in sparse spatial channel recovery for
beamforming in mmWave outdoor systems”, in Proc. IEEE Wireless Communications and Networking Conference (WCNC),
Apr. 2014, pp. 887-892.
[11] J. Lee, G. Gye-Tae, and Y. H. Lee, ”Exploiting spatial sparsity for estimating channel sof hybrid MIMO systems in
millimeter wave communications”, in Proc. IEEE Globecom, 2014.
[12] R. Méndez-Rial, C. Rusu, N. González-Prelcic, A. Alkhateeb, and R. W. Heath Jr., ”Hybrid MIMO architectures for
millimeter wave communications”, in Proc. IEEE Globecom, 2014.
[13] K. Venugopal, A. Alkhateeb, N. González-Prelcic, and R. W. Heath Jr., ”Time-domain channel estimation for wideband
millimeter wave systems with hybrid architecture”, in submitted to Int. Conf. Acoust. Speech and Sig. Proc. (ICASSP), Sept.
2016, pp. 1-5.
[14] K. Venugopal, A. Alkhateeb, N. González-Prelcic, and R. W. Heath Jr., ”Channel estimation for hybrid architecture based
wideband millimeter wave systems”, IEEE Journal on Selected Areas in Communications, to appear, 2017.
[15] J. Rodrı́guez-Fernández, N. González-Prelcic, K. Venugopal, and R. W. Heath Jr., ”A frequency-domain approach to
wideband channel estimation in millimeter wave systems”, in Proc. IEEE Int. Conf. Commun. (ICC), 2017.
[16] Z. Gao, L. Dai, and Z. Wang, ”Channel estimation for mmwave massive MIMO based access and backhaul in ultra-dense
network”, in Proc. IEEE Int. Conf. on Commun. (ICC), May 2016, pp. 1-6.
[17] Z. Gao, C. Hu, L. Dai, and Z. Wang, ”Channel estimation for millimeter-wave massive MIMO with hybrid precoding over
frequency-selective fading channels”, IEEE Communications Letters, vol. 20, no 6, pp. 1259-1262, June 2016.
[18] P. Schniter and A. Sayeed, ”Channel estimation and precoder design for millimeter-wave communications: The sparse
way”, in Proc. Asilomar Conf. Signals, Syst., Comput., Nov. 2014, pp. 273-277.
[19] A. Alkhateeb and R. W. Heath Jr., ”Freuqnecy selective hybrid precoding for limited feedback millimeter wave systems”,
IEEE Trans. Commun., vol. 64, no. 5, pp. 1801-1818, May 2016.
[20] S. G. Larew, T. A. Thomas, M. Cudak, and A. Ghosh, ”Air interface design and ray tracing study for 5g millimeter wave
communicaitons”, in Proc. IEEE Globecom Workshops (GC Wkshps), Dec. 2013.
[21] T. S. Rappaport, G. R. MacCartney, M. K. Samimi, and S. Sun, ”Wideband millimeter-wave propagation measurements
and channel models for future wireless communication system design”, IEEE Transactions on Communications, vol. 63,
no. 9, pp. 3029-3056, Sept. 2015.
[22] J. A. Tropp, A. C. Gilbert, and M. J. Strauss, ”Simultaneous sparse approximation via greedy pursuit”, in Proc. IEEE
Conf. Acous. Speech, and Signal Processing (ICASSP), 2005.
[23] S. M. Kay, Fundamentals of Statistical Signal Processing, Volume I: Estimation Theory. Prentice Hall PTR, 1993.

View publication stats