1

An Approximate ML Detector for MIMO
Channels Corrupted by Phase Noise

arXiv:1708.02321v1 [] 7 Aug 2017

Richard Combes and Sheng Yang

Abstract
We consider the multiple-input multiple-output (MIMO) communication channel impaired by phase
noises at both the transmitter and receiver. We focus on the maximum likelihood (ML) detection problem
for uncoded single-carrier transmission. We derive an approximation of the likelihood function, based
on which we propose an efficient detection algorithm. The proposed algorithm, named self-interference
whitening (SIW), consists in 1) estimating the self-interference caused by the phase noise perturbation,
then 2) whitening the said interference, and finally 3) detecting the transmitted vector. While the exact
ML solution is computationally intractable, we construct a simulation-based lower bound on the error
probability of ML detection. Leveraging this lower bound, we perform extensive numerical experiments
demonstrating that SIW is, in most cases of interest, very close to optimal with moderate phase noise.
More importantly and perhaps surprisingly, such near-ML performance can be achieved by applying
only twice the nearest neighbor detection algorithm. In this sense, our results reveal a striking fact: nearML detection of phase noise corrupted MIMO channels can be done as efficiently as for conventional
MIMO channels without phase noise.

Index Terms
MIMO systems, phase noise, maximum likelihood detection, probability of error.

I. I NTRODUCTION
We consider the signal detection problem for the following discrete-time multiple-input multipleoutput (MIMO) channel


y = diag ejθr,1 , . . . , ejθr,nr H diag ejθt,1 , . . . , ejθt,nt x + z ,

(1)

The authors are with the Laboratory of Signals and Systems (L2S), CentraleSupélec, 3 rue Joliot-Curie, 91190 Gif-sur-Yvette,
France. (e-mail: richard.combes,sheng.yang@centralesupelec.fr)

August 9, 2017

DRAFT

2

where H ∈ Cnr ×nt is the channel matrix known to the receiver; z ∈ Cnr ×1 represents a realization
of the additive noise whereas θt,l and θr,k are the phase noises at the l th transmit antenna and
the k th receive antenna, respectively; the input vector x ∈ Cnt ×1 is assumed to be carved
from a quadratic amplitude modulation (QAM). The goal is to estimate x from the observation
y ∈ Cnr ×1 , with only statistical knowledge on the additive noise and the phase noises.
In the case where the phase noise is absent, the problem is well understood, and the maximum likelihood (ML) solution can be found using any nearest neighbor detection (NND)
algorithm (see [1] and the references therein). In particular, the sphere decoder [2] has been
shown to be very efficient [3], so that its expected complexity (averaged over channel realizations)
is polynomial in the problem dimension nt . Furthermore, there exist approximate NND algorithms (e.g., based on lattice reduction) that achieve near-ML performance when applied for
MIMO detection [4].
The presence of phase noise in (1) is both a practical and long-standing problem in communication. In their seminal paper [5] back in the 70’s, Foschini et al. used this model to
capture the residual phase jitter at the phase-locked loop of the receiver side, and investigated
both the performance of decoding algorithms as well as constellation design in the scalar
case (nt = nr = 1). As a matter of fact, in most wireless communication systems, phase noise
is present due to the phase and frequency instabilities in the radio frequency oscillators used at
both the transmitter and the receiver [6]. The channel (1) can be seen as a valid mathematical
model when the phase noise varies slowly as compared to the symbol duration.1 While phase
noise can be practically ignored in conventional MIMO systems, its impact becomes prominent
at higher carrier frequencies since it can be shown that phase noise power increases quadratically
with carrier frequency [6], [8]. The performance degradation due to phase noise becomes even
more severe with the use of higher order modulations for which the angular separation between
constellation points can be small. At medium to high SNR, phase noise dominates additive noise,
becoming the capacity bottleneck [9], [10]. As for signal detection, finding the ML solution for
the MIMO phase noise channel (1) is hard in general. Indeed, unlike for conventional MIMO
channels, the likelihood function of the transmitted signal cannot be obtained in closed form.
1

As pointed out in [7] and the references therein, an effective discrete-time channel is usually obtained from a waveform

phase noise channel after filtering. When the continuous-time phase noise varies rapidly during the symbol period, the filtered
output also suffers from amplitude perturbation. More discussion is provided in Section VI-D.

August 9, 2017

DRAFT

3

Our Contribution. In this work, we propose an efficient MIMO detection algorithm which
finds an approximate ML solution in the presence of phase noise. The main contributions of this
work are summarized as follows.
(i) We derive a tractable approximation of the likelihood function of the transmitted signal.
While the exact likelihood does not have a close-form expression, the proposed approximation has a simple form and turns out to be accurate for weak to medium phase noises.
(ii) Since maximizing the approximate likelihood function over a discrete signal set is still
hard, we propose a heuristic method that finds an approximate solution by applying twice
the nearest neighbor detection algorithm. The proposed algorithm, called self-interference
whitening (SIW), has a simple geometric interpretation. Intuitively, the phase noise perturbation generates self-interference that depends on the transmitted signal through the
covariance matrix. The main idea is to first estimate the covariance of the self-interference
with a potentially inaccurate initial signal solution, then perform the whitening with the
estimated covariance, followed by a second detection. From the optimization point of view,
our algorithm can be seen as a (well-chosen) concave approximation to a non-concave
objective function.
(iii) We assess the performance of SIW and competing algorithms in different communication
scenarios. Since the error probability of ML decoding is unknown, we propose a simulationbased lower bound which we use as a benchmark. Simulation results show that SIW
achieves near ML performance in most scenarios. In this sense, our work reveals that
near optimal MIMO detection with phase noise can be done as efficiently as without phase
noise. Although the likelihood approximation is derived using the assumption that the phase
noise has a Gaussian distribution, our numerical experiments show that SIW works well
even with non-Gaussian phase noises.
Related Work. Receiver design with phase noise mitigation has been extensively investigated
in the past years (see [11], [12] and references therein). More complex channel models, including
multi-carrier systems (e.g., OFDM) and time-correlated phase noises (e.g., the Wiener process)
have also been considered. In particular, joint data detection and phase noise estimation algorithms have been proposed in [13], [11]. A phase noise estimation based scheme to improve the
system performance for smaller alphabets has been proposed in [12]. However, the challenging
problem of signal detection in MIMO phase noise channels using higher order modulation, where

August 9, 2017

DRAFT

4

performance is extremely sensitive to phase noise, has not been addressed adequately before.
The remainder of the paper is organized as follows. We start with a formal description of
the problem in the next section. The approximation of the likelihood function is derived in
Section III-A, followed by the proposed algorithm described in Section III-B. The hardness of
finding the exact ML solution is investigated in Section IV. We present the numerical experiments
in Section V. Further discussion on the proposed algorithm and relevance of the considered
channel model is provided in Section VI. Section VII concludes the paper.
II. A SSUMPTIONS AND P ROBLEM F ORMULATION
Notation
Throughout the paper, we use the following notation. For random quantities, we use upper
case letters, e.g., X, for scalars, upper case letters with bold and non-italic fonts, e.g., V , for
vectors, and upper case letter with bold and sans serif fonts, e.g., M , for matrices. Deterministic
quantities are denoted in a rather conventional way with italic letters, e.g., a scalar x, a vector
v , and a matrix M . The Euclidean norm of a vector v is denoted by kvv k, respectively. The
transpose and conjugated transpose of M are M T and M H , respectively.
A. System model
In the following, we describe formally the channel model (1) presented in the previous section.
We assume a MIMO channel with nt transmit and nr receive antennas. Let H denote the
channel matrix, where the (k, l)-th element of H , denoted as hk,l , represents the channel gain
between the l th transmit antenna and k th receive antenna. The transmitted vector is denoted
by x = [x1 , . . . , xnt ]T , where xl ∈ X , l = 1, . . . , nt , X being typically a QAM constellation
P
with normalized average energy, i.e., |X1 | x∈X |x|2 = 1. For a given transmitted vector x , the
received vector in base-band can be written as the following random vector
Y = Λ RH Λ T x + Z ,

(2)



where the diagonal matrices Λ R := diag ejΘr ,1 , . . . , ejΘr ,nr and Λ T := diag ejΘt,1 , . . . , ejΘt,nt
capture the phase perturbation at the receiver and transmitter, respectively; Z is the additive
white Gaussian noise (AWGN) vector with Z ∼ CN (0, γ −1I ), where γ is the nominal signalto-noise ratio (SNR). We assume that the phase noise Θ := [Θt,1 · · · Θt,nt Θr,1 · · · Θr,nr ]T is
jointly Gaussian with Θ ∼ N (0, Q θ ) where the covariance matrix Q θ can be arbitrary. Note
August 9, 2017

DRAFT

5

that this model includes as a special case the uplink channel in which nt is the number of
single-antenna users. In such a case, the transmit phase noises are independent. For simplicity,
we consider uncoded transmission in which each symbol xl can take any value from X with
equal probability.
Further, we assume that the channel matrix can be random but is perfectly known at the
receiver, whereas such knowledge at the transmitter side is irrelevant in uncoded transmission.
We also define H Θ := Λ RH Λ T and accordingly H θ for some realization of Θ = θ . By definition,
we have H 0 = H . Finally, we ignore the temporal correlation of the phase noise process and
the channel process, and focus on the spatial aspect of the signal detection problem.
B. Problem formulation
With AWGN, we have the following conditional probability density function (pdf)
p(yy | x , θ , H ) =

γ nr −γkyy −H
Hθ x k2
e
,
π nr

(3)

from which we obtain the likelihood function by integrating over Θ


p(yy | x , H ) = EΘ p(yy | x , Θ , H )
i
 h
γ nr
y −H
HΘ x k2
−γky
+ ln nr .
= ln EΘ e
π

(4)
(5)

The ML detector finds an input vector from the alphabet X nt such that the likelihood function
is maximized. In practice, it is often more convenient to use the log-likelihood function as the
objective function, i.e., after removing a constant term
 h
i
y −H
HΘ x k2
−γky
x
y
H
Q
f (x , , , γ, θ ) := ln EΘ e
,

(6)

where the arguments γ and Q θ can be omitted whenever confusion is not likely. Thus,
xML (yy , H ) := arg maxn f (x
x, y , H ).
x̂
x ∈X

t

(7)

From (7) we see two main challenges to compute the optimal solution:
1) The expectation in (7) cannot be obtained in closed form. A numerical implementation is
equivalent to finding the numerical integral in nt + nr dimensions. This can be extremely
hard in high dimensions.
2) The size of the optimization space, |X |nt , can be prohibitively large when the modulation
size |X | and the input dimension nt become large.

August 9, 2017

DRAFT

6

In Section IV, we examine in more details why both of these issues are indeed challenging.
In a conventional MIMO channel without phase noise, finding the ML solution is reduced to
solving the following problem
x0ML (yy , H ) := arg minn kyy − H 0 x k2 ,
x̂
x ∈X

t

(8)

which is also called the minimum Euclidean distance detection or nearest neighbor detection (NND). Although the search space in (8) remains large, the expectation is gone. Furthermore,
since the objective function is the Euclidean distance, efficient algorithms (e.g., sphere decoder [2]
or lattice decoder [1]) exploiting the geometric structure of the problem can be applied without
searching over the whole space X nt . It is shown in [3] that the sphere decoder has a polynomial
average complexity with respect to the input dimension nt when the channel matrix is drawn
i.i.d. from a Rayleigh distribution.
In practice, one may simply ignore the existence of phase noise and still apply (8) to obtain
x0ML which we refer to as the naive ML solution in our work. While this can work relatively
x̂
well when the phase noise is close to 0, it becomes highly suboptimal with stronger phase noise
which is usually the case in high frequency bands with imperfect oscillators. In this paper, we
provide a near ML solution by circumventing the two challenges mentioned earlier. We first
propose an approximation of the likelihood function. Then we propose an algorithm to solve
approximately the optimization problem (7).
III. P ROPOSED S CHEME
A. Proposed Approximation of the Likelihood Function
In this section, we propose to approximate the likelihood function with the implicit assumption
that the phase noise is not large. Indeed, in practice, the standard deviation of the phase noise
is typically smaller than 10 degrees ≈ 0.174 rad. For stronger phase noises, it is no longer
reasonable to use QAM and the problem should be addressed differently.
ΛRH Λ Tx k = kΛ
ΛH
H Λ Tx k,
The likelihood function (3) depends on the Euclidean norm kyy −Λ
R y −H
in which the difference vector ΛH
R y − H Λ Tx can be rewritten and approximated as follows


jθθ t
e

HDx Dy] 
ΛH
(9)
R y − H Λ Tx = [−H
e−jθθ r
HDx Dy] θ,
≈ (yy − H x ) − j[H

August 9, 2017

(10)

DRAFT

7


where we define D x := diag(x1 , . . . , xnt ), D y := diag(y1 , . . . , ynr ), and recall that θ := θ Tt

θ Tr

T

(10) is from the linear approximation2 ejθθ = 1 + jθθ + o(θθ ). Thus the Euclidean norm has the
corresponding real approximation:
Aθ + b k2 ,
kyy − Λ RH Λ Tx k2 ≈ kA
where A ∈ R2nr ×(nt +nr ) and b ∈ R2nr ×1 are defined as




H D x)
Dy )
Im (H
Im (D
Re (yy − H x )
 , b := 
.
A := 
H D x ) − Re (D
Dy )
− Re (H
Im (yy − H x )

(11)

(12)

With the above approximation, we can derive the approximation of the log-likelihood function.
Proposition 1. Let A and b be defined as in (12). Then we have the following approximation of
 h
i
y −H
H θ xk2
−γky
x, y , H , γ, Q θ ) with
the log-likelihood function ln EΘ e
≈ fˆ(x
1
x, y , H ) := −γ b TWx−1b − ln det (W
Wx ) ,
fˆ(x
2

(13)

AQ θA T .
Wx := I + 2γA

(14)

where Wx is defined as

Hence, the proposed approximate ML (aML) solution is


1
T
−1
xaML (yy , H ) := arg minn γ b Wx b + ln det (W
Wx ) .
x̂
x ∈X t
2

(15)

Proof. Since we assume that Θ ∼ N (0, Q θ ), we have


Z
h
i
1
1 −1
T
T
T
AΘ +bbk2
−γkA
2
A A + Q θ ) θ − 2γbb Aθ − γkbbk
EΘ e
=p
d θ exp −θθ (γA
(16)
Qθ ) θ
det(2πQ
{z 2
}
|
1 −1
Q
2

s


Q)
det(Q
Q)A
ATb
exp −γkbbk2 + γ 2b TA (2Q
Qθ )
det(Q


Z
1
1 −1
T
T
T
Q)A
A b ) Q (θθ + γ(2Q
Q)A
A b)
· dθ p
exp −(θθ + γ(2Q
2
Q)
det(2πQ
θ
s

Q)
det(Q
Q)A
ATb
=
exp −γkbbk2 + γ 2b TA (2Q
Qθ )
det(Q
s

Q)
det(Q
A((2Q
Qθ )−1 + γA
ATA )−1A T )bb
=
exp −γbbT (II − γA
Qθ )
det(Q

=

2

(17)

(18)

(19)

Here we use, with a slight abuse of notation, ejθθ to denote the vector obtained from the element-wise complex exponential

operation. Similarly, the little-o Landau notation o(θθ ) is element-wise.

August 9, 2017

DRAFT

;

8

(a) 4-QAM, f = fˆ = −10.

(b) 16-QAM, f = fˆ = −4.

(c) 64-QAM, f = fˆ = −1.6.

Fig. 1: The proposed approximation of the likelihood function in the scalar case. Solid line is
the actual likelihood level set, dashed line is the approximation. Here γ = 30dB and phase noise
has standard deviation 3◦ at the transmitter and at the receiver.

s
=


Q)
det(Q
A(2Q
Qθ )A
AT )−1b
exp −γbbT (II + γA
Qθ )
det(Q


1
T
T −1
b
I
A
Q
A
b
=p
exp
−γb
(I
+
2γA
)
,
θ
AQ θA T )
det(II + 2γA

(20)
(21)

where (18) holds since the integrand in (17) is a pdf with respect to θ ; (20) is from the Woodbury
−1
V . Taking the logarithm and we obtain
matrix identity (II + U C V )−1 = I − U C −1 + V U
the approximated log-likelihood function (13).
In Figure 1, we illustrate the proposed approximation for 4-, 16-, and 64-QAM. In all three
cases, we plot for each constellation point a level set of the likelihood function with respect to “yy ”
in solid line. The level sets of the approximated likelihood function are plotted similarly in dashed
line. While the likelihood function is evaluated using numerical integration, the approximation is
in closed form given by (13). In this figure, we observe that the approximation is quite accurate,
especially for signal points with smaller amplitude. Further, the resemblance of the level sets for
the approximate likelihood to ellipsoids suggests that the main contribution in the right hand side
of (13) comes from the first term −γ b TWx−1b . We shall exploit this feature later on to construct
the proposed algorithm.
Qθ → 0, i.e., when the phase noise vanishes faster than
Remark 1. We can check that when γQ
the AWGN does, the above solution minimizes kbbk2 = kyy − H x k2 , which corresponds to the
optimal NND solution in the conventional MIMO case. Nevertheless, when the phase noise does

August 9, 2017

DRAFT

9

not vanish, the approximate log-likelihood function (13) depends on x in a rather complex way,
due to the presence of the matrix Wx . Focusing on the term b TWx−1b , we can think of Wx as
the covariance matrix of some equivalent noise. Indeed, if we approximate the multiplicative
phase noise as an additive perturbation, then the perturbation is a self-interference that depends
on the input vector x . This perturbation is not isotropic nor circularly symmetric, and can be
captured by the covariance matrix Wx . Some more discussions in this regard will be given in
the following subsection.
While the proposed approximation simplifies significantly the objective function, the optimization problem (15) remains hard when the search space is large. For instance, with 64-QAM
and 4 × 4 MIMO, the number of points in Xtn is more than 107 ! Therefore, we need further
simplification by exploiting the structure of the problem.
B. The Self-Interference Whitening Algorithm
As mentioned above, the difficulty of the optimization (15) is mainly due to the presence of
the matrix Wx that depends on x . Let us first assume that the Wx corresponding to the optimal
xaML were somehow known, and is denoted by Wx̂x . Then the optimization problem (15)
solution x̂
would be equivalent to


1
T
−1
Wx̂x )
xaML (yy , H ) = arg minn γ b Wx̂x b + ln det (W
x̂
x ∈X t
2
= arg minn b TWx̂x−1b
x ∈X

(23)

t

−1

−1

Wx̂x 2 ỹy − Wx̂x 2 H̃
H x̃
x k2 ,
= arg minn kW
x ∈X

− 21

where Wx̂x

(22)

t



−1 H
−1
x, ỹy , and H̃
H are defined by
is any matrix such that Wx̂x 2 Wx̂x 2 = Wx̂x−1 ; x̃






x)
H ) − Im (H
H)
Re (x
Re (yy )
Re (H
 , ỹy := 
 , H̃
.
x := 
H := 
x̃
x)
H ) Re (H
H)
Im (x
Im (yy )
Im (H

(24)

(25)

Note that for a given Wx̂x , (24) can be solved efficiently with any NND algorithm. Unfortunately,
xaML , the exact Wx̂x cannot be found. Therefore, the idea is
without knowing the optimal solution x̂
x, and then solve the optimization
to first estimate the matrix Wx̂x with some suboptimal solution x̂
problem (24) with a NND. We call this two-step procedure self-interference whitening (SIW).

August 9, 2017

DRAFT

10

x0ML as the initial estimate to obtain Wx̂x , and
For instance, we can use the naive ML solution x̂
have
−1

−1

ML

ML

x0aML (yy , H ) = arg minn kW
Wx̂x0 2 ỹy − Wx̂x0 2 H̃
H x̃
x k2 .
x̂
x ∈X

t

(26)

Remark 2. The intuition behind the SIW scheme is as follows. From the definition of Wx in
(14) and A in (12), we see that Wx depends on x only through H D x . First, the column space
of H D x does not vary with x since D x is diagonal. Second, a small perturbation of x does
x0ML is close to the
not perturb Wx too much in the Euclidean space. Since the naive ML point x̂
actual point x in the column space of H , it provides an accurate estimate of Wx . This can also
be observed on Figure 1c, where we see that the ellipsoid-like dashed lines have similar sizes
and orientations for constellation points that are close to each other.
Remark 3. Another possible initial estimate is the naive linear minimum mean square error (LMMSE) solution. As the naive ML, the naive LMMSE ignores the phase noise and returns
H H (γ −1I + H H H )−1y − x k2 .
x0LMMSE (yy , H ) := arg minn kH
x̂
x ∈X

t

(27)

It is worth mentioning that in the presence of phase noise the naive LMMSE is not necessarily
dominated by the naive ML solution, as will shown in the numerical experiments of Section V.
The main algorithm of this work is described in Algorithm 1. In the algorithm, the complex
function NND(yy , H , X ) finds among the points from the alphabet X the closest one to y in
H , X̃ ) is the real counterpart of NND. The
the column space of H ; the function realNND(ỹy , H̃
x0 )” embeds the real vector x̃
x0 to the complex space by taking the upper
function “complex(x̃
half as the real part and the lower half as the imaginary part. It is worth noting that the newly
obtained point is accepted only when it has a higher approximate likelihood value than the naive
ML point does. An example of the scalar case is provided in Figure 2 where 256-QAM is used.
The transmitted point is x and the received point is y. The solid line is the level set of the
likelihood function. If the likelihood function was computed for each point in the constellation,
then one would recover x from y successfully. But this would be hard computationally. With
the Euclidean detection, x̂ that is closer to y than x is would be found instead, which would
cause an erroneous detection. The SIW algorithm can “correct” the error as follows. First, to
estimate the unknown matrix Wx , we compute the matrix Wx̂x which is represented by the red
dashed ellipse around x̂. We can see that the estimate Wx̂x is very close to the correct value Wx ,

August 9, 2017

DRAFT

11

Algorithm 1 Self-interference whitening
Input: y , H , γ, Q θ
x0LMMSE from (27)
Find x̂
x0ML ← NND(yy , H , X )
Find x̂
x0ML , y , H , γ, Q θ ) then
x0LMMSE , y , H , γ, Q θ ) > fˆ(x̂
if fˆ(x̂
x ← x̂
x0LMMSE
x̂
else
x ← x̂
x0ML
x̂
end if
x using (12) and (14)
Generate Wx̂x from x̂
1

Find Wx̂x2 using the Cholesky decomposition
H according to (25)
Generate ỹy and H̃
−1

−1

x0 ← realNND(W
Wx̂x 2 ỹy , Wx̂x 2 H̃
H x̃
x, X̃ )
x̃
x0 ← complex(x̃
x0 )
x̂
x0 , y , H , γ, Q θ ) > fˆ(x̂
x, y , H , γ, Q θ ) then
if fˆ(x̂
x0aML ← x̂
x0
x̂
else
x0aML ← x̂
x
x̂
end if
x(2)
Output: x̂
aML

given by the actual x (blue dashed line). Then, we generate the coordinate system with Wx̂x and
search for the closest constellation point to y in this coordinate system. In this example, x can
be recovered successfully. More importantly, computationally efficient NND algorithms can be
used to perform the search.
Remark 4. The complexity of the SIW algorithm is essentially twice that of the NND algorithm
used, since the other operations including the LMMSE detection have at most cubic complexity
with respect to the dimension of the channel. The complexity of the NND algorithm depends directly on the conditioning of the given matrix. If the columns are close to orthogonal, then channel
inversion is almost optimal. However, in the worse case, when the matrix is ill-conditioned, the
NND algorithm can be slow and its complexity is exponential in the problem dimension. As

August 9, 2017

DRAFT

12

x

y

x̂ = x̂0ML

2

3

Wx̂

x

1

y

x̂ = x̂0ML

Fig. 2: Illustration of the proposed detection in the scalar case. An example with 256-QAM,
PN 2◦ . The dashed lines represent the ellipse defined by the matrix Wx̂ (in red) and Wx (in blue).

mentioned earlier, there exist approximate NND algorithms, e.g., based on lattice reduction, that
can achieve near optimal performance with much lower complexity.

IV. H ARDNESS OF ML D ECODING
In the section we explore how ML decoding may be implemented. While in most cases
of interest the SIW algorithm gives near-optimal performance, and ML decoding is too costly
computationally, it is useful to simply provide lower bounds on the performance of any decoding
algorithm. We comment on the difficulty of implementing ML decoding, in particular when the
dimensions nr , nt are large. For simplicity, we assume that {Θt,k }k=1,...,nt are i.i.d. N (0, σt2 ) and
{Θr,k }k=1,...,nr are i.i.d. N (0, σr2 ).
A. Hardness of computing the likelihood
h
i
HΘ x k2
To compute the likelihood one needs to compute EΘ e−γkyy −H
. For large dimensions,
this seems impossible to do in closed form or using numerical integration. However a natural

August 9, 2017

DRAFT

13

alternative is to use the Monte-Carlo method, using the following estimate:
!
s
X
1
2
H (t) x k
Θ
x, y , H , γ, Q θ ) ≈ F̂s := ln
fˆ(x
,
e−γkyy −H
s t=1

(28)

where Θ (1) , . . . , Θ (s) are s i.i.d. copies of Θ . Using the delta method [14] and the central limit
theorem, the asymptotic variance of this estimator is

HΘ x k2
1 Var e−γkyy −H
Var(F̂s ) ∼ 
2 ,
s
HΘ x k2 ]
E [e−γkyy −H

s → ∞.

(29)

Let us compute this quantity in the (arguably easiest) case where H is the identity matrix with
nt = nr = n:
2

kyy − HΘ x k =

n
X

2

yk − xk ej(Θr,k +Θt,k ) .

(30)

k=1

Define the one-dimensional likelihood
h
i
j(Θr +Θt ) |2
g(x, y, γ, σ) := E e−γ|y−xe
,
and using independence we calculate the moments:
n
h
i Y
y −H
HΘ xk2
−γky
E e
=
g(xk , yk , γ, σ),
Var(e

y −H
HΘ xk2
−γky

)=

k=1
n
Y

g(xk , yk , 2γ, σ) −

k=1

(31)

(32)
n
Y

g(xk , yk , γ, σ)2 .

(33)

k=1

From (29), the asmyptotic variance is hence:
Var(F̂s ) =
where v(x, y, γ, σ) :=

g(x,y,2γ,σ)
.
g(x,y,γ,σ)2
j(Θr +Θt ) |2

random variable e−γ|y−xe

1
s

n
Y

!
v(xk , yk , γ, σ) − 1 .

(34)

k=1

It is noted that, unless σr2 + σt2 = 0 or x = 0, or y = 0, the
cannot be a constant. Thus, we have v(xk , yk , γ, σ) > 1 for

k = 1, . . . , n. As a result, the asymptotic error (34) grows exponentially with the dimension n,
so that the Monte-Carlo method is infeasible in high dimensions, since the number of samples
must also scale exponentially with n to maintain a constant error.
For instance, consider xk = yk = 1 for all k, γ = 40dB, nr = nt = 20, and σr = σt = 3
degrees. Then v ≈ 7, so that, to obtain an error smaller than 0.1, one would require s ≈ 1018
samples, which is clearly not feasible in practice.

August 9, 2017

DRAFT

14

B. Hardness of maximum likelihood search
Assume that one is able to estimate the value of the likelihood with high accuracy (as seen
above this is typically hard), and denote by f¯ this value. We may then consider the following
algorithm. Given the received symbol y , and a radius ρ, one first computes the set of points
x ∈ X nt : kyy − H x k2 ≤ ρ2 }, then one computes the value of f¯ for each of those
S := {x
points and returns the point maximizing f¯. In the large system limit, the following concentration
phenomenon occurs.
Proposition 2. Assume that H has i.i.d. entries with distribution CN (0, 1) and X is chosen
Y − HX k2 . Then we have:
uniformly at random from X nt . Define the radius R2 := kY

2 +σ 2 
σr
 
t
+ γ −1 nr ,
E R2 = 2nt nr 1 − e− 2

(35)

and for any η > 0 we have:
n
 
 o
P (1 − η)E R2 ≤ R2 ≤ (1 + η)E R2 → 1,

(36)

nt , nr → ∞.

Note that the above result is general and does not impose that the number of antennas nr , nt
scale at the same speed. We draw two conclusions from this result: (i) Any such algorithm applied
p
with radius (1 + η) E [R2 ] for any η > 0 is guaranteed to inspect the optimal point with high
probability. (ii) Any such algorithm
which
r has a large success

 probability needs to inspect every

2
√
σ 2 +σ
point in a sphere of radius O
nt nr 1 − exp − r 2 t
, and therefore typically has a very
high complexity.
From the above analysis we see that the second difficulty of the decoding problem for phase
noise channels, even when the likeihood can be computed, lies in the number of points to be
inspected which is exponentially large in nr , nt . The problem is that computing the likelihood at
any given point does not give us any information about the structure of f , and does not help in
maximizing f efficiently.
C. Non-concavity in the high SNR regime
Indeed maximizing f is difficult, and it seems that even performing zero-forcing, i.e., maximizing f over x ∈ Cnt rather than over x ∈ X nt , is difficult since f is non-concave, at least
in the high SNR regime. We first show that, in the high SNR regime, the log-likelihood can be
approximated by a function of the minimal value of kyy − Hθ xk2 , where the minimum is taken
over all possible phases θ .

August 9, 2017

DRAFT

15

Lemma 1. For all x, y , H , Qθ , we have the following high SNR behavior
1
x, y , H , γ, Q θ ) = m(x
x, y , H ) := min kyy − Hθ x k2 .
lim − f (x
θ ∈Rnt +nr
γ

(37)

γ→∞

Proof. Consider  > 0, we have that

1 kyy − HΘ x k2 ≤ m(xx, y , H ) +  e−γ(m(xx,yy ,HH )+) ≤ e−γkyy −HHΘ xk ≤ e−γm(xx,yy ,HH ) .
2



(38)

Taking expectations and then logarithms:
n
o 1
1
2
x, y , H )+)+ ln P kyy −H
HΘx k ≤ m(x
x, y , H )+ ≤ f (x
x, y , H , γ, Q θ ) ≤ −m(x
x, y , H ).
−(m(x
γ
γ
Since the mapping θ 7→ kyy − Hθ x k2 is continuous, for any given  > 0, the probability

x, y , H ) +  is strictly positive. Letting γ → ∞, we have
P kyy − HΘ x k2 ≤ m(x
1
1
x, y , H , γ, Q θ ) ≤ lim sup f (x
x, y , H , γ, Q θ ) ≤ −m(x
x, y , H ).
x, y , H ) + ) ≤ lim inf f (x
−(m(x
γ→∞
γ→∞
γ
γ
x, y , H , γ, Q θ ) = m(x
x, y , H ).
Since the above holds for all  > 0, we have limγ→∞ − γ1 f (x
We can now show that in general the log-likelihood is not concave, hence maximizing it is
not straightforward.
x, y , H , γ, Q θ )
Proposition 3. For γ large enough and H 6= 0 , there exists y such that x 7→ f (x
is a non-concave function.
Proof. Assume that f is concave, then for all x we must have:



1
x + x∗
∗
x, y , H , γ, Q θ ) + f (x
x , y , H , γ, Q θ ) ≤ f
, y , H , γ, Q θ .
f (x
2
2

(39)

From Lemma 1, the above inequality implies that



1
x + x∗
∗
x, y , H ) + m(x
x , y, H ) ≥ m
m(x
, y, H .
2
2

(40)

We shall construct an example to show that the above does not hold in general. Consider any z
x. Then x and z are
such that H z 6= 0 , and let y = H z and x = j(|z1 |, ..., |znt |), so that x ∗ = −x
x, y , H ) = m(zz , y , H ) = 0. Similarly m(x
x∗ , y , H ) = 0.
equal up to a phase transformation, so m(x
∗

x
By definition m( x +x
, y , H ) = m(00, y , H ) = kyy k2 . In this example, (40) would imply 0 ≥ kyy k2 ,
2

which is clearly a contradiction since y 6= 0 . Hence f cannot be concave for γ large enough.
This fact that the log-likelihood is in general non-concave gives another important insight:
SIW can in fact be seen as a (well chosen) concave approximation of a non-concave function.
To circumvent the problem of non-concavity of the log-likelihood SIW approximates it by a

August 9, 2017

DRAFT

16

function which, when Wx is fixed, is concave. As non-concavity appears mainly for high SNRs,
the discrepancy between the performance of ML and SIW (if any) should be more visible in the
high SNR regime, and this will be confirmed by our numerical experiments in the next section.
V. N UMERICAL E XPERIMENTS
In this section, we look at different communication scenarios in which we compare the
proposed scheme to some baseline schemes, including two schemes that ignore the phase noise:
•

The naive LMMSE solution given by (27)

•

The naive ML solution (8)

and a scheme that takes the phase noise into account:
•

Selection between naive LMMSE and ML: the receiver first finds the naive LMMSE and
naive ML solutions, then computes the proposed approximate likelihood function and selects
the one with higher value.

Note that we focus on the vector detection error rate3 as our performance metric.
A. Simulation-based lower bounds
In order to appreciate the performance of the proposed algorithm, we need to compare it
not only with the existing schemes, but also to the fundamental limit given by ML detection,
which is optimal. Let us recall that the proposed SIW algorithm may suffer from two levels of
suboptimality. First, the approximate likelihood function (13) may be inaccurate in some cases.
Second, even if (13) is accurate, the SIW algorithm is not guaranteed to find the optimal solution
(15). Therefore, in order to identify the source of the potential suboptimality, it would be useful
to compare the performance of the SIW scheme with the performance given by (13) and with
that given by (7).
Unfortunately, as pointed out earlier, finding (13) requires an exhaustive search with complexity
growing as |X |nt . Finding (7) is even harder because of the numerical multi-dimensional integration. As such, we resort to lower bounds on the detection error for (13) and (7), respectively,
which are enough for our purpose of benchmarking. To that end, we write


aML
X, Y , H ) < maxn fˆ(x
x, Y , H )
Pe
≥ P fˆ(X
x ∈X

3

t

(41)

The detection is considered successful only when all the symbols in x are recovered correctly. Otherwise an error is declared.

August 9, 2017

DRAFT

17


x, Y , H ) ,
X, Y , H ) < max fˆ(x
≥ P fˆ(X


x ∈L

∀ L ⊆ X nt ,

(42)

where the first lower bound is from the definition of the detection criterion, namely, error occurs
if there exists at least one input vector that has a strictly higher approximate likelihood value.
Note that although the second inequality is valid for all L, it becomes equality if L contains all
the points in X nt that have a higher approximate likelihood value than X does. In this work, we
only take a large set around X to obtain the lower bound (42), without any theoretical guarantee
of tightness of (42). Similarly, for ML detection, we have


ML
X, Y , H ) < maxn f (x
x, Y , H )
Pe ≥ P f (X
x ∈X

(43)

t

X, Y , H ) < f (X
X0 , Y , H )}
≥ P {f (X
X 6= X 0 , f (X
X, Y , H ) < f (X
X0 , Y , H )} ,
= P {X

(44)
∀ X 0 ∈ X nt ,

(45)

where the first lower bound is from the definition of the ML detection criterion, namely, error
occurs if there exists at least one input vector that has a strictly higher likelihood value; the
X, Y , H ) < f (X
X0 , Y , H ); Note that the
equality (45) holds since X 6= X 0 is a consequence of f (X
second lower bound (44) holds for any vector X 0 from the alphabet X nt and with equality when
X 0 is the exact ML solution. Since the ML solution is unknown, we can use any suboptimal
solution instead and still obtain a valid lower bound. Now we can see that (45) is much easier
to evaluate than (43) is since there is no need to perform the maximization over X nt . Intuitively,
if X 0 is a near ML solution, then the lower bound should be tight enough. We shall have some
more discussions on this with the upcoming numerical examples. The lower bound (45) can be
obtained by simulation:
1) For a given observation y and channel H , find a suboptimal solution x 0 .
x, y , H ) and f (x
x0 , y , H ), count an ML error only when f (x
x, y , H ) <
2) If x 0 6= x , compute f (x
x0 , y , H ); otherwise the counter remains unchanged.
f (x
With the proposed method, we need to perform twice the numerical integration (e.g., MonteCarlo integration) only when x 0 6= x . If the latter event happens with small probability, then the
average complexity to evaluate (45) is low. In other words, using a x 0 from a better reference
scheme not only makes the lower bound tighter but also makes it easier to evaluate.

August 9, 2017

DRAFT

18

(a) 64-QAM, PN 3◦ .

(b) 256-QAM, PN 2◦ .

(c) 1024-QAM, PN 1◦ .

Fig. 3: SISO Rayleigh fading with i.i.d. phase noise.

B. Scenario 1: Point-to-point SISO channel
The first scenario focuses on the point-to-point Rayleigh fading single-antenna channels, also
known as single-input single-output (SISO) channels. We consider three different modulation
orders (64, 256, and 1024) with correspondingly three values of phase noise strength (3◦ , 2◦ ,
and 1◦ ) in terms of the standard deviation at both the transmitter and receiver sides. The idea is
to assess the performance of the proposed algorithm in different phase noise limited regimes. In
the SISO case, we compare the proposed scheme with the naive ML scheme which consists in a
simple threshold detection for the real and imaginary parts. Several remarks are in order. First,
we see that ignoring the existance of phase noise incurs a significant performance loss. Second,
if exhaustive search is done with the proposed likelihood approximation, then it achieves the ML
performance. This can be seen from the fact that the proposed simulation-based lower bound
overlaps with the curve with exhaustive search. This confirms the accuracy of the closed-form
approximation (13) at least in the SISO case. It is worth mentioning that the exact likelihood
in this case can also be derived via the Tikhonov distribution as shown in [5], [15] where the
analytic expression involving the Bessel function has been provided. Finally, more remarkably,
the SIW algorithm almost achieves the ML performance without exhaustive search.
C. Scenario 2: Point-to-point LoS-MIMO channel
The second scenario is the point-to-point line-of-sight (LoS) MIMO system commonly deployed as microwave backhaul links [16], [9], [17]. We assume that the channel is constant over
time but each antenna is driven by its own oscillator. This is the worst-case assumption but
also often motivated by the fact that the communication distance is large and thus the distance

August 9, 2017

DRAFT

19

(a) 0.33 Opt. distance, 64-QAM.

(b) 0.7 Opt. distance, 256-QAM.

(c) Opt. distance, 1024-QAM.

Fig. 4: 4 × 4 LoS MIMO. Each antenna has i.i.d. phase noise of standard deviation 1◦ .

between antenna elements is increased accordingly to make sure that the channel matrix is well
conditioned [16], [17]. Here, we adopt the model with two transmit and two receive antennas
each one with dual polarizations. This is effectively a 4×4 MIMO channel. The optimal distance
between the antenna elements at each side can be derived as a function of the communication
distance [16]. However, it may not always be possible to install the antennas with the optimal
spacing due to practical constraints. The condition number of the channel matrix increases when
the antenna spacing decreases away from the optimal distance. In Figure 4, we consider three
configurations with distances 0.33, 0.7, and 1 of the optimal value, generated using the model
from [17]. Accordingly, we use 64-, 256-, and 1024-QAM. For simplicity, we do not consider any
precoding although it may further improve the performance as shown in [17]. We assume that
the phase noises are i.i.d. with standard deviation of 1◦ . We make the following observations.
First, as in the SISO case, phase noise mitigation substantially improves performance. Also,
the proposed likelihood approximation remains accurate as shown by the comparison between
the exhaustive search (15) and the lower bound on ML detection. Further, the proposed SIW
algorithm follows closely the exhaustive search curve and hence achieves near ML performance.
Finally, although in the considered scenario the naive LMMSE is outperformed by the naive ML
scheme, the selection between them can provide a non-negligible gain as shown in Figure 4b.
D. Scenario 3: Uplink SIMO channel with centralized receiver oscillator
The third scenario is the uplink cellular communication channel with four single-antenna
users and one multi-antenna base station receiver. It is assumed that the phase noises at the
users’ side are i.i.d., whereas there is no phase noise at the receiver side. This is a reasonable

August 9, 2017

DRAFT

20

(a) 4 × 4, 64-QAM, PN 4◦ .

(b) 4 × 4, 256-QAM, PN 2◦ .

(c) 4 × 10, 256-QAM, PN 2◦ .

Fig. 5: Uplink SIMO channel. Four single-antenna users with i.i.d. phase noise, multi-antenna
receiver without phase noise.

assumption since the oscillators at the base station are usually of higher quality than those used
by mobile devices. We assume i.i.d. Rayleigh fading in this scenario where three configurations
are considered, as shown in Figure 5. Unlike in the previous scenarios, the naive ML is dominated
by the naive LMMSE at high SNR. This somewhat counter-intuitive observation can be explained
as follows. Without receiver phase noise, the channel can be inverted and we obtain a spatial
parallel channel. Although channel inversion incurs some power loss when the channel is not
orthogonal, each of the resulting parallel subchannels sees an independent phase noise. Therefore,
the demodulation only suffers from a scalar self-interference. On the other hand, with naive ML
the receiver tries to find the closest vector in the image space of H to the received vector y .
Since the linear map H mixes the perturbation of the different transmit phase noises, the naive
ML detection, ignoring the presence of phase noise, suffers from the aggregated perturbation
from all the phase noises. That is why the naive LMMSE can be better than the naive ML
scheme in the high SNR regime where phase noise dominates the additive noise. In the case
when both the transmitter and the receiver have comparable phase noises, such a phenomenon
is rarely observed since the channel inversion also increases the perturbation with the presence
of receiver phase noises.
From Figure 5, we remark that as before the proposed SIW scheme is superior to the other
schemes. With a relatively strong phase noise of 4◦ , the error rate of SIW is 3 to 4 times lower
than that of the naive schemes and is at most twice that of the ML lower bound. With a smaller
phase noise of 2◦ , the SIW scheme can support 256-QAM with four receive antennas, achieving
an error rate 5 times lower than that of the naive schemes. In Figure 5c, we increase the number

August 9, 2017

DRAFT

21

Fig. 6: Uniform phase noise at both transmitter and receiver, Rayleigh fading.

of receive antennas to 10, we see that the gap between the naive schemes is decreased due to the
increased orthogonality of the channel. Nevertheless, the gap between the naive schemes and the
proposed scheme does not decrease since the orthogonality between the users does not reduce
the impact of the phase noise from the transmitter side. We could expect the same observation
even with massive MIMO. Nevertheless, with massive MIMO uplink, the receiver phase noise
can be mitigated substantially due to the asymptotic orthogonality [18].
VI. F URTHER D ISCUSSION AND E XPERIMENTS
A. Robustness to the phase noise distribution
One of the main assumptions of our work is that the phase noise follows a Gaussian distribution. Indeed, our derivation of the closed-form approximation (13) depends on this assumption.
We have seen that this approximation is very accurate in various practical scenarios with Gaussian
phase noise. In practice, however, phase noises may not be Gaussian, which leads to the following
natural question on the robustness: Does the proposed algorithm still work when the phase noise
is not Gaussian? The answer turns out to be positive when we let the phase noise be uniformly
distributed. In Figure 6, we consider three previous configurations (shown in Figure 3c, 5a, and
5b, respectively) but with uniform phase noises. From the results, we see that the proposed
algorithm works well as in the Gaussian phase noise, especially when the phase noise is small.
In fact, we believe that the phase noise distribution with certain regularity (e.g., continuous and

August 9, 2017

DRAFT

22

Fig. 7: Impact of iterations.

bounded density) should not have great impact on the performance of the proposed algorithm
when the phase noise is not too strong.
B. Further improvement with more iterations
Although the proposed scheme achieves near ML performance in many cases, there is still
room for improvement in some situations. As shown in Figure 5 for the 4 × 4 channel with 64QAM, the error floor of the proposed scheme is three times higher than the one of the exhaustive
search. In order to reduce the error floor, we propose to introduce a list of candidate solutions
by iteration. We recall that in the SIW scheme we start from the best solution between naive ZF
and naive ML, and we use this solution as a starting point to estimate the matrix W x followed
by a nearest neighbor detection. The SIW algorithm replaces the starting point with the newly
found one if the latter has a higher approximated likelihood. We can extend the SIW algorithm
with more iterations. Specifically, we can fix a maximum number of iterations. As long as the
newly found point is not inside the list, it should be added to the list and we continue to iterate.
The procedure stops either when we hit the maximum number of iterations or we find a point
already inside the list. At the end, we select the point with the highest likelihood value from the
list.

August 9, 2017

DRAFT

23

C. Potential issues in the very high SNR regime
At very high SNR, the proposed approximation (13) may become less accurate and the perforHΘ x k
mance loss increases with SNR. It can be seen from (6) that the approximation error of kyy −H
has greater impact as γ grows. Numerical experiments show that such performance loss may
be apparent for SNR higher than 40dB depending on the constellation size. A straightforward
workaround is to impose a ceiling value of the SNR γmax for the decoder. In other words, for
SNR higher than the threshold γmax , we let the decoder work as if it were at γmax . Intuitively,
the probability of error cannot be larger than the one at γmax , since we are feeding less noisy
observations to the decoder than what they are supposed to be. With the decoding function, the
observation space can be partitioned into |X |nt regions, each one of which corresponds to a
vector in X nt . The probability of decoding error is the probability that the observation is outside
of the decoding region corresponding to the actual input. With a smaller variance of AWGN, the
observation has a higher probability to be inside the region, hence a lower probability of error.
Nevertheless, the formal proof of this argument is not trivial, and is outside of the scope of the
current paper.4
Another issue at very high SNR is that the likelihood may be too small as compared to the
finite numerical precision. Therefore, it becomes impossible to obtain the simulation-based lower
bound in a reliable way. Furthermore, the number of Monte-Carlo samples required to reach any
given accuracy grows with the SNR.
D. On the practical validity of the adopted discrete-time model
The discrete-time channel model (1) that we adopt in this work is a simplification of the
waveform phase noise channel. Indeed, the discrete-time output sequence is obtained from
filtering followed by sampling in most communication systems. Filtering a waveform corrupted
by phase noise results in not only phase perturbation but also amplitude variation [19]. In the
following, we shall show that the amplitude variation is negligible in the practical regime of
interest. For simplicity, we focus on the single-antenna case with a rectangular filter (i.e., an
integrator). We adhere to the commonly accepted Wiener model for the continuous-time phase
4

We can always add an artificial noise to the observation in order to reduce the SNR to γmax before the detection. In this

way we have an error floor for SNR beyond γmax .

August 9, 2017

DRAFT

24

noise process {Θ(t)}. In particular, Θ(t) ∼ N (0, βt). The equivalent filtered channel gain for
the k th symbol interval of duration T is
Z
Z
1 (k+1)T jΘ(t) d jΘ(kT ) 1 T jΘ(t)
e
dt = e
e
dt
T kT
T 0

(46)

d

from the property of a Wiener process where = means equality in distribution. Assuming that
we can somehow track the past state Θ(kT ) perfectly, we now focus on the following random
variable due to the residual phase noise corresponding to the innovation part:
Z
1 T jΘ(t)
e
dt
B(β, T ) :=
T 0
Z βT
d 1
=
ej Θ̃(t) dt
βT 0
d

= B(1, βT )

(47)
(48)
(49)

d

where in (48) follows from Θ(t) = Θ̃(βt) for some normalized Wiener process {Θ̃(t)} with
Θ̃(t) ∼ N (0, t). We notice that the random variable B(β, T ) depends on the parameters β
and T only through the product βT . The distribution of B(β, T ) has been characterized both
approximately [19] and exactly [20]. In particular, an approximate moment-generating function
has been derived in [19] for small phase noise. We can use the results therein to obtain the
following characterization.
Proposition 4. Let us consider the polar representation B(β, T ) = GejΦ with G ≥ 0 and
Φ ∈ [−π, π). When S := βT is small, we have
S
S2
,
Var (Φ) ≈ , Var (G) ≈
3
180
and thus,
2
1
Var (G) ≈
Var (Φ) .
20

(50)

(51)

In Figure 8, we compare the approximation given by (51) to the correct value obtained
numerically. The approximation (51) is surprisingly accurate even for a standard deviation of
20◦ for the phase. More importantly, such results show that for a small phase perturbation up
to 5 degrees – this is the regime of interest in the present work – the amplitude variation is
less than −55 dB. Therefore, the interference caused by the amplitude variation is dominated
by the AWGN and can be treated as noise without any performance loss. In other words, the
discrete-time model adopted in this work – ignoring the amplitude variation – is indeed valid
for phase noise with moderate variance.

August 9, 2017

DRAFT

25

Fig. 8: Amplitude vs phase perturbation for filtered channel gain.

E. Constellation design
The purpose of the current work is to design a detector that takes into account the phase
perturbation for existing systems in which typical QAM signaling is used. As we observe from the
numerical experiments, in some of the scenarios, even the lower bound on the ML detection error
exhibits an error floor. The error floor is however not a fundamental limit of either the detector or
the channel. With a carefully designed signaling scheme, the probability of ML detection error
can be arbitrarily small with an increasing SNR. For instance, if we use amplitude modulation,
then when the SNR grows, the amplitude ambiguity decreases and the detection error vanishes.
The cost of such an extreme scheme is a reduced spectral efficiency. Algorithms to design a
constellation based on the statistics of both the additive and phase noises have been proposed in
the literature (see [5], [21] and the references therein). To the best of the authors’ knowledge,
only SISO has been considered so far. In the MIMO case, such constellations can still be used
and should provide improvements over QAM. The difficulty with non-QAM constellation lies
in the MIMO detection part, since efficient NND algorithms cannot be applied directly. The
constellation design problem for MIMO phase noise channels is a challenging and interesting
problem in its own right, which is however out of the scope of the current work.
VII. C ONCLUSIONS
In this work, we have studied the ML detection problem for uncoded MIMO phase noise
channels. We have proposed an approximation of the likelihood function that has been shown

August 9, 2017

DRAFT

26

to be accurate in the regimes of practical interest. More importantly, thanks to the geometric
interpretation of the approximate likelihood function, we have designed a simple algorithm that
can solve approximately the optimization problem with only two nearest neighbor detections.
Numerical experiments show that the proposed algorithm can greatly mitigate the impact of
phase noises in different communication scenarios.
R EFERENCES
[1] E. Agrell, T. Eriksson, A. Vardy, and K. Zeger, “Closest point search in lattices,” IEEE Trans. Inf. Theory, vol. 48, no. 8,
pp. 2201–2214, Aug. 2002.
[2] E. Viterbo and J. Boutros, “A universal lattice code decoder for fading channel,” IEEE Trans. Inf. Theory, vol. 45, no. 10,
pp. 1639–1642, Jul. 1999.
[3] B. Hassibi and H. Vikalo, “On the sphere-decoding algorithm I. Expected complexity,” IEEE Trans. Signal Process.,
vol. 53, no. 8, pp. 2806– 2818, Aug. 2005.
[4] J. Jaldén and P. Elia, “DMT optimality of LR-aided linear decoders for a general class of channels, lattice designs, and
system models,” IEEE Trans. Inf. Theory, vol. 56, no. 10, pp. 4765–4780, Oct. 2010.
[5] G. J. Foschini, R. D. Gitlin, and S. B. Weinstein, “On the selection of a two-dimensional signal constellation in the presence
of phase jitter and gaussian noise,” Bell Labs Technical Journal, vol. 52, no. 6, pp. 927–965, 1973.
[6] A. Hajimiri and T. H. Lee, “A general theory of phase noise in electrical oscillators,” IEEE Journal of Solid-State Circuits,
vol. 33, no. 2, pp. 179–194, Feb. 1998.
[7] H. Ghozlan and G. Kramer, “Models and information rates for wiener phase noise channels,” IEEE Trans. Inf. Theory,
vol. 63, no. 4, pp. 2376–2393, Apr. 2017.
[8] A. Demir, A. Mehrotra, and J. Roychowdhury, “Phase noise in oscillators: A unifying theory and numerical methods for
characterization,” IEEE Transactions on Circuits and Systems I: Fundamental Theory and Applications, vol. 47, no. 5, pp.
655–674, May 2000.
[9] G. Durisi, A. Tarable, C. Camarda, R. Devassy, and G. Montorsi, “Capacity bounds for MIMO microwave backhaul links
affected by phase noise,” IEEE Trans. Commun., vol. 62, no. 3, pp. 920–929, Mar. 2014.
[10] S. Yang and S. Shamai (Shitz), “On the multiplexing gain of discrete-time MIMO phase noise channels,” IEEE Trans. Inf.
Theory, vol. 63, no. 4, pp. 2394–2408, Apr. 2017.
[11] R. Krishnan, G. Colavolpe, A. G. i Amat, and T. Eriksson, “Algorithms for joint phase estimation and decoding for MIMO
systems in the presence of phase noise and quasi-static fading channels,” IEEE Trans. Signal Process., vol. 63, no. 13, pp.
3360–3375, July 2015.
[12] T. Datta and S. Yang, “Improving MIMO detection performance in presence of phase noise using norm difference criterion,”
in 53rd Annual Allerton Conference on Communication, Control, and Computing (Allerton), Sep. 2015, pp. 286–292.
[13] H. Mehrpouyan, A. A. Nasir, S. D. Blostein, T. Eriksson, G. K. Karagiannidis, and T. Svensson, “Joint estimation of
channel and oscillator phase noise in MIMO systems,” IEEE Trans. Signal Process., vol. 60, no. 9, pp. 4790–4807, Sep.
2012.
[14] G. W. Oehlert, “A note on the delta method,” The American Statistician, vol. 46, no. 1, pp. 27–29, 1992.
[15] G. Colavolpe, A. Barbieri, and G. Caire, “Algorithms for iterative decoding in the presence of strong phase noise,” IEEE
J. Select. Areas Commun., vol. 23, no. 9, pp. 1748–1757, Sep. 2005.

August 9, 2017

DRAFT

27

[16] F. Bøhagen, P. Orten, and G. E. Øien, “Design of optimal high-rank line-of-sight MIMO channels,” IEEE Trans. Wireless
Commun., vol. 6, no. 4, pp. 1420–1425, Apr. 2007.
[17] P. Ferrand and S. Yang, “Blind precoding in line-of-sight MIMO channels,” in IEEE SPAWC, Edinburgh, UK, 2016.
[18] E. Björnson, J. Hoydis, M. Kountouris, and M. Debbah, “Massive MIMO systems with non-ideal hardware: Energy
efficiency, estimation, and capacity limits,” IEEE Trans. Inf. Theory, vol. 60, no. 11, pp. 7112–7139, Nov. 2014.
[19] G. J. Foschini and G. Vannucci, “Characterizing filtered light waves corrupted by phase noise,” IEEE Trans. Inf. Theory,
vol. 34, no. 6, pp. 1437–1448, Jun. 1988.
[20] Y. Wang, Y. Zhou, D. K. Maslen, and G. S. Chirikjian, “Solving phase-noise Fokker-Planck equations using the motiongroup Fourier transform,” IEEE Trans. Commun., vol. 54, no. 5, pp. 868–877, May 2006.
[21] R. Krishnan, A. G. i Amat, T. Eriksson, and G. Colavolpe, “Constellation optimization in the presence of strong phase
noise,” IEEE Trans. Commun., vol. 61, no. 12, pp. 5056–5066, Dec. 2013.

A PPENDIX
A. Proof of Proposition 2
ΛRHΛ T − H )X
X so that R2 = kV
V + Z k2 . We have:
Define the vector V := (Λ
V k2 + kZ
Z k2 + Z H V + V H Z
R2 = kV

(52)

V k2 + kZ
Zk2 )2 + (Z
ZH V + V H Z )2 + 2(kV
V k2 + kZ
Zk2 )2 (Z
ZH V + V H Z ).
R4 = (kV

(53)

Using the fact that V and Z are independent and Z has i.i.d. CN (0, γ −1 ) entries, we verify that
 




V k2 + E kZ
Z k2 ,
E R2 = E kV
 





 

V k4 + E kZ
Zk4 + 2E kZ
Zk2 E kV
V k2 + 2γ −1 kV
V k2 ,
E R4 = E kV
V k2 ) + Var(kZ
Zk2 ) + 2γ −1 EkV
V k2 .
Var(R2 ) = Var(kV

(54)
(55)
(56)

Zk2 ] = nr γ −1 and Var(kZ
Zk2 ) = 2nr γ −2 , it follows that
Since E [kZ
 


V k2 + γ −1 nr ,
E R2 = E kV

(57)



V k2 ) + 2γ −2 nr + 2γ −1 E kV
V k2 .
Var(R4 ) = Var(kV

(58)

V k2 .
We now calculate the moments of kV
Expectation over H : Let us define Ak,l := Xl (ej(Θt,l +Θr,k ) − 1) and write, conditional on A ,
2

Vk =
kV

nt
nr X
X
k=1

August 9, 2017

l=1

2

Hk,l Ak,l

d

=

nr
X

Ek ,

(59)

k=1

DRAFT

28

Ak k2 ) with Ak := [Ak,l ]l=1,...,nt , k = 1, . . . , nr ; we used the fact that H has
where Ek ∼ Exp(1/kA
V k2 for a given
i.i.d. CN (0, 1) entries. Then we can calculate the first and second moments of kV
A = [Aj,k ]j,k :
2

Vk ] =
EH |AA [kV

nr
X

EH |AA [Ek ] =

k=1

nt
nr X
X

|Ak,l |2 ,

"


V k4 = EH |AA
EH |AA kV

nr
X



#


E2k + EH |AA 

k=1

=

nr
X

(60)

k=1 l=1



A k k4 + 
2kA

k=1

nt
nt X
X
k0 =1
k0 6=k

nt X
nt
X
k0 =1
k0 6=k



Ek Ek0 

(61)

k=1



Ak k2 kA
A k 0 k2 
kA

(62)

k=1

nr X
nr
X
Ak k2 kA
Ak0 k2 .
(1 + 1{k = k 0 })kA
=

(63)

k0 =1 k=1



Moments of A : We recall that |Ak,l |2 = 2|Xl |2 (1 − cos(Θr,k + Θt,l )). We have that E |Xl |2 =


1 and define E |Xl |4 = P̄ 2 , for l = 1, . . . , nt . Using the independence and the identity

for zero-mean Gaussian Θ, we obtain
E [cos(Θ)] = exp − Var(Θ)
2


σ 2 +σ 2




2
2
− r2 t
E |Ak,l | = 2E |Xl | E [1 − cos(Θt,l + Θr,k )] = 2 1 − e
.
(64)
We now calculate the correlation between the entries of A
E(|Ak,l |2 |Ak0 ,l0 |2 ) = 4E(|Xl |2 |Xl0 |2 ) ρk,k,0 l,l0

(65)

with
ρk,k,0 l,l0 := E((1 − cos(Θt,l + Θr,k ))(1 − cos(Θt,l0 + Θr,k0 )))
σ 2 +σ 2
− r2 t

= 1 − 2e

+ e−σr −σt cosh(σr2 (1{k = k 0 } + σt2 1{l = l0 })),
2

2

(66)
(67)

where to obtain the last equality we use the trigonometric identities and again apply the identity

x
−x
E [cos(Θ)] = exp − Var(Θ)
; we recall that cosh(x) = e +e
.
2
2
V k2 : From (60) and (64), we have the first moment
Moments of kV
nt
nr X


σ 2 +σ 2

 X


2
2
− r2 t
Vk =
E kV
E |Ak,l | = nt nr 2 1 − e
.

(68)

k=1 l=1

August 9, 2017

DRAFT

29

For the variance, we apply (63) and (68)
2




V k2
V k4 − E kV
V k2 = E kV
(69)
Var kV
nt 
nt X
nr X
nr X
X



 

2
2
0
2
2
0
0
0
0
E |Ak,l | |Ak ,l | (1 + 1{k = k }) − E |Ak,l | E |Ak ,l | .
=
k0 =1 k=1 l0 =1 l=1

(70)
Noting that E [|Ak,l |2 |Ak0 ,l0 |2 ] = E [|Ak,l |2 ] E [|Ak0 ,l0 |2 ] if k 6= k 0 and l 6= l0 , we obtain the variance


V k2 ) = 4nt nr w1 + w2 (nt − 1) + w3 (nr − 1) ,
Var(kV
(71)
where w1 , w2 , w3 > 0 do not depend on nt , nr and correspond to the cases (k = k 0 , l = l0 ),
(k = k 0 , l 6= l0 ), (k 6= k 0 , l = l0 ), respectively, in the summation (70),
2

2

2

w1 := 2P̄ 2 (1 − 2e−σ + e−2σ cosh(2σ 2 )) − (1 − e−σ )2 ,
2

2

2

w2 := 2(1 − 2e−σ + e−2σ cosh(σr2 )) − (1 − e−σ )2 ,
2

2

2

w3 := P̄ 2 (1 − 2e−σ + e−2σ cosh(σt2 )) − (1 − e−σ )2 ,
where we define σ 2 :=

(72)
(73)
(74)

σt2 +σr2
.
2

Putting it together: From (54) and (68), we have

2 +σ 2 
σr
 


t
V k2 + γ −1 nr = nt nr 2 1 − e− 2
+ γ −1 nr ,
E R2 = E kV

(75)

which yields the first result. Note that E [R2 ] ≥ c1 nt nr for some constant c1 > 0 with respect to
(nt , nr ). From (56) and (71), we have proven that there exists a constant c2 > 0 such that


V k2 ) + 2γ −2 nr + 2γ −1 E kV
V k2 ≤ c2 nr nt (nt + nr ).
Var(R2 ) = Var(kV

(76)

Hence
Var(R2 )
c2
≤
(E [R2 ])2
c21



1
1
+
nr nt


→ 0,

nt , nr → ∞,

(77)

and applying Chebychev’s inequality yields the second result:

 
 
P (1 − η)E R2 ≤ R2 ≤ (1 + η)E R2 → 1, nt , nr → ∞.

August 9, 2017

(78)

DRAFT

30

B. Proof of Proposition 4
From [19, eq.6]5 , we can derive the moment-generating function (MGF) of (G, Φ)
1

MG,Φ (ξ, η) = eξ sinhc− 2

p 
Sξ
"
##
"
r !  −1
p 
Sξ
Sξ
Sη 2
cothc
−
+ tanhc
Sξ
, (79)
· exp
8
4
4

where we define sinhc(x) :=

sinh(x)
,
x

coth(x)
,
x

cothc(x) :=

tanhc(x) :=

tanh(x)
,
x

and recall that

S := βT . It follows that the MGF of G and Φ are
p 
MG (ξ) = MG,Φ (ξ, 0) = e sinhc
Sξ ,


1 2
Sη ,
MΦ (η) = MG,Φ (0, η) = exp
6
ξ

− 12

(80)
(81)

where we used the fact that, when x → 0, sinhc(x) = 1 + O(x2 ), tanhc(x) = 1 + O(x2 ), and
cothc(x) −

1
1
= + O(x2 ).
2
x
3

(82)

After finding the first and second derivatives of both MGF (80) and (81) with some elementary
manipulations, we obtain the desired variances

S2
2
Var (G) = E G2 − (E (G))2 = MG00 (0) − (MG0 (0)) =
,
180

S
2
Var (Φ) = E Φ2 − (E (Φ))2 = MΦ00 (0) − (MΦ0 (0)) = .
3

(83)
(84)

Note that we use approximate equality in (50) and (51) since the MGF derived in [19] is indeed
approximative with the assumption of small S.

5

Note that the random variable B in (47) differs with the one in [19, eq.1] in a normalization factor T . The MGF has been

scaled accordingly.

August 9, 2017

DRAFT

