1

Approximate Gram-Matrix Interpolation
for Wideband Massive MU-MIMO Systems

arXiv:1610.00227v2 [] 19 Jun 2017

Charles Jeon, Zequn Li, and Christoph Studer

Abstract
Numerous linear and non-linear data-detection and precoding algorithms for wideband massive multiuser (MU) multiple-input multiple-output (MIMO) wireless systems that rely on orthogonal frequencydivision multiplexing (OFDM) or single-carrier frequency-division multiple access (SC-FDMA) require
the computation of the Gram matrix for each active subcarrier. Computing the Gram matrix for each
active subcarrier, however, results in excessively high computational complexity. In this paper, we propose
novel, approximate algorithms that significantly reduce the complexity of Gram-matrix computation by
simultaneously exploiting correlation across subcarriers and channel hardening. We show analytically that
a small fraction of Gram-matrix computations in combination with approximate interpolation schemes
are sufficient to achieve near-optimal error-rate performance at low computational complexity in massive
MU-MIMO systems. We also demonstrate that the proposed methods exhibit improved robustness against
channel-estimation errors compared to exact Gram-matrix interpolation algorithms that typically require
high computational complexity.
Index Terms
Equalization, interpolation, massive MU-MIMO, orthogonal frequency-division multiplexing (OFDM),
precoding, single-carrier frequency-division multiple access (SC-FDMA).

I. I NTRODUCTION

M

ASSIVE multi-user (MU) multiple-input multiple-output (MIMO) is widely believed to be a key
technology for fifth-generation (5G) wireless systems [1], [2]. The idea of massive MU-MIMO

is to equip the base-station (BS) with hundreds or thousands of antenna elements while serving tens
C. Jeon and C. Studer are with the School of Electrical and Computer Engineering, Cornell University, Ithaca, NY, USA
(e-mail: jeon@csl.cornell.edu, studer@cornell.edu). Website: vip.ece.cornell.edu
Z. Li was a visitor at the School of Electrical and Computer Engineering, Cornell University, Ithaca, NY and is now with
Department of Automation, Tsinghua University, Beijing, China (e-mail: li-zq13@mails.tsinghua.edu.cn).
A simulation framework to reproduce our results will be made available on GitHub after (possible) acceptance of this paper.

2

of users in the same time-frequency resource. Such large antenna arrays enable extremely fine-grained
beamforming in the uplink (users transmit to the BS) and in the downlink (BS transmits to the users),
which offers improved spectral efficiency compared to traditional, small-scale MIMO technology.
In the uplink, linear data-detection algorithms that rely on minimum-mean square error (MMSE)
equalization or zero-forcing (ZF) equalization are known to achieve near-optimal error-rate performance
in realistic massive MU-MIMO systems with a finite number of transmit antennas [3]–[5]. Non-linear
data-detection algorithms [6]–[8] have recently been shown to outperform linear methods in massive
MU-MIMO systems in which the number of users is comparable to the number of BS antennas. Most
of these linear and non-linear data-detection algorithms entail high computational complexity, often
dominated by the computation of the so-called Gram matrix G = HH H [9], where H ∈ CB×U is the

(uplink) channel matrix, B is the number of BS antennas, and U is the number of (single-antenna) users.
The computational complexity is orders-of-magnitudes higher in wideband systems that use orthogonal
frequency-division multiplexing (OFDM) or single-carrier frequency-division multiple access (SC-FDMA),
in which a Gram matrix must be computed for each active subcarrier (i.e., subcarriers used for pilots
or data transmission) [9]. In the massive MU-MIMO downlink, precoding algorithms are necessary to
avoid multi-user interference [1]. In wideband systems, the complexity of most precoding algorithms
is—analogously to the uplink—dominated by computation of the Gram matrices on all active subcarriers.
While some data-detection and precoding algorithms have been proposed that avoid the computation of
the Gram matrix altogether (see, e.g., [10]–[12]), these methods do not allow the re-use of intermediate
results in time-division duplexing (TDD) systems. Specifically, the Gram matrix and its inverse cannot be
re-used in the uplink (for equalization) and downlink (for precoding), which would significantly lower
the computational complexity. Hence, such algorithms inevitably perform redundant computations during
data-detection and precoding, which leads to inefficient transceiver designs.
A. Interpolation-Based Matrix Computations
In practical wideband communication systems, e.g., building upon IEEE 802.11n [17] and 3GPPLTE [18], the channel’s delay spread is often significantly smaller than the number of active subcarriers.
Hence, the channel coefficients are strongly correlated across subcarriers, which can be exploited to
reduce the computational complexity of commonly-used matrix computations required in multi-antenna
systems. More specifically, the papers [19]–[22] avoid a brute-force approach in traditional, small-scale,
and point-to-point MIMO-OFDM systems by using exact, interpolation-based algorithms for matrix
inversion and QR factorization. While a few hardware designs [23], [24] have demonstrated the efficacy
of these exact interpolation methods in small-scale MIMO systems, they do not scale well to wideband

3

massive MU-MIMO systems with hundreds or thousands of BS antennas, tens of users, and thousands of
subcarriers. In addition, the impact of imperfect channel-state information (CSI) and antenna correlation
on such exact, interpolation-based matrix computation algorithms is routinely ignored, but significantly
affects their performance in practical scenarios (see Section VI for a detailed discussion).
B. Contributions
Inspired by the exact, interpolation-based matrix computation algorithms put forward in [19]–[22]
for traditional, small-scale MIMO systems, we propose novel algorithms for approximate Gram matrix
computation in wideband massive MU-MIMO systems. We start by establishing the minimum number of
Gram matrix base-points that are required for exact interpolation. We then show that channel-hardening
in massive MU-MIMO enables approximate interpolation schemes that achieve near-exact error-rate
performance, even with strong undersampling in the frequency domain. In particular, we provide analytical
results that characterize the approximation errors of the proposed interpolation methods depending on
the delay spread and the antenna configuration. We furthermore derive exact mean-squared error (MSE)
expressions of our approximate interpolation algorithms for imperfect CSI and BS-antenna correlation.
We characterize the trade-offs between computational complexity and error-rate performance in realistic
massive MU-MIMO-OFDM systems, and we demonstrate the robustness of our approximate interpolation
methods for realistic scenarios with imperfect CSI and BS-antenna correlation.
C. Relevant Prior Art
Data detection and precoding for small-scale, single- and multi-carrier MIMO systems is a well studied
topic; see, e.g., [4], [13]–[16] and the references therein. However, all of these results focus on small-scale,
point-to-point MIMO systems and do (i) not exploit specific properties of massive MU-MIMO systems
and (ii) ignore the fact that time-division duplexing (TDD)-based systems must perform data-detection
and precoding, and hence, can re-use intermediate results (such as the Gram matrix) to reduce the
computational complexity. In contrast, our results exploit the specifics of massive MU-MIMO systems,
namely channel hardening, and enable a re-use of the computations carried out in the uplink for downlink
precoding.
The recent report [25] proposed an approximate interpolation-based ZF-based equalizer for a wideband
massive MU-MIMO testbed. In contrast to our work, the authors interpolate the inverse of the Gram
matrix. While simulation results in [25] show that the method works well in practice, no theoretical
results have been provided. In contrast, we use approximate methods to interpolate the Gram matrix, and

4

we provide exact analytical results that provide a solid foundation of approximate interpolation methods
in massive MU-MIMO systems.
D. Notation
Lowercase and uppercase boldface letters stand column vectors and matrices, respectively. The transpose,
Hermitian, and pseudo-inverse of the matrix A are denoted by AT , AH , and A† , respectively. We use
[A]m,n to represent the mth row and nth column entry of the matrix A. Sets are designated by uppercase

calligraphic letters, and |A| denotes the cardinality of the set A. The complex conjugate of a ∈ C is

a∗ . The indicator function is defined as I(a), where I(a) = 1 if a is true and 0 otherwise. CN (µ, Σ)

designates a complex Gaussian distribution with mean vector µ and covariance matrix Σ.
E. Paper Outline

The rest of the paper is organized as follows. Section II introduces the necessary prerequisites. Section III
proposes our exact and approximate interpolation-based Gram matrix computation algorithms. Section IV
and Section V provide an approximation error and complexity analysis of the proposed approximate
interpolation methods, respectively. Section VI shows numerical simulation results. We conclude in
Section VII.
II. P REREQUISITES
We start by summarizing the considered wideband massive MU-MIMO system and channel model. We
then outline computationally-efficient ways for linear data detection and precoding that make use of the
Gram matrix.
A. System Model
Without loss of generality, we focus on the uplink1 of a wideband massive MU-MIMO system with B
base-station antennas, U single-antenna user terminals (with U  B ), and W subcarriers. For each active

subcarrier ω ∈ Ω with set Ω containing the indices of the active (data and pilot) subcarriers, we model
the received frequency-domain (FD) signal as follows:

yω = Hω sω + nω .
1

(1)

By assuming channel reciprocity [1], the proposed methods directly apply to the downlink, in which HT is the downlink
channel matrix. Once all of the Gram matrices have been computed, they can be re-used in the downlink for linear (e.g.,
zero-forcing or minimum-mean square error) precoding. See Section II-C for the details.

5

Here, Hω ∈ CB×U is the FD channel matrix, sω ∈ CU is the transmit vector, and nω ∈ CB models

additive noise. The FD input-output relation in (1) is able to model both OFDM and SC-FDMA systems.
For OFDM systems, the entries of the transmit vector sω are taken from a discrete constellation set O

(e.g., 16-QAM); in SC-FDMA systems, the constellation points are assigned in the time-domain and the
resulting vectors are transformed into the FD to obtain the transmit vectors sω . See, e.g., [9], for more
details on SC-FDMA transmission.
B. Wideband Channel Model
In wideband MIMO multicarrier systems, the FD channel matrices Hω , ω = 0, . . . , W − 1, are directly

b ` ∈ CB×U , where ` = 0, . . . , W − 1 are the channel “taps” in
related to the time-domain (TD) matrices H
the TD. We first introduce the model used to characterize the presence of antenna correlation at the BS

side2 which occurs in the TD. Specifically, we use the standard correlation model from [4] and express
the `th TD channel matrix as follows:
b ` = R1/2 H
b uncor .
H
`
`

(2)

b uncor represents an uncorrelated TD channel matrix and R` ∈ CB×B is a correlation matrix that
Here, H
`
contains ones on the main diagonal and δ ∈ R on the off-diagonals. We allow δ to be either real positive

or negative but it must satisfy δ 2 ≤ 1. We rewrite the correlation matrix as R` = (1 − δ)IB + δ1B ,

where IB and 1B is the B × B identity and all-ones matrix, respectively. We note that in the absence of
b` = H
b uncor .
receive-side correlation, i.e., δ = 0, we have R = IB and H
`

In order to take into account the practically-relevant case of imperfect CSI at the BS, we assume
b ` ∈ CB×U ,
that the FD channel matrices Hω , ω = 0, . . . , W − 1, are obtained from the TD matrices H
` = 0, . . . , W − 1, via the discrete Fourier transform [26] as follows:3
Hω =

L−1
X
`=0



b ` exp − j2πω`
H
W


+ σEω ,

(3)

where the matrix Eω ∈ CB×U models channel estimation error on subcarrier ω and the parameter σ ∈ R+
determines the intensity of channel-estimation errors; σ = 0 corresponds to the case for perfect CSI. We
assume that the entries of the matrix Eω are i.i.d. (across entries and subcarriers) circularly-symmetric
2

In massive MU-MIMO systems, the users are likely to be uncorrelated as they are spatially well-separated over potentially
large cells; in contrast, the antennas at the BS are typically confined to a comparably small area, which increases the potential
for receive-side correlation.
3

Although one could extract refined channel estimates Hω by exploiting the fact that only the first L taps of Hω are active in
the TD, we do not further pursue this direction. An analysis of such channel estimation algorithms is left for future work.

6

complex Gaussian with unit variance. Equation (3) relies on the assumption that at most L ≤ W of

the first channel taps are non-zero (or dominant) and the remaining ones are zero (or insignificant), i.e.,
b ` = 0B×U for ` = L, . . . , W − 1. In practical OFDM and SC-FDMA systems, the maximum number
H

of non-zero channel taps should not exceed the cyclic prefix length (assuming perfect synchronization).
Hence, we can safely assume that L  W in practical scenarios and for most standards, such as IEEE
802.11n [17] or 3GPP-LTE [18].

C. Linear Data Detection and Precoding
In the uplink, linear minimum mean-square error (MMSE) equalization was shown to achieve nearoptimal error-rate performance in massive MU-MIMO systems [5]. The most common approach to linear
MMSE equalization first forms the Gram matrix Gω = HH
ω Hω and then, computes an estimate of
H
0 −1
the transmit vector as ŝω = (Gω + I N
Es ) Hω yω , where N0 and Es stand for the noise variance and

average energy per transmit symbol, respectively. Directly computing the Gram matrix Gω for every
subcarrier results in excessively high complexity. In fact, even by exploiting symmetries, 2BU 2 real-valued
multiplications are required, which is more than 16 k multiplications per subcarrier for a system with 128
BS antennas and 8 users. The hardware design for massive MU-MIMO data detection in [9] confirms
this observation and shows that computing the Gram matrix dominates the overall hardware complexity
and power consumption.
In the downlink, ZF precoding is the most common linear precoding method [1], which computes
xω = H∗ G−1
ω sω , where xω is the B -dimensional transmit signal and sω the data vector. If the Gram

matrix Gω has been precomputed for equalization in the uplink phase, then it can be re-used for ZF
precoding in the downlink phase to minimize recurrent operations. Hence, in order to minimize the overall
complexity of equalization and precoding, efficient ways to compute the Gram matrix Gω on all active
subcarriers ω ∈ Ω are required.
III. I NTERPOLATION - BASED G RAM M ATRIX C OMPUTATION
We now discuss exact and approximate interpolation-based methods for low-complexity Gram matrix
computation. We note that the exact Gram matrix interpolation assumes that we have perfect CSI, so
throughout this section, we will assume that σ = 0. In Section IV, however, we will relax the perfect-CSI
assumption and study the performance of exact and approximate interpolation methods with imperfect
CSI. As a consequence of (3), the Gram matrices in the FD are given by
Gω =

L−1
X L−1
X
0

`=0 ` =0



j2πω(` − `0 )
Hb
b
H` H`0 exp
W

(4)

7

for w = 0, 1, . . . , W − 1. Given the FD channel matrices Hω for all active subcarriers ω ∈ Ω, a

straightforward “brute-force” approach simply computes Gω = HH
ω Hω for each active subcarrier ω ∈ Ω.
In order to reduce the complexity of such a brute-force approach, we next discuss exact and approximate
Gram-matrix interpolation methods that take advantage of the facts that (i) the channel matrices (and
hence, the Gram matrices) are “smooth” (or correlated) across subcarriers if L < W and (ii) massive
MU-MIMO causes an effect known as channel hardening [1], [2].
A. Exact Gram-Matrix Interpolation
The Gram matrix Gω in (4) is a Laurent polynomial matrix in the variable xω = exp(j2πω/W ); we
refer the reader to [22] for more details on Laurent polynomial matrices. Hence, we can establish the
following result for exact Gram-matrix interpolation; a short proof is given in Appendix A.
Lemma 1. The Gram matrices Gω in (4) for all subcarriers ω = 0, . . . , W are fully determined from
2L − 1 distinct and non-zero Gram-matrix base-points.

As an immediate consequence of this lemma, one can interpolate all of the Gram matrices exactly from
only 2L − 1 distinct and non-zero Gram-matrix base-points that have been computed explicitly.

In order to perform exact interpolation, we first define a set of base points P ⊂ Ω that contains

|P| ≥ 2L−1 distinct subcarrier indices. We denote the k th base-point index as pk , where k = 0, . . . , |P|−1,

and the set of all base-point indices as P = {p0 , . . . , p|P|−1 }. For each subcarrier index in the base-point

set P , we then explicitly compute |P| ≥ 2L − 1 Gram matrices Gω = HH
ω Hω , ω ∈ P , and perform
entry-wise interpolation for the gram matrices Gω on all remaining active subcarriers ω ∈ Ω\P .

The exact interpolation procedure for each entry is as follows. For a fixed entry (m, n), we define the

vector gP ∈ C|P| , which is constructed from the (m, n) entries [Gω ]m,n taken from base-points ω ∈ P ,

T
i.e., gP = [Gp0 ]m,n · · · [Gp|P|−1 ]m,n . Then, the vector gΩ\P ∈ C|Ω|−|P| that contains the entry (m, n)
for all remaining Gram matrices Gω , ω ∈ Ω\P , is given by

gΩ\P = FΩ\P,L (F†P,L gP ).

(5)

Here, FP,L represents a |P| × (2L − 1) matrix where we take the |P| rows indexed by P and the first

L and last L − 1 columns from the W -point discrete Fourier transform (DFT) matrix; the entries of


the DFT matrix are defined as [F]m,n = √1W exp − j2π
(m
−
1)(n
−
1)
. Similarly, FΩ\P,L represents
W

a (|Ω| − |P|) × (2L − 1) matrix where we take |Ω| − |P| rows indexed by Ω\P and the first L and
last L − 1 columns from the W -point DFT matrix. In words, the exact interpolation method in (5) first

computes the 2L − 1 TD Gram-matrix entries and then, transforms these elements into the frequency

8

10

[Gω ]1,1

9.5

9

8.5

8
300

exact Gram matrices
0th order interpolation
1st order interpolation

310

320

subcarrier index ω

330

340

Fig. 1. Illustration of 0th and 1st order interpolation for the entry [Gω ]1,1 across subcarriers. We explicitly compute the Gram
matrix for every fifth subcarrier index ω and interpolate the remaining matrices. We assume a 128 BS antenna, 8 user massive
MU-MIMO system with W = 2048 subcarriers, a delay spread of L = 144, and perfect CSI without BS-antenna correlation.

domain via the DFT. See [19]–[22] for additional details on other exact interpolation methods developed
for small-scale, point-to-point MIMO systems.
Although the method in (5) is able to exactly interpolate the Gram matrix across all W subcarriers, it is
in many situations not practical due to the high complexity of the matrix inversion required in F†P,L gP . If,
however, one can sample the base points uniformly over all W tones, the complexity of matrix inversion
can be reduced significantly. Unfortunately, this approach is often infeasible in practice due to the presence
of guard-band constraints in OFDM-based or SC-FDMA-based standards [17], [18]. Another issue of
exact interpolation methods, such as the ones in [19]–[22] and ours in (5), is that they generally assume
perfect CSI and no BS-antenna correlation. As we will show in Section VI-C, imperfect CSI results in
poor interpolation performance—this is due to the fact that the matrix FP,L is typically ill conditioned,
especially when sampling Gram-matrices close to the minimum number of 2L − 1 required base points.
We next propose two approximate interpolation schemes that not only require (often significantly) lower

complexity than a brute-force approach or exact interpolation in (5), but also approach the performance
of a brute-force approach in massive MU-MIMO systems and are robust to channel-estimation errors.

9

B. Approximate Gram-Matrix Interpolation
We consider the following two approximate Gram-matrix interpolation methods illustrated in Fig. 1.
0th Order Interpolation: We select a set of |P| distinct base-points with P = {p0 , . . . , p|P|−1 } ⊂ Ω.

We explicitly compute Gp = HH
p Hp on these base points and perform 0th order (or nearest-neighbor)
interpolation for the remaining subcarriers in the set Ω\P according to:
e ω = Gp ,
G

p = arg min |p̃ − ω|,
p̃∈P

∀ω ∈ Ω\P.

(6)

e ω equal to the nearest Gram matrix that has been
In words, we set the interpolated Gram matrix G
computed explicitly for one of the neighboring base points.
1st Order Interpolation: Analogously to the 0th order interpolation method, we explicitly compute
Gp = HH
p Hp on a selected set of base-points p ∈ P . Then, for each target subcarrier ω ∈ Ω\P we pick

two nearest base-points pk and pk+1 , i.e., pk ≤ ω ≤ pk+1 , and perform entry-wise linear interpolation
according to

e ω = λω Gpk + (1 − λω )Gpk+1 ,
G

ω ∈ Ω\P,

(7)

where λω = (pk+1 − ω)/(pk+1 − pk ) and pk < pk+1 .
IV. A PPROXIMATION E RROR A NALYSIS
We now analyze the approximation error associated with the approximate interpolation schemes proposed
eω
in Section III-B. We use Gω to represent the Gram matrices that have been computed exactly and G
to represent the Gram matrices that have been obtained via approximate interpolation. We note that the
exact interpolation scheme in Section III-A exhibits no approximation error.
A. Mean-Square-Error of Approximate Interpolation
We study the mean-squared error (MSE) on each entry (m, n) for the ω -th subcarrier, which we define
as follows:
h
i
2
e
o-MSE(m,n)
,
E
[
G
]
−
[G
]
.
ω m,n
ω m,n
ω

(8)

Here, o represents the order of interpolation, i.e., we have either o = 0 or o = 1. Our results make
extensive use of the scaled Fejér kernel [27] given by
fL (φ) = L−2

1 − cos(Lφ)
sin2 (Lφ/2)
= L−2
1 − cos(φ)
sin2 (φ/2)

(9)

10

and rely on the following key properties; the proof is given in Appendix B.
Lemma 2. The scaled Fejér kernel (9) is non-negative, bounded from above by one, and monotonically
decreasing in φ for φ ∈ [0, 2π/L] with L > 1.
B. MSE of 0th Order Interpolation
The following result precisely characterizes the MSE of 0th order interpolation for imperfect CSI as
modeled in (3) and BS-antenna correlation as modeled in (2); the proof is given in Appendix C.
b uncor , ` = 1, . . . , L, be distributed CN (0, 1/(BL)) per
Theorem 1. Let the entries of the TD matrices H
`

complex entry. Assume that the off-diagonal of the receive correlation matrix be δ , the variance of the
channel estimation error to be σ , and p ∈ P is the closest base point to the target subcarrier ω . Then,
for any (m, n) entry of the Gram matrix Gω , the MSE for the 0th order interpolation method in (6) is
given by



2π
2
(p − ω) ,
0-MSEω = εCSI + (1 + εcor ) 1 − fL
B
W

(10)

where we use the definitions
εCSI = 2σ 2 (2 + Bσ 2 )

and εcor = δ 2 (B − 1).

From this result, we observe that, as the number of BS antennas B increases, εCSI increases quadratically
with respect to σ 2 . In addition, for perfect CSI, i.e., σ = 0 so εCSI = 0, the MSE for 0th order interpolation
decreases with an increasing number of BS antennas B as

∂
∂B 0-MSEω

< 0, if δ 2 < 1. Also, we note in

the case for non-zero correlation, i.e., δ 6= 0, the MSE for 0th order interpolation is amplified (compared

to that with no correlation) by a factor of 1 + δ 2 (B − 1). Furthermore, we observe that the MSE is

independent of the entry of the Gram matrix (i.e., the MSE is identical for diagonal as well as off-diagonal
b uncor .
entries); this is a result of the i.i.d. assumption of the TD channel matrices H
`
To gain insight into the behavior of 0th order interpolation in the large BS-antenna limit, we have the
following result.
Corollary 1. Assume the conditions in Theorem 1 and εCSI → 0 as B → ∞. Then, as B → ∞, MSE of
0th order interpolation is given by

lim 0-MSEω = 2δ

B→∞

2




1 − fL

2π
(p − ω)
W


.

11

Corollary 1 demonstrates that in the large-BS antenna limit, the MSE of 0th order interpolation is zero
across all subcarriers if and only if the BS antennas are uncorrelated, i.e., δ = 0. For δ 6= 0, the MSE
depends on the distance between the nearest base-point and the target subcarrier.

While the MSE expression in (10) is exact, it does not provide much intuitive insights. We define the
following quantity that enables us to further analyze the MSE in (10).
Definition 1. The maximum distance between any subcarrier ω and the nearest base point is given by:


pk+1 − pk
dmax = max
.
(11)
k∈P
2
With the maximum distance dmax for a given set of base points P , Corollary 2 shows that the 0th

order approximation can be bounded from above using simple analytic expressions; the proof is given in
Appendix D.

Corollary 2. Let dmax be the maximum distance in (11) and assume the conditions in Theorem 1 hold.
Then, the maximum MSE of 0th order interpolation over all active subcarriers ω ∈ Ω is bounded by:
max{0-MSEω }
ω∈Ω

2


εCSI + B (1 + εcor ),



≤
2(1 + εcor )
2π


1 − fL
dmax ,
εCSI +
B
W

W
L
W
<
.
L

dmax ≥
dmax

Corollary 2 implies that regardless of small or large maximum distance dmax , the MSE given by the 0th
order approximation always decreases with the number of BS antennas B only if δ = 0 and σ 2 = 0 (see
also Theorem 1). In addition, if the distance between the interpolated subcarrier index ω and its closest
base point is sufficiently small, i.e., dmax < W/L, then we obtain a sharper upper bound on the MSE
than εCSI +

2
B (1

+ εcor ). In a scenario with a large delay spread L, Corollary 2 reveals that one requires

finer-spaced base points for 0th order interpolation in order to keep the approximation error strictly smaller
than εCSI +

2
B (1

+ εcor ). In addition, since the maximum error is mainly determined by dmax , a good

strategy for selecting base points with 0th order approximation is uniformly spacing them in the set of
active subcarriers Ω.
C. MSE of 1st Order Interpolation
We now present the approximation error analysis of 1st order interpolation. The following result
precisely characterizes the MSE of 1st order interpolation; the proof is given in Appendix E.

12

b uncor , ` = 1, . . . , L, be distributed CN (0, 1/(BL)) per
Theorem 2. Let the entries of the TD matrices H
`

complex entry. Assume that the off-diagonal of the receive correlation matrix be δ , the variance of the
channel estimation error to be σ across all subcarriers ω = 0, . . . ,, and p ∈ P is the closest base point

to the target subcarrier ω . Then, for any (m, n)-th entry of the Gram matrix Gω , the MSE for the 1st
order interpolation method in (6) is given by
1-MSEω = εCSI (1 − λω (1 − λω ))

2
(1 + εcor ) 1 − λω (1 − λω ) + λω (1 − λω )fL (θ)
B

− (1 − λω )fL (λω θ) − λω fL ((1 − λω )θ) ,

+

where θ =

2π
W (pk+1

(12)

− pk ) and λω = (pk+1 − ω)/(pk+1 − pk ).

Analogously to 0th order interpolation, we observe that the MSE of 1st order interpolation is independent
of the entry (m, n) and impacted by CSI errors and receive correlation (see Section IV-B for detailed
discussion). The result shown next in Corollary 3 reveals that if the spacing between the two base-points
pk and pk+1 defined as dk = pk+1 − pk is sufficiently small, then the 1st order interpolation strictly

outperforms 0th order interpolation, i.e., 1-MSEω < 0-MSEω for all ω ∈ (pk , pk+1 ); the proof is given in

Appendix F.

Corollary 3. Let dk denote the spacing between two base-points pk and pk+1 , and assume the conditions
in Theorem 1 and Theorem 2 hold. If dk ≤ W/(3L), then
1-MSEω ≤ 0-MSEω ,

for all ω ∈ (pk , pk+1 ),

(13)

which holds with equality if and only if L = 1 and εCSI = 0.
We note that the condition dk < W/(3L) is not sharp; Appendix F outlines the details on how it can
be sharpened. Furthermore, given that dk is significantly larger than W/(3L), we can construct situations
for which 0th order interpolation outperforms 1st order interpolation. Note that for L = 1, the FD channel
is flat (i.e., Gω is constant for all ω ) and hence, 1st and 0th order interpolation have the same MSE.
In summary, we see that for both approximate interpolation methods, we can lower the MSE by
increasing the number of BS antennas B assuming that the channel estimation error εCSI decreases with
respect to B . In the large-antenna limit B → ∞ with perfect CSI and no BS-antenna correlation, the MSE

vanishes, which is an immediate consequence of channel hardening in massive MU-MIMO systems [1].
Furthermore, 1st order interpolation generally outperforms 0th order interpolation for a sufficiently small

13

minimum spacing between adjacent base points, i.e., for dk ≤ W/(3L).
V. C OMPLEXITY A NALYSIS
We next compare the computational complexity of the four studied Gram-matrix computation algorithms:
brute-force computation, exact interpolation, 0th order interpolation, and 1st order interpolation. We
measure the computational complexity by counting the number of real-valued multiplications.4
A. Brute-Force Computation
We start by deriving the total computational complexity required by the brute-force (BF) method. We
only compute the upper triangular part of Gω (since the matrix is Hermitian). Each off-diagonal entry
requires B complex-valued multiplications, which corresponds to 4B real-valued multiplications; each
diagonal entry requires only 2B real-valued multiplications. Hence, the computational complexity of
computing Gω using the BF method is


U (U − 1)
CBF = |Ω| 4B
+ 2BU = 2|Ω|BU 2
2

(14)

for a total number of |Ω| active subcarriers.
B. Exact Interpolation
We now derive the computational complexity of exact interpolation as discussed in Section III-A. Exact
interpolation requires a BF computation of the Gram matrix at each of the |P| base points. We will use
the |P| precomputed base points of Gω to interpolate the remaining |Ω| − |P| Gram matrices.

We will assume that the base points and the assumed channel delay spread L are fixed a-priori

so that FΩ\P,L F†P,L in (5) can be precomputed and stored. We emphasize that this approach does
not include the computational complexity of computing the interpolation matrix itself, which favors
this particular interpolation scheme from a complexity perspective. In fact, we only need to multiply
the precomputed interpolation matrix FΩ\P,L F†P,L with the vector gP , which requires 4(|Ω| − |P|)|P|
real-valued multiplications. Hence, the total computational complexity of exact interpolation is:
CExact =

|P|
U (U + 1)
CBF + 4(|Ω| − |P|)|P|
|Ω|
2

= 2|P|(|Ω| − |P| + B)U 2 + 2|P|(|Ω| − |P|)U.

(15)

4
We assume that a complex-valued multiplication requires four real-valued multiplications; computation of the squared
magnitude of a complex number is assumed to require two real-valued multiplications.

14

We note that if the number of users U is large and the number of base points is similar to the number of
BS antennas, i.e., |P| ' B , then the BF method in (14) and exact interpolation (15) exhibit a similar
complexity. We also observe that the complexity of exact interpolation (15) is lower than that given by the
BF method (14) if |P| < (1 + U )−1 BU . Since the use of |P| ≥ 2L − 1 distinct base points guarantees

exact interpolation (assuming perfect CSI), we observe that exact interpolation has lower complexity than
the BF method if L is (approximately) less than B/2.
C. 0th Order Interpolation
The computational complexity of the 0th order interpolation method is given by
C0th = 2|P|BU 2 ,

(16)

as we only need to compute the Gram matrices on all the base points. We note that since typically
|P|  |Ω| the savings (in terms of real-valued multiplications) are significant compared to the BF approach

and exact interpolation, but does so at the cost of approximation errors (cf. Section VI-D).
D. 1st Order Interpolation
The computational complexity of the 1st order interpolation is given by
C1st = C0th + 4(|Ω| − |P|)

U (U + 1)
2

= 2|P|BU 2 + 2(|Ω| − |P|)U (U + 1),

(17)

where we assume that the interpolation weight λω was precomputed. We note that the linear interpolation
stage for each subcarrier ω ∈ Ω\P requires four real-valued multiplications. By comparing (16) to (17),

we observe that the complexity of 1st order interpolation always exceeds the complexity of the 0th
order method, but the complexity is significantly lower than that of the BF method as we generally have
|P|  |Ω|.

VI. N UMERICAL R ESULTS
In this section, we study the MSE, the error-rate performance, and the computational complexity of the
proposed Gram-matrix interpolation schemes. We consider a MU-MIMO-OFDM system with 128 BS
antennas and with 8 single-antenna users. We assume a total of W = 2048 subcarriers, with |Ω| = 1200
active subcarriers, similar to that used in 3GPP LTE [18]. Unless stated otherwise, we assume that the

entries of the TD channel matrices are i.i.d. circularly-symmetric complex Gaussian with variance 1/(BL)

15

and we consider 16-QAM transmission (with Gray mapping). We use the linear MMSE equalizer as
outlined in Section II-C for data detection. For situations with imperfect CSI, we consider training-based
maximum-likelihood (ML) channel estimation, where we use a single orthogonal training sequence of
length U with the same transmit and noise power as for the data symbols.
A. Complexity Comparison
We now assess the complexity of the various Gram-matrix computation methods in comparison to
the overall complexity required for linear MMSE-based data detection, which includes Gram-matrix and
matched-filter computation as well as matrix inversion for each active subcarrier. The results shown here
are for a 128 × 8 (the notation represents B × U ) massive MU-OFDM-MIMO system with 1200 active
subcarriers and a delay spread of L = 144.

Fig. 2 compares the complexity of Gram matrix computation for four different methods, brute-force,
exact, 0th-, and 1st-order interpolation methods for |P| = {0.25|Ω| , 0.5|Ω| , 0.75|Ω| , |Ω|}. The solid part

of the bar plot shows the complexity of Gram matrix computation; the fenced part corresponds to the
remaining complexity required for data detection (including matched-filter computation and a matrix
inversion for each active subcarrier). The percentage values indicate the relative complexity of Gram-matrix
computation compared to the total complexity required for data detection.
We see that the exact interpolation method incurs significant complexity in the simulated system (see
Section V-B for exact details when exact interpolation achieves lower complexity than a BF approach). We
also see that the proposed 0th and 1st order approximation methods both achieve significant complexity
reductions. For |P| = 0.25|Ω|, the proposed methods requires less than half the complexity of a BF
approach. As we will show in Section VI-C, the proposed approximate interpolation methods will exhibit

similar error-rate performance as that the BF approach (see Figs. 4 and 5), but does so at fraction of the
computational complexity.
B. MSE of Approximate Interpolation
Fig. 3 compares the MSE of 0th and 1st order interpolation as proposed in Section III-B. Note that
the BF method and exact interpolation have an MSE of zero and hence, we exclude these results. We
select two base points at subcarriers 500 and 600 and one target point at subcarrier 512, and compare the
MSE for different numbers of BS antennas and under ideal and non-ideal scenarios. In the ideal scenario,
we assume perfect CSI and no BS-antenna correlation, whereas in the non-ideal scenario we assume
channel-estimation at SNR = 25 dB across all subcarriers with the signal-to-noise ratio (SNR) defined by
SNR = U/(Bσ 2 ) and a BS-antenna correlation of δ = 0.1. In order to assess the approximation error

number of real-valued multiplications

16

8

·107

7

brute-force
exact
0th order
1st order

92%
91%

6
89%

5
4
3

78%

78%

0

78% 78% 78% 78%
73% 72%

2
1

78%
65% 64%

49% 47%

0.25

0.5

0.75

1

fraction of |Ω| base-points
Fig. 2. We compare the complexity of various Gram matrix computation methods in comparison to the complexity of
data detection for a 128 × 8 MU-OFDM-MIMO system with |Ω| = 1200, L = 144, and four different sets of base-points
|P| = {0.25|Ω| , 0.5|Ω| , 0.75|Ω| , |Ω|}. The percentage values shown in the above bar plots show the relative percentage of
Gram-matrix computation compared to the total complexity required for data detection.

with respect to different channel delay spreads, we set L ∈ {36, 72, 144}. The resulting MSE is shown in

Fig. 3. Note that the MSE for both 0th and 1st are independent of the entry (as predicted by Theorems 1
and 2); hence, we consider the average MSE across all entries.
We observe that the 1st order interpolation method achieves a lower MSE than that given by 0th order
interpolation, where the performance gap increases with larger delay spreads L. This observation can
intuitively be explained by the fact that for small delay spreads L, the channel is more smooth across
subcarriers. For larger delay spreads L, however, 1st order interpolation captures the faster-changing
behavior of the Gram matrix, whereas the 0th order interpolation ignores such changes. We also see
that the MSE degrades in the non-ideal scenario, even if we increase the number of BS antennas; this
behavior can also be observed in our analytical results. Finally, we observe that the simulated MSE
matches perfectly our theoretical results provided in Theorems 1 and 2.

17

mean-squared error (MSE)

100
L = 144
10−1

0th order, ideal (sim)
1st order, ideal (sim)
0th order, non-ideal (sim)
1st order, non-ideal (sim)

10−2
L = 72
10−3 3
2

24

L = 36
25

26

base-station antennas B

27

28

Fig. 3. MSE of 0th and 1st order interpolation for an entry of G512 using two base points G500 and G600 for three different
delay spreads L ∈ {36, 72, 144} for the ideal and non-ideal scenarios. The markers represent simulation results whereas the
lines represent our approximation-error analysis. Evidently, our theory matches perfectly with the simulated values.

C. Error-rate Performance
We now compare the error-rate performance of the proposed Gram-matrix computation schemes. We
simulate the bit-error rate (BER) for a MU-MIMO-OFDM system for a different number of base-points |P|

and for perfect as well as imperfect CSI. We also investigate the impact of a more realistic channel
model. For all results, we simulate three different numbers of base-points |P| = L/4, |P| = 2L − 1, and

|P| = 4L, and select equally-spaced base points. Figures 4(a) and 4(b) show BER simulation results for

an i.i.d. Rayleigh fading scenario with perfect and imperfect CSI, respectively. Figure 4(c) shows BER
simulation results for the QuaDRiGa channel model5 with imperfect CSI [28].

Figure 4(a) shows that exact interpolation for |P| ≥ 2L − 1 base points provides identical results as

the BF method (up to machine precision) for a system with perfect CSI. For |P| = L/4 base points, the
proposed 0th and 1st order interpolation exhibit an error floor; this performance loss can be mitigated

substantially by increasing the number of base points to |P| = 2L − 1. By setting |P| = 4L < 0.5|Ω|,
the both the 0th and 1st order interpolation methods exhibit virtually no BER performance loss.

5
We simulate a square antenna array with a non-line-of-sight scenario with a 2 GHz carrier frequency, 20 MHz bandwidth, and
200 m distance between BS antenna and the users. Our algorithms assume L = 144 but the true delay spread is slightly smaller.

18

bit-error rate (BER)

100

10−1

10−2

brute-force
0th order, L/4
1st order, L/4
0th order, 2L−1
1st order, 2L−1
0th order, 4L
1st order, 4L

10−3

10−4 0

4

8

12

16

average SNR per receive antenna [dB]

20

(a) Perfect CSI (i.i.d. Rayleigh fading).

bit-error rate (BER)

100

10−1

10−2

10−3

10−4 0

brute-force
0th order, L/4
1st order, L/4
exact, L/4
0th order, 2L−1
1st order, 2L−1
exact, 2L−1
0th order, 4L
1st order, 4L
exact, 4L

4

8

12

16

average SNR per receive antenna [dB]

20

(b) Imperfect CSI (i.i.d. Rayleigh fading).

bit-error rate (BER)

100

10−1

10−2

10−3

10−4 0

brute-force
0th order, L/4
1st order, L/4
exact, L/4
0th order, 2L−1
1st order, 2L−1
exact, 2L−1
0th order, 4L
1st order, 4L
exact, 4L

4

8

12

16

average SNR per receive antenna [dB]

20

(c) Imperfect CSI (QuaDRiGa channel model).
Fig. 4. Uncoded bit error-rate (BER) comparison in a 128 BS antenna, 8 (single-antenna) user, wideband massive MU-MIMOOFDM system. The values next to the legend entries correspond to the number of base points |P|. The proposed approximate
interpolation schemes (0th and 1st order interpolation) outperform exact interpolation for scenarios with imperfect CSI and
approach the performance of the exact brute-force method for a small number of base points.

19

fraction of brute-force complexity CBF

1

brute-force
0th order interpolation
1st order interpolation

0.8
0.6

SNR = 32.56 dB

1 dB

0.4
0.2
0

9

10

SNR = 9.51 dB
11
12
13

14

15

minimum SNR [dB] for 10−3 BER
Fig. 5. Trade-off between SNR performance and computational complexity for 0th and 1st order interpolation in an i.i.d. Rayleigh
fading channel with imperfect CSI. Both approximate interpolation methods approach the performance of an exact brute-force
approach at a fraction of the complexity.

Figure 4(b) shows the situation for imperfect CSI (with channel estimation). We observe that the
performance of the BF method and that of exact interpolation are no longer equal. In fact, for |P| = 2L − 1
and |P| = 4L base points, exact interpolation exhibits a significant error floor. The reason is due to the fact
that the interpolation matrix is extremely ill-conditioned, which results in significant noise enhancement

artifacts. Although the error floor is decreased for |P| = 4L, a floor remains at 10−3 BER. In contrast,
the error floor of 0th order and 1st order interpolation for |P| = 2L − 1 and |P| = 4L base points is

well-below 10−4 BER and hence, the proposed approximate interpolation schemes are more resilient to
scenarios with imperfect CSI than exact interpolation.
Figure 4(c) shows the BER performance for the QuaDRiGa channel model with imperfect CSI. We

observe that all considered interpolation methods achieve a lower error floor than that given in Fig. 4(b)
for |P| = L/4; this is due to the fact that the effective delay spread for the considered channel is smaller
than L = 144 (which is assumed in the algorithms). Once again, we observe a strong BER floor of exact

interpolation for all considered numbers of base points. In summary, we see that the proposed approximate
interpolation methods are more robust in practical scenarios than the exact interpolation method.

20

D. Performance/Complexity Trade-off
We now investigate the BER performance vs. computational complexity trade-off for the proposed
approximate interpolation methods with imperfect CSI. We use the complexity CBF of the BF method in
(14) as our baseline, and we compare it to that of the proposed 0th and 1st order interpolation methods
in (16) and (17), respectively. We vary the number of base points |P| from L to D and simulate the
minimum SNR required for the linear MMSE equalizer to achieve 10−3 BER.

Figure 5 shows the trade-off results for 0th and 1st order interpolation. For a fixed fraction of the
complexity of CBF , we observe that the 1st order interpolation method always outperforms the 0th order
interpolation method. Hence, Fig. 5 clearly reveals that the additional complexity required by linear
interpolation is beneficial from a performance and complexity perspective. In addition, we see that the 1st
order interpolation method is able to approach the SNR performance of the BF method by 1 dB with only
45% of the complexity.

VII. C ONCLUSIONS
We have studied the performance of exact and approximate interpolation-based Gram matrix computation
for wideband massive MU-MIMO-OFDM systems. Instead of performing a brute-force (BF) computation
of the Gram matrix for all subcarriers or using exact interpolation schemes, we have proposed two simple,
yet efficient approximate interpolation methods. We have demonstrated that channel hardening in massive
MU-MIMO enables the proposed 0th and 1st order interpolation schemes to perform close to that of
an exact BF computation at only a fraction of the computational complexity. In addition, the proposed
approximate interpolation methods are more robust to channel-estimation errors and receive-side antenna
correlation than exact interpolation methods.
There are many avenues for future work. We expect the use of higher-order approximate Gram-matrix
interpolation schemes to perform better at an increase in complexity. An analysis of such methods is left
for future work. Our results also indicate that the broad range of existing exact interpolation schemes for
small-scale, point-to-point MIMO systems, e.g., matrix inversion and QR decomposition [19], [20] can
be made more robust and less complex in massive MU-MIMO systems if combined with approximate,
low-order interpolation schemes. In fact, the recent result in [25] for approximate interpolation of matrix
inversion demonstrates this claim in a massive MU-MIMO scenario via simulations. In addition, a
performance analysis of the proposed or other approximate interpolation methods in more realistic systems
that suffer from frequency and timing offset is left for future work.

21

ACKNOWLEDGMENTS
The authors would like to thank M. Wu and C. Dick from Xilinx, Inc. for insightful discussions on
interpolation-based matrix computations. The authors would also like to thank the anonymous reviewers
for their suggestions. The work of C. Jeon and C. Studer was supported in part by Xilinx Inc., and by
the US National Science Foundation (NSF) under grants ECCS-1408006, CCF-1535897, and CAREER
CCF-1652065.
A PPENDIX A
P ROOF OF L EMMA 1
Since all of the possible values in the exponent of (4), i.e., ` − `0 , are integers ranging from −(L − 1)

to L − 1, the Gram matrix Gω in (4) is a polynomial with degree no larger than 2L − 1. Consequently,
the Gram matrices Gω in (4) for all subcarriers ω = 0, . . . , W are fully determined from 2L − 1 distinct
and non-zero Gram-matrix base-points.

A PPENDIX B
P ROOF OF L EMMA 2
Evidently, (9) is non-negative, i.e., fL (φ) ≥ 0. To show that fL (φ) ≤ 1, we use the fact that the

Fejér kernel LfL (φ) = L−1 1−cos(Lφ)
1−cos(φ) is upper bounded by L [27, Eq. 1.2.24]. Now, we show that fL (φ)
is monotonically decreasing in φ ∈ [0, 2π/L] for L > 1. We start by defining an auxiliary function

g(φ) = sin(Lφ/2)/ sin(φ/2), so that fL (φ) = gL2 (φ)/L2 . For L > 1, φ ∈ (π/L, 2π/L], sin(Lφ/2)

and sin(φ/2) are monotonically decreasing and increasing respectively, and hence, g(φ) and fL (φ), are
monotonically decreasing. For φ ∈ [0, π/L], the derivative of g(φ) with respect to φ is given by:
dg(φ)
(L cos(Lφ/2) sin(φ/2) − cos(φ/2) sin(Lφ/2))
=
,
dφ
2 sin2 (φ/2)

and we have that
L cos(Lφ/2) sin(φ/2)
L tan(φ/2)
=
< 1.
cos(φ/2) sin(Lφ/2)
tan(Lφ/2)

Hence, g 0 (φ) < 0, and therefore, fL (φ) is monotonically decreasing in φ ∈ [0, 2π/L].

22

A PPENDIX C
P ROOF OF T HEOREM 1
Suppose we use Gpk at base point pk to approximate Gw at the target subcarrier index ω . Hence, the
MSE in (8) is given by the following expression:
h
i
e ω ]m,n − [Gω ]m,n 2 .
0-MSE(m,n)
= E [G
ω

(18)

We will obtain an analytical expression for (18) with imperfect CSI and the BS-antenna correlation model
introduced in (3) and (2), respectively. We start with expressing the received channel Hω from (3) as
(a)

Hω =

W
−1
X
`=0




σ
j2πω`
b
H` I(` < L) + √ Ẽ` exp −
,
W
W

(19)

where (a) follows from noting that the DFT is orthogonal (with normalization constant) with each entries
of Ẽ` distributed CN (0, 1). Hence, the received TD channel matrix under imperfect CSI is given as
bσ = H
b ` I(` < L) +
H
`

√σ Ẽ` .
W

Now, we start with the correlation model introduced in (2). We note that the B × B BS correlation
√
matrix R` = (1 − δ)IB + δ1B can be expressed by (αIB + β1B )2 = R` where α = 1 − δ , and
p
b ` is
β = B1 ( 1 + (B − 1)δ − α). Since the matrix (αIB + β1B ) is symmetric, we note that the H
b ` = (αIB + β1B )H
b uncor so that
expressed as H
`

b ` ]m,n = α[H
b uncor ]m,n + β
[H
`

B
X
b uncor ]b,n
[H
`

(20)

b=1

For our derivation of the MSE, we will utilize the following auxiliary function:




2πω
2πpk
0
0
(` − ` ) − exp j
(` − ` ) .
R``0 = exp j
W
W

(21)

e ω ]m,n − [Gω ]m,n in (18), we obtain the following expression for the 0th order
By substituting R``0 into [G
MSE:

0-MSE(m,n)
ω

 W
−1 X
B
2
X
σ ∗
σ
b
b
=E
[H` ]b,m [H`0 ]b,n R``0
`,`0 =0 b=1

=

W
−1
X

W
−1
X

B 
X

R`∗1 `2 R`3 `4

`1 ,`2 =0 `3 ,`4 =0 b1 ,b2 =1


 σ

σ ∗
σ ∗
σ
b
b
b
b
E [H`1 ]b1 ,m [H`2 ]b1 ,n [H`3 ]b2 ,m [H`4 ]b2 ,n

23

(a)

=

W
−1
X
`1 ,`2 =0
`1 6=`2

B
X

|R`1 `2 |2

 σ ∗
  σ ∗

b ]
bσ
b
bσ
E [H
`1 b2 ,m [H`1 ]b1 ,m E [H`2 ]b1 ,n [H`2 ]b2 ,n

(22)

b1 ,b2 =1

where (a) follows from R`1 `2 = 0 if `1 = `2 , and independence and the zero-mean assumption on the
TD channel for `1 6= `2 , which enforces `1 = `3 and `2 = `4 for (22). By inspection of (22), we observe
(m,n)

that 0-MSEω

is independent of m and n and thus, the MSE of the off-diagonal and diagonal entries

are equal. We simplify (22) for imperfect CSI and the BS-antenna correlation model. We first note that
 σ ∗
 σ2
b ]
bσ
E [H
I(b1 = b2 )
`1 b2 ,m [H`1 ]b1 ,m =
W
1
(α2 I(b1 = b2 ) + 2αβ + β 2 B)I(`1 < L)
+
BL
σ2
1
=
I(b1 = b2 ) +
(α2 I(b1 = b2 ) + δ)I(`1 < L),
W
BL
where the last step is obtained by noting that 2αβ + β 2 B = δ . Therefore, the inner sum is evaluated by
B
X

 σ ∗
  σ ∗

b ]
bσ
b
bσ
E [H
`1 b2 ,m [H`1 ]b1 ,m E [H`2 ]b1 ,n [H`2 ]b2 ,n

b1 ,b2 =1

Bσ 4
σ2
+
(α2 + δ)(I(`1 < L) + I(`2 < L))
W2
LW
1
+
(α2 (α2 + 2δ) + Bδ 2 )I(`1 < L)I(`2 < L).
(23)
BL2
√
P −1
2
Now, we simplify (22) using the results from (23) with the fact that α = 1 − δ and W
`1 =0 |R`1 `2 | = 2W
P −1
2
2
so that W
`1 ,`2 =0 |R`1 `2 | = 2W . Hence,
=

0-MSE(m,n)
ω

L−1
1 + δ 2 (B − 1) X
|R`1 `2 |2
= 2Bσ + 4σ +
BL2
4

2

`1 ,`2 =0
`1 6=`2

= 2Bσ 4 + 4σ 2 +

L−1
`
X X
2(1 + δ 2 (B − 1)) 
ejθ∆` 
1−
B
L2
`=0 ∆`=−`

2
= 2σ 2 (2 + Bσ 2 ) + (1 + δ 2 (B − 1))(1 − fL (θ)),
B

(a)

where (a) comes from the definition of Fejér kernel [27]. Note that we defined the shorthand variable
θ=

2π
W (pk

− ω).

24

A PPENDIX D
P ROOF OF C OROLLARY 2
Since fL (φ) is non-negative, it is obvious that maxω∈Ω {0-MSEω } ≤ εCSI + 2(1 + εcor )/B for all ω ∈ Ω.


b2π
The equality is satisfied if (p−ω)/W = b/L for some integer b > 0 so that fL 2π
= 0.
W (p−ω) = fL L
We note that this can only happen if dmax ≥ W/L, where dmax is the maximum distance between any
target subcarrier ω point and its nearest base point; this is due to the fact that

2π
W (p − ω)

≤

2π
W dmax

<

2π
L.

Assume dmax < W/L. Then, by Lemma 2, the maximum MSE of 0th order interpolation is given by:
max{0-MSEω }
ω∈Ω




2
2π
= εCSI + (1 + εcor ) 1 − min fL
(pk − ω)
B
W
ω∈Ω\P



2π
2
= εCSI + (1 + εcor ) 1 − fL
dmax
.
B
W
A PPENDIX E
P ROOF OF T HEOREM 2
The proof is similar to that of Theorem 1 in Appendix C. We start by defining the following auxiliary
function




2πpk+1
2πpk
∆` + (1 − λω ) exp j
∆`
Q``0 =λω exp j
W
W


2πω
− exp j
∆` ,
W
where we introduced the variable ∆` = ` − `0 . The result is obtained by substituting Q``0 in place
P −1
2
of R``0 at (22) in Appendix C. Note that W
`1 =0 |Q`1 `2 | = 2(1 − λω (1 − λω ))W which shows that
PW −1
2
2
`1 ,`2 =0 |Q`1 `2 | = 2(1 − λω (1 − λω ))W .
A PPENDIX F
P ROOF OF C OROLLARY 3
Without loss of generality, we will assume that the target subcarrier index ω is closer to pk+1 so that
λω ∈ [0, 0.5]. We will assume that pk < ω < pk+1 so that dk = pk+1 − pk > 0. Using the results from

Appendix C and Appendix E, the difference of 1-MSEω and 0-MSEω is given by:
1-MSEω − 0-MSEω = −εCSI λω (1 − λω ))
+

L−1
1 + εcor X
(|Q`1 `2 |2 − |R`1 `2 |2 ).
BL2
`1 ,`2 =0

25

Without loss of generality, we assume that `1 − `2 > 0 since |Q`1 `2 |2 − |R`1 `2 |2 is even and is 0 if `1 = `2 .
We simplify the term |Q`1 `2 |2 − |R`1 `2 |2 by denoting φ = θ(`1 − `2 ) > θ, where θ =

2π
W dk

and expand

the expression |Q`1 `2 |2 − |R`1 `2 |2 as follows:
(a)

|Q`1 `2 |2 − |R`1 `2 |2 = 2λω (1 − λω )(cos(φ) − 1)

− (cos((1 − λω )φ) − cos(λω φ))
(b)

(24)

= 2λω sin(φ/2)(sin((1 − 2λω )φ/2) − (1 − λω ) sin(φ/2)).

Here, (a) follows from the definition of Q`1 `2 and R`1 `2 and (b) is a results from simplifying the expression
cos((1 − λω )φ) − cos(λω φ)) = −2 sin((1 − 2λω )φ/2) sin(φ/2).

We first note that since εCSI = 2σ 2 (2 + Bσ 2 ) ≥ 0 and by (24), 1-MSEω = 0-MSEω if εCSI = 0, and

L = 1 or λω = 0. This behavior can be explained intuitively because when εCSI = 0 and L = 1, then the

channel is flat across all subcarriers, and hence 1-MSEω = 0-MSEω
Hence, we now show that 1-MSEω < 0-MSEω for L > 1 and λω 6= 0 by showing that |Q`1 `2 |2 −

|R`1 `2 |2 < 0. First note that `1 − `2 < L and if dk < W/(3L), then θ < φ < Lθ =

2π
W Ldk

<

2π
3 .

Hence,

sin(φ/2) > 0. Therefore, showing that (24) is negative is equivalent to:
g(λω ) =

sin((1 − 2λω )φ/2)
< sin(φ/2).
1 − λω

(25)

We now prove (25) by noting that g(λω ) = sin(φ/2) if λω = 0 and g 0 (λω ) < 0 for all λω ∈ (0, 0.5] so
g(λω ) is monotonically decreasing in (0, 0.5]. The proof is straightforward by:
g 0 (λω ) =

−φ cos((1 − 2λω )φ/2) sin((1 − 2λω )φ/2)
,
+
1 − λω
(1 − λω )2

(26)

and, hence, g 0 (λω ) < 0 in (26) can be expressed as:
tan((1 − 2λω )φ/2) < (1 − λω )φ.

(27)

To show (27), we introduce the shorthand notation γ = 1 − 2λω ∈ [0, 1). With the new notation γ , the
proof is straightforward by:

(a)

(b)

tan((1 − 2λω )φ/2) ≤ γ tan(φ/2) ≤ γφ < (γ + λω )φ
= (1 − λω )φ,

where (a) follows from the convexity of tan(x) in x ∈ [0, π/2) and (b) follows from tan(φ/2) < φ for

all φ ∈ (0, 2π/3]. Since g 0 (λω ) < 0 for all λω ∈ (0, 0.5], g(λω ) is monotonically decreasing and thus,

26

from (24), it follows that |Q`1 `2 |2 − |R`1 `2 |2 < 0 for L > 1.

We conclude by noting that a sharper upper bound on dk can be obtained by directly computing the

bounds for 1-MSEω − 0-MSEω , i.e.,
2λω (1 + εcor )
1-MSEω − 0-MSEω = −εCSI λω (1 − λω )) +
B


× fL (λω θ) − fL ((1 − λω )θ) − (1 − λω )(1 − fL (θ)) < 0,

for all λω ∈ (0, 0.5], but we leave an analysis of such refined bounds for future work.
R EFERENCES
[1] E. G. Larsson, F. Tufvesson, O. Edfors, and T. L. Marzetta, “Massive MIMO for next generation wireless systems,” IEEE
Commun. Mag., pp. 186–195, Feb. 2014.
[2] T. L. Marzetta, “Non-cooperative cellular wireless with unlimited numbers of base station antennas,” IEEE Trans. Wireless
Comm., vol. 9, no. 11, pp. 3590–3600, Nov. 2010.
[3] J. Hoydis, S. ten Brink, and M. Debbah, “Massive MIMO: How many antennas do we need?” in Proc. IEEE 49th Ann.
Allerton Conf. on Comm. Control, and Computing (Allerton), Monticello, IL, USA, Sept. 2011, pp. 545–550.
[4] A. Paulraj, R. Nabar, and D. Gore, Introduction to Space-Time Wireless Communications. Cambridge Univ. Press, 2003.
[5] E. Larsson, O. Edfors, F. Tufvesson, and T. Marzetta, “Massive MIMO for next generation wireless systems,” IEEE Commun.
Mag., vol. 52, no. 2, pp. 186–195, Feb. 2014.
[6] T. Narasimhan and A. Chockalingam, “Channel hardening-exploiting message passing (CHEMP) receiver in large-scale
MIMO systems,” IEEE J. Sel. Topics Signal Process., vol. 8, no. 5, pp. 847–860, Oct. 2014.
[7] C. Jeon, R. Ghods, A. Maleki, and C. Studer, “Optimality of large MIMO detection via approximate message passing,” in
Proc. IEEE Int. Symp. Inf. Theory (ISIT), Jun. 2015, pp. 1227–1231.
[8] C. Jeon, A. Maleki, and C. Studer, “On the performance of mismatched data detection in large MIMO systems,” in Proc.
IEEE Int. Symp. Inf. Theory (ISIT), Jul. 2016, pp. 180–184.
[9] M. Wu, B. Yin, G. Wang, C. Dick, J. Cavallaro, and C. Studer, “Large-scale MIMO detection for 3GPP LTE: Algorithm
and FPGA implementation,” IEEE J. Sel. Topics Signal Process., vol. 8, no. 5, pp. 916–929, Oct. 2014.
[10] B. Yin, M. Wu, J. R. Cavallaro, and C. Studer, “Conjugate gradient-based soft-output detection and precoding in massive
MIMO systems,” in Proc. IEEE Global Telecommun. Conf. (GLOBECOM), Dec. 2014, pp. 3696–3701.
[11] Z. Wu, C. Zhang, Y. Xue, S. Xu, and X. You, “Efficient architecture for soft-output massive MIMO detection with
Gauss-Seidel method,” in Proc. IEEE Int. Symp. Circuits and Syst. (ISCAS), May 2016, pp. 1886–1889.
[12] M. Wu, C. Dick, J. R. Cavallaro, and C. Studer, “FPGA design of a coordinate descent data detector for large-scale
MU-MIMO,” in Proc. IEEE Int. Symp. Circuits and Syst. (ISCAS), May 2016, pp. 1894–1897.
[13] H. Sampath and A. Paulraj, “Linear precoding for space-time coded systems with known fading correlations,” IEEE Commun.
Lett., vol. 6, no. 6, pp. 239–241, Jun. 2002.
[14] A. Scaglione, P. Stoica, S. Barbarossa, G. B. Giannakis, and H. Sampath, “Optimal designs for space-time linear precoders
and decoders,” IEEE Trans. Signal Process., vol. 50, no. 5, pp. 1051–1064, May 2002.
[15] H. Sampath, P. Stoica, and A. Paulraj, “Generalized linear precoder and decoder design for MIMO channels using the
weighted MMSE criterion,” IEEE Trans. Commun., vol. 49, no. 12, pp. 2198–2206, Dec. 2001.

27

[16] C. Studer, S. Fateh, and D. Seethaler, “ASIC implementation of soft-input soft-output MIMO detection using MMSE
parallel interference cancellation,” IEEE J. Solid-State Circuits, vol. 46, no. 7, pp. 1754–1765, Jul. 2011.
[17] IEEE Draft Standard; Part 11: Wireless LAN Medium Access Control (MAC) and Physical Layer (PHY) specifications;
Amendment 4: Enhancements for Higher Throughput, P802.11n/D3.0, Sep. 2007.
[18] 3GPP, “Evolved Universal Terrestrial Radio Access (E-UTRA); Physical channels and modulation,” 3rd Generation
Partnership Project (3GPP), TS 36.211, Jan. 2016. [Online]. Available: http://www.3gpp.org/ftp/Specs/html-info/36211.htm
[19] D. Cescato, M. Borgmann, H. Bölcskei, J. Hansen, and A. Burg, “Interpolation-based QR decomposition in MIMO-OFDM
systems,” in Proc. IEEE Int. Workshop Signal Process. Advances Wireless Commun. (SPAWC), June 2005, pp. 945–949.
[20] M. Borgmann and H. Bölcskei, “Interpolation-based efficient matrix inversion for MIMO-OFDM receivers,” in Proc.
Asilomar Conf. Signals, Syst., Comput., vol. 2, Nov. 2004, pp. 1941–1947 Vol.2.
[21] D. Cescato and H. Bölcskei, “Algorithms for interpolation-based QR decomposition in MIMO-OFDM systems,” IEEE
Trans. Signal Process., vol. 59, no. 4, pp. 1719–1733, Apr. 2011.
[22] ——, “QR decomposition of Laurent polynomial matrices sampled on the unit circle,” IEEE Trans. Inf. Theory, vol. 56,
no. 9, pp. 4754–4761, Sep. 2010.
[23] L. W. Chai, P. L. Chiu, and Y. H. Huang, “Reduced-complexity interpolation-based QR decomposition using partial layer
mapping,” in Proc. IEEE Int. Symp. Circuits and Syst. (ISCAS), May 2011, pp. 2381–2384.
[24] P. L. Chiu, L. Z. Huang, L. W. Chai, and Y. H. Huang, “Interpolation-based QR decomposition and channel estimation
processor for MIMO-OFDM system,” IEEE Trans. Circuits Syst. I, vol. 58, no. 5, pp. 1129–1141, May 2011.
[25] C. Desset, E. Björnson, S. Kashyap, E. G. Larsson, C. Mollén, L. Liu, S. Malkowsky, H. Prabhu, Y. Liu, J. Vieira, O. Edfors,
E. Karipidis, F. Dielacher, and D. V. Pop, “Distributed and centralized baseband processing algorithms, architectures, and
platforms,” MAMMOET, Tech. Rep. ICT-619086-D3.2, 2016.
[26] H. Bölcskei, D. Gesbert, and A. Paulraj, “On the capacity of OFDM-based spatial multiplexing systems,” IEEE Trans.
Commun., vol. 50, no. 2, pp. 225–234, Feb. 2002.
[27] P. L. Butzer and R. J. Nessel, Fourier analysis and approximation.

Academic Press New York, 1971.

[28] S. Jaeckel, L. Raschkowski, K. Börner, and L. Thiele, “QuaDRiGa: A 3-D multi-cell channel model with time evolution
for enabling virtual field trials,” IEEE Trans. Antennas Propag., vol. 62, no. 6, pp. 3242–3256, Jun. 2014.

