1

Message passing-based joint CFO and
channel estimation in millimeter wave systems
with one-bit ADCs
arXiv:1803.09012v1 [] 23 Mar 2018

Nitin Jonathan Myers, Student Member, IEEE, and
Robert W. Heath Jr., Fellow, IEEE.

Abstract
Channel estimation at millimeter wave (mmWave) is challenging when large antenna arrays are
used. Prior work has leveraged the sparse nature of mmWave channels via compressed sensing based
algorithms for channel estimation. Most of these algorithms, though, assume perfect synchronization and
are vulnerable to phase errors that arise due to carrier frequency offset (CFO) and phase noise. Recently
sparsity-aware, non-coherent beamforming algorithms that are robust to phase errors were proposed
for narrowband phased array systems with full resolution analog-to-digital converters (ADCs). Such
energy based algorithms, however, are not robust to heavy quantization at the receiver. In this paper,
we develop a joint CFO and wideband channel estimation algorithm that is scalable across different
mmWave architectures. Our method exploits the sparsity of mmWave MIMO channel in the angle-delay
domain, in addition to compressibility of the phase error vector. We formulate the joint estimation as
a sparse bilinear optimization problem and then use message passing for recovery. We also give an
efficient implementation of a generalized bilinear message passing algorithm for the joint estimation in
mmWave systems with one-bit ADCs. Simulation results show that our method is able to recover the
CFO and the channel compressively, even in the presence of phase noise.

Index Terms
Millimeter wave communication, wideband channel estimation, synchronization, compressed sensing, message passing, one-bit receivers
N. J. Myers (nitinjmyers@utexas.edu) and R. W. Heath Jr. (rheath@utexas.edu) are with the Wireless Networking and
Communications Group, The University of Texas at Austin, Austin, TX 78712 USA. This material is based upon work supported
in part by the National Science Foundation under Grant No. NSF-CCF-1527079, and by a gift from Huawei Technologies, Inc.

2

I. I NTRODUCTION
Millimeter wave communication introduces new challenges in the design of MIMO communication systems [1]. For instance, large antenna arrays at the transmitter (TX) and the receiver
(RX) are necessary to meet the link budget requirements [2]. As a result, the channel has a
higher dimension compared to what is typical in lower frequency MIMO systems, and must be
estimated more frequently thanks to the smaller coherence time [3]. Furthermore, cost and power
consumption are major issues at the larger bandwidths that accompany mmWave, primarily due
to high resolution ADCs [4]. Typical mmWave hardwares that limit power consumption at large
bandwidths introduce compression in the channel measurements. For example, the one-bit ADC
architecture [4] allows access to the output of every antenna at the expense of heavy quantization.
The compression of channel measurements and the use of large antenna arrays complicate signal
processing at mmWave.
Compressed sensing (CS) [5] [6] is an efficient technique to recover sparse high-dimensional
signals with few projections. As MIMO channel matrices at mmWave are sufficiently sparse when
expressed in an appropriate dictionary, applying tools from CS to mmWave channel estimation
can potentially reduce the training overhead. CS-based sparse channel estimation algorithms have
been proposed for various hardware architectures [7]–[9]. Recent developments in approximate
message passing [10] [11] have enabled channel estimation algorithms in low resolution receivers
[12]. Most CS-based channel estimation algorithms, however, assume perfect synchronization and
fail in practice because of the CFO and phase noise [13].
CFO and phase noise are hardware impairments that corrupt the phase of the channel measurements. The mismatch between the carrier frequencies of the local oscillators at the TX and
the RX results in CFO. Phase noise in the system arises due to short-term random fluctuations
in the frequency of the oscillators. Both these non-idealities are larger at mmWave due to the
high carrier frequency and ignoring them can result in significant channel estimation error [14].
Correcting for the CFO and then performing channel estimation seems like a possible solution.
The disadvantage, however, is that prior to beamforming or channel estimation, mmWave systems
operate at very low SNR, which can result in significant error in the CFO estimate. Prior work
has considered joint CFO and channel estimation [15] [16] in lower frequency systems. These
joint estimation algorithms, however, cannot be applied to typical mmWave systems due to
differences in the hardware architectures. Furthermore, they are not designed to incorporate the

3

sparse nature of mmWave channels. Therefore, there is a need to design either phase error robust
channel estimation algorithms, or joint CFO and channel estimation algorithms that can exploit
the sparsity of mmWave channel.
Recent work on phase error robust channel estimation is limited to narrowband systems and has
focussed on specific mmWave hardware. In [13], [17] and [18], phase error robust compressive
beamforming algorithms were proposed for the analog beamforming architecture. In [13], phase
tracking followed by phase error compensated compressive beamforming was proposed. The
compensation, however, was done prior beamforming and therefore suffers from low SNR. The
non-coherent algorithms in [17] and [18] are not robust to heavy quantization at the receiver.
These algorithms cannot be used in one-bit receivers because the energy information of the
channel measurements is completely lost due to one-bit quantization.
Recent joint CFO and sparse narrowband channel estimation algorithms in [14] and [19]
require high computational complexity when extended to wideband systems. In [14], we proposed
a joint CFO and narrowband channel estimation algorithm using third-order tensors [20]. We
also developed a sparsity-aware joint estimation algorithm for the one-bit ADC architecture in
[19]. The main idea underlying our approach in [19] was to use the lifting technique [21] along
with message passing for the joint estimation with one-bit channel measurements. Extending
the narrowband solutions in [14] or [19] to typical wideband mmWave systems would require
convex optimization over millions of variables, which may be prohibitive in a practical setting.
In this paper, we propose a sparse bilinear formulation of the joint CFO and wideband channel
estimation problem, and solve it using message passing. We assume that all the RF chains at
the TX or the RX are driven by the same reference oscillator. Hence, there is a unique CFO
and a phase noise process in our MIMO system model, that corrupt the channel measurements.
We also assume that there is perfect frame timing synchronization between the TX and the RX.
We summarize the main contributions of our work as follows.
•

We formulate the joint CFO and wideband channel estimation problem as a noisy quantized
sparse bilinear optimization problem. Our framework leverages the sparse nature of the wideband channel in the angle and delay domains, and also exploits the compressibility of the
phase error vector in the frequency domain.

•

To solve the non-convex problem at hand, we use the vector variance version of the Parametric
Bilinear Generalized Approximate Message Passing (PBiGAMP) algorithm [22] and optimize

4

it for fast joint estimation. The parameters of the sparse priors corresponding to the wideband
channel and the phase error vector are learned using an Expectation Maximization (EM)
algorithm.
•

We provide insights into the design of training matrices for joint CFO and channel estimation
using PBiGAMP. Specifically, we show that shifted Zadoff-Chu training proposed in [12] to
accelerate message passing cannot be used for joint estimation as it results in a continuum
of optimal solutions for the bilinear optimization problem. We explain the “CFO propagation
effect” to highlight the trade-off between fast message passing and identifiability in the sparse
bilinear problem.

•

We evaluate the performance of our joint estimation algorithm assuming a digital receiver
architecture with one-bit ADCs and compare it with the hypothetical full resolution case.
Simulation results show that the proposed approach is able to recover both the channel and
the CFO compressively with IID Gaussian and IID QPSK training matrices, even in the
presence of phase noise uncertainity.

Our algorithm is advantageous over the existing sparsity-aware methods for joint estimation or
phase error robust channel estimation in terms of the capability to efficiently handle frequency
selective channels and scalability to other mmWave architectures.
Notation: A is a matrix, a is a column vector and a, A denote scalars. Using this notation
AT , A and A∗ represent the transpose, conjugate and conjugate transpose of A. The matrices
|A| and |A|2 contain the element-wise magnitude and squared magnitude of the entries of A.
We use A(i) and A(j) to denote the ith row and j th column of A. We use diag (a) to denote
a diagonal matrix with entries of a on its diagonal. The scalar am denotes the mth element
of a. The symbol ⊗ is used to denote the kronecker product. vec (A) is a vector obtained by
stacking all the columns of A and vecm (A) denotes the mth element of vec (A). We define

Ai,j = veci A(j) . We use IN to denote the set {1, 2, 3, ..N }. The matrix UN ∈ CN ×N denotes
the unitary discrete Fourier transform matrix. N (m, R) is the probability density function of
complex Gaussian random vector with mean m and covariance R. We define e`,N ∈ RN ×1 as
the N dimensional canonical basis vector with its `th coordinate as 1.

5

II. S YSTEM AND C HANNEL M ODELS
In this section, we describe the underlying hardware architecture, CFO and phase noise model,
and the wideband mmWave channel model used for our simulations. In particular, we focus on the
digital receiver architecture with one-bit ADCs, to highlight the differences with the existing noncoherent algorithms. Nevertheless, our algorithm can be extended to other mmWave architectures
as the underlying joint estimation problem is bilinear in nature.
A. System Model

RF

Re()

Im()

Transmitter

Ntx
RF
Block

H

ADC
Digital
Baseband
Processing

Nrx

RF

Re()
Im()

f1

q bit

q bit

ADC

f2

Fig. 1. A MIMO system with local oscillators operating at f1 and f2 , and q-bit ADCs at the receiver. Each antenna is associated
with an RF chain and a pair of q-bit ADCs. In this work, we consider the extreme cases of q = 1 and q = ∞.

We consider a MIMO system with uniform linear array of Ntx antennas at the TX and Nrx
antennas at the RX, as shown in Fig. 1. We use linear arrays for a concise representation of the
simplifications involved in PBiGAMP; our framework can be extended to other array geometries
using appropriate array response vectors in the formulation. We do not impose constraints on the
number of RF chains or the resolution of the digital-to-analog converters (DACs) at the TX. The
resolution of the Nrx ADCs at the RX, however, is assumed to be limited. The baseband signal
at the TX is upconverted to the mmWave band, using a local oscillator at a carrier frequency
f1 . The transmitted RF signal propagates through the wireless channel and is downconverted at
the RX using a carrier frequency f2 , that slightly differs from f1 . Although the MIMO system
can have multiple RF chains, we assume that all the RF chains at a given end are driven by
the same reference oscillator. Even if the RF chains at a given end were driven by different

6

oscillators, achieving carrier synchronization locally is feasible [23]. After downconversion at
the RX, the output at each antenna is sampled using a pair of q bit ADCs, one each for the
in-phase and the quadrature phase components. We use Qq (.) to represent the q bit quantization
function corresponding to the ADCs. In this work, we consider the extreme case of q = 1 and
provide a performance comparison relative to q = ∞. The quantization functions for the two
cases are Q1 (x) = sign (real (x)) + j sign (imag (x)) and Q∞ (x) = x. Note that the functions
sign (·) , real (·) and imag (·) are applied element-wise on the vector.
The impact of CFO on channel estimation algorithms is more significant in one-bit receivers
than the full resolution ones. The mismatch in the carrier frequencies, i.e., |f2 − f1 | is typically
in the order of several parts per millions (ppms) of f1 or f2 . Due to the high carrier frequencies
at mmWave, even such small differences can significantly perturb the channel estimate when
ignored [13]. For a symbol duration of T seconds, we define the digital domain CFO as  =
2π (f1 − f2 ) T . CFO results in unknown phase errors in the received samples that linearly increase
with time. Hence, the impact of CFO on standard channel estimation algorithms is determined
by the length of training. As one-bit receivers relatively need a longer training for channel
estimation when compared to the full resolution ones, channel estimation algorithms that ignore
phase errors are more vulnerable to the CFO in one-bit systems than the full resolution ones.
Phase noise in wireless systems arises due to jitter in the frequency of the oscillators. For
2
2
a phase noise variance of βtx
at the TX and βrx
at the RX, the phase noise variance in the
2
2
received samples can be approximated as β 2 = βrx
+ βtx
. The approximation is valid when

the 3 dB bandwidth of the phase noise power spectral density is significantly smaller than the
channel coherence bandwidth [24]. Let φk denote the phase error introduced in the k th received
sample, due to phase noise at the TX and the RX. As is common in prior work, we model
the phase errors using a Wiener process [25] in which the increments, i.e., φk − φk−1 , are IID
Gaussian random variables with zero mean and variance β 2 . As β 2 is proportional to f12 [13],
phase noise is higher at mmWave carrier frequencies for a given quality of oscillator.
Now, we describe the received signal model in the digital receiver architecture. Let t [n] ∈
CNtx ×1 be the nth transmit symbol satisfying the power constraint E [t∗ [n] t [n]] = P . The
discrete time baseband representation of the MIMO channel is assumed to be limited to L
taps. Let H [`] ∈ CNrx ×Ntx be the `th tap of the equivalent baseband channel, where ` ∈
{0, 1, 2, ..., L − 1}. We assume that perfect frame timing synchronization can be achieved using
the control channel. Our assumption can be justified in situations where the mmWave system

7

co-exists with a lower frequency system [26] that can perform timing synchronization. Frequency
synchronization, however, may not be achieved due to different offsets and phase noise processes
for each of these systems. With the timing synchronization assumption, the sampled baseband
vector in the nth symbol duration can be given by
Y(n) = Qq

ej(n+φn )

L−1
X

!
H [`] t [n − `] + V(n) ,

(1)

`=0
2

where V(n) ∼ N (0, σ INrx ) is additive white Gaussian noise. In this work, we develop an
algorithm to estimate  and {H [`]}L−1
`=0 from the series of observations Y(n) . Our joint estimation
algorithm can be extended to any q-bit ADC architecture by defining appropriate output likelihood
functions in message passing.
B. Channel Model
We consider a clustered channel model for the frequency selective mmWave MIMO channel.
The channel consists of Ncs clusters with Mn rays in the nth cluster. Let γn,m , τn,m , θr,n,m and
θt,n,m denote the complex gain, delay, angle-of-arrival (AoA) and angle-of-departure (AoD) of
the mth ray in the nth cluster. We assume that the transmitted signal is bandlimited to 1/T Hz.
With ωr,n,m = π sin θr,n,m , ωt,n,m = π sin θt,n,m and the Vandermonde vector

T
aN (∆) = 1 , ej∆ , ej2∆ , · · · , ej(N −1)∆ ,

(2)

the `th tap of the wideband MIMO channel for a half wavelength spaced uniform linear array
is given by
Ncs X
Mn
X

τn,m 
.
H [`] =
γn,m aNrx (ωr,n,m ) aNtx (ωt,n,m ) sinc ` −
T
n=1 m=1
∗



(3)

The wideband channel can be represented using Nrx Ntx L complex entries, and the matrix in (3)
is large in typical mmWave systems. The channel impulse response in (3) is represented using a
linear combination of bandlimited sinc(·) functions. Notice that each of these sinc(·) functions
is delayed by the normalized delay spread, i.e, τn,m /T and evaluated at periodic time instants
to obtain the discrete time representation in (3). Other filtering functions could also be used to
incorporate the effect of pulse shaping at the TX or filtering at the RX [27].
The mmWave MIMO channel is aproximately sparse in an appropriate dictionary due to the
propagation characteristics of the environment at mmWave frequencies. Compared to the lower
frequency channels, mmWave channels comprise of fewer clusters [4]. Each of the channel taps

8

H [`], is approximately sparse in the spatial Fourier basis at mmWave [12]. Furthermore, the
channel is approximately sparse along the time dimension as the delays of the propagation rays
are heavily clustered within the delay spread. As the delays τn,m may not necessarily be an
integer multiple of T , there is a leakage effect along the time dimension. Let C [`] ∈ CNrx ×Ntx
be the 2-D Fourier transform of H [`], such that
H [`] = UNrx C [`] U∗Ntx ,

∀` ∈ {0, 1, 2, ..., L − 1} .

(4)

The approximate sparsity of the mmWave MIMO channel along the angle and delay domains
[28] is directly reflected in the matrices {C [`]}L−1
`=0 . A higher resolution dictionary can be used
for the angle and delay dimensions to increase sparsity of the mmWave channel at the expense
of higher dimensionality and higher frame coherence [29]. Our GAMP based approach will be
robust to leakage effects that arise due to approximate sparsity.
III. D EMYSTIFYING JOINT CFO

AND CHANNEL ESTIMATION AT MM WAVE

In this section, we propose a sparse bilinear formulation for the joint estimation problem. We
also identify existing techniques to solve the problem and describe their limitations in terms of
scalability to other mmWave architectures and computational complexity.
A. Bilinear formulation
We derive a compact form for the received signal model in (1) for a SC-FDE system [3]. Let
T ∈ CNtx ×Np be a training block of length Np such that T(k) = t [k]. The TX transmits a training


sequence with a cyclic prefix of length L − 1, i.e., T(Np −L+2) , T(Np −L+3) , ..., T(Np ) , T . The
cyclic prefix padded transmission gets convolved with the frequency selective MIMO channel
before sampling at the receiver. As usual, the first L − 1 samples of the received block that
experience interference from the previous transmit block are discarded [3]. Let Y ∈ CNrx ×Np
be the received block obtained after discarding the first L − 1 received vectors. We define
an ` circulant delay matrix J` ∈ CNp ×Np , such that its first column is the canonical basis
vector e1+`,Np . We define a vector d (, β) ∈ CNp ×1 such that its nth entry has the phase error
corresponding to the received vector Y(n) , i.e., dn (, β) = ej(n+φn ) . Note that β is the standard
deviation of the incremental phase errors in the Wiener phase noise process. The received samples
corresponding to the transmit block defined by T can be expressed using (1) as
!
L−1
X
Y = Qq
H [`] TJ` diag (d (, β)) + V .
`=0

(5)

9

The phase errors in (5) are invariant along any column of the unquantized received block as
there is a unique CFO and a phase noise process in the system. The representation in (5) can
be further simplified to capture the structure in the phase errors and the channel.
We exploit the structure in the joint estimation problem using sparsity of the channel and the
phase error vector in appropriate dictionaries. A compact representation of (5) can be obtained
by following the same steps in [12], except for the diagonal matrix containing the phase errors.
We use b to denote the DFT of the phase error vector d (, β). With Z used to denote the
noiseless unquantized version of Y in (5) such that Y = Qq (Z + V), we have
Z=

L−1
X

H [`] TJ` diag (d (, β))

`=0

=

L−1
X

UNrx C [`] U∗Ntx TJ`

diag





U∗Np b

`=0



U∗Ntx TJ0


 U∗ TJ
1

Ntx
= UNrx [C [0] C [1] C [2] . . . C [L − 1]] 
.
..
{z
}
|

∆
=C
U∗Ntx TJL−1
|
{z
∆







 diag U∗Np b .



(6)

}

=F
Nrx ×Ntx L

The matrix C ∈ C

is just a concatenation of the angle domain representation in (4)

corresponding to each tap of the MIMO channel and is approximately sparse. Furthermore, the
vector b, i.e., the DFT of d (, β) can be considered to be approximately sparse. The sparse
representation of the phase error vector is valid in practice as the spread of the oscillator’s
spectrum about the center frequency is relatively small compared to the bandwidth of the signal.
Now, we derive a sparse bilinear formulation for the joint estimation problem. Using (6), the
quantized received block Y in (5) can be expressed as




Y = Qq UNrx CFdiag U∗Np b + V .

(7)

We define the Np Nrx ×1 vectors y = vec (Y), z = vec (Z) and v = vec (V), and the Nrx Ntx L×1
vector c = vec (C). In vector notation, (6) can be written as


z = diag U∗Np b ⊗ aNrx (0) vec (UNrx CF) .

(8)

10

Notice that aNrx (0) is just the all ones vector in Nrx dimension. Using the property vec (PQR) =

RT ⊗ P vec (Q), the received vector y is expressed as





(9)
y = Qq diag U∗Np b ⊗ aNrx (0) FT ⊗ UNrx c + v .
We define the matrices G = U∗Np ⊗ aNrx (0) and A = (F ⊗ UNrx )T to rewrite (9) as
y = Qq (diag (Gb) Ac + v) .

(10)

Estimating the phase errors and the channel is equivalent to estimating b and c from y in (10).
The joint estimation problem in (10) can be observed to be a noisy quantized bilinear problem
in b and c, subject to the sparsity of b and c.
B. Limitations of existing techniques
1) CFO robust methods: Existing sparse channel estimation methods that are robust to CFO
discard the phase of the channel measurements. A hashing technique based non-coherent beam
alignment algorithm was proposed in [17] for analog beamforming systems. This method, however, assumes fine control over the phase shifters, which is not necessarily the case with mmWave
systems. The received signal strength (RSS) based method in [18] accounts for the limited phase
control and uses pseudo-random phase shifts for compressive beam-training. The solutions in
[17] and [18] assume a narrowband mmWave system and perform beam-alignment with just the
magnitude of the channel measurements. These methods, however, cannot be used in one-bit
receivers as energy detection with undithered one-bit ADCs is not feasible unless additional
circuit components are used. For example, Q1 (r) and Q1 (α r) are the same for any α > 0.
As the only information provided by undithered one-bit ADCs is phase quantized to 4 levels,
discarding it due to phase errors leaves no information.
2) Joint estimation using lifting: Lifting [21] [30] is a convex relaxation technique that
transforms a bilinear problem to a higher dimensional one and then recovers the original vectors
by solving the higher dimensional problem. We describe the lifting technique applied to the joint
estimation problem in (10). With M denoting the number of entries in y or z, i.e., M = Nrx Np ,
the mth entry of z can be given as
zm = vecm (diag (Gb) Ac)
= G(m) bA(m) c

(11)
(12)

11

T
= G(m) bcT A(m)
(13)


∀m ∈ IM .
(14)
= A(m) ⊗ G(m) vec bcT ,

We define a lifted variable x = vec bcT and a measurement matrix Φ ∈ CM ×Np Nrx Ntx L , such
that Φ(m) = A(m) ⊗ G(m) . Hence, the quantized measurements in (10) can be expressed as
y = Qq (Φx + n) .

(15)

The lifted vector x in (15) is sparse as it is just an outer product [20] of the sparse vectors b
and c. Several CS-based algorithms [6] [31] can be used to recover x from the possibly underdetermined noisy quantized system in (15). Using the SVD of the higher dimensional matrix
estimate, the vectors in the bilinear problem can be estimated upto a scale factor.
Lifting followed by the SVD was applied to joint CFO and narrowband channel estimation
for one-bit receivers in our previous work [19]. The main issue in extending our method in [19]
to wideband systems arises due to the large dimensionality of the lifted problem. For instance,
the dimension of x to perform joint CFO and channel estimation in wideband systems would
be Nrx Ntx Np L. Using lifting necessarily implies solving for millions of variables for typical
wideband mmWave systems, due to the large number of antennas and the need for additional
pilots to compensate for the heavy quantization in low resolution systems.
The limitations of the existing phase error robust and joint estimation solutions in terms
of architectural scalability and computational complexity, motivate the need to develop new
low complexity joint estimation algorithms that can be applied to wideband systems and low
resolution receivers.
IV. M ESSAGE PASSING BASED JOINT CFO

AND CHANNEL ESTIMATION

In this section, we give a brief introduction to PBiGAMP [22] and discuss its application
to the joint estimation in (10). We exploit the inherent structure in our problem to derive a
low complexity and memory efficient implementation of PBiGAMP for the joint estimation.
Furthermore, we explain the CFO propagation effect induced by special training matrices that
prevents further reduction in the computational complexity.
A. Introduction to PBiGAMP
The joint estimation problem in (10) can be solved using PBiGAMP [22] by considering b, c,
z and y of (8) and (10) as realizations of random vectors, say b, c, z and y. Let bi , ck , zm and

12

ym be the elements of these random vectors. In general, deriving the closed form Minimum
Mean-Squared Error (MMSE) estimates [32] of b and c is difficult as it requires marginalizing
the joint PDF of b and c conditioned on y = y. Using ideas from message passing, PBiGAMP
can obtain the MMSE estimates of both the vectors in the bilinear problem.
We explain message passing using the factor graph [33] in Fig. 2, that shows the dependency
between y, the random variables (b, c) and their prior distributions. The circular nodes in the
factor graph are called as variable nodes as they represent the random variables. The rectangular
nodes in the factor graph are called as factor nodes and they contain the prior distribution
of a random variable or the likelihood function associated with an observation. The messages
in message passing are essentially probability distributions, also called as beliefs. The idea
underlying message passing is to perform belief flows iteratively between the factors and the
variables until all the variable nodes reach a consensus on their marginal probability distributions.
As the factor graph for joint estimation is strongly connected, standard message passing can be
computationally intractable.
Using ideas from Approximate Message Passing [10], PBiGAMP simplifies the messages by
assuming a large number of variable nodes. Simulation results in [22] that show that PBiGAMP
outperforms lifting techniques in Section III-B2, for IID Gaussian measurement matrices, motivate applying it to our problem. Furthermore, PBiGAMP performs optimization over the same
number of variables in the problem, unlike lifting [21] that solves the problem in a higher
dimensional space. For the joint estimation problem in typical wideband mmWave systems,
PBiGAMP is memory efficient over lifting by several orders of magnitude.
B. PBiGAMP for joint estimation
In this section, we explicitly state PBiGAMP [22] for joint estimation in (10) and describe
the information contained in the factor nodes of Fig. 2. To be consistent with the notation used
in [22], we rewrite the random variable dependency corresponding to (11) in the tensor notation
as
zm =

Nb X
Nc
X

(i,k)
zm
bi ck ,

(16)

i=1 k=1
(i,k)

where Nb = Np , Nc = Nrx Ntx L, and zm

is an element of a third order tensor given by

(i,k)
zm
= Gm,i Am,k .

(17)

13

Messages

Likelihood
functions

Messages

Channel
variables

Priors

Priors

Synchronization
variables

Fig. 2. The factor graph of bilinear message passing for joint CFO and channel estimation. The rectangular nodes, called as
factors, contain the likelihood functions corresponding to the received samples or the sparse priors. Messages are sent between
the factor nodes and the variable nodes until the marginal probability distributions of the variables converge.

The output likelihood function in Fig. 2, denoted by pym |zm (ym | z) is given by
 √
 √


2sign(Re{ym })Re{z}
2sign(Im{ym })Im{z}
F
F
q=1
σ
σ
pym |zm (ym | z) =
,
ky −zk2

 1 2 e− mσ2
q=∞
πσ

(18)

where F (·) is the cumulative distribution function of the standard normal distribution. The
sparsity of the vectors b and c is incorporated by assuming parametrized Bernoulli-Gaussian
distributions for their priors pb (b) and pc (c). For simplicity, it is assumed that each entry of
b is independent of the other and identically distributed as pb . Similarly, the entries of c are
assumed to be IID, with pc as the distribution. Furthermore, the vectors b and c are assumed to
be independent of each other. Let λb and λc denote the sparsity fraction of b and c. Let σb2 and
σc2 be the variances of the coefficients corresponding to the non-zero support of the vectors b
and c. With δ (x) used to represent the Dirac-delta function, the Bernoulli-Gaussian distributions
pb and pc can be given as
pb (x) = λb δ (x) + (1 − λb ) N (0, σb2 ),

(19)

pc (x) = λc δ (x) + (1 − λc ) N (0, σc2 ).

(20)

14

Definitions:

pzm |pm z | pb; ν p ,
pck |rk(c | rb; ν r ) ,
(b | qb; ν q ) ,

pbi |qi

pym |zm(ym | z) N (z;b
p,ν p )
pym |zm(ym | z 0 ) N (z 0 ;b
p,ν p )
r ,ν r )
R pc(c) N (c;b
0
0 r ,ν r )
c0 pc(c ) N (c ;b
q ,ν q )
R pb(b) N (b;b
0
0 q ,ν q )
b0 pb(b ) N (b ;b

(D1)

R

z0

(D2)
(D3)

Initializations:
∀m : sbm (0) = 0
∀i, k : choose b
bi (1), ν b (1), b
ck (1), ν c (1)
i

(I1)
(I2)

k

for t = 1, . . . Tmax
(i,∗)

∀m, i : zbm

(i,k)
b
ck (t)
k=1 zm
P Nb b
(i,k)
= i=1 bi (t)zm
P Nb b
(i,∗)
bi (t)b
zm (t)
= i=1

(t) =

∀m, k :

(∗,k)
zbm (t)

∀m :

(∗,∗)
zbm (t)

∀m : ν pm (t) =

P Nc

(R1)
(R2)
or

(∗,k)
ck (t)b
zm (t)
k=1 b
P Nc c
(∗,k)
zm (t)|2
k=1 νk (t)|b

P Nc

(i,∗)

P Nb

b
zm (t)|2 +
i=1 νi (t)|b
P Nb b P Nc c
(i,k)
p
p
∀m : νm (t) = ν m (t) + i=1 νi (t) k=1 νk (t)|zm |2

∀m :

(∗,∗)
pbm (t) = zbm (t)

∀m :

z (t) = var{z
νm
m

−

sbm (t−1)ν pm (t)

(R4)
(R5)
(R6)

p
| pm = pbm (t); νm
(t)}

(R7)

p
∀m : zbm (t) = E{zm | pm = pbm (t); νm
(t)}

(R8)

s (t) = (1 − ν z (t)/ν p (t))/ν p (t)
∀m : νm
m
m
m

(R9)

p
pbm (t))/νm
(t)

∀m : sbm (t) = (b
zm (t) −
P
−1
(∗,k)
M
s
r
zm (t)|2
∀k : νk (t) =
m=1 νm (t)|b
P
(∗,k)
∀k : rbk (t) = b
ck (t) + νkr (t) M
bm (t)b
zm (t)∗
m=1 s
P
P Nb b
(i,k) 2
M
s (t)
− νkr (t)b
ck (t) m=1 νm
|
i=1 νi (t)|zm
P
−1
(i,∗)
q
M
s
∀i : νi (t) =
zm (t)|2
m=1 νm (t)|b
P
(i,∗)
q
M
∀i : qbi (t) = b
bi (t) + νi (t) m=1 sbm (t)b
zm (t)∗
P
P
(i,k)
− ν q (t)b
bi (t) M ν s (t) Nc ν c (t)|zm |2
m=1

i

m

k=1

k

∀k : νkc (t+1) = var{ck | rk = rbk (t); νkr (t)}

if

(R3)

∀k : b
ck (t+1) = E{ck | rk = rbk (t); νkr (t)}
∀i : νib (t+1) = var{bi | qi = qbi (t); νiq (t)}
∀i : b
bi (t+1) = E{bi | qi = qbi (t); νiq (t)}
PM
P
(∗,∗)
(∗,∗)
(∗,∗)
zm (t) − zbm (t−1)|2 ≤ τstop M
zm (t)|2 ,
m=1 |b
m=1 |b

(R10)
(R11)

(R12)
(R13)

(R14)
(R15)
(R16)
(R17)
(R18)

stop

(R19)

end
b=b
b (t) , b
Output: b
c=b
c (t) .

TABLE I
T HE PB I GAMP A LGORITHM FROM [22]

The parameters governing pb and pc , however, are not known apriori and can be learned by
embedding PBiGAMP within the Expectation Maximization (EM) algorithm [11]. For a given
set of likelihood functions and prior distributions, the vector variance PBiGAMP algorithm [22]
to obtain the MMSE estimates of b and c in (10) is summarized in Table I.
The channel and the phase error vector can be derived using appropriate transformations
over the PBiGAMP estimates. The vector b
c obtained from PBiGAMP is just an estimate of

15

the vectorized version of C, the angle-delay domain representation of the wideband channel in
b ∈ CNrx ×Ntx L denote the angle-delay domain estimate of the wideband channel, such
(7). Let C
 
b =b
that vec C
c. It can be seen from (6) that C is just a concatenation of the angle domain
representation of all the L taps of the MIMO channel. Therefore, the `th tap of the antenna
domain MIMO channel can be derived as
h
i
b
b
b
b
H [`] = UNrx C(`Ntx +1) , C(`Ntx +2) , . . ., C(`Ntx +Ntx ) U∗ Ntx .

(21)

b derived from PBiGAMP is an estimate of the DFT of the phase error vector
The vector b
b is just a single tone frequency estimation problem. As the
d (, β). Estimating the CFO from b
phase noise is modelled as a Wiener process, we apply the Extended Kalman Filter (EKF) [34]
b to estimate .
over the time domain samples, i.e., U∗ b,
Np

In this work, we exploit the compressibility of the phase error vector in the DFT basis for
joint estimation. It is possible, however, to incorporate the statistics of the phase noise by
replacing the nodes corresponding to b in Fig. 2 with the phase error variables and adding
factors corresponding to the phase noise process. In such case, the message passing algorithm
must handle the non-linear dependence of z on the phase errors.
C. PBiGAMP: From theory to practice
As seen from Table I, the generic implementation of PBiGAMP is memory and computationally intensive, as it involves repeated operations over a third order tensor in (17). In this
section, we provide insights into the key equations of Table I and describe our low complexity
implementation for joint estimation.
The equations in PBiGAMP are essentially determined by the belief flows from the factor
nodes to the variable nodes, and vice versa. Without loss of generality, we describe the message
flow between the variable node b1 and the factor node y1 in Fig. 2. In standard message passing,
the message sent from a factor node to b1 essentially represents the PDF of b1 presumed by
that factor node. In a fully connected factor graph, the node b1 receives messages from all the
M factor nodes ({yi }M
i=1 ) in addition to its prior distribution pb . As b1 receives several beliefs
from different factors, the belief sent by b1 to y1 is just a normalized product of all the beliefs
received by b1 except the one from y1 . For generic priors and likelihood functions in the factor
graph, performing standard message passing can be difficult as the belief flows are flows of
PDFs that are functions.

16

PBiGAMP simplifies standard message passing using the central limit theorem (CLT) and the
Taylor series approximation. Notice that y1 receives beliefs from all its neighbouring variable
nodes and contains the likelihood function in itself. The belief sent from y1 to b1 is computed
by multiplying the likelihood function, with all the incoming beliefs to y1 except the one from
b1 , and then integrating over all the random variables except b1 . The multiplication followed
by integration essentially yields the PDF of b1 presumed by the factor node y1 . The integration
is often multidimensional and can be difficult to compute. If y1 depends on a large number
of independent variable nodes through a linear function, then the linear combination can be
approximated as a Gaussian random variable using the CLT [10]. In such case, the variable nodes
can send just the mean and variances of the PDFs to the factors and this information is sufficient
to compute the mean and variance of the Gaussian random variable, as seen in (R3) and (R5)
of Table I. The message sent from y1 to b1 can now be computed by multiplying the likelihood
at y1 with the compound Gaussian PDF, and marginalizing the product with respect to b1 . In
general, the likelihood function can be non-linear in nature and can yield a complicated PDF
of b1 . Using a second order Taylor series expansion for the log PDF, PBiGAMP [10] simplifies
the belief sent from y1 to b1 to a Gaussian whose mean and variance can be computed from
(R7)-(R9) of Table I. For q-bit ADCs, the closed form expressions for the conditional mean and
variance in (R7) and (R8) can be found in [35, Appendix A]. Unlike standard message passing,
PBiGAMP is computationally tractable as the messages contain just the mean and variances of
the PDFs.
After several approximate message flows between the factor nodes and the variable nodes,
PBiGAMP is expected to converge. The MMSE of b1 is computed as the expectation of the
effective marginal, i.e., the normalized product of all the M Gaussian PDFs received by b1 from
the factor nodes and the prior distribution on b1 . Notice that the normalization has to be done to
ensure that the PDF integrates to 1. The expectation step for the MMSE of b1 is given in (R18)
using (R11) and (R12) as the intermediate steps. Similarly, the expectation for the MMSE of c1
is given in (R16) using (R13) and (R14) as the intermediate steps. Thus, PBiGAMP provides
estimates of the vectorized sparse channel and the DFT of the phase error vector.
The generic implementation of PBiGAMP is computationally expensive primarily due to (R1),
(R2), (R5), (R12) and (R14) in Table I. It can be verified that each of these operations have a

2
complexity of O Nrx
Np2 Ntx L for a single PBiGAMP iteration to perform joint estimation. For
(i,∗)

instance, (R1) requires computing the scalar zbm (t) for every m ∈ IM and i ∈ INb , Therefore,

17

(R1) demands O(M Nb Nc ) computations, as each scalar computation requires Nc multiplications.
As M = Nrx Np , Nb = Np and Nc = Nrx Ntx L for the joint estimation, (R1) has a complexity of

2
O Nrx
Np2 Ntx L for every PBiGAMP iteration. By exploiting the structure in the joint estimation
problem, the complexity of PBiGAMP can be significantly reduced.
We describe our fast implementation of PBiGAMP for the joint estimation in the following
b (t) ∈ CNp ×1 , C
b (t) ∈
sub-sections. The tth iteration variables of PBiGAMP are defined as b
b (t) ∈ CNrx ×Np , such that bbi (t) , b
CNrx ×Ntx L and Z
ck (t) and zbm (t) are the ith , k th and mth entries
b (t) , vec(C(t))
b
b
of b
and vec(Z(t)).
1) Operations (R1)-(R3): For an io ∈ INp , we have
(io ,∗)
zbm
(t)

=

Nc
X

(io ,k)
zm
b
ck (t)

(22)

k=1

=

Nb X
Nc
X

i ,Np

(i,k) o
zm
ei

b
ck (t)

(23)

i=1 k=1



∗
io ,Np
b
= vecm UNrx C(t)Fdiag(U
e
)
Np


∗
b
= vecm UNrx C(t)Fdiag(U
)
,
Np (i )
o

(24)
(25)

where the compact form in (25) is obtained by going back from the tensor formulation to the
original bilinear model in (6). Likewise, for a ko ∈ INc , we have
(∗,ko )
zbm
(t) =

Nb X
Nc
X

(i,k)b
zm
bi (t)ekko ,Nc

(26)

i=1 k=1



b
= vecm UNrx ero ,Nrx (eco ,Ntx L )T Fdiag(U∗Np b(t))


b
= vecm UNrx (ro ) F(co ) diag(U∗Np b(t))
,

(27)
(28)

where (ro , co ) correspond to the row and column of the Nrx × Ntx L matrix version of eko ,Nc ,
with
ko = (co − 1)Nrx + ro .

(29)



(∗,∗)
∗ b
b
zbm
(t) = vecm UNrx C(t)Fdiag(U
b(t))
.
Np

(30)

Similarly, (R3) can be computed as

b
To evaluate (30), the product C(t)F
can be found using Nrx Ntx LNp computations. Using the
b
FFT over the resultant product, UNrx C(t)F
can be computed with an additional complexity

18

of O(Np Nrx logNrx ). Finally, the complexity to multiply the resultant Nrx × Np matrix with
b
the IFFT of b(t)
is O(Np logNp ) + Nrx Np . Therefore, the computational complexity of (30) is
2
O(Nrx Ntx Np L), unlike O(Nrx
Np2 Ntx L) of the generic implementation using (R1)-(R3).
(i,∗)

2) Operations (R4) , (R5): It can be noticed from (25) that zbm (t) is invariant with respect
to i and is given by
(i,∗)
(t)
zbm



1
b
p
vecm UNrx C(t)F .
=
Np

(31)

Using the invariance property in (31), a compact version of the first summand in (R4) can be
expressed as
Np
X

vib

2
(i,∗)
zbm
(t)

PNp

i=1

=

Np

i=1

vib


vecm

2

b
UNrx C(t)F


.

(32)

For the second summand in (R4), it can be shown from (28) that




1
(∗,ko )
b
zbm
(t) = √
⊗ aNrx (0) .
vecm F(co ) diag U∗Np b(t)
Nrx

(33)
(∗,k)

Furthermore, as co denotes the column number corresponding to ko (see (29)), zbm

(t) is

invariant ∀k ∈ IaNrx \I(a−1)Nrx , where a ∈ INtx L . To use this invariance for efficient computation
of the second summand in (R4), we define a row vector µc (t) ∈ R1×Ntx L containing the columnwise mean corresponding to Nrx × Ntx L matrix version of {vkc (t)}k∈INc as
X
1
µcn (t) =
vkc (t) .
Nrx

(34)

k∈InNrx \I(n−1)Nrx

With some algebraic manipulation, the second summand in (R4) can be simplified as



Nc
X
2
2
c
(∗,k)
∗ b
c
vk (t) zbm (t) = vecm µ (t) Fdiag(UNp b(t))
⊗ aNrx (0) .

(35)

k=1

To simplify the computations involved in (R5), we expand the summand using (17) as
Nb X
Nc
X

vib (t)vkc (t)

(i,k) 2
zm

=

Nb X
Nc
X

i=1 k=1

(36)

i=1 k=1

PNp
=
where (37) follows from (36) as |Gm,i | =

vib (t)vkc (t) |Gm,i Am,k |2

√1 ,
Np

i=1

N

c
vib (t) X

Np

2

vkc (t) ATk,m ,

(37)

k=1

∀m, i. Besides, as AT = F⊗UNrx , the entries of

AT are invariant within blocks of size Nrx ×Nrx . With arguments similar to the simplifications
involved in the second summand of (R4), (R5) can be efficiently evaluated as
PNp b
Nb X
Nc
X



v (t)
c
(i,k) 2
b
vi (t)vk (t) zm
= i=1 i vecm µc (t) |F|2 ⊗ aNrx (0) .
Np
i=1 k=1

(38)

19

(∗,k)

2

3) Operations (R11) , (R13): From (33), it can be observed that zbm (t)

is fixed for m ∈

IaNrx \ I(a−1)Nrx and k ∈ IbNrx \ I(b−1)Nrx , where a ∈ INp and b ∈ INtx L . The invariance of
(∗,k)

zbm (t)

2

with respect to m and k arise because of the kronecker product with aNrx (0) and the

column invariance in (29). To exploit this property in computing vkr (t) of (R11), we construct a
vector µz (t) ∈ RNp ×1 to contain the column-wise mean corresponding to the Nrx × Np matrix
s
version of {vm
(t)}m∈INp Nrx , i.e.,

µzk (t) =

1
Nrx

X

s
vm
(t)

.

(39)

m∈IkNrx \I(k−1)Nrx

−1
With the above definitions, a simplified version of vjr (t)
in (R11) can be given as



M
X
2
s
(∗,k) 2
∗ b
z
vm (t) zbm
Fdiag(UNp b) µ (t) ⊗ aNrx (0) .
= veck

(40)

m=1

From (31), we rewrite (viq (t))−1 in (R13) as
(viq (t))−1

=

M
X

s
vm
vecm



2



b
UNrx C(t)F

m=1
s



= v vec

2

b
UNrx C(t)F


,

∀i ∈ INp .



b
b = Fdiag U∗ b(t)
and consider the term in
4) Operations (R12) , (R14): We define F(t)
Np
the second summand of (R12) for k = ko . From (28) and (29), we have


M
M
X
X
(co )
(∗,ko )
∗
b
sbm (t)vecm UNrx (ro ) F(t)
sbm (t)b
zm (t) =
.

(41)

m=1

m=1



b
b ∈ CNrx ×Np defined such that sbm (t) = vecm S(t)
, (41) can be expressed as,
With S(t)
M
X

(∗,ko )
sbm (t)b
zm
(t)∗

D
E
(co )
b
b
= S(t), UNrx (ro ) F(t)

(42)

m=1


∗
∗
b
b (co )
= UNrx (ro ) S(t)
F(t)


d∗ .
b F(t)
= vecko U∗Nrx S(t)
The term in the third summand of (R12) can be given by
PNp b
Np
M
M
X
X
X
v (t) |Am,k |2
s
b
(i,k) 2
s
vm (t)
vi (t) zm
=
vm (t) i=1 i
Np
m=1
m=1
i=1
!
PNp b


i=1 vi (t)
T 2 s
veck A v .
=
Np

(43)
(44)

20

Once again exploiting invariance within AT , we efficiently compute the third summand in
(R12) as
M
X

s
vm
(t)

m=1

Np
X

vib (t)

(i,k) 2
zm

PNp

b
i=1 vi (t)
Np

=

i=1

!
veck





|F|2 µz (t) ⊗ aNrx (0) .

(45)

The term in the second summand of (R14) can be rewritten using (25) as
M
X

D
E
(i,∗)∗
∗
b
b
= S(t),
)
.
sbm (t)b
zm
UNrx C(t)Fdiag(U
Np (i)

(46)

m=1

b
For fast implementation of (46), we first compute the column wise inner product between S(t)
b
then perform a Fast Fourier
Transform (FFT). We construct g ∈ CNp ×1 ,
and UNrx C(t)F
and
 



b
b
. Hence, the second summand in (R14) simplifies
such that gk =
S(t)
, UNrx C(t)F
(k)

(k)

to

M
X


(i,∗)∗
sbm (t)b
zm
= veci UNp g .

(47)

m=1

The term in the third summand of (R14) can be expressed as
M X
Nc
X

s
(i,k)
vm
(t)vkc (t) zm

2

=

m=1 k=1

M X
Nc
X

s
vm
(t)vkc (t) |Gm,i Am,k |2

m=1 k=1
Nc
M
X
1 X s
2
=
vkc (t) ATk,m .
vm (t)
Np m=1
k=1

Using the compact form of
M X
Nc
X

(48)

PNc

k=1

vkc (t) ATk,m

s
(i,k)
vm
(t)vkc (t) zm

m=1 k=1

2

=

2

(49)

from (38), we rewrite (49) as

M



1 X s
vm (t)vecm µc (t) |F|2 ⊗ aNrx (0) .
Np m=1

(50)

Furthermore, exploiting the kroenecker product with aNrx (0), we have
M X
Nc
X
m=1 k=1

(i,k)
s
vm
(t)vkc (t) zm

2

=

Nrx c
µ (t) |F|2 µz (t).
Np

(51)

5) Complexity of the simplified operations: Using the complexity of matrix multiplications
and FFTs, the complexity of our simplifications is summarized in Table II. It can be noticed
that the overall complexity of PBiGAMP using our implementation is O (Nrx Np Ntx L), thereby
achieving a speedup factor of Np Nrx compared to the generic implementation. After all possible
simplifications of PBiGAMP for joint estimation, training design is possibly the only frontier
that can be exploited to further reduce the computational complexity.

21

Operation

Complexity

b (t) F
UNrx C

O (Nrx Ntx LNp )

b (t)
U∗Np b

O (Np logNp )

(R3), (R13)

O (Nrx Np )

(R4), (R5), (R11)

O (Ntx LNp )

(R12), (R14)

O (Nrx Ntx LNp )

TABLE II
C OMPLEXITY OF A SINGLE PB I GAMP

ITERATION USING A FAST IMPLEMENTATION .

D. Insights into training design
In this section, we explain how a reasonable training solution that allows a low complexity
implementation does not permit joint estimation. Structured training blocks that aid fast transforms can reduce the complexity of PBiGAMP operations involving multiplications with the
training matrix. Furthermore, such blocks also occupy a lower memory footprint relative to
unstructured ones of the same dimension. For example, training blocks that contain circulantly
shifted rows of a fixed Zadoff-Chu (ZC) sequence were proposed in [12] for fast channel
estimation using EM-GAMP. Similar to EM-GAMP algorithm, circulant training matrices also
b
aid fast matrix multiplications in PBiGAMP. For instance, the computation of C(t)F
in Table II
can be accelerated in every PBiGAMP iteration for the joint estimation. Shifted ZC training
blocks, however, result in a continuum of optimal solutions for bilinear optimization problem.
We illustrate the “CFO propagation effect” to show the trade-off between fast message passing
and identifiability in the joint estimation problem.
To explain the CFO propagation effect, we consider Np = Ntx L pilots and a circulant training
matrix T as per [12], for a narrowband system, i.e., L = 1. Because T∗ is also circulant, the
eigenvectors of T∗ are the columns of the DFT matrix UNtx . Let ΛT be a diagonal matrix
containing the conjugated eigenvalues of T∗ , such that T∗ UNtx = UNtx Λ∗T . With Λ defined as
diag (d (, β)), the noiseless unquantized received block in (6) can be given as
Z = UNrx CU∗ Ntx TΛ
= UNrx CΛT U∗ Ntx Λ .

(52)
(53)

For our analysis, we assume that there is no phase noise, i.e., β = 0 and choose the CFO () to
be an integer multiple of 2π/Np . We define the integer d = Np /2π and interpret Λ as a matrix

22

containing the eigenvalues of the d circulant delay matrix Jd ∈ CNp ×Np , i.e., Jd UNtx = UNtx Λ .
As Λ∗ = Λ− , it can be shown that JNtx −d UNtx = UNtx Λ∗ . Furthermore, as JNtx −d is a real
matrix and UNtx = U∗Ntx , we have U∗Ntx Λ = JNtx −d U∗Ntx . It follows from (53) that
Z = UNrx CΛT JNtx −d U∗ Ntx .

(54)

Comparing (54) with (53), it can be observed that the beamspace matrix C and a CFO of
2πd/Np , result in the same received samples as the beamspace matrix CΛT JNtx −d and zero
CFO. We call this as the CFO propagation effect and study its impact on the identifiability of
the channel and the CFO in joint estimation.
Now, we show that training blocks consisting of circulantly shifted ZC sequences propagate
the CFO into the channel in an inseparable manner. From the perfect autocorrelation property
of ZC sequences, it can be concluded that the diagonal entries of ΛT have constant modulus.
In such case, Λ−1
T is well defined and (54) can be rewritten as
∗
Z = UNrx CΛT JNtx −d Λ−1
T ΛT U Ntx

(55)

∗
= UNrx CΛT JNtx −d Λ−1
T U Ntx T
|
{z
}

(56)

∆

=C()

= UNrx C()U∗ Ntx TΛ0 .

(57)

Notice that Λ0 in (57) is just an identity matrix, also interpreted as a zero CFO perturbation. We
compare (57) with (53) and define a beamspace matrix matrix C() = CΛT JNtx −d Λ−1
T . It can
be noticed that C() has the same sparsity as that of C, as permutation and scaling operations
determined by JNtx −d and ΛT preserve sparsity. Therefore, if (C, ) is a solution to the joint
estimation problem, (C(), 0) is also a solution. In fact, any frequency that is an integer multiple
of 2π/Np can be propagated into the channel matrix so that (C(2πk/Np ),  − 2πk/Np ) is a
solution for every k ∈ INp . Therefore, the circulantly shifted ZC training that aids fast message
passing results in a continuum of optimal solutions for the joint estimation problem.
V. S IMULATIONS
In this section, we provide simulation results for the proposed bilinear message passing based
joint CFO and channel estimation algorithm. We consider a hardware architecture in Fig. 1 that
uses a uniform linear array of antennas and choose Ntx = Nrx = 32. In this work, we assume that
the receiver is equipped with one-bit ADCs, i.e., q = 1. To obtain a performance benchmark for

23

one-bit receivers, we evaluate our algorithm for full resolution receivers, i.e., q = ∞, although
they may not be practical at large bandwidths due to high power consumption. We consider
a mmWave carrier frequency of 38 GHz and an operating bandwidth of W = 100 MHz [36],
which corresponds to a symbol duration of T = 10 ns.
We describe the simulation parameters of the clustered mmWave channel model in (3). We
consider L = 16 taps to model the mmWave channel in the 100 MHz bandwidth. Such model
is good enough for RMS delay spreads upto 33 ns, assuming L = 5W τrms . Our assumption
is reasonable as 95% of the measured RMS delay spreads at 38 GHz were found to be less
than 37.8 ns [37]. The measurements in [37] were made in a UMi-LoS environment using a
resolution of 2 ns. We assume Ncs = 4 clusters, each comprising of 10 rays and the complex ray
gains as IID standard normal random variables. Furthermore, the AoAs and AoDs of the rays
within a cluster are chosen from a laplacian distribution corresponding to an angle spread of

P
2
2
15◦ . The wideband channel is scaled so that E
` kH [`]kF = E kCkF = Nrx Ntx , where the
expectation is taken across several channel realizations. The channel matrix generated with the
aforementioned parameters is practical at mmWave and the angle-delay domain representation,
i.e., C, can be verified to be approximately sparse. Our joint estimation algorithm does not
require any knowledge about the sparsity order of C and learns it using the EM algorithm.
Our PBiGAMP based approach exploits compressibility of the phase error vector in the Fourier
basis. To evaluate the worst case performance of our algorithm, we choose a CFO that is
maximally off grid and within the practical limits. As the resolution of the DFT grid is 2π/Np ,
we choose  = 45π/1024 for a training block of Np = 1024 pilots. The corresponding analog
domain CFO can be verified to be 2.2 MHz, which is about 58 ppm of the carrier frequency.
We set βtx = βrx = 0.047 rad so that the standard deviation of the effective Wiener phase noise
process is β = 0.067 rad. These parameters translate to a phase noise level of about −85 dBc/Hz
at 1 MHz offset for each of the TX and RX oscillators [25], and meet the specifications of a
38 GHz oscillator [38].
The performance of our joint estimation algorithm is evaluated using the Normalised Mean
Square Error (NMSE) of the channel estimate and the mean square error (MSE) of the CFO
estimate . For a given SNR, the variance of the IID Gaussian noise, i.e., σ 2 in (1), is chosen
such that
SNR = 10log10

kTk2F
Np σ 2

!
.

(58)

24

For our simulations, we consider training blocks comprising of IID QPSK entries, IID Gaussian
entries, and shifted ZC sequences proposed in [12]. For mmWave systems, the IID QPSK training
is more practical than the IID Gaussian one, as it can be generated using a TX architecture that
is as simple as analog-beamforming with 2-bit phase shifters.
A. NMSE of the channel estimate
Due to the bilinear nature of the problem, we can only estimate the wideband channel or
equivalently C upto a scale factor. Furthermore, any positive amplification of C results in the
same received block in one-bit receivers at high SNR. Therefore, the NMSE of the channel
estimate is defined as



NMSE = E 

2

b
C − γC
kCk2F

b
where γ is a scalar such that γ = arg min C − aC
a

F



F

,

(59)

b The matrix γ C
b
for a given C and C.

can be considered as the normalized angle-delay domain estimate of the wideband channel. The
function E[.] denotes the empirical expectation and is taken across several channel and training
realizations. For a given realization of the channel and its estimate, we define the Normalized
2
b /kCk2 .
Squared Error (NSE) as C − γ C
F

F

For a sequence of Np = 1024 pilots, our joint estimation algorithm recovered the channel
within acceptable limits with a probability greater than 0.95, at a SNR of 0 dB using IID QPSK
training. The number of outliers in this case is determined by the phase-transition region of
PBiGAMP [22]. It can be observed from Fig. 3 that the probability of successful recovery
monotonically increases as a function of the training length and quickly approaches 1. To ignore
the effect of outliers, the channel NMSE and CFO MSE results we report are averages over
95% of the realizations for Np = 1024. In practice, the failure probability can be lowered by
increasing the training length or by designing a retransmission protocol that accounts for the
failure. As seen in Fig. 4, the NMSE monotonically decreases with the number of pilots. In
practical wireless systems, the choice of the number of pilots is determined by the channel
coherence time [3].
In Fig. 5, we plot the NMSE as a function of the SNR for various training sequences. It can
be observed that the reconstruction error is approximately the same for IID QPSK training and
IID Gaussian training matrices. For the one-bit case, the NMSE saturates at high SNR because

25

1
0.9
0.8

Pr(NSE 5 x)

0.7
0.6
0.5
0.4
0.3

Np
Np
Np
Np

0.2
0.1

= 512
= 1024
= 2048
= 4096

0
-15

-12.5

-10

-7.5

-5

-2.5

0

x(dB)

Fig. 3. The empirical CDF of the Normalized Squared Error of the channel estimate obtained using an IID QPSK training,
at a SNR of 0 dB in a one-bit receiver. The reconstruction performance in terms of the recovery probability and the mean
monotonically improve with the number of pilots.
-5

IID QPSK, 1-bit
IID QPSK, 1-bit

-6
-7

NMSE (dB)

-8
-9
-10
-11
-12
-13
-14
-15
512

1024

2048

3072

4096

5120

Number of pilots (Np )

Fig. 4. The NMSE of the channel estimate as a function of the training length, for an IID QPSK training at a SNR of 0 dB.
The NMSE monotonically decreases with the training length for one-bit and full resolution receivers. Due to the quantization
noise in one-bit receivers, NMSE in the one-bit case is higher than the full resolution one.

the recovery performance is limited by the quantization noise. As shown in Fig. 4, the channel
reconstruction error in one-bit receivers can be further decreased by using a higher number of
pilots for the training. It can be noticed from Fig. 5 that joint estimation with circulantly shifted
ZC sequences proposed in [12] performs poorly. The failure due to such structured matrices can
be attributed to the CFO propagation effect and confirms with our analysis in Section IV-D.

26

0

NMSE (dB)

-5

-10

-15

-20

-25
-10

IID QPSK, 1-bit
IID Gaussian , 1-bit
Shifted ZC, 1-bit
IID QPSK , 1-bit
IID Gaussian , 1-bit
Shifted ZC, 1-bit
-5

0

5

10

15

20

SNR(dB)

Fig. 5. NMSE of the channel estimate obtained with one-bit and full resolution channel measurements for Np = 1024. Joint
estimation is possible with IID Gaussian and IID QPSK training, but not with shifted ZC training due to the CFO propagation
effect discussed in Section IV-D. It can be observed that the NMSE for the one-bit case saturates at high SNR due to quantization
noise.

B. MSE of the CFO estimate
The CFO in our algorithm is obtained using an Extended Kalman Filter [34] on the inverse
b If b
DFT of b.
 is the estimate of the CFO in the digital domain, the MSE of the CFO estimate


is given by E ( − b
)2 . For a sequence of 1024 pilots, the MSE of b
 is shown as a function of
the SNR in Fig. 6, for IID QPSK and IID Gaussian training matrices. It can be noticed from
Fig. 6 that the MSE of the CFO estimate saturates even for the full resolution case because the
performance of the EKF is limited by the phase noise at high SNR [39]. Similar to the NMSE
of the channel estimate, the MSE of b
 is expected to saturate at high SNR for one-bit receivers
due to quantization noise. The CFO estimation error, however, is determined by the phase noise
[39] because the CFO MSE for the one-bit case approaches that of the full resolution one. As
expected, the CFO MSE decreases with the number of pilots and is shown in Fig. 7.
C. Performance invariance with the CFO
In this section, we show that the performance of our joint estimation approach is invariant
to the CFO within practical limits. As an example, we evaluate our algorithm for different
values of the CFO in the range [−40 ppm, 40 ppm] of f1 . These limits were chosen according
to the IEEE 802.11ad specifications. From Fig. 8, we see that the reconstruction errors in the

27

10-1

IID
IID
IID
IID

MSE of CFO estimate

10-2

QPSK , 1-bit
Gaussian , 1-bit
QPSK , 1-bit
Gaussian , 1-bit

10-3

10-4

10-5

10-6
-10

-8

-6

-4

-2

0

5

10

SNR(dB)

Fig. 6. MSE of the CFO estimate as a function of SNR for a training length of 1024. A maximally off-grid CFO of  = 45π/1024
was chosen to evaluate the performance. At high SNR, the MSE saturates for both the one-bit case and the infinite resolution
case due to phase noise [39].

MSE of the CFO estimate

10-2

IID QPSK, 1-bit
IID QPSK, 1-bit

10-3

10

-4

10-5

10-6
512

1024

2048

3072

4096

Number of pilots (Np )

Fig. 7. MSE of the CFO estimate as a function of the training length for a SNR of 0 dB. Here, the CFO in the system was fixed
to  = 45π/1024, and the phase noise standard deviation to β = 0.067 rad. The MSE decreases with the number of pilots and
the performance gap between the one-bit case and the full-resolution case is negligible. The sharp decrease in the CFO MSE
can be attributed to the phase-transition effect of PBiGAMP [22].

channel is constant across the practical range of the CFO. The invariance arises due to the use
of Bernoulli-Gaussian prior for the synchronization variables, that achieves robustness against
off-grid leakage effects.

28

-6

IID QPSK, 1-bit
IID QPSK, 1-bit

-7

NMSE (dB)

-8

-9

-10

-11

-12
-40

-30

-20

-10

0

10

20

30

40

CFO in ppm of f1

Fig. 8. For a SNR of 0 dB, the plot shows the invariance of the channel NMSE over the practical range of the CFO. The phase
noise standard deviation was set to 0.067 rad, and 1024 pilots were used for the joint estimation. It can be observed from the
plot that our joint estimation algorithm is robust to leakage effects that arise due to an off-grid CFO.

VI. C ONCLUSIONS
Most sparsity-aware channel estimation algorithms for mmWave systems assume perfect synchronization and perform poorly in the presence of such errors. In this work, we propose a
message passing based algorithm that can leverage the sparsity of the wideband mmWave channel
in addition to the compressibility of the phase error vector. We exploit the structure in the joint
estimation problem to provide a low complexity implementation of a bilinear message passing
algorithm. Unlike the existing methods that are specific to certain hardware architectures or use
non-coherent techniques, our technique can be adapted to perform joint estimation with other
mmWave architectures.
Simulation results show that it is possible to perform joint CFO and channel estimation
using IID QPSK and IID Gaussian training matrices. The CFO propagation effect proposed
in this work shows that there is a trade-off between fast message passing using structured
training and identifiability in joint estimation. In our future work, we will also consider frame
synchronization while ensuring low complexity and scalability of our algorithm to different
mmWave architectures.
ACKNOWLEDGMENT
The authors would like to thank Jianhua Mo and Philip Schniter for helping them in getting
started with the GAMP algorithm.

29

R EFERENCES
[1] S. Rangan, T. S. Rappaport, and E. Erkip, “Millimeter-wave cellular wireless networks: Potentials and challenges,” in Proc.
of the IEEE, vol. 102, no. 3, pp. 366–385, 2014.
[2] A. Alkhateeb, J. Mo, N. Gonzalez 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, 2014.
[3] R. W. Heath, Introduction to Wireless Digital Communication: A Signal Processing Perspective. Pearson Education, 2017.
[4] R. W. Heath, N. Gonzalez-Prelcic, S. Rangan, W. Roh, and A. M. Sayeed, “An overview of signal processing techniques
for millimeter wave MIMO systems,” IEEE J. Sel. Topics Signal Process., vol. 10, no. 3, pp. 436–453, 2016.
[5] E. J. Candès and M. B. Wakin, “An introduction to compressive sampling,” IEEE Signal Process. Mag., vol. 25, no. 2,
pp. 21–30, 2008.
[6] P. T. Boufounos and R. G. Baraniuk, “One bit compressive sensing,” in Proc. of the IEEE 42nd Annual Conf. on Inform.
Sciences and Systems, 2008., 2008, pp. 16–21.
[7] Z. Marzi, D. Ramasamy, and U. Madhow, “Compressive channel estimation and tracking for large arrays in mm-wave
picocells,” IEEE J. Sel. Topics Signal Process., vol. 10, no. 3, pp. 514–527, 2016.
[8] A. Alkhateeb, O. El Ayach, G. Leus, and R. W. Heath, “Channel estimation and hybrid precoding for millimeter wave
cellular systems,” IEEE J. Sel. Topics Signal Process., vol. 8, no. 5, pp. 831–846, 2014.
[9] J. Mo, P. Schniter, N. G. Prelcic, and R. W. Heath, “Channel estimation in millimeter wave MIMO systems with one-bit
quantization,” in Proc. of the Asilomar Conf. on Signals, Systems and Computers, 2014, pp. 957–961.
[10] D. L. Donoho, A. Maleki, and A. Montanari, “Message-passing algorithms for compressed sensing,” Proc. of the National
Academy of Sciences, vol. 106, no. 45, pp. 18 914–18 919, 2009.
[11] J. P. Vila and P. Schniter, “Expectation-maximization Gaussian-mixture approximate message passing,” IEEE Trans. on
Signal Process., vol. 61, no. 19, pp. 4658–4672, 2013.
[12] J. Mo, P. Schniter, and R. W. Heath Jr., “Channel estimation in broadband millimeter wave MIMO systems with few-bit
ADCs,” IEEE Trans. on Signal Process., 2017.
[13] H. Yan and D. Cabric, “Compressive sensing based initial beamforming training for massive MIMO millimeter-wave
systems,” in Proc. of the IEEE Global Conf. on Signal and Info. Process. (GlobalSIP), 2016, pp. 620–624.
[14] N. J. Myers and R. W. Heath Jr., “A compressive channel estimation technique robust to synchronization impairments,” in
Proc. of the 18th IEEE Intl. Workshop on Signal Process. Adv. in Wireless Commun. (SPAWC), 2017, pp. 1–5.
[15] J. Chen, Y.-C. Wu, S. Ma, and T.-S. Ng, “Joint CFO and channel estimation for multiuser MIMO-OFDM systems with
optimal training sequences,” IEEE Trans. on Signal Process., vol. 56, no. 8, pp. 4008–4019, 2008.
[16] W. Zhang, Q. Yin, W. Wang, and F. Gao, “One-shot blind CFO and channel estimation for OFDM with multi-antenna
receiver,” IEEE Trans. on Signal Process., vol. 62, no. 15, pp. 3799–3808, 2014.
[17] O. Abari, H. Hassanieh, M. Rodreguez, and D. Katabi, “Millimeter wave communications: From point-to-point links to
agile network connections,” in Proc. of the ACM Workshop on Hot Topics in Networks, 2016, pp. 169–175.
[18] M. E. Rasekh, Z. Marzi, Y. Zhu, U. Madhow, and H. Zheng, “Noncoherent mmwave path tracking,” in Proc. of the 18th
ACM Intl. Workshop on Mobile Computing Systems and Applications, 2017, pp. 13–18.
[19] N. J. Myers and R. W. Heath Jr., “Joint CFO and channel estimation in millimeter wave systems with one-bit ADCs,” in
Proc. of the 7th IEEE Intl. Workshop on Comput. Adv. in Multi-Sensor Adapt. Process. (CAMSAP), 2017, pp. 1–5.
[20] N. D. Sidiropoulos, L. De Lathauwer, X. Fu, K. Huang, E. E. Papalexakis, and C. Faloutsos, “Tensor decomposition for
signal processing and machine learning,” IEEE Trans. on Signal Process., vol. 65, no. 13, pp. 3551–3582.

30

[21] E. J. Candes, T. Strohmer, and V. Voroninski, “Phaselift: Exact and stable signal recovery from magnitude measurements
via convex programming,” Commun. on Pure and Applied Mathematics, vol. 66, no. 8, pp. 1241–1274, 2013.
[22] J. T. Parker and P. Schniter, “Parametric bilinear generalized approximate message passing,” IEEE J. Sel. Topics Signal
Process., vol. 10, no. 4, pp. 795–808, 2016.
[23] Y. Fu, C. Tellambura, and W. A. Krzymien, “Limited-feedback precoding for closed-loop multiuser MIMO OFDM systems
with frequency offsets,” IEEE Trans. on Wireless Commun., vol. 7, no. 11, 2008.
[24] V. Syrjala and M. Valkama, “Receiver DSP for OFDM systems impaired by transmitter and receiver phase noise,” in Proc.
of the IEEE Intl. Conf. on Commun. (ICC), 2011, pp. 1–6.
[25] A. Demir, A. Mehrotra, and J. Roychowdhury, “Phase noise in oscillators: A unifying theory and numerical methods
for characterization,” IEEE Trans. on Circuits and Systems I: Fundamental Theory and Applications, vol. 47, no. 5, pp.
655–674, 2000.
[26] N. G. Prelcic, A. Ali, V. Va, and R. W. Heath Jr., “Millimeter-wave communication with out-of-band information,” IEEE
Commun. Magazine, vol. 55, no. 12, pp. 140–146, 2017.
[27] K. Venugopal, A. Alkhateeb, N. G. Prelcic, and R. W. Heath, “Channel estimation for hybrid architecture-based wideband
millimeter wave systems,” IEEE Journal on Selected Areas in Commun., vol. 35, no. 9, pp. 1996–2009, 2017.
[28] J. Wang, “Beam codebook based beamforming protocol for multi-Gbps millimeter-wave WPAN systems,” IEEE Journal
on Selected Areas in Commun., vol. 27, no. 8, 2009.
[29] M. F. Duarte and R. G. Baraniuk, “Spectral compressive sensing,” Applied and Computational Harmonic Analysis, vol. 35,
no. 1, pp. 111–129, 2013.
[30] S. Ling and T. Strohmer, “Self-calibration and biconvex compressive sensing,” Inverse Problems, vol. 31, no. 11, p. 115002,
2015.
[31] S. Rangan, “Generalized approximate message passing for estimation with random linear mixing,” in Proc. of the IEEE
Intl. Symposium on Inform. Theory, 2011, pp. 2168–2172.
[32] S. M. Kay, Fundamentals of statistical signal processing.

Prentice Hall PTR, 1993.

[33] F. R. Kschischang, B. J. Frey, and H.-A. Loeliger, “Factor graphs and the sum-product algorithm,” IEEE Trans. on Inform.
Theory, vol. 47, no. 2, pp. 498–519, 2001.
[34] P. Zhao, L. Kuang, and J. Lu, “Carrier frequency offset estimation using extended Kalman filter in uplink OFDMA systems,”
in Proc. of the IEEE Intl. Conf. on Commun., 2006., vol. 6, pp. 2870–2874.
[35] C.-K. Wen, C.-J. Wang, S. Jin, K.-K. Wong, and P. Ting, “Bayes-optimal joint channel-and-data estimation for massive
MIMO with low-precision ADCs,” IEEE Trans. on Signal Process., vol. 64, no. 10, pp. 2541–2556.
[36] F. B. Mismar and B. L. Evans, “Machine learning approach to estimating mmWave signal measurements during handover,”
arXiv preprint arXiv:1710.01879, 2017.
[37] Aalto Univ., AT & T, BUPT, CMCC, Ericsson, Huawei, Intel, KT Corporation, Nokia, NTT DOCOMO, New York Univ.,
Qualcomm, Samsung, Univ. of Bristol, and Univ. of Southern California, “5G channel model for bands up to 100 GHz
(v2.3), Tech. Rep.” http://www.5gworkshops.com/5gcm.html, 2016, [Online].
[38] “MMIC Voltage Controlled Oscillator,” http://www.tlcprecision.com/Data%20Sheets/TLCO01981.pdf, 2013, [Online].
[39] O. H. Salim, A. A. Nasir, H. Mehrpouyan, W. Xiang, S. Durrani, and R. A. Kennedy, “Channel, phase noise, and frequency
offset in OFDM systems: Joint estimation, data detection, and hybrid Cramer-Rao lower bound,” IEEE Trans. on Commun.,
vol. 62, no. 9, pp. 3311–3325, 2014.

