# Phase Noise Detection via Expectation Propagation and Related Algorithms
††thanks: This research was granted by University of Parma through the action
Bando di Ateneo 2022 per la ricerca co-funded by MUR-Italian Ministry of
Universities and Research - D.M. 737/2021 - PNR - PNRR - NextGenerationEU.

Elisa Conti, Armando Vannucci, Amina Piemontese, and Giulio Colavolpe

###### Abstract

In the context of signal detection in the presence of an unknown time-varying
channel parameter, receivers based on the Expectation Propagation (EP)
framework appear to be very promising. EP is a message-passing algorithm based
on factor graphs with an inherent ability to combine prior knowledge of system
variables with channel observations. This suggests that an effective
estimation of random channel parameters can be achieved even with a very
limited number of pilot symbols, thus increasing the payload efficiency.
However, achieving satisfactory performance often requires ad-hoc adjustments
in the way the probability distributions of latent variables - both data and
channel parameters - are combined and projected. Here, we apply EP to a
classical problem of coded transmission on a strong Wiener phase noise
channel, employing soft-input soft-output decoding. We identify its
limitations and propose new strategies which reach the performance benchmark
while maintaining low complexity, with a primary focus on challenging
scenarios where the state-of-the-art algorithms fail.

###### Index Terms:

expectation propagation (EP), phase noise, iterative detection/decoding,
factor graphs (FGs), sum-product algorithm (SPA), low-density parity-check
(LDPC) codes.

##  I Introduction

Detection algorithms for channels affected by a time-varying phase noise have
received a lot of attention in the literature considering either linear or
continuous-phase modulations and different scenarios (see, e.g., [1, 2, 3, 4,
5, 6, 7] and references therein). This is because in many communication links
phase noise must be considered one of the major impairments. Examples are
represented by scenarios employing a high carrier frequency, such as
(coherent) optical communications, communications from geostationary
satellites, etc.

One of the most effective algorithms is that proposed in [1]. It is designed
using the framework based on factor graphs (FGs) and the sum-product algorithm
(SPA). In a scenario like the one at hand, where continuous random variables
(the channel phase) appear in the FG, a common approach to implement the SPA
is to resort to the use of canonical distributions [8]. In particular, in [1],
the messages representing the a-posteriori probability density functions
(pdfs) of the channel phase are represented through Tikhonov distributions,
which can be described by a single complex parameter. Although suboptimal, the
resulting algorithm exhibits an excellent trade-off between performance and
complexity. The suboptimality is related to the joint presence of discrete
(the code symbols) and continuous random variables in the graph, which brings
up mixture pdfs with exponential proliferation, approximated in [1] with
unimodal distributions. The presence of distributed pilot symbols is thus
required to make the algorithm bootstrap. In fact, the algorithm in [1] based
on the Tikhonov parametrization is very sensitive to the positioning of pilot
symbols, so that increasing the pilot spacing causes a performance
degradation. In particular, when pilot symbols are concentrated at both ends
of the codeword, as a preamble/postamble, or, more in general, when a maximum
distance between consecutive pilot blocks is exceeded, it is known to fail [1,
9]. This occurs since this kind of algorithm cannot achieve any effective
phase detection without the aid of extrinsic information, provided by the
decoding part of the FG during turbo iterations, or without distributed pilots
(with sufficiently close spacing). Therefore, other alternative solutions have
to be analyzed to this purpose.

As an extension of [1], the results of Raphaeli and coworkers [10, 6] are
obtained by letting mixture messages propagate one step further into the FG;
their exponential proliferation being limited by an appropriate pruning of the
mixture components, performed at the level of the Markov chain governing phase
noise. This way, the phase uncertainty inherent in the observation of channel
output is free to interact with the provisional estimation of previous phase
samples (or with the following ones, in backward block processing). Such a
mixture message reduction approach, however, is characterized by a
considerable complexity.

The same objective of letting channel observations interact with provisional
estimates of adjacent phase samples is achieved by the expectation propagation
(EP) algorithm [11, 12, 13, 14], where similar mixture reduction techniques
are applied to the marginal distributions of the variables of interest, rather
than to individual messages. Despite the projection of messages or marginals
do have some features in common, and can even reduce to the same algorithm in
some cases [15, 16, 17], there is a profound conceptual difference between
them. In fact, only in the latter case, a message coming from a channel
observation is merged, i.e., multiplied, with the (provisional) prior belief
on the destination variable that comes from the rest of the FG. This is the
key to the potential success of EP in many similar applications, including,
e.g., transmission over fading channels [18, 19].

The EP framework has already been successfully applied to phase noise
channels, with either distributed or concentrated pilots [9, 20], resorting to
joint detection and decoding. This is not, however, a practical solution,
since a separate (and sequential) detection of channel parameters and decoding
is desirable in the design of digital receivers. On the contrary, the
challenge of making EP work with separate phase detection and decoding has not
been solved. In fact, previous analyses have shown that, despite its ability
to bootstrap the phase estimation process from a proper block of
preamble/postamble pilots, the native EP algorithm is not able to effectively
refine the phase estimates beyond the first phase-detecting iteration [21];
and this is true even when pilot blocks are close enough, a situation where
the algorithm in [1] shows a better performance.

In this work, we first examine the EP framework applied to phase noise
channels, identify its weaknesses and the underlying reasons. Preliminary
results of this analysis can be found in [21]. We then propose novel
strategies to improve its performance, resulting in a flexible algorithm
capable of adapting to various pilot positioning and iterations scheme
(whether separate or joint detection and decoding) always outperforming the
state-of-the-art solutions.

The rest of this paper is organized as follows. In Section II, we introduce
the system model and the FG representation of the joint posterior pdf of the
information symbols and channel parameters. Section III is devoted to the
description of the EP algorithm while, in Section IV, we present the proposed
strategies to overcome the EP limitations without increasing its complexity.
The message scheduling is outlined in Section V. A complexity comparison among
the analyzed algorithms is carried out in Section VI. Finally, in Section VII
we show and discuss the numerical results and, in Section VIII, we draw some
concluding remarks.

Notation. The following notational convention will be used throughout the
paper. Bold letters denote vectors, the superscript $(\cdot)^{*}$ applied to a
scalar term denotes its complex conjugate. We indicate with $\Re[\cdot]$ and
$\Im[\cdot]$ the real and imaginary parts of a complex quantity, respectively.
A complex circularly symmetric (respectively, real) Gaussian random variable
${x}$ with mean $\eta_{x}$ and variance $\sigma_{x}^{2}$ is denoted by
$x\sim{\cal N}_{\mathbb{C}}(\eta_{x},\sigma^{2}_{x})$ [respectively, by
$x\sim{\cal N}(\eta_{x},\sigma^{2}_{x})$]. We denote the complex
(respectively, real) Gaussian pdf with argument $x$, mean $\eta_{x}$, variance
$\sigma_{x}^{2}$ by $g_{\mathbb{C}}({x}-\eta_{x};\sigma_{x}^{2})$
[respectively, by $g({x}-\eta_{x};\sigma_{x}^{2})$]. We denote by $p(\cdot)$
the pdf of continuous variables (vectors) and with a capital $P(\cdot)$ the
probability mass function (pmf) of discrete variables as well as the pdf of
mixed random vectors. We use $E_{p}[\cdot]$ to indicate the expectation under
the distribution $p(\cdot)$. When $p(\cdot)$ is not specified, it is implicit
from the context.

##  II System Model and Related Factor Graph

Figure 1: Factor graph representing the joint probability distribution of
phase noise samples and coded symbols.

In the system we wish to investigate, a sequence of $K$ complex coded symbols
$\boldsymbol{c}=[c_{0},c_{1},\cdots,c_{K-1}]$, modulated over an $M$-ary
constellation (with alphabet $\mathcal{A}$), is transmitted over an additive
white Gaussian noise (AWGN) channel also affected by Wiener phase noise, so
that the received samples are

$r_{k}=c_{k}{\rm e}^{j\theta_{k}}+n_{k}\,,$ (1)

where $\\{n_{k}\\}$ is a sequence of independent (complex and circularly
symmetric) Gaussian noise samples, i.e., $n_{k}\sim{\cal
N}_{\mathbb{C}}(0,2\sigma^{2})$, with zero mean and given variance per
component $\sigma^{2}$, while the phase noise sequence
$\boldsymbol{\theta}=[\theta_{0},\theta_{1},\cdots,\theta_{K-1}]$ follows the
Wiener model, where each sample

$\theta_{k}=\theta_{k-1}+\Delta_{k}$ (2)

results from the previous one plus a zero-mean (real) Gaussian increment
$\Delta_{k}\sim{\cal N}(0,\sigma_{\Delta}^{2})$ whose variance
$\sigma_{\Delta}^{2}$ dictates the severity of phase noise.

Since the system model is the same as the one discussed in [1], the FG
representing the joint distribution of the latent variables
$\boldsymbol{\theta}$ and $\boldsymbol{c}$, conditioned on the value of the
observed variables $\boldsymbol{r}=[r_{0},r_{1},\cdots,r_{K-1}]$, is the same
as the one shown therein and is reported in Fig. 1. The subscripts of messages
identify their direction in the FG ($u$p, $d$own, $f$orward, $b$ackward),
whereas the subscript of variable and factor nodes coincides with a time
index.

Before observing the channel output, the phase noise sequence is independent
of that of coded symbols, so that the joint posterior distribution of all
latent system variables is

$\begin{split}P(\boldsymbol{c,\theta\\!\\!\mid\\!r})&\\!\propto
P(\boldsymbol{c})p(\boldsymbol{\theta})p(\boldsymbol{r\mid c,\theta})\\\
&\propto\chi(\boldsymbol{c})p(\theta_{0})\prod_{k=1}^{K-1}p(\theta_{k}\mid\theta_{k-1})\prod_{k=0}^{K-1}f_{k}(c_{k},\theta_{k})\end{split}$
(3)

where $\propto$ is the proportionality symbol, $\chi(\boldsymbol{c})$ is an
indicator function implementing the code constraints and we assume that all
the allowed codewords in the employed codebook are equally likely.

The factor node $f_{k}$ implements the observation of the received samples in
(1), hence

$\begin{split}f_{k}(c_{k},\theta_{k})&=p(r_{k}\mid
c_{k},\theta_{k})=g_{\mathbb{C}}(r_{k}-c_{k}{\rm
e}^{j\theta_{k}};\sigma^{2})\\\
&\propto\exp\left(-\frac{1}{2\sigma^{2}}\left|r_{k}-c_{k}{\rm
e}^{j\theta_{k}}\right|^{2}\right)\end{split}$ (4)

is a complex Gaussian distribution.

Besides $f_{k}$, the only other factor nodes in the FG of Fig. 1 (if we
disregard the check nodes included in the code sub-graph) are the ones related
to the conditional pdfs in (3) that implement the Markov chain governing the
Wiener phase noise,

$p(\theta_{k}\mid\theta_{k-1})\triangleq
p_{\Delta}(\theta_{k}-\theta_{k-1})=g(\theta_{k}-\theta_{k-1};\sigma_{\Delta}^{2})\,,$
(5)

which are Gaussian too, as per (2).

Owing to the fact that the additive sources of randomness in (1) and (2) are
Gaussian, the FG in Fig. 1 seems to represent a Gaussian belief network, where
all messages are Gaussian [22] and the message passing procedure resembles the
simple operations of a Kalman smoother. This is clearly not true, for two
reasons. First, the problem at hand entails a “mixed model”, where continuous
($\boldsymbol{\theta}$) and discrete ($\boldsymbol{c}$) variables coexist,
which produces probabilistic mixtures in the FG, specifically arising from the
messages $p_{d}(\theta_{k})$ which are linear combinations of simpler
component distributions. In addition, (4) is a Gaussian pdf only if seen as a
function of $r_{k}$ whereas the factor nodes $f_{k}$ send downward messages
$p_{d}(\theta_{k})$ to the variable nodes $\theta_{k}$, hence (4) must be seen
as a function of $\theta_{k}$ and is thus proportional to a Tikhonov
distribution [1],

$\begin{split}f_{k}(c_{k},\theta_{k})&\propto\exp\left(-\frac{|c_{k}|^{2}}{2\sigma^{2}}+\frac{1}{\sigma^{2}}\Re\left[r_{k}c_{k}^{*}{\rm
e}^{-j\theta_{k}}\right]\right)\\\
&=\exp\left(-\frac{|c_{k}|^{2}}{2\sigma^{2}}\right)2\pi
I_{0}(|z_{k}|)t(\theta_{k};z_{k})\end{split}$ (6)

where $z_{k}=r_{k}c_{k}^{*}/\sigma^{2}$ is the complex parameter
characterizing the general Tikhonov distribution:

$t(\theta_{k};z_{k})=\frac{1}{2\pi I_{0}(|z_{k}|)}\exp\left(\Re\left[z_{k}{\rm
e}^{-j\theta_{k}}\right]\right)\,.$ (7)

In directional statistics [23], the phase $\angle z_{k}$ corresponds to the
circular mean of the circular random variable $\theta_{k}$ in (7), while the
magnitude $|z_{k}|$ is a measure of its precision, i.e., of the inverse of the
variance
$\sigma^{2}_{\theta_{k}}=1-\frac{I_{1}(|z_{k}|)}{I_{0}(|z_{k}|)}\in[0;1]$,
where $I_{p}(x)$ is the modified Bessel function of the first kind of order
$p$.111More in general, the $p$-th trigonometric moment of the Tikhonov
variable $\theta_{k}$ is
$E_{t}[\exp{(jp\theta_{k})}]=\exp{(jp\angle{z_{k}})}\frac{I_{p}(|z_{k}|)}{I_{0}(|z_{k}|)}$,
so that the phase of the first trigonometric moment is the circular mean.

The Bayesian inference problem for the estimation of symbols and channel
parameters cannot be solved by the plain SPA since the message

$p_{d}(\theta_{k})=\sum_{m=0}^{M-1}P_{d}(c_{k}^{m})f_{k}(c_{k}^{m},\theta_{k})\propto\sum_{m=0}^{M-1}\alpha_{k}^{m}t(\theta_{k};z_{k}^{m})$
(8)

is a Tikhonov mixture. Its coefficients $\alpha_{k}^{m}$, defined as

$\alpha_{k}^{m}=P_{d}(c_{k}^{m})\exp\left(-\frac{|c_{k}^{m}|^{2}}{2\sigma^{2}}\right)2\pi
I_{0}(|z_{k}^{m}|)\,,$ (9)

depend on the extrinsic information on each symbol (i.e., on the values of its
pmf for each of the possible $M$ symbols $c_{k}^{m}$), as provided by the
soft-input soft-output decoding part of the FG (top of Fig. 1), as well as on
the magnitude of the symbols. In the SPA, the mixture messages (8) would
propagate through the bottom half of the FG and eventually proliferate to
produce untractable mixtures with an exponentially increasing number of
components. Approximate inference is thus demanded, which can be performed in
a variety of ways.

##  III Transparent Propagation Algorithms and Expectation Propagation

The most straightforward approach to approximate message passing is to
discretize the distribution of $\theta_{k}$, with a given number of samples
$N_{\theta}$, so that all the latent system variables appear to be discrete
and modelled by their pmf. This produces an algorithm called discretized-phase
BCJR (dp-BCJR) in [4], since the forward-backward message passing procedure,
along the Markov chain in the bottom half of the FG, is the same as that of
the celebrated BCJR algorithm. Its computational complexity is very high, thus
we shall implement this algorithm, denoted as dp-BCJR in the simulation
results, as a practical benchmark, assuming that its performance gets close to
that of an ideal SPA when $N_{\theta}$ is sufficiently large.

A totally different approach relies on projecting messages and/or marginals
onto a selected family of approximating distributions, usually in the
exponential form (so that the family is closed under the multiplication
operation),

$q(x)=\exp\left(\sum_{i}\eta_{i}g_{i}(x)\right)\,,$ (10)

where $\eta_{i}$ are called natural parameters and the functions $g_{i}(x)$
are the features of the family. This way, only the natural parameters
$\eta_{i}$ are updated, for each distribution, resulting in a parametric
message passing procedure. In [1], the Tikhonov approximating family is
selected quite naturally, since it is known to be the marginal distribution of
$\theta_{k}$, conditioned on the values taken by the corresponding symbol and
by the corresponding channel output, i.e., $p(\theta_{k}\mid
c_{k},r_{k})\propto f_{k}(c_{k},\theta_{k})$, as per (4). Curiously, the way
in which message $p_{d}(\theta_{k})$ has been projected onto a Tikhonov pdf in
[1] relies on the Gaussian expression (4) of the factor node, so that (8) is
approximated by the Gaussian pdf with minimum Kullback-Leibler (KL) divergence
from the mixture, which is further interpreted as a Tikhonov message towards
$\theta_{k}$.

A more natural solution would have been to project the Tikhonov mixture (8)
onto a Tikhonov pdf with minimum KL divergence, which is achieved (for this as
well as for any other exponential approximating family) by matching the
expectations of the features [12]. We refer to this approach, along with the
Tikhonov parameterization of [1], as Transparent Propagation (TP) algorithms
(the reason will be clear later). The features of a Tikhonov pdf like (7) are
$\cos(\theta_{k})$ and $\sin(\theta_{k})$, associated with the natural
parameters $\Re[z_{k}]$ and $\Im[z_{k}]$ respectively, so that their
expectations can be jointly computed as the $\Re/\Im$ parts of the first
trigonometric moment $E[\exp(j\theta_{k})]$. Exploiting known results [23] and
the linearity of expectation, the mixture in (8) is approximated by the
Tikhonov pdf $p_{d}^{TP}(\theta_{k})=t(\theta_{k};z_{k}^{TP})$, which achieves
the minimum KL divergence when $z_{k}^{TP}$ obeys the following moment
matching equation

$\frac{I_{1}(|z_{k}^{TP}|)}{I_{0}(|z_{k}^{TP}|)}{\rm e}^{j\angle
z_{k}^{TP}}=\sum_{m=0}^{M-1}\overline{\alpha}_{k}^{m}\frac{I_{1}(|z_{k}^{m}|)}{I_{0}(|z_{k}^{m}|)}{\rm
e}^{j\angle z_{k}^{m}}$ (11)

where $\overline{\alpha}_{k}^{m}$ stands for the normalized version of the
coefficients in (8), i.e.,

$\sum_{m=0}^{M-1}\overline{\alpha}_{k}^{m}=1\,.$ (12)

The complex equation (11) yields the circular mean $\angle z_{k}^{TP}$ as well
as the variance of phase noise, related to the Bessel ratio of the magnitude
$|z_{k}^{TP}|$ [10]. Despite the different analytical approach, the two ways
of projecting the mixture message $p_{d}(\theta_{k})$ described above, i.e.,
that of [1] and the one in (11), yield extremely similar results, as verified
by numerical simulation in all phase noise scenarios that we analyzed. For
this reason, their equivalent performance is collectively shown under the TP
label in the simulation results in Section VII.

In the EP algorithmic framework [11], it is the entire marginal of each latent
variable that is approximated/projected; in our case,

$p(\theta_{k})=p_{d}(\theta_{k})p_{f}(\theta_{k})p_{b}(\theta_{k})=p_{d}(\theta_{k})p_{u}(\theta_{k})$
(13)

of which the message $p_{d}(\theta_{k})$ only represents one factor. In (13),
message $p_{u}(\theta_{k})\triangleq p_{f}(\theta_{k})p_{b}(\theta_{k})$ is
defined as the product of forward and backward messages, and is sent upwards,
opposite to message $p_{d}(\theta_{k})$, although not reported in the FG in
Fig. 1. Indeed, the conceptual meaning of the two factors of the marginal in
(13) is the following: $p_{d}(\theta_{k})$ is an observation message, that
carries information on $\theta_{k}$ after the channel observation entailed in
the factor node (6), while $p_{u}(\theta_{k})$ plays the role of a temporary
prior belief for the phase noise variable $\theta_{k}$, as provided by the
rest of the FG. Denoting the projection operation as222We denote with
$\text{KL}[p(\cdot)\,||\,q(\cdot)]$ the KL divergence of an approximated
distribution $q(\cdot)$ from the exact distribution $p(\cdot)$.

$\textrm{proj}[p(\cdot)]=\arg\min_{q(\cdot)\in{\cal
F}}\text{KL}\left[p(\cdot)\parallel q(\cdot)\right]\,,$ (14)

where $\cal F$ is the approximating family, the approximating message

$p_{d}^{EP}(\theta_{k})=\frac{\textrm{proj}[p_{d}(\theta_{k})p_{u}(\theta_{k})]}{p_{u}(\theta_{k})}$
(15)

is computed and sent in place of the mixture $p_{d}(\theta_{k})$. After the
observation $f_{k}(c_{k},\theta_{k})$, the estimated marginal for the phase
noise variable $\theta_{k}$ is thus updated in a different way, compared to
the plain SPA rule (13), and expressed as

$p^{EP}(\theta_{k})=p_{d}^{EP}(\theta_{k})p_{u}(\theta_{k})\,,$ (16)

which is equal to the numerator of (15) and thus belongs, by construction, to
the approximating family (Tikhonov, in our case).

The projection operation in (15) still amounts to a moment matching operation
like (11), despite a marginal is projected, in EP, rather than a message, as
in TP. Of course, in the case of EP, the temporary prior
$p_{f}(\theta_{k})p_{b}(\theta_{k})$ must be accounted for in the right hand
side of (11), by simply adding the (complex) parameters $z_{k,f}$ and
$z_{k,b}$ of the forward and backward messages to those of each mixture
component, i.e., to $z_{k}^{m}$ (hence the coefficients $\alpha_{k}^{m}$ must
be calculated accordingly). Therefore, from a conceptual standpoint, the
mixture to be projected is now the result of the interaction between the prior
belief, deriving from the Markov chain at the bottom half of the FG in Fig. 1,
and the current channel observation, represented by message
$p_{d}(\theta_{k})$. The expression for the marginal mixture parameters thus
becomes

$z_{k,MIX}^{m}=z_{k,u}+z^{m}_{k}$ (17)

where we denoted with $z_{k,u}=z_{k,f}+z_{k,b}$ the temporary prior parameter.

The general approach to approximate message passing can be either a projection
of the marginal pdf of each variable node (onto the selected approximating
family), as prescribed by EP, or a projection of individual mixture messages
onto the same family, as described above for TP algorithms. If, for the sake
of argument, the EP algorithm did not produce any information about the
temporary prior, then we would assume $p_{f}(\theta_{k})p_{b}(\theta_{k})$ to
be uniform in (15) and hence it would transparently shift out of the
projection operation, to be simplified with the denominator, so that (15)
would reduce to the simple projection of message $p_{d}(\theta_{k})$. This is
the reason for which we denote the message projection approach as TP [16]. It
can be demonstrated that in some specific problems, the TP and EP approaches
lead to the same solution [17], although in general they differ remarkably.

No matter if one of the TP algorithms or if EP is adopted, the forward and
backward messages $p_{f}(\theta_{k})$ and $p_{b}(\theta_{k})$ are assumed to
belong to the Tikhonov family too, despite the factor nodes in the Markov
chain of the FG, e.g., $p(\theta_{k+1}\mid\theta_{k})$, introduce a
convolution between the Tikhonov message
$p_{d}^{TP/EP}(\theta_{k})p_{f}(\theta_{k})$ and a Gaussian like (5), for the
computation of $p_{f}(\theta_{k+1})$. As it is shown in the Appendix of [1],
if the phase noise variance $\sigma_{\Delta}^{2}$ is not exaggeratedly large
then $p_{f}(\theta_{k+1})$ is very well approximated by a Tikhonov pdf whose
parameter $|z_{f,k+1}|$ is less than that of the incoming message,
$|z_{f,k}+z_{k}^{TP/EP}|$, while the circular mean is the same (see (36)-(38)
in [1], for details). This corresponds conceptually to the fact that the
Wiener phase noise update in (2) does not bias the temporary estimate for
$\theta_{k+1}$ but decreases its precision, according to its Gaussian
variance. A similar conclusion holds for the backward messages too.

##  IV The EP Modification

The EP algorithm has been already applied to transmission over Wiener phase
noise channels in [9] (and later in [20]), in the absence of distributed pilot
symbols, achieving a good performance when turbo iterations333We call turbo
iterations the exchange of information between the detector and decoder parts
of the FG, implemented by the vertical messages $p_{d}(\cdot)$ and
$P_{u}(\cdot)$ in Fig. 1. are employed. However, existing literature does not
demonstrate successful operation of EP with separate phase detection and
decoding.

As it is known, the performance of TP algorithms, such as [1], significantly
deteriorates by increasing the pilot spacing, leading to complete failure when
they are concentrated at both ends of the codeword. In contrast, in these
scenarios, EP represents a viable alternative, due to its ability to self-
sustain the process of phase estimation across a long block of payload coded
symbols.

However, the native EP algorithm is not able to effectively refine the phase
estimates beyond the first phase-detecting iteration [21], resulting in
unsatisfying performance both with distributed pilots (where it is
outperformed by TP algorithms) and with more concentrated pilots (where,
unable to benefit from inner detector iterations, its performance is limited
as shown in Section VII).

Through a detailed numerical analysis of EP, we found that its failure is
confined to some critical data packets that, due to a long sequence of noisy
observations, bring the sequence of precision values $|z^{EP}_{k,f/b}|$ to a
collapse, while propagating forward or backward messages, $p_{f}(\theta_{k})$
and $p_{b}(\theta_{k})$. Once a critical lower threshold is exceeded, the
precision approaches zero and is rarely able to recover, even if the channel
in (1) outputs reliable observations, i.e., samples $r_{k}$ with little noise.
In order to make EP overcome these limiting situations, we introduce the
modifications described hereafter, which improve its robustness even in very
challenging conditions, as further verified in Section VII.

###  IV-A Bessel Ratio Approximation

A first countermeasure that we introduce in the classical EP algorithm
concerns the way in which the Bessel Ratio (BR) $I_{1}(x)/I_{0}(x)$ is
inverted. Since $|z^{EP}_{k,f/b}|$ are calculated after the EP message (15),
which is in turn computed by applying the moment matching (11), we verified
that the occasional collapse of their value is largely due to a numerical
instability entailed in the approximation of the BR in (11) whose inverse, not
available in closed form, is necessary to obtain $|z^{EP}_{k}|$. In fact,
despite the approximation

$I_{1}(|z|)/I_{0}(|z|)\simeq\exp(-0.5/|z|)\,,$ (18)

which is commonly adopted in the literature [9, 20], is very close to the
exact function when $|z|\gg 0$, it is not equally accurate for lower values of
$|z|$, e.g, in the range $[0,5]$ as shown in Fig. 2. When this traditional
approximation (18) is adopted, the moment matching equation (11) becomes

${\rm e}^{-0.5/|z^{EP}_{k}|+j\angle
z_{k}^{EP}}=\sum_{m=0}^{M-1}\overline{\alpha}_{k,MIX}^{m}{\rm
e}^{-0.5/|z_{k,MIX}^{m}|+j\angle z_{k,MIX}^{m}},$ (19)

where, differently from [9], we adopted BR approximation also in the right-
hand side. In fact, since the exponential approximation is not perfectly
overlapped with the true BR, as seen in Fig. 2, using the same mapping
function in both sides of the equation yields better results. In [21], we
followed a different approach, by replacing the BR straightforwardly with its
argument $|z|$, so that the resulting moment matching equation,

$z_{k}^{EP}=\sum_{m=0}^{M-1}\overline{\alpha}_{k}^{m}z_{k}^{m}\,,$ (20)

gives the Tikhonov parameter for the marginal of $\theta_{k}$ as a simple
linear combination of the mixture components’ parameters. In [21], our aim was
not to accurately approximate the BR but rather to avoid the compression of
the $|z_{k}^{EP}|$ values that, as stated, seems to be inherent to the native
EP message passing, due to the shape of the BR (and of its inverse), and tends
to depress the level of confidence in the transmitted messages. Although this
strategy showed satisfactory performance, in [21], the resulting algorithm
requires an additional monitoring of the precision values $|z_{k,f/b}|$, along
the forward and backward propagation of $p_{f/b}(\theta_{k})$, in order to
identify sudden drops of the precision values, which would lead inevitably to
errors in phase detection. Such a monitoring is implemented in [21] as a post-
processing of the computed precision values, followed by the rejection of
messages whose parameters are critical. This allows the algorithm to withstand
challenging scenarios, e.g, a gap of $4000$ symbols among preamble and
postamble without the aid of turbo iterations.

In this paper, we follow a different approach focusing on the accuracy of the
inverse BR approximation. In fact, because of the very steep shape of the BR
in this range, even approximations that are apparently quite similar may
produce very different results. We analyzed various kinds of approximations
and we found that a good compromise between simplicity and performance, in
terms of achievable results, is a piecewise function composed by a parabola,
for low values of $|z|$, followed by an exponential profile in the saturating
range. Therefore, the sought magnitude of the Tikhonov parameter is
$|z|=f^{-1}(BR)$, where $f^{-1}(\cdot)$ denotes the following inverse BR
approximation,

$\\!\\!\\!\\!f^{-1}(BR)\\!=\\!\left\\{\begin{array}[]{cl}\\!\\!\\!2.55-3.02\cdot\sqrt{0.71-BR},&\\!BR\leq
0.59\\\ -0.5/\log(BR)+0.55,&\text{elsewhere}\end{array}\right.\\!\\!\\!.$ (21)

In Section VII, simulation results demonstrate the significant advantages
related to this approximation which lead to a higher accuracy in the
estimation of the Tikhonov precision $|z|$ obtained by solving the moment
matching equation.

Figure 2: On the left, the comparison among the Bessel ratio, its commonly
adopted approximation and the proposed one. On the right, the comparison
between the two approximations for the inverse BR.

###  IV-B Rejection Strategy

The rejection of critical messages is, indeed, a well assessed technique
within the EP framework, whenever the projection and division operations (see
(15)) make improper distributions arise [22, 19]. These distributions are
easily recognizable when dealing with Gaussian mixtures, since the resulting
message is characterized by a negative variance. When, on the contrary, the
approximating family is that of Tikhonov distributions, as in the present
case, any value for the parameter $z_{k}$ in (7) characterizes a valid
distribution. Although improper distributions can sometimes be allowed to
propagate along the graph, when the EP algorithms faces critical conditions it
is usually advantageous to discard the corresponding updates [24] and specific
rejection criteria [22] become useful, and often necessary, for ensuring the
success of EP.

Typically, improper messages arise whenever there is an inconsistency between
the two factors of the distribution to be projected in (15), i.e, the
temporary prior, represented by the parameter $z_{k,u}$ in (17), and the
observation message $p_{d}(\theta_{k})$, which, being a mixture, is
represented by the set of values $z^{m}_{k}$ in (17), one for each possible
transmitted symbol.

The modified moment matching equation, including the approximation discussed
in Section IV-A, is

$\\!\\!\\!f(|z^{EP}_{k}|){\rm e}^{j\angle
z_{k}^{EP}}\\!\\!=\\!\\!\sum_{m=0}^{M-1}\overline{\alpha}_{k,MIX}^{m}\frac{I_{1}(|z_{k,MIX}^{m}|)}{I_{0}(|z_{k,MIX}^{m}|)}{\rm
e}^{j\angle z_{k,MIX}^{m}}\,,$ (22)

where the definition of ${\alpha}_{k,MIX}^{m}$ is analogous to that in (9) and
$\overline{\alpha}_{k,MIX}^{m}$ are normalized coefficients, as in (12) .
Before the calculation of (22), a real-time monitoring is introduced on the
mixture parameters $z^{m}_{k,MIX}$ that result from (17). We define the angle

$\gamma^{m}_{k}=\arg\left[z^{m}_{k,MIX}z^{*}_{k,u}\right]\,,$ (23)

where the superscript $m$ refers to the $m$-th mixture mode, which quantifies
the inconsistency between the phase estimate brought about by each mode and
the one provided by the prior distribution. Considering the PN evolution until
a given time epoch $k$, the argument of $z_{k,u}$ typically yields a stable
estimate of phase noise, which should be refined by the current channel
observation. When this is very noisy, or there is little agreement between the
phase angles, i.e., the absolute value of (23) is large, the inclusion of the
observation message can introduce spurious phase slips in the PN estimate. We
verified that this problem occurs when the estimated marginal distribution of
$\theta_{k}$ before the EP projection, i.e.,
$p_{d}(\theta_{k})p_{u}(\theta_{k})$, is a mixture of modes (associated to the
complex parameter $z^{m}_{k,MIX}$) which is spread around the (temporary)
prior distribution of phase noise, $p_{u}(\theta_{k})$ (associated to the
complex parameter $z_{k,u}$). In this case, the projection produces a rather
flat Tikhonov distribution, i.e., characterized by a small precision
$\left|z_{k}^{EP}\right|$ and a phase estimate $\arg\left[z_{k}^{EP}\right]$
which gets far from both the prior and the mixture modes, thus creating an
uncertainty condition in the estimation process from which it is sometimes
difficult to recover, in the absence of pilot symbols.

In order to avoid such situations, the specific rejection strategy that we
adopt is thus based on the inconsistency angles (23). We assume that the
$m$-th mixture mode is inconsistent if
$\left|\gamma^{m}_{k}\right|>\Gamma_{th}$, i.e, if its absolute angular
deviation from the prior phase angle $\arg\left[z_{k,u}\right]$ exceeds a
given threshold $\Gamma_{th}$. Before performing the projection (22), at each
time epoch $k$, the algorithm checks the rejection condition,

$\sum^{M-1}_{m=0}\mathds{1}(\left|\gamma^{m}_{k}\right|>\Gamma_{th})>\overline{M}\,,$
(24)

where $\mathds{1}(\cdot)$ is an indicator function equal to $1$ if the
condition is true and $0$ otherwise, to verify if the number of inconsistent
modes is above a given limit $\overline{M}<M$. When (24) is satisfied, the
observation message $p_{d}(\theta_{k})$ is considered unreliable and the
information contained in it is rejected by setting $z^{m}_{k}$ to zero in
(17). Hence, both $p_{d}(\theta_{k})$ and $p^{EP}_{d}(\theta_{k})$ are treated
as uniform distributions, which do not carry any information on the phase,
whose marginal equals that of the prior, i.e., a Tikhonov distribution with
parameter to $z_{k,u}$.

The values of $\Gamma_{th}$ and $\overline{M}$ (24) depend on the modulation
format ($M$) and on the severity of phase noise ($\sigma_{\Delta}$).
Especially at high signal-to-noise ratios (SNRs), there can be more than one
suitable pair and, possibly, multiple rejection conditions can be adopted by
logical disjunction ($\vee$) as, e.g.,

$\left(\sum^{M-1}_{m=0}\\!\\!\mathds{1}(\left|\gamma^{m}_{k}\right|>\\!\Gamma_{th,1})>\overline{M}_{1}\\!\right)\vee\left(\sum^{M-1}_{m=0}\\!\\!\mathds{1}(\left|\gamma^{m}_{k}\right|>\\!\Gamma_{th,2})>\overline{M}_{2}\\!\right)$
(25)

in the case of two message rejection conditions, whose parameters should be
set such that $\Gamma_{th,1}<\Gamma_{th,2}$ and
$\overline{M}_{1}>\overline{M}_{2}$. In any case, as stated, the rejection
condition is evaluated before the EP projection operation, thus reducing the
number of unnecessary operations.

###  IV-C Message Damping

Another solution to overcome the convergence problems of EP, which is known in
the literature, is the use of “damped” updates [14], [22]. In iterative
parametric message passing, the damping technique consists in propagating a
convex combination of the new parameter value and the old one. In our case,
instead of computing from (15)

${z}^{EP}_{k,d}=z^{EP}_{k}-z_{k,u}\,,$ (26)

which is the parameter of $p_{d}^{EP}(\theta_{k})$ at the current $n$-th
iteration, the linear combination

$z^{(n)}_{k,d}=\xi\cdot(z^{EP}_{k}-z_{k,u})+(1-\xi)\cdot z^{(n-1)}_{k,d}$ (27)

is used, which partly recovers the value computed at the previous, $(n-1)$-th,
inner detector iteration. The value for the damping parameter $\xi$ must be
selected, as discussed in Section VII, where we apply damping and demonstrate
its advantages through numerical results.

##  V Message Scheduling

Regarding message scheduling, for any algorithm, we take advantage of the
structure of the FG, which evidences the logical separation of the decoding
part from the phase detection part, corresponding to the upper and lower
halves of Fig. 1. Indeed, $P_{d}(c_{k})$ is the extrinsic information on the
$k$-th code symbol, sent from the decoder to the phase detector, as well as
$p_{u}(\theta_{k})$ which, as stated, represents the prior message in (13) and
is thus seen as a temporary estimate of phase noise, sent towards the decoder.

We shall assume that the phase-detector subgraph operates first, by exchanging
horizontal messages along the Markov chain in the lower part of Fig. 1, so as
to update $p_{f}(\theta_{k})$ and $p_{b}(\theta_{k})$ by message passing (see
[1], to recall the rules of computation). This horizontal message passing
could be iterated, with multiple passes, before sending vertical messages to
the decoder subgraph so that the phase estimate could be refined. A wide
variety of message scheduling can be adopted, such as a flooding schedule
[25], a forward filtering followed by either a backward filtering or a
backward smoothing and so on. After a detailed comparison, we concluded that
the best performing scheduling is a forward/backward filtering, composed by
two independent passes, one for each direction. Therefore, during the first
inner-detector iteration the temporary prior parameter $z_{k,u}$ is equal to
either $z_{k,f}$ or $z_{k,b}$, depending on the considered direction. In other
words, the opposite message has no impact on the observation message
approximation. On the other hand, from the $2$nd iteration onwards, the
temporary prior is determined by the current message in the analyzed direction
entering the variable node $\theta_{k}$ and the previous iteration’s opposite
message, i.e.,

$z^{(n)}_{k,u}=z^{(n)}_{k,f}+z^{(n-1)}_{k,b}$ (28)

during the forward pass and vice versa for the backward pass. The number of
“inner” detector iterations is referred to as $N_{D}$.

After $N_{D}$ phase-detection iterations, vertical messages
$p_{u}(\theta_{k})$ are sent from the phase detector to the upper part of the
FG, implementing the decoder, where an inner message passing procedure takes
place, until convergence or until a maximum number of decoding iterations is
reached.

The exchange of information, through vertical messages, between detector and
decoder can be iterated for $N_{T}$ times,444In order to indicate the number
of inner detector, inner decoder and turbo iterations, the notation
$N_{D}$-$N_{C}$-$N_{T}$ is used. so as to yield a refinement of both the
symbols’ and the phase samples’ estimates, from which the other receiver half
can benefit. On the other hand, detector-decoder iterations are considered
impractical and they are usually avoided in the implementation of digital
receivers, especially at high baud rates.

Since phase detection is performed before symbol decoding, the
forward/backward messages $p_{f/b}(\theta_{k})$ propagated at the first turbo
iteration cannot benefit from any symbol estimate, hence phase estimation is
accomplished (locally) in a non data-aided fashion, while a data-aided phase
estimation occurs, instead, at time epochs where pilot symbols are inserted.
In the case of turbo iterations, on the contrary, a (provisional) soft-
decision-directed strategy is implemented, thanks to the extrinsic information
$P_{d}(c_{k})$, even for the payload symbols.

##  VI Complexity Analysis

| Two-terms additions | Two-terms multiplications | LUT accesses
---|---|---|---
TP [1] | $7M+12$ | $11M+22$ | $2M+2$
EP | $16M+18$ | $34M+25$ | $12M+8$
dp-BCJR | $5N_{\theta}^{2}+(18M-6)N_{\theta}-(M+2)$ | $N_{\theta}(14M+1)+1$ | $2N_{\theta}^{2}+N_{\theta}(3M-1)-M$
TABLE I:  Complexity comparison per code symbol per iteration

We compared the computational complexity of the proposed algorithm based on
the EP framework with the TP algorithm proposed in [1] and the dp-BCJR.
Referring to the $k$-th time epoch, we considered the operations needed for
the computation of the observation message $p_{d}(\theta_{k})$ approximation,
the forward $p_{f}(\theta_{k})$ and backward $p_{b}(\theta_{k})$ messages and,
finally, the message $P_{u}(c_{k})$. In order to compute the algorithmic
complexity, we counted the number of required operations distinguishing among
two real terms multiplications, additions and accesses to the lookup table
(LUT) for the computation of nonlinear functions. The total number of required
operations is reported in Table I. As it can be seen, the complexity of all
the algorithms increases linearly with the constellation cardinality $M$.
However, for the dp-BCJR algorithm, it is much higher. In fact, a critical
parameter which has a strong impact on the computational load is the number of
discretization levels $N_{\theta}$. An empirical rule for the choice of this
parameter is provided in [1], which suggests that it is linearly related to
$M$ according to $N_{\theta}=8M$. We found that the value obtained through
this rule of thumb is not sufficient for an accurate quantization of the PN
variables, which requires at least $512$ levels in the QPSK case.

Figure 3: Comparison among the complexity of the analyzed algorithms.

In order to visualize the complexity gap among the analyzed algorithms, the
number of required operations is plotted versus $M$ in Fig. 3. It is possible
to notice that the EP and TP algorithms have a comparable complexity whereas
the dp-BCJR is characterized by a much higher computational load and hence
represents a totally impractical solution. Because of this reason and since,
with a sufficiently large $N_{\theta}$, achieves the performance of the plain
SPA, it is considered here just as a performance benchmark.

##  VII Simulation Results

We analyzed the performance of different receivers for a coded transmission,
employing a low-density-parity-check (LDPC) encoder, through the system under
investigation. Bit error rate (BER) is plotted versus the SNR $E_{b}/N_{0}$,
where $N_{0}=\sigma^{2}$ is the variance per component of noise samples in (1)
and $E_{b}$ is the average energy per information (payload) bit.

In all simulation results, in order to allow the algorithms convergence, pilot
symbols, which are known at the receiver, were inserted in the transmitted
codeword. The presence of pilots involves a slight decrease of the effective
information rate, resulting in an increase in the required average energy per
bit $E_{b}$. In order to make a fair comparison, this increase has been
artificially introduced in the curve referring to the Known Phase scenario,
even though pilot symbols were useless in this case since there was no need
for phase estimation. Thus, the Known Phase curve shows the performance of the
adopted LDPC code. Furthermore, we included the All Pilots curve, which
represents an ideal bound since a “genie-aided” receiver is assumed to perform
the estimation of PN samples resorting to a perfect knowledge of the entire
transmitted codeword. As discussed in Section VI, the dp-BCJR algorithm is
regarded as a performance benchmark. The number of phase discretization levels
is set to $N_{\theta}=512$ for all the analyzed modulation formats.

The performance of the algorithms analyzed in Figs. 4 and 5 refers to the
practical case of separate detection and decoding, i.e., the absence of turbo
iterations ($N_{T}=1$). The maximum number of LDPC decoding iterations $N_{C}$
is set to $200$.

In Fig. 4, we considered the transmission of a block of $2000$ coded symbols
punctured by distributed pilot symbols. The transmitted sequence started with
a pilot, then one pilot symbol was regularly alternated with $19$ payload
symbols, amounting to a $5.3\%$ overhead555We use the notation $\\{$number of
pilots in a block$\\}/\\{$payload symbols among two consecutive pilot
blocks$\\}$ to describe a pilot distribution.. The payload symbols were the
output of a $(3,6)$-regular LDPC code with rate-$1/2$ and length $4000$ [26].
The quadrature phase-shift keying (QPSK) modulation format ($M=4$) was
adopted, since it is robust to the strong phase noise that was assumed, with
$\sigma_{\Delta}=6^{\circ}$. Although representing an ideal situation, the
performance of the All Pilots algorithm shows a loss of about $0.2\,$dB,
compared to the Known Phase curve, since PN still has to be estimated in the
presence of AWGN.

Figure 4: Performance for an LDPC-coded [26] QPSK modulation in the separate
detection and decoding scenario ($N_{T}=1$). The phase noise standard
deviation is $\sigma_{\Delta}=6\degree$ and $(1/19)$ pilot symbols are
distributed among the payload (for other parameters, see text).

The performance of the dp-BCJR algorithm is more than $1\,$ dB away from the
theoretical All Pilots reference. The TP algorithm of [1] shows
$0.4\,\text{dB}$ of penalty, compared to the practical dp-BCJR limit.
Therefore, thanks to the help of distributed pilots among the payload, a
satisfactory performance is achieved by this low-complexity algorithm.

Moving to the primary object of our investigation, we compared the native
implementation of the EP algorithm (labelled EP in figures), which employs the
exponential BR approximation as in (19), with two improved variants. The first
one only exploits the damping technique described in Section IV-C and is
labelled EP Damping in figures. The second one relies, in addition to damping,
on the inverse BR approximation (21) and on the message rejection strategies
discussed in Section IV-B (labelled EP Mod in figures). From Fig. 4, we notice
that the native EP algorithm loses about $2\,\text{dB}$, compared to the dp-
BCJR benchmark, a penalty that increases towards higher SNRs, where a slight
change of slope can be noticed. The introduction of damping allows to
partially solve the convergence problems of EP, reducing the loss to
$1\,\text{dB}$. By numerical simulations, we found that the best value for the
damping parameter is $\xi=0.4$, hence we used this value in all simulations,
if not otherwise stated. Since the native EP and its damped version do not
benefit from inner detector iterations, we set $N_{D}=1$, so that the
application of damping actually reduces to an attenuation of the precision
values, i.e., the propagation of message parameter (27) is nothing but (26)
multiplied by $\xi$.

The performance of the Modified EP algorithm (EP Mod) shows that the adoption
of a more accurate moment matching equation (22) combined with a suitable
message rejection strategy allows to reach the performance benchmark, with a
remarkably lower complexity, compared to dp-BCJR. Therefore, the Modified EP
algorithm outperforms not only the native and damped versions of EP analyzed
above but also the simpler TP algorithm. In this algorithm, inner detector
iterations were indeed useful to refine the phase estimate and $N_{D}=2$ was
sufficient for achieving the convergence of phase estimation. We investigated
by simulation the parameters adopted for the rejection condition (24) and
found that $\Gamma_{th}=\pi/2$ and $\overline{M}=0$ achieve the optimal
performance reported in Fig. 4.

Figure 5: Performance for QPSK and $16$-QAM modulations in the separate
detection and decoding scenario ($N_{T}=1$) employing the DVB-S2 LDPC code
(rate $1/2$) and pilots positioning ($36/1440$) [27]. The phase noise standard
deviation is $\sigma_{\Delta}=1\degree$ (for other parameters, see text).

In Fig. 4, however, the successful operation of the proposed solution occurs
in a context where the placement of pilots is quite favorable. Nevertheless,
in many communications standards the position of pilots is predefined.
Therefore, with the aim of designing a receiver able to cope with sparser
pilots’ distributions, in Fig. 5 we tested the same algorithms analyzed above
in the case of the DVB-S2 (Digital Video Broadcasting Satellite Second
Generation) system, where pilot symbols are organized into bursts of $36$
symbols every $1440$ information symbols with the exception of a longer
sequence of $90$ symbols in the preamble and postamble [27]. We considered the
standardized LDPC code with codeword length $64800$ and rate $1/2$ whereas the
analyzed modulation formats are QPSK and $16$-QAM. For this system, we assumed
$\sigma_{\Delta}=1\degree$ for the phase noise and we still focus on receivers
employing separate detection and decoding, i.e., $N_{T}=1$.

As seen in Fig. 5, the algorithm that is most vulnerable to the sparser
distribution of pilots is TP, as expected. In fact, a loss of $1\,\text{dB}$
with respect to the dp-BCJR benchmark is registered for the QPSK modulation
while, moving to a higher constellation cardinality, it considerably increases
reaching almost $3\,\text{dB}$ for $16$-QAM. The penalty is reduced to about
$0.7\,\text{dB}$ for both modulation formats, using the native EP algorithm,
which however shows the onset of a floor around $10^{-4}$ for the $16$-QAM and
around $10^{-5}$ for the QPSK. The adoption of damping improves the
performance but an error floor is still present for the $16$-QAM while for the
QPSK there is still a residual loss with respect to the benchmark. It is only
with the proposed Modified EP algorithm that the BER curves in Fig. 5
practically coincide with those of the dp-BCJR benchmark. The refinement of
phase estimates converges after $N_{D}=2$ inner detector iterations. The
employed rejection condition is in the form (25) with
$\left(\Gamma_{th,1},\Gamma_{th,2}\right)=\left(\frac{\pi}{12},\frac{\pi}{6}\right)$
for both modulations,
$\left(\overline{M}_{1},\overline{M}_{2}\right)=\left(1,0\right)$ for the QPSK
and $\left(\overline{M}_{1},\overline{M}_{2}\right)=(12,7)$ for the $16$-QAM.

Figure 6: Performance for QPSK modulation in the joint detection and decoding
scenario ($1$-$1$-$200$ iterations) employing the DVB-S2 LDPC code (rate
$1/2$) and concentrated pilots, positioned in a preamble and postamble each
composed by $90$ symbols and separated by $32400$ information symbols. The
phase noise standard deviation is $\sigma_{\Delta}=1\degree$.

In Fig. 6, we analyzed the most challenging scenario where pilot symbols are
concentrated in the preamble and postamble of the data packet, both consisting
of $90$ symbols. We still employed the DVB-S2 LDPC-code with rate $1/2$ and
the QPSK modulation, so that a block of $32400$ information symbols is
transmitted without any pilot in-between. We assumed
$\sigma_{\Delta}=1\degree$ and, since none of the tested receivers could
achieve a satisfying performance when $N_{T}=1$, we considered the joint
detection and decoding scenario, setting the maximum number of allowed
iterations to $1$-$1$-$200$.

Due to the absence of distributed pilots, the TP algorithm fails completely,
as expected. In fact, in the presence of a (non pilot) payload symbol $c_{k}$,
the observation message $p_{d}(\theta_{k})$ is a uniform pdf, that does not
carry any information, since the decoder has no knowledge of the transmitted
bits. Hence, the TP algorithm cannot bootstrap, even with multiple turbo
iterations. On the contrary, the algorithms based on the EP framework benefit
from the exchange of messages between detector and decoder because the PN
estimates approximation is carried out after combining the observation message
$p_{d}(\theta_{k})$ with the prior distribution $p_{u}(\theta_{k})$.
Nevertheless, as evidenced by Fig. 6, the native implementation of EP yields
unsatisfactory performance and the use of damping allows to gain about
$1\,\text{dB}$ but the gap from the ideal All Pilots curve is still
significant.

The adoption of the proposed improvements leads to an excellent performance
achieving an almost perfect overlap to the All Pilots curve. For this scenario
we set $\xi=0.5$ and we exploited, also inside the rejection strategy, the
knowledge of the current estimate of the symbols pmf $P_{d}(c_{k})$, provided
by the decoder. In fact, we set666We denote with $[\cdot]$ the rounding
function to the closest integer.

$\displaystyle{(\Gamma_{th,1},\Gamma_{th,2})=\left(\frac{\pi}{6},\frac{\pi}{4}\right)\,,}$
(29) $\displaystyle{(\overline{M}_{1},\overline{M}_{2})=([2\max
P_{d}(c_{k})],[\max P_{d}(c_{k})])}$ (30)

so that the clearer the decoder’s belief regarding the transmitted symbol, the
more stringent the condition to be met for rejection. Moreover, with reference
to the $n$-th turbo iteration, we apply message rejection only if (25) and

${\max_{c_{k}\in\mathcal{A}}P^{(n)}_{d}(c_{k})\leq\max_{c_{k}\in\mathcal{A}}P^{(n-1)}_{d}(c_{k})\,,}$
(31)

are both satisfied. In fact, condition (31) identifies the cases where the
iterations are not properly refining the symbols and PN estimates. Finally, at
high SNRs, the value of $N_{0}$ that is passed to the detector is slightly
larger than the actual one to avoid the overestimation of the message
reliability as commonly done in the literature when dealing with suboptimal
iterative algorithms.

##  VIII Conclusions

We analyzed the expectation propagation algorithm operating in a digital
receiver for transmissions over channels affected by Wiener phase noise. We
identified and discussed the weaknesses of the native implementation of EP and
examined various supplementary techniques that could be applied to enhance its
performance. Our investigation highlighted the challenges that arise in
applying the EP framework to the considered context and, while solutions taken
from the existing literature proved beneficial, they still presented some
relevant limitations.

By carefully analyzing the message passing procedure, we found that the EP
rules of computation for the messages can sometimes induce a collapse of the
level of confidence (precision) in the resulting estimates, hence lead to
phase slips and subsequent errors.

To address this, we introduced a more accurate approximation of the message
projection rule of EP and proposed a novel message rejection strategy aimed at
identifying and discarding misleading messages. These countermeasures prevent
the appearance of spurious phase detection errors so as to reach the
performance of the reference benchmark even in the most challenging scenarios
where either the separation of phase detection and decoding is required, or
the absence of distributed pilots is imposed by the adopted standard.

## References

  * [1] 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.
  * [2] A. Barbieri and G. Colavolpe, “Simplified soft-output detection of CPM signals over coherent and phase noise channels,” _IEEE Trans. Wireless Commun._ , vol. 6, no. 7, pp. 2486–2496, Jul. 2007.
  * [3] G. Colavolpe, “Communications over phase noise channels: a tutorial review,” _International Journal of Satellite Communications and Networking_ , vol. 32, pp. 167–185, May/June 2014, article first published online: July 2013\.
  * [4] M. Peleg, S. Shamai (Shitz), and S. Galán, “Iterative decoding for coded noncoherent MPSK communications over phase-noisy AWGN channel,” _IEE Proc. Commun._ , vol. 147, pp. 87–95, Apr. 2000.
  * [5] R. Krishnan, M. R. Khanzadi, T. Eriksson, and T. Svensson, “Soft metrics and their performance analysis for optimal data detection in the presence of strong oscillator phase noise,” _IEEE Transactions on Communications_ , vol. 61, no. 6, pp. 2385–2395, 2013.
  * [6] A. Kreimer and D. Raphaeli, “Efficient low-complexity phase noise resistant iterative joint phase estimation and decoding algorithm,” _IEEE Transactions on Communications_ , vol. 66, no. 9, pp. 4199–4210, 2018.
  * [7] 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 Transactions on Signal Processing_ , vol. 60, no. 9, pp. 4790–4807, 2012.
  * [8] B. D. Woerner and W. E. Stark, “Trellis-coded direct-sequence spread-spectrum communications,” _IEEE Trans. Commun._ , vol. 42, pp. 3163–3170, Dec. 1994\.
  * [9] G. Colavolpe and A. Modenini, “Iterative carrier synchronization in the absence of distributed pilots for low SNR applications,” in _6th ESA International Workshop on Tracking, Telemetry and Command Systems for Space Applications (TTC’13)_. European Space Agency, 2013.
  * [10] S. Shayovitz and D. Raphaeli, “Message passing algorithms for phase noise tracking using tikhonov mixtures,” _IEEE Transactions on Communications_ , vol. 64, no. 1, pp. 387–401, 2016.
  * [11] T. P. Minka, “Expectation propagation for approximate Bayesian inference,” in _Proc. Conf. on Uncertainty in Artificial Intelligence (UAI)_ , Seattle, Washington, USA, Aug. 2001.
  * [12] ——, “Divergence measures and message passing,” Microsoft Research Cambridge, Tech. Rep., 2005, mSR-TR-2005-173 (2005).
  * [13] T. Heskes and O. Zoeter, “Expectation propagation for approximate inference in dynamic Bayesian networks,” in _Proc. Conf. on Uncertainty in Artificial Intelligence (UAI)_ , Edmonton, Canada, Jul. 2002.
  * [14] A. Vehtari, A. Gelman, T. Sivula, P. Jylänki, D. Tran, S. Sahai, P. Blomstedt, J. P. Cunningham, D. Schiminovich, and C. Robert, “Expectation propagation as a way of life: A framework for Bayesian inference on partitioned data,” 2014\. [Online]. Available: https://arxiv.org/abs/1412.4869
  * [15] A. Vannucci, G. Colavolpe, and L. Veltri, “Estimation of correlated gaussian samples in impulsive noise,” _IEEE Communications Letters_ , vol. 24, no. 1, pp. 103–107, 2020.
  * [16] A. Mirbadin, A. Vannucci, and G. Colavolpe, “Expectation propagation and transparent propagation in iterative signal estimation in the presence of impulsive noise,” in _2020 43rd International Conference on Telecommunications and Signal Processing (TSP)_ , 2020, pp. 175–179.
  * [17] A. Mirbadin, A. Vannucci, G. Colavolpe, R. Pecori, and L. Veltri, “Iterative receiver design for the estimation of gaussian samples in impulsive noise,” _Applied Sciences_ , vol. 11, no. 2, 2021. [Online]. Available: https://www.mdpi.com/2076-3417/11/2/557
  * [18] Y. Qi and T. P. Minka, “Window-based expectation propagation for adaptive signal detection in flat-fading channels,” _IEEE Trans. Wireless Commun._ , vol. 6, no. 1, pp. 348–355, 2007.
  * [19] E. Conti, A. Piemontese, G. Colavolpe, and A. Vannucci, “Expectation propagation for flat-fading channels,” _IEEE Wireless Communications Letters_ , vol. 12, no. 11, pp. 1846–1850, 2023.
  * [20] L. Szczecinski, H. Bouazizi, and A. Aharony, “Parametric phase tracking via expectation propagation,” 2020. [Online]. Available: https://arxiv.org/abs/2005.01844
  * [21] G. Colavolpe, E. Conti, A. Piemontese, and A. Vannucci, “The difficult road of expectation propagation towards phase noise detection,” in _ICC 2023 - IEEE International Conference on Communications_ , Rome (Italy), 2023, pp. 4640–4645.
  * [22] M. Seeger, “Bayesian Gaussian process models: PAC-bayesian generalisation error bounds and sparse approximations,” _PhD thesis, University of Edinburgh_ , 2003.
  * [23] K. V. Mardia and P. E. Jupp, _Directional statistics_. Wiley Online Library, 2000, vol. 2.
  * [24] E. Conti, A. Piemontese, G. Colavolpe, and A. Vannucci, “Detection for flat-fading channels based on expectation propagation,” 2023, available at https://arxiv.org/abs/2305.19439.
  * [25] F. R. Kschischang, B. J. Frey, and H.-A. Loeliger, “Factor graphs and the sum-product algorithm,” _IEEE Trans. Inform. Theory_ , vol. 47, pp. 498–519, Feb. 2001.
  * [26] D. MacKay, “http://www.inference.org.uk/mackay/,” web page.
  * [27] ETSI, “ETSI - DVBS2 74r13, Digital Video Broadcasting (DVB): Second generation framing structure, channel coding and modulation systems for Braoadcasting, Interactive Services, News Gathering and other broadband satellite applications,” 2003.