# Joint CFO, Gridless Channel Estimation and Data Detection for Underwater
Acoustic OFDM Systems

Lei Wan, Jiang Zhu, En Cheng and Zhiwei Xu Lei Wan and En Cheng are with
Department of Information and Communication Engineering and Key Laboratory of
Underwater Acoustic Communication and Marine Information Technology (Xiamen
University), Ministry of Education, China (email: {leiwan,
chengen}@xmu.edu.cn). Jiang Zhu and Zhiwei Xu are with the engineering
research center of oceanic sensing technology and equipment, Ministry of
Education, Ocean College, Zhejiang University, and also with the key
laboratory of ocean observation-imaging tested of Zhejiang Province, No.1
Zheda Road, Zhoushan, 316021, China (email: {jiangzhu16, xuzw}@zju.edu.cn).
The corresponding author is Jiang Zhu (email: jiangzhu16@zju.edu.cn).

###### Abstract

In this paper, we propose an iterative receiver based on gridless variational
Bayesian line spectra estimation (VALSE) named JCCD-VALSE that _j_ ointly
estimates the _c_ arrier frequency offset (CFO), the _c_ hannel with high
resolution and carries out _d_ ata decoding. Based on a modularized point of
view and motivated by the high resolution and low complexity gridless VALSE
algorithm, three modules named the VALSE module, the minimum mean squared
error (MMSE) module and the decoder module are built. Soft information is
exchanged between the modules to progressively improve the channel estimation
and data decoding accuracy. Since the delays of multipaths of the channel are
treated as continuous parameters, instead of on a grid, the leakage effect is
avoided. Besides, the proposed approach is a more complete Bayesian approach
as all the nuisance parameters such as the noise variance, the parameters of
the prior distribution of the channel, the number of paths are automatically
estimated. Numerical simulations and sea test data are utilized to demonstrate
that the proposed approach performs significantly better than the existing
grid-based generalized approximate message passing (GAMP) based _j_ oint _c_
hannel and _d_ ata decoding approach (JCD-GAMP). Furthermore, it is also
verified that joint processing including CFO estimation provides performance
gain.

keywords: Iterative receivers, message passing, expectation propagation, grid-
less, OFDM channel estimation

##  I Introduction

Due to the propagation characteristics, underwater acoustic communications
typically suffer from strong multipath and Doppler effect [1]. Since the
propagation speed is low (around 1500 m/s), the reflection of acoustic waves
from surface and bottom of the ocean, along with the refraction caused by
nonuniformity of sound speed in water, cause significant multipath delay
spread. The typical delay spread of underwater acoustic communication channels
is from several milliseconds to tens of milliseconds, and could be hundreds of
milliseconds and more in extreme cases. On the other hand, due to the movement
of transmitter and receiver platforms, as well as water media, Doppler effect
is almost inevitable in underwater acoustic communications. Considering the
wide band nature of underwater acoustic communication systems, the Doppler
effect usually exhibits itself as signal dilation or compression [2].

With the advantages of low complexity frequency domain multipath channel
equalization and high bandwidth efficiency, orthogonal frequency division
multiplexing (OFDM) technology was introduced into underwater acoustic
communications, and soon became a popular choice especially for high speed
communications [3, 4, 5]. However, Doppler effect in underwater acoustic
channels could destroy the orthogonality between the subcarriers in OFDM and
hence deteriorates its performance. In [4], a two-step approach was proposed
to compensate the Doppler distortion, in which a resampling operation is
carried out first to convert the wide band Doppler distortion into narrow band
carrier frequency offset (CFO). In the second step, the CFO is estimated and
compensated. Experimental results in [4] demonstrated that the proposed scheme
can effectively mitigate the Doppler distortion.

In coherent underwater acoustic communication systems, accurate channel
estimation is essential. Since compressed sensing based channel estimation
methods exploit the structure of the channel and require less number of
measurements, they have been the most popular choices for underwater acoustic
OFDM systems [6, 7]. Among compressed sensing algorithms, orthogonal matching
pursuit (OMP) algorithm [8] features lower computational complexity than basis
pursuit methods, and it is widely adopted in channel estimations [9, 10]. On
the other hand, as an approximation for sum-product algorithm (SPA),
approximate message passing (AMP) has been utilized for sparse linear
estimation problems [11], and hence is also employed for sparse channel
estimations [12].

To further improve the performance of wireless communication systems in
challenging channel conditions, joint receivers in which modules such as
channel estimation, channel equalization interactively work together are
proposed, both for radio and underwater acoustic communication systems. In
[13], a joint channel estimation, equalization and data detection system based
on expectation-maximization (EM) algorithm is proposed for OFDM systems with
high mobility. In the proposed scheme, the time-varying channel is represented
by discrete cosine orthogonal basis functions and the equalized symbols are
quantized into nearest data symbol constellation point in each iteration. A
low complexity generalized AMP (GAMP) based joint channel estimation and data
decoding approach for OFDM system is proposed [14]. Later, a joint channel
estimation and symbol detection scheme based on parametric bilinear GAMP
algorithm is proposed for single-carrier systems [15], which supports the
sparse characteristic of wireless communication channels and soft-input soft-
output (SISO) equalizer. Recently, based on Markov chain Monte Carlo (MC)
algorithm with Gibbs sampling and convolutional coding plus differential
coding, a joint phase shift keying (PSK) data detection and channel estimation
scheme utilizing channel sparsity is proposed in [16].

For underwater acoustic communications, the concept of joint processing has
been widely pursued. In the early work [17], an iterative receiver based on
message passing techniques is proposed for single carrier underwater acoustic
communication systems, in which equalization and decoding are performed
multiple times in a turbo manner. However, channel estimation is carried out
based on training symbols, utilizing maximum-likelihood approach. In [18, 19,
20], iterative receivers were developed for single carrier multiple-input
multiple-output (MIMO) underwater acoustic communication systems, for both
time-domain and frequency domain equalization. In the proposed iterative
receivers, soft equalization outputs are feedback for interference cancelation
and channel estimation. For OFDM underwater acoustic communication systems in
doubly spread channels, an iterative sparse channel estimation and data
detection receiver using partial interval demodulation is proposed in [21]. In
[22], an iterative receiver is proposed for inter-frequency interference (IFI)
cancelation for single-carrier underwater acoustic communication systems
adopting frequency domain equalization. Channel impulse responses, explicit
and residual IFI are updated in an iterative fashion based on the soft symbol
estimates in the proposed receiver. Considering the impulse noise in
underwater acoustic communication channels, a joint channel estimation and
impulsive noise mitigation receiver utilizing sparse Bayesian learning (SBL)
algorithm is proposed [23], in which the information on data symbols are
employed to carry out joint time-varying channel estimation and data
detection, to further improve the system performance. The performance of the
proposed algorithm is verified through both numerical simulations and real
data. More recently, a joint channel estimation and equalization scheme with
superimposed training is proposed for single carrier underwater acoustic
communication systems [24]. A message-passing-based bidirectional channel
estimation algorithm is adopted in the proposed scheme, utilizing the
correlation between consecutive channels.

It is worth noting that all the above compressed sensing based approaches
utilize the fixed dictionary and suffer grid mismatch [25], and the gridless
parametric channel estimation approach indeed improves both the channel
estimation and data decoding performance over the grid based approaches [26,
27]. Moreover, existing iterative receivers seldom incorporate the CFO
estimation and compensation into the design. In this paper, according to
message passing and from the modularized point of view [28], an iterative
receiver based on gridless variational Bayesian line spectra estimation
(VALSE) [29] which _j_ ointly estimates the residual _C_ FO and _c_ hannel,
meanwhile also achieves _d_ ata decoding named JCCD-VALSE is proposed for
underwater acoustic OFDM systems. The JCCD-VALSE consists of three modules
named as high resolution gridless channel estimation, minimum mean squared
error (MMSE) estimation, and low-density parity-check code (LDPC) decoder.
Through the exchange of soft information between modules, the channel
estimation and data decoding accuracy progressively improves. Compared with
existing receivers, the joint processing including residual CFO estimation and
compensation, the gridless channel estimation, and the feature of automatical
estimation of nuisance parameters are the advantages of the proposed scheme.
Hence, the contributions of this paper include:

  1. 1.

Targeting at underwater acoustic OFDM systems in channels with an
approximately common Doppler scaling factor, the system model for received
signal after resampling is derived, and the corresponding receiver JCCD-VALSE
which jointly estimates residual CFO, channel, and finishes data decoding is
proposed.

  2. 2.

In JCCD-VALSE, all the nuisance parameters such as the number of paths, the
noise variance are automatically estimated. This reduces the reliance on
priori information, and also improves the estimation accuracy.

  3. 3.

In JCCD-VALSE, gridless channel estimation based on VALSE is adopted, which
improves the channel estimation and data decoding performance, compared with
the existing channel estimation and data decoding methods. Besides, it is also
shown that taking residual CFO into consideration benefits both the channel
estimation and data decoding performances.

  4. 4.

Sea trial data of mobile acoustic communication experiment (MACE10) is
utilized to further demonstrate the effectiveness of the proposed scheme.

In addition, it is also worth noting that our proposed approach is very
flexible and can deal with both uniform and nonuniform pilots.

The rest of this paper is organized as follow: Section II introduces the
system model. The probabilistic formulation of the channel, the data symbol
and the normalized residual CFO are presented in Section III. The ensuing
Section IV proposes the joint CFO, channel estimation and data decoding
algorithm. Section V presents the performance metrics in terms of the coded
bit error rate (BER) and the normalized mean squared error (NMSE) of the
various algorithms. The sea trial data decoding results are introduced in
Section VI. Finally, we conclude the paper in Section VII.

Notation: For a matrix $\mathbf{A}$, let $[{\mathbf{A}}]_{i,j}$ or $A_{ij}$
denote the $(i,j)$th element of ${\mathbf{A}}$, and ${\rm diag}(\mathbf{A})$
returns a vector with elements being the diagonal elements of $\mathbf{A}$.
Let $\Re\\{\cdot\\}$ and $\Im\\{\cdot\\}$ denote the real and imaginary part
operator, respectively. For random vectors $\mathbf{x}$ and $\mathbf{y}$ with
joint probability density function (PDF) $p({\mathbf{x}},{\mathbf{y}})$, the
marginal PDF $p({\mathbf{x}})=\int p({\mathbf{x}},{\mathbf{y}}){\rm
d}{\mathbf{y}}$ with ${\rm d}{\mathbf{y}}$ being ${\rm
d}\Re\\{\mathbf{y}\\}{\rm d}\Im\\{\mathbf{y}\\}$. Let $\odot$ denote the
Hadamard product operator. For a vector ${\mathbf{a}}$, let ${\rm
diag}({\mathbf{a}})$ return a diagonal matrix with diagonal elements being
$\mathbf{a}$. Let ${\mathcal{S}}\subset\\{1,\cdots,N\\}$ be a subset of
indices. For a square matrix ${\mathbf{A}}\in{\mathbb{C}}^{N\times N}$, let
$\mathbf{A}_{{\mathcal{S}},{\mathcal{S}}}$ denote the submatrix by choosing
both the rows and columns of $\mathbf{A}$ indexed by $\mathcal{S}$. Let
${(\cdot)}^{*}_{\mathcal{S}}$, ${(\cdot)}^{\rm T}_{\mathcal{S}}$ and
${(\cdot)}^{\rm H}_{\mathcal{S}}$ be the conjugate, transpose and Hermitian
transpose operator of ${(\cdot)}_{\mathcal{S}}$, respectively. For two vectors
${\mathbf{a}}\in{\mathbb{C}}^{n}$ and ${\mathbf{b}}\in{\mathbb{C}}^{n}$,
${\mathbf{a}}/{\mathbf{b}}$ denotes the componentwise division.

##  II System model

This section establishes the OFDM model similar to [4], and the assumption
that all paths have a similar Doppler scaling factor is adopted, which seems
to be justified as long as the dominant Doppler effect is caused by the direct
transmitter/receiver motion. For completeness, the detailed derivations are
presented.

Consider the cyclic prefix (CP) OFDM system, and let $T$ and $T_{\rm cp}$
denote the OFDM symbol duration and the CP duration. The number of subcarriers
is $N$ and the frequency spacing is $\Delta f=1/T$. With carrier frequency
$f_{c}$, the frequency of the $k$th subcarrier is

$\displaystyle f_{k}=f_{c}+k\Delta f,~{}k=-N/2,\dots,N/2-1,$

and the bandwidth is $B=N\Delta f$.

Let $d[k]$ denote the information symbol to be transmitted on the $k$th
subcarrier. The transmitted continuous time OFDM block in passband is

$\displaystyle
s(t)=\frac{2}{\sqrt{N}}\Re\left\\{\sum\limits_{k=-N/2}^{N/2-1}d[k]{\rm
e}^{{\rm j}2\pi f_{k}t}g(t)\right\\},$

where $g(t)$ is

$\displaystyle g(t)=\begin{cases}1,\quad t\in[-T_{\rm cp},T],\\\ 0,\quad{\rm
otherwise}.\end{cases}$

A multipath underwater channel where all paths have a similar Doppler scaling
factor $a$ is considered, and the channel impulse response is

$\displaystyle
c(\tau,t)=\sum\limits_{p=1}^{L}A_{p}\delta(\tau-\tau_{p}(t))=\sum\limits_{p=1}^{L}A_{p}\delta(\tau-(\tau_{p}-at)),$

where $A_{p}$ and $\tau_{p}$ denote the gain and delay of the $p$th path, $L$
is the total number of paths. The received continuous time signal in passband
$\tilde{y}(t)$ is

$\displaystyle\tilde{y}(t)=s(t)\ast c(\tau,t)$
$\displaystyle=\sum\limits_{p=1}^{L}A_{p}s(t-\tau_{p}+at)$
$\displaystyle=\frac{2}{\sqrt{N}}\Re\left\\{\sum\limits_{k=-N/2}^{N/2-1}d[k]\sum\limits_{p=1}^{L}A_{p}{\rm
e}^{{\rm j}2\pi f_{k}(t-(\tau_{p}-at))}g((1+a)t-\tau_{p})\right\\}.$

After downconverting and low pass filtering, the baseband signal ${y}(t)$ can
be obtained as

$\displaystyle{y}(t)$ $\displaystyle={\rm LPF}\left[\tilde{y}(t){\rm e}^{-{\rm
j}2\pi f_{c}t}\right]$
$\displaystyle\approx\frac{1}{\sqrt{K}}\sum\limits_{k=-N/2}^{N/2-1}d[k]{\rm
e}^{{\rm j}2\pi k\Delta ft}{\rm e}^{{\rm j}2\pi
af_{k}t}\left[\sum\limits_{p=1}^{L}A_{p}{\rm e}^{-{\rm j}2\pi
f_{k}\tau_{p}}g((1+a)t-\tau_{p})\right]+n(t),$

where $n(t)$ is the additive noise in baseband. Note that each subcarrier
experiences a Doppler-induced frequency shift ${\rm e}^{{\rm j}2\pi af_{k}t}$.
Provided that the bandwidth of the OFDM signal is comparable to the center
frequency, the Doppler-induced frequency shift differs considerably on
different subcarriers, which introduces the intercarrier interference.

In order to mitigate the Doppler effect and convert the wideband OFDM system
into a narrowband system, a resampling approach is carried out first [4].
Suppose the Doppler factor $a$ is estimated as $\hat{a}$. Resample ${y}(t){\rm
e}^{-{\rm j}2\pi\hat{a}f_{c}t}$ with a resampling factor $\hat{a}$ to yield

$\displaystyle{z^{\prime}}(t)$ $\displaystyle={y}(\frac{t}{1+\hat{a}}){\rm
e}^{-{\rm j}2\pi\hat{a}f_{c}\frac{t}{1+\hat{a}}}$
$\displaystyle=\frac{1}{\sqrt{N}}\sum\limits_{k=-N/2}^{N/2-1}d[k]{\rm e}^{{\rm
j}2\pi\frac{1+a}{1+\hat{a}}k\Delta ft}{\rm e}^{{\rm
j}2\pi\frac{a-\hat{a}}{1+\hat{a}}f_{c}t}\left[\sum\limits_{p=1}^{L}A_{p}{\rm
e}^{-{\rm j}2\pi
f_{k}\tau_{p}}g(\frac{1+a}{1+\hat{a}}t-\tau_{p})\right]+n(t){\rm e}^{{\rm
j}2\pi\frac{\hat{a}}{1+\hat{a}}f_{c}t},$ (1)

The resampling factor $\hat{a}$ is chosen to make $\frac{1+a}{1+\hat{a}}$ as
close to one as possible. As a result, (II) can be approximated as

$\displaystyle{z^{\prime}}(t)\approx\frac{1}{\sqrt{N}}{\rm e}^{{\rm
j}2\pi\frac{a-\hat{a}}{1+\hat{a}}f_{c}t}\sum\limits_{k=-N/2}^{N/2-1}d[k]{\rm
e}^{{\rm j}2\pi k\Delta ft}\left[\sum\limits_{p=1}^{L}A_{p}{\rm e}^{-{\rm
j}2\pi f_{k}\tau_{p}}g\left(t-\tau_{p}\right)\right]+v(t).$ (2)

From (2), usually CFO estimation and compensation are carried out as the
second step to further remove the impact of Doppler [4]. Let ${r}(t)$ denote
the output of the CFO compensation, sampling $r(t)$ at $(n-1)T_{s}=(n-1)/B$,
$n=1,\dots,N$ yields 111In engineering application, fftshift operator ${\rm
fftshift}(\cdot)$ is usually introduced and model (3) can also be formulated
as ${\mathbf{r}}={\mathbf{e}}(\omega)\odot\left({\mathbf{F}}^{\rm H}{\rm
fftshift}({\mathbf{h}}\odot{\mathbf{d}})\right)+{\mathbf{w}}$.

$\displaystyle{\rm e}^{{\rm
j}\pi(n-1)}{\mathbf{r}}={\mathbf{e}}(\omega)\odot\left({\mathbf{F}}^{\rm
H}({\mathbf{h}}\odot{\mathbf{d}})\right)+{\mathbf{w}},$ (3)

where ${r}_{n}=r(t)|_{t=(n-1)T_{s}}$; the $n$th entry $d_{n}$ of
${\mathbf{d}}$ is $d_{n}=d[n-1-N/2]$; ${\mathbf{F}}$ denotes the normalized
DFT matrix with the $(m,n)$th entry being $\frac{1}{\sqrt{N}}{\rm e}^{-{\rm
j}2\pi(n-1)(m-1)/N}$; $\omega$ denotes the normalized residual CFO which is
usually very small and close to $0$, ${\mathbf{e}}(\omega)=[1,{\rm e}^{{\rm
j}\omega},\cdots,{\rm e}^{{\rm j}(N-1)\omega}]^{\rm T}$; $\mathbf{h}$ is the
frequency-domain channel response vector and can be viewed as a line spectra,
which can be represented as

$\displaystyle{\mathbf{h}}={\mathbf{A}}({\bm{\theta}}){\bm{\beta}},$ (4)

where ${{\bm{\beta}}\in{\mathbb{C}}^{L}}$ denotes the complex coefficient
vector related to the path gains and delays of the channel and

$\displaystyle\beta_{p}=A_{p}{\rm e}^{-{\rm j}2\pi f_{-\frac{N}{2}}\tau_{p}},$

${{\bm{\theta}}\in{\mathbb{R}}^{L}}$ denotes the frequency component related
to the path time delay of the channel and

$\displaystyle\theta_{p}=-2\pi\Delta f\tau_{p},$

$L$ denotes the number of paths, ${\mathbf{A}}({\bm{\theta}})$ is

$\displaystyle{\mathbf{A}}({\bm{\theta}})=[{\mathbf{a}}(\theta_{1}),{\mathbf{a}}(\theta_{2}),\cdots,{\mathbf{a}}(\theta_{L})],$

and ${\mathbf{a}}(\theta)=[1,{\rm e}^{{\rm j}\theta},\cdots,{\rm e}^{{\rm
j}(N-1)\theta}]^{\rm T}$; $\mathbf{w}$ denotes the noise and approximation
error which is assumed to follow the Gaussian distribution, i.e.,
${\mathbf{w}}\sim{\mathcal{CN}}({\mathbf{0}},\sigma^{2}{\mathbf{I}}_{N})$ with
$\sigma^{2}$ being the unknown variance.

Define

$\displaystyle{y}_{n}={\rm e}^{{\rm j}\pi(n-1)}{z}_{n}.$

Then, (3) is simplified as

$\displaystyle{\mathbf{y}}={\mathbf{e}}(\omega)\odot\left({\mathbf{F}}^{\rm
H}({\mathbf{h}}\odot{\mathbf{d}})\right)+{\mathbf{w}}.$ (5)

Assume that the symbol vector $\mathbf{d}$ are partitioned into three parts
corresponding to the pilot, null, and data, which are represented as
${\mathbf{d}}_{{\mathcal{P}}}$, ${\mathbf{d}}_{{\mathcal{N}}}$ and
${\mathbf{d}}_{{\mathcal{D}}}$ and
${\mathcal{P}}\cup{\mathcal{N}}\cup{\mathcal{D}}=\\{1,2,\cdots,N\\}$. Then the
goal is to jointly estimate the data symbol ${\mathbf{d}}_{{\mathcal{D}}}$,
the channel vector $\mathbf{h}$, the normalized residual CFO $\omega$, and the
nuisance parameters such as $\sigma^{2}$ by fully exploiting the channel
structure $\mathbf{h}$. In the following, we design an approximate Bayesian
algorithm.

##  III Probabilistic Formulation

This section describes the probabilistic formulation of the channel
$\mathbf{h}$, the data symbol $\mathbf{d}$ and the normalized residual CFO
$\omega$.

###  III-A Parametric Channel Model

For the channel $\mathbf{h}$, the number of path $L$ is usually unknown. Here
the number of paths is assumed to be $L_{\rm max}$ and $L<L_{\rm max}$ [29],
i.e.,

$\displaystyle{\mathbf{h}}=\sum\limits_{l=1}^{L_{\rm
max}}{\beta}_{l}{\mathbf{a}}({\theta}_{l})\triangleq{\mathbf{A}}({\bm{\theta}}){\bm{\beta}},$
(6)

where
${\mathbf{A}}({\bm{\theta}})=[{\mathbf{a}}({\theta}_{1}),\cdots,{\mathbf{a}}({\theta}_{L_{\rm
max}})]$. Since the number of path is supposed to be $L_{\rm max}$, the binary
hidden variables ${\mathbf{s}}=[s_{1},...,s_{L_{\rm max}}]^{\rm T}$ are
introduced, where $s_{l}=1$ means that the $l$th path is active, otherwise
deactive ($\beta_{l}=0$). Let $\rho$ denote the prior probability of each path
being active, i.e.,

$\displaystyle p(s_{l})=\rho^{s_{l}}(1-\rho)^{(1-s_{l})},\quad
s_{l}\in\\{0,1\\}.$ (7)

Then the probability mass function (PMF) of ${\mathbf{s}}$ is
$p({\mathbf{s}})=\prod\limits_{l=1}^{L_{\rm max}}p(s_{l})$. Given that
$s_{l}=1$, we assume that ${\beta}_{l}\sim{\mathcal{CN}}({\beta}_{l};0,\nu)$,
where $\nu$ denotes the priori variance of the path being active. Thus
$(s_{l},{\beta}_{l})$ follows a Bernoulli-Gaussian distribution, that is

$\displaystyle
p({\beta}_{l}|s_{l};\nu)=(1-s_{l}){\delta}({\beta}_{l})+s_{l}{\mathcal{CN}}({\beta}_{l};0,\nu).$
(8)

In addition, the distribution of $\bm{\beta}$ conditioned on $\mathbf{s}$ is
$p({\bm{\beta}}|{\mathbf{s}})=\prod\limits_{l=1}^{L_{\rm
max}}p({\beta}_{l}|s_{l};\nu)$.

From (7) and (8), it can be seen that the parameter $\rho$ denotes the
probability of the $l$th component being active and $\nu$ is a variance
parameter. The variable ${\bm{\theta}}=[\theta_{1},...,\theta_{L_{\rm
max}}]^{\rm T}$ has the prior PDF
$p({\bm{\theta}})=\begin{matrix}\prod_{l=1}^{L_{\rm
max}}p(\theta_{l})\end{matrix}$. Generally, $p(\theta_{l})$ is encoded through
the von Mises distribution [30, p. 36]

$\displaystyle
p(\theta_{l})={\mathcal{VM}}(\theta_{l};\mu_{0,l},\kappa_{0,l})=\frac{1}{2\pi{I_{0}}(\kappa_{0,l})}{\rm
e}^{\kappa_{0,l}{\cos(\theta-\mu_{0,l})}},$ (9)

where $\mu_{0,l}$ and $\kappa_{0,l}$ are the mean direction and concentration
parameters of the prior of the $l$th frequency $\theta_{l}$, $I_{p}(\cdot)$ is
the modified Bessel function of the first kind and the order $p$ [30, p. 348].
For a given von Mises distribution, ${\rm arg}({\rm
E}_{{\mathcal{VM}}(\theta;\mu,\kappa)}[{\rm e}^{{\rm j}\theta}])={\rm
arg}\left({\rm e}^{{\rm
j}\mu}\frac{I_{1}(\kappa)}{I_{0}(\kappa)}\right)=\mu={\rm
E}_{{{\mathcal{VM}}(\theta;\mu,\kappa)}}[\theta]$. In addition, ${\rm E}[{\rm
e}^{{\rm j}m\theta}]={\rm e}^{{\rm j}m\mu}(I_{m}(\kappa)/I_{0}(\kappa))$ [30,
pp. 26]. Without any knowledge of the frequency $\theta_{l}$, the
uninformative prior distribution $p(\theta_{l})={1}/({2\pi})$ is used. For
more details please refer to [29].

###  III-B OFDM System

A single-input single-output OFDM system with $N$ subcarriers is considered.
Since the structure of the consecutive OFDM symbols is not exploited, a single
OFDM system is modeled 222As you see, the proposed approach outputs the PDFs
of the amplitudes and delays of all the paths, which is very suitable for
channel tracking and will be left for future work.. Of the $N$ subcarriers,
$N_{\mathcal{P}}$ and $N_{\mathcal{N}}$ are dedicated as pilots and nulls,
respectively, and the remaining
$N_{\mathcal{D}}=N-N_{\mathcal{P}}-N_{\mathcal{N}}$ are used to transmit the
coded/interleaved data bits. The sets $\mathcal{N}$, $\mathcal{P}$,
$\mathcal{D}$ denote the indices of the null, pilot and data subcarriers,
respectively. It follows that
${\mathcal{P}}\cup{\mathcal{N}}\cup{\mathcal{D}}=\\{1,2,\cdots,N\\}$ and
${\mathcal{P}}\cap{\mathcal{N}}=\emptyset$,
${\mathcal{P}}\cap{\mathcal{D}}=\emptyset$ and
${\mathcal{N}}\cap{\mathcal{D}}=\emptyset$. Besides, the null subcarriers are
not necessary to delicate and can be set as the empty set
${\mathcal{N}}=\emptyset$ to improve the data transmission rate.

Let the vector ${\mathbf{u}}\in\\{0,1\\}^{K}$ be the transmitted (equi-
probable) information bits. These bits are coded by a rate-$R$ encoder and
interleaved to get the length-$K/R$ vector $c={\mathcal{C}}({\mathbf{u}})$,
where ${\mathcal{C}}:\\{0,1\\}^{K}\rightarrow\\{0,1\\}^{K/R}$ is the
interleaving and coding function such as a turbo code and a LDPC code. The
coded/interleaved bits are partioned into subvectors
${\mathbf{c}}^{(i)}\in\\{0,1\\}^{Q}$, $i\in{\mathcal{D}}$, which are then
mapped to the $i$th subcarrier. Let
${\mathcal{M}}:\\{0,1\\}^{Q}\rightarrow{\mathbb{A}}_{D}\subset{\mathbb{C}}$
denote the $2^{Q}$-ary mapping, where ${\mathbb{A}}_{D}$ is the symbol
alphabet. The pilots are drawn from the pilot symbol alphabet
${\mathbb{A}}_{P}\subset{\mathbb{C}}$. In OFDM, ${\mathbb{A}}_{D}$ is
typically a $2^{Q}$-ary quadrature amplitude modulation (QAM) alphabet and
${\mathbb{A}}_{P}$ is a quadrature phase shift keying (QPSK) alphabet. The
nulls, pilots and data symbols are stacked in vector $\mathbf{d}$. Vectors
${\mathbf{d}}_{\mathcal{N}}$, ${\mathbf{d}}_{\mathcal{P}}$ and
${\mathbf{d}}_{\mathcal{D}}$ contain the nulls, pilot symbol and data symbols,
respectively.

For the OFDM symbol, the information bits ${\mathbf{u}}$ has the following PDF

$\displaystyle
p({\mathbf{u}})=\prod\limits_{k\in{\mathcal{K}}}p(u_{k})=\prod\limits_{k\in{\mathcal{K}}}\left(0.5\times{1}_{[u_{k}\in\\{0,1\\}]}\right),$

where ${1}_{[\cdot]}$ denotes the indicator function, ${\mathcal{K}}$ denotes
the index set of the information bits. For the interleaving and coding
function, the PMFs $p({\mathbf{c}}|{\mathbf{u}})$ of subvectors ${\mathbf{c}}$
conditioned on $\mathbf{u}$ are

$\displaystyle
p({\mathbf{c}}|{\mathbf{u}})={1}_{[{\mathbf{c}}=C({\mathbf{u}})]}.$

Note that $p({\mathbf{c}}|{\mathbf{u}})$ describe the structure of the channel
code and interleaver, whose details are not provided. Conditioned on
$\mathbf{c}$, the PMF of $\mathbf{d}$ is $p({\mathbf{d}}|{\mathbf{c}})$.

The normalized residual CFO $\omega$ is treated as an unknown deterministic
parameter and $\omega\in[-\pi,\pi]$.

Figure 1: Factor graph of the joint PDF (10) describing the parametric channel
and OFDM system. Here the circle denotes the variable (vector) node, and the
square denotes the factor node. The nulls are not plotted here for simplicity,
as it can only be used for noise variance estimation.

According to the above probabilistic model, the factor graph representation
describing the channel and OFDM system model is given in Fig. 1. The joint PDF
is

$\displaystyle
p({\mathbf{y}},{\mathbf{d}}_{{\mathcal{D}}},{\mathbf{c}},{\mathbf{u}},{\bm{\beta}},{\bm{\theta}},{\mathbf{s}},;\omega,\nu,\rho,\sigma^{2})=p({\mathbf{s}};\rho)p({\bm{\beta}}|{\mathbf{s}};\nu)p({\bm{\theta}})p({\mathbf{u}})p({\mathbf{c}}|{\mathbf{u}})p({\mathbf{d}}|{\mathbf{c}})p({\mathbf{y}}|{\bm{\theta}},{\bm{\beta}},{\mathbf{d}}_{{\mathcal{D}}};\omega,\sigma^{2}).$
(10)

Firstly, the maximum likelihood estimation of the nuisance parameters are

$\displaystyle(\omega_{\rm{ML}},\nu_{\rm{ML}},\rho_{\rm{ML}},\sigma_{\rm{ML}}^{2})$
$\displaystyle=\underset{\omega,\nu,\rho,\sigma^{2}}{\operatorname{argmax}}\quad
p({\mathbf{y}};\omega,\nu,\rho,\sigma^{2})$ $\displaystyle=\int
p({\mathbf{y}},{\mathbf{d}}_{{\mathcal{D}}},{\mathbf{c}},{\mathbf{u}},{\bm{\beta}},{\bm{\theta}},{\mathbf{s}},;\omega,\nu,\rho,\sigma^{2}){\rm
d}{\mathbf{d}}_{{\mathcal{D}}}{\rm d}{\bm{\beta}}{\rm d}{\mathbf{u}}{\rm
d}{\mathbf{c}}{\rm d}{\bm{\theta}}{\rm d}{\mathbf{s}}.$ (11)

Then the BER optimal receiver, i.e., maximum a posterior (MAP) estimate is

$\displaystyle\hat{u}_{k}=\underset{u_{k}\in\\{0,1\\}}{\operatorname{argmax}}~{}p(u_{k}|{\mathbf{y}};\omega_{\rm{ML}},\nu_{\rm{ML}},\rho_{\rm{ML}},\sigma_{\rm{ML}}^{2}),$
(12)

where
$p(u_{k}|{\mathbf{y}};\omega_{\rm{ML}},\nu_{\rm{ML}},\rho_{\rm{ML}},\sigma_{\rm{ML}}^{2})\propto
p(u_{k},{\mathbf{y}};\omega_{\rm{ML}},\nu_{\rm{ML}},\rho_{\rm{ML}},\sigma_{\rm{ML}}^{2})$
can be obtained via marginalizing all the random variables but $u_{k}$ in the
joint PDF
$p({\mathbf{y}},{\mathbf{d}}_{{\mathcal{D}}},{\mathbf{c}},{\mathbf{u}},{\bm{\beta}},{\bm{\theta}},{\mathbf{s}},;\omega_{\rm{ML}},\nu_{\rm{ML}},\rho_{\rm{ML}},\sigma_{\rm{ML}}^{2})$.
Note that solving either (III-B) or (12) is computationally intractable.
Therefore, we resort to approximate Bayesian methods in the ensuing section.

##  IV Algorithm

This section utilizes expectation propagation (EP) [34] to develop the JCCD-
VALSE algorithm through novelly combining the gridless VALSE, the MMSE module
and the LDPC decoder module.

Figure 2: Factor graph of the joint PDF (10) with an additional hidden
variable node $\mathbf{h}$ and the module of the JCCD-VALSE algorithm. This
factor graph is equivalent to the factor graph shown in Fig. 1.

It is worth noting that in conventional underwater OFDM systems, null
subcarriers are beneficial for both Doppler and noise estimation. While in
this work, the pilot and data can be used for both Doppler and noise
estimation. As a consequence, the null subcarriers can be replaced with the
data subcarriers, which improves the spectrum efficiency. Define the ordered
set ${\mathcal{M}}={\mathcal{P}}\cup{\mathcal{D}}$.

Before designing the algorithm, an important point is emphasized. Note that

$\displaystyle p({\mathbf{y}}|{\mathbf{h}},{\mathbf{d}};\omega,\sigma^{2})$
$\displaystyle=\frac{1}{(\pi\sigma^{2})^{N}}{\rm
e}^{-\frac{\left({\mathbf{y}}-{\mathbf{e}}(\omega)\odot\left({\mathbf{F}}^{\rm
H}({\mathbf{h}}\odot{\mathbf{d}})\right)\right)^{\rm
H}\left({\mathbf{y}}-{\mathbf{e}}(\omega)\odot\left({\mathbf{F}}^{\rm
H}({\mathbf{h}}\odot{\mathbf{d}})\right)\right)}{\sigma^{2}}}$
$\displaystyle=\frac{1}{(\pi\sigma^{2})^{N}}{\rm
e}^{-\frac{\|\tilde{\mathbf{y}}(\omega)-{\mathbf{h}}\odot{\mathbf{d}}\|_{2}^{2}}{\sigma^{2}}}$
$\displaystyle=\prod\limits_{i=1}^{N}p(\tilde{y}_{i}(\omega)|{h}_{i},{d}_{i};\omega,\sigma^{2}),$
(13)

where $\tilde{\mathbf{y}}(\omega)$ is

$\displaystyle\tilde{\mathbf{y}}(\omega)={\mathbf{F}}\left({\mathbf{y}}\odot{\mathbf{e}}(-\omega)\right).$
(14)

This means that

$\displaystyle\tilde{\mathbf{y}}(\omega)={\mathbf{h}}\odot{\mathbf{d}}+\tilde{\mathbf{w}}$
(15)

is equivalent to the original observation model (5), where
$\tilde{\mathbf{w}}\sim{\mathcal{CN}}({\mathbf{0}},\sigma^{2}{\mathbf{I}})$.

It is beneficial to introduce an additional hidden variable $\mathbf{h}$
defined in (4), and an equivalent factor shown in Fig. 2 (a) can be obtained.
Compared to the original factor graph in Fig. 1, a delta factor node [28] is
introduced, which is the key to decompose the original problem into
subproblems to be solved by designing the respective modules. According to
Fig. 2 (b), the JCCD-VALSE can be summarized as follows: Firstly, initialize
the extrinsic message $m_{{\rm A}\rightarrow{\rm
B}}({\mathbf{h}}_{\mathcal{M}})$ from module A to module B as

$\displaystyle m_{{\rm A}\rightarrow{\rm
B}}({\mathbf{h}}_{{\mathcal{M}}})={\mathcal{CN}}({\mathbf{h}}_{{\mathcal{M}}};{\mathbf{h}}_{{\rm
A},{\mathcal{M}}}^{\rm ext},{\rm diag}({\mathbf{v}}_{{\rm
A},{\mathbf{h}}_{{\mathcal{M}}}}^{\rm ext})),$ (16)

and calculate the extrinsic message $m_{{\rm B}\rightarrow{\rm
C}}({\mathbf{d}}_{\mathcal{D}})$ from module B to module C. Secondly, running
the LDPC decoder algorithm in module C to obtain the posterior means and
variances of ${\mathbf{d}}_{\mathcal{D}}$. Thirdly, update the residual CFO
and noise variance estimate in module B. Fourthly, calculate the extrinsic
message $m_{{\rm B}\rightarrow{\rm A}}({\mathbf{h}}_{\mathcal{M}})$ from
module B to module A. Finally, calculate the extrinsic message $m_{{\rm
A}\rightarrow{\rm B}}({\mathbf{h}}_{\mathcal{M}})$ from module A to module B,
which closes the algorithm. The details of the above steps are described
below.

###  IV-A Calculate the Message $m_{{\rm B}\rightarrow{\rm
C}}({\mathbf{d}}_{\mathcal{D}})$

According to the belief propagation (BP), the message $m_{{\rm
B}\rightarrow{\rm C}}({\mathbf{d}}_{\mathcal{D}})$ is

$\displaystyle m_{{\rm B}\rightarrow{\rm C}}({\mathbf{d}}_{\mathcal{D}})$
$\displaystyle\propto\int m_{{\rm A}\rightarrow{\rm
B}}({\mathbf{h}})p({\mathbf{y}}|{\mathbf{h}},{\mathbf{d}}_{\mathcal{D}};\omega,\sigma^{2}){\rm
d}{\mathbf{h}}$
$\displaystyle=p({\mathbf{y}}|{\mathbf{d}}_{\mathcal{D}};\omega,\sigma^{2})=\prod\limits_{i\in{\mathcal{D}}}{\mathcal{CN}}(\tilde{y}_{i}(\omega);h_{{\rm
A},i}^{\rm ext}d_{i},|d_{i}|^{2}v_{{\rm A},{h}_{i}}^{\rm ext}).$ (17)

For each $i\in\mathcal{D}$, $d_{i}$ is drawn from a finite alphabet set, and
the PMF $m_{{\rm B}\rightarrow{\rm C}}({\mathbf{d}}_{\mathcal{D}})$ can be
obtained by evaluating the above Gaussian PDF at the points of the symbol
alphabet ${\mathbb{A}}_{D}$ followed by normalization.

###  IV-B Calculate the Posterior Means and Variances of
${\mathbf{d}}_{\mathcal{D}}$ in LDPC Module

The message $m_{{\rm B}\rightarrow{\rm C}}({\mathbf{d}}_{\mathcal{D}})$ is
input to the LDPC module, where the SPA is applied. After several iterations
where the change of the loglikelihood ratio of the bits is less than a
threshold or the number of iterations exceeds the setting maximum, the beliefs
of the data symbols are obtained. These beliefs are further used to calculate
the posterior means $d_{{\rm C},i}^{\rm post}$ and variances $v_{{\rm
C},d_{i}}^{\rm post}$ of the data.

###  IV-C Residual CFO and Noise Variance Estimation

The type II approach [31] and the EM algorithm are used to estimate the
normalized residual CFO $\omega$ and noise variance $\sigma^{2}$,
respectively. The type-II method is also termed as type-II maximum likelihood
(ML) method [32], which aims to to maximize the marginal likelihood in
Bayesian models. To estimate the residual CFO via type II approach, the
marginal likelihood $p(\tilde{\mathbf{y}}(\omega);\omega)$ is obtained. Define

$\displaystyle{\mathbf{z}}={\mathbf{h}}\odot{\mathbf{d}}.$ (18)

and model (19) can be formulated as

$\displaystyle\tilde{\mathbf{y}}(\omega)={\mathbf{z}}+\tilde{\mathbf{w}},$
(19)

where $\tilde{\mathbf{w}}\triangleq{\mathbf{F}}^{\rm
H}({\mathbf{w}}\odot{\mathbf{e}}(-\omega))\sim{\mathcal{CN}}({\mathbf{0}},\sigma^{2}{\mathbf{I}})$.
Straightforward calculation shows that the posterior means ${\mathbf{z}}_{\rm
B}^{\rm post}$ and variances ${\mathbf{v}}_{{\rm B},{\mathbf{z}}}^{\rm post}$
of $\mathbf{z}$ are

$\displaystyle{\mathbf{z}}_{\rm B}^{\rm post}$
$\displaystyle={\mathbf{h}}_{\rm A}^{\rm post}\odot{\mathbf{d}}_{\rm C}^{\rm
post},$ (20a) $\displaystyle{\mathbf{v}}_{{\rm B},{\mathbf{z}}}^{\rm post}$
$\displaystyle=|{\mathbf{d}}_{\rm C}^{\rm post}|^{2}\odot{\mathbf{v}}_{{\rm
A},{\mathbf{h}}}^{\rm post}+|{\mathbf{h}}_{\rm A}^{\rm
post}|^{2}\odot{\mathbf{v}}_{{\rm C},{\mathbf{d}}}^{\rm
post}+{\mathbf{v}}_{{\rm A},{\mathbf{h}}}^{\rm post}\odot{\mathbf{v}}_{{\rm
C},{\mathbf{d}}}^{\rm post},$ (20b)

where the posterior means and variances of $\mathbf{h}$ are obtained in module
A, i.e., the VALSE module, see Subsection IV-E (eq. (53a) and eq. (53b)).
Suppose that $\mathbf{z}$ follows Gaussian distribution with diagonal
covariance matrix and independent of $\tilde{\mathbf{w}}$, i.e.,
${\mathbf{z}}\sim{\mathcal{CN}}({\mathbf{z}};{\mathbf{z}}_{\rm B}^{\rm
post},{\rm diag}({\mathbf{v}}_{{\rm B},{\mathbf{z}}}^{\rm post}))$. For the
pseudo measurement $\tilde{\mathbf{y}}(\omega)$, it follows
$\tilde{\mathbf{y}}(\omega)\sim{\mathcal{CN}}(\tilde{\mathbf{y}}(\omega);{\mathbf{z}}_{\rm
B}^{\rm post},{\rm diag}({\mathbf{v}}_{{\rm B},{\mathbf{z}}}^{\rm
post})+\sigma^{2}{\mathbf{I}})$. The type II ML estimation problem can be
formulated as

$\displaystyle\hat{\omega}=\underset{\omega}{\operatorname{argmin}}~{}(\tilde{\mathbf{y}}(\omega)-{\mathbf{z}}_{\rm
B}^{\rm post})^{\rm H}({\rm diag}({\mathbf{v}}_{{\rm B},{\mathbf{z}}}^{\rm
post})+\sigma^{2}{\mathbf{I}})^{-1}(\tilde{\mathbf{y}}(\omega)-{\mathbf{z}}_{\rm
B}^{\rm post})\triangleq g(\omega).$ (21)

One can apply the Newston step to refine the previous estimate
$\hat{\omega}_{g}^{{}^{\prime}}$ as

$\displaystyle\hat{\omega}_{g}=\hat{\omega}_{g}^{{}^{\prime}}-\frac{\partial
g(\omega)}{\partial\omega}/\frac{\partial^{2}g(\omega)}{\partial\omega^{2}}\Big{|}_{\omega=\hat{\omega}_{g}^{{}^{\prime}}},$
(22)

where

$\displaystyle\frac{\partial
g(\omega)}{\partial\omega}=2\Re\left\\{\left(\frac{\partial\tilde{\mathbf{y}}(\omega)}{\partial\omega}\right)^{\rm
H}({\rm diag}({\mathbf{v}}_{{\rm B},{\mathbf{z}}}^{\rm
post})+\sigma^{2}{\mathbf{I}})^{-1}(\tilde{\mathbf{y}}(\omega)-{\mathbf{z}}_{\rm
B}^{\rm post})\right\\}$ (23)

and

$\displaystyle\frac{\partial^{2}g(\omega)}{\partial\omega^{2}}$
$\displaystyle=2\Re\left\\{\left(\frac{\partial^{2}\tilde{\mathbf{y}}(\omega)}{\partial\omega^{2}}\right)^{\rm
H}({\rm diag}({\mathbf{v}}_{{\rm B},{\mathbf{z}}}^{\rm
post})+\sigma^{2}{\mathbf{I}})^{-1}(\tilde{\mathbf{y}}(\omega)-{\mathbf{z}}_{\rm
B}^{\rm post})\right\\}$
$\displaystyle+2\left(\frac{\partial\tilde{\mathbf{y}}(\omega)}{\partial\omega}\right)^{\rm
H}({\rm diag}({\mathbf{v}}_{{\rm B},{\mathbf{z}}}^{\rm
post})+\sigma^{2}{\mathbf{I}})^{-1}\frac{\partial\tilde{\mathbf{y}}(\omega)}{\partial\omega},$
(24)

where

$\displaystyle\frac{\partial\tilde{\mathbf{y}}(\omega)}{\partial\omega}={\mathbf{F}}\left({\mathbf{y}}\odot{\bm{\chi}}\odot{\mathbf{e}}(-\omega)\right)$
(25)

and

$\displaystyle\frac{\partial^{2}\tilde{\mathbf{y}}(\omega)}{\partial\omega^{2}}={\mathbf{F}}\left({\mathbf{y}}\odot{\bm{\chi}}^{2}\odot{\mathbf{e}}(-\omega)\right),$
(26)

where ${\bm{\chi}}=-{\rm j}[0,1,\cdots,N-1]^{\rm T}$. It is worth noting that
the residual CFO estimation approach (21) is very general and generalize the
null subcarriers based approach in [4]. By letting ${\mathbf{v}}_{{\rm
B},z_{i}}^{\rm post}\rightarrow\infty$, $i\in{\mathcal{P}}\cup{\mathcal{D}}$,
i.e., the posterior variances of $\mathbf{z}$ corresponding to pilot and data
subcarriers tend to infinity, (21) reduces to the null subcarriers based
approach [4]. Since we utilize more data to estimate the residual CFO, it
makes sense that our approach estimates the residual CFO more accurately,
leading to better channel estimation and data decoding performance, as
illustrated in the numerical experiments.

Once the residual CFO $\omega$ is updated, the pseudo measurements
$\tilde{\mathbf{y}}(\omega)$ are refined according to (14).

The EM approach is adopted to estimate the noise variance, which is

$\displaystyle\hat{\sigma}^{2}=\underset{{\sigma}^{2}}{\operatorname{argmax}}~{}{\rm
E}_{{\mathbf{z}}}\left[\ln
p(\tilde{\mathbf{y}}(\omega);{\mathbf{z}},\sigma^{2})\right]=-\frac{{\rm
E}_{{\mathbf{z}}}\|\tilde{\mathbf{y}}(\omega)-{\mathbf{z}}\|_{2}^{2}}{\sigma^{2}}-N\ln\sigma^{2}+{\rm
const},$ (27)

where the expectation is taken with respect to the posterior PDF of
${\mathbf{z}}$.

Substituting (20) in (27), one obtains

$\displaystyle\hat{\sigma}^{2}=\underset{{\sigma}^{2}}{\operatorname{argmax}}~{}-\frac{\|\tilde{\mathbf{y}}(\omega)-{\mathbf{z}}_{\rm
B}^{\rm post}\|_{2}^{2}+{\mathbf{1}}^{\rm T}{\mathbf{v}}_{{\rm
B},{\mathbf{z}}}^{\rm post}}{\sigma^{2}}-N\ln\sigma^{2}.$ (28)

Setting the objective function of (28) with respect to $\sigma^{2}$ to zero,
the noise variance is estimated as

$\displaystyle\hat{\sigma}^{2}=\frac{\|\tilde{\mathbf{y}}(\hat{\omega}_{g})-{\mathbf{z}}_{\rm
B}^{\rm post}\|^{2}+{\mathbf{1}}^{\rm T}{\mathbf{v}}_{{\rm
B},{\mathbf{z}}}^{\rm post}}{N}.$ (29)

###  IV-D Calculate the Extrinsic Message $m_{{\rm B}\rightarrow{\rm
A}}({\mathbf{h}}_{\mathcal{M}})$ from Module B to Module A

According to variational message passing [33], the message $m_{{\rm
B}\rightarrow{\rm A}}({\mathbf{h}}_{\mathcal{M}})$ can be calculated as

$\displaystyle\ln m_{{\rm B}\rightarrow{\rm A}}({\mathbf{h}}_{\mathcal{M}})$
$\displaystyle={\rm E}_{\mathbf{d}}\left[\ln
p({\mathbf{y}}|{\mathbf{h}},{\mathbf{d}})\right]+{\rm const}$
$\displaystyle={\rm E}_{\mathbf{d}}\left[\ln
p(\tilde{\mathbf{y}}(\omega)|{\mathbf{h}},{\mathbf{d}})\right]+{\rm const}$
$\displaystyle={\rm
E}_{\mathbf{d}}\left[-\frac{|\tilde{\mathbf{y}}(\omega)-{\mathbf{h}}\odot{\mathbf{d}}|^{2}}{\sigma^{2}}\right]+{\rm
const},$ (30)

where ${\rm const}$ is to ensure that the PDF $m_{{\rm B}\rightarrow{\rm
A}}({\mathbf{h}}_{\mathcal{M}})$ is normalized and ${\rm
E}_{\mathbf{d}}[\cdot]$ denotes the expectation of $\mathbf{d}$ with respect
to the posterior PDF. It can be seen that $m_{{\rm B}\rightarrow{\rm
A}}({\mathbf{h}}_{\mathcal{M}})$ is Gaussian distributed.

From (IV-D), one has $m_{{\rm B}\rightarrow{\rm
A}}({\mathbf{h}}_{\mathcal{M}})=\prod\limits_{i\in{\mathcal{M}}}m_{{\rm
B}\rightarrow{\rm A}}({h}_{i})$. Now we calculate the message $m_{{\rm
B}\rightarrow{\rm A}}({h}_{i})$ corresponding to the pilots and data
subcarriers, respectively.

For $i\in{\mathcal{P}}$, straightforward calculation shows that

$\displaystyle m_{{\rm B}\rightarrow{\rm
A}}({h}_{i})={\mathcal{CN}}(h_{i};h_{{\rm B},i}^{\rm ext},v_{{\rm
B},h_{i}}^{\rm ext}),$ (31)

where

$\displaystyle h_{{\rm B},i}^{\rm ext}=\tilde{y}_{i}(\omega)/d_{i},\quad
v_{{\rm B},h_{i}}^{\rm ext}=\sigma^{2}/|{d}_{i}|^{2},i\in{\mathcal{P}}.$ (32)

For the data symbols where $i\in{\mathcal{D}}$,

$\displaystyle\ln m_{{\rm B}\rightarrow{\rm A}}({h}_{i})={\rm
E}_{d_{i}}\left[\ln
p(\tilde{\mathbf{y}}(\omega)|{\mathbf{h}},{\mathbf{d}})\right]+{\rm const}$
$\displaystyle={\rm
E}_{\mathbf{d}}\left[-\frac{|\tilde{y}_{i}(\omega)-{h_{i}}d_{i}|^{2}}{\sigma^{2}}\right]+{\rm
const}=-{\rm
E}[|d_{i}|^{2}]\frac{|h|_{i}^{2}}{\sigma^{2}}+2\frac{\Re\\{\tilde{y}_{i}^{*}(\omega)h_{i}{\rm
E}[d_{i}]\\}}{\sigma^{2}}+{\rm const}$ $\displaystyle=-\frac{(|d_{{\rm
C},i}^{\rm post}|^{2}+v_{{\rm C},d_{i}}^{\rm
post})|h|_{i}^{2}}{\sigma^{2}}+2\frac{\Re\\{\tilde{y}_{i}^{*}(\omega)h_{i}d_{{\rm
C},i}^{\rm post}\\}}{\sigma^{2}}+{\rm const}.$ (33)

As a consequence, one has

$\displaystyle h_{{\rm B},i}^{\rm ext}$
$\displaystyle=\frac{\tilde{y}_{i}(\omega)d_{i}^{*}}{|d_{{\rm C},i}^{\rm
post}|^{2}+v_{{\rm C},d_{i}}^{\rm post}},$ (34a) $\displaystyle v_{{\rm
B},h_{i}}^{\rm ext}$ $\displaystyle=\sigma^{2}/(|d_{{\rm C},i}^{\rm
post}|^{2}+v_{{\rm C},d_{i}}^{\rm post}),i\in{\mathcal{D}}.$ (34b)

Note that if $v_{{\rm C},d_{i}}^{\rm post}=0$, (34) is consistent with (32).

###  IV-E Calculate the message $m_{{\rm A}\rightarrow{\rm
B}}({\mathbf{h}}_{\mathcal{M}})$

Define ${\mathbf{A}}_{\mathcal{M}}$ as the matrix chosen from the rows of
$\mathbf{A}$ indexed by ${\mathcal{M}}$. The pseudo linear measurement model
under heterogenous noise is

$\displaystyle\widetilde{\mathbf{h}}_{\mathcal{M}}={\mathbf{A}}_{\mathcal{M}}({\bm{\theta}}){\bm{\beta}}+\widetilde{\bm{\epsilon}}_{\mathcal{M}},$
(35)

where $\widetilde{\mathbf{h}}_{\mathcal{M}}={\mathbf{h}}_{{\rm
B},{\mathcal{M}}}^{\rm ext}$,
$\widetilde{\bm{\epsilon}}_{\mathcal{M}}\sim{\mathcal{CN}}({\mathbf{0}},{\rm
diag}(\tilde{\bm{\sigma}}_{\mathcal{M}}^{2}))$ and
$\widetilde{\bm{\sigma}}_{\mathcal{M}}^{2}={\mathbf{v}}_{{\rm
B},{\mathbf{h}}_{\mathcal{M}}}^{\rm ext}$.

The VALSE algorithm tries to construct a structured PDF
$q({\bm{\theta}},{\bm{\beta}},{\mathbf{s}}|\widetilde{\mathbf{h}}_{\mathcal{M}})$
to approximate the true PDF
$p({\bm{\theta}},{\bm{\beta}},{\mathbf{s}}|\widetilde{\mathbf{h}}_{\mathcal{M}})$
by minimizing their Kullback-Leibler (KL) divergence
${\rm{KL}}\left(q({\bm{\theta}},{\bm{\beta}},{\mathbf{s}}|\widetilde{\mathbf{h}}_{\mathcal{M}})||p({\bm{\theta}},{\bm{\beta}},{\mathbf{s}}|\widetilde{\mathbf{h}}_{\mathcal{M}})\right)$.
For
$q({\bm{\theta}},{\bm{\beta}},{\mathbf{s}}|\widetilde{\mathbf{h}}_{\mathcal{M}})$,
it is supposed to be factored as

$\displaystyle
q({\bm{\theta}},{\bm{\beta}},{\mathbf{s}}|\widetilde{\mathbf{h}}_{\mathcal{M}})=\prod_{i=1}^{L_{\rm
max}}q(\theta_{i}|\widetilde{\mathbf{h}}_{\mathcal{M}})q({\bm{\beta}}|\widetilde{\mathbf{h}}_{\mathcal{M}},{\mathbf{s}})q({\mathbf{s}}|\widetilde{\mathbf{h}}_{\mathcal{M}}),$
(36)

where
$q({\mathbf{s}}|\widetilde{\mathbf{h}}_{\mathcal{M}})=\delta({\mathbf{s}}-\hat{\mathbf{s}})$
is restricted to a point estimate. Given that VALSE outputs the approximated
PDF
$q({\bm{\theta}},{\bm{\beta}},{\mathbf{s}}|\widetilde{\mathbf{h}}_{\mathcal{M}})$,
according to EP [34], the extrinsic message $m_{{\rm A}\rightarrow{\rm
B}}({\mathbf{h}}_{\mathcal{M}})={\mathcal{CN}}({\mathbf{h}}_{\mathcal{M}};{\mathbf{h}}_{{\rm
A},{\mathcal{M}}}^{\rm ext},{\rm diag}({\mathbf{v}}_{{\rm
A},{\mathbf{h}}_{\mathcal{M}}}^{\rm ext}))$ of ${\mathbf{h}}_{\mathcal{M}}$
from module A to module B is

$\displaystyle m_{{\rm A}\rightarrow{\rm B}}({\mathbf{h}}_{\mathcal{M}})$
$\displaystyle\propto\frac{{\rm Proj}[\int
q({\bm{\beta}}_{\hat{S}}|\widetilde{\mathbf{h}}_{\mathcal{M}})\delta({\mathbf{h}}_{\mathcal{M}}-{\mathbf{A}}_{{\mathcal{M}},\hat{S}}({\bm{\theta}}){\bm{\beta}}_{\hat{S}})q({\bm{\theta}}|\widetilde{\mathbf{h}}_{\mathcal{M}}){\rm
d}{\mathbf{w}}_{\hat{S}}{\rm d}{\bm{\theta}}]}{m_{{\rm B}\rightarrow{\rm
A}}({\mathbf{h}}_{\mathcal{M}})}$ $\displaystyle\triangleq\frac{{\rm
Proj}[q_{\rm A}({\mathbf{h}}_{\mathcal{M}})]}{m_{{\rm B}\rightarrow{\rm
A}}({\mathbf{h}}_{\mathcal{M}})},$ (37)

where ${\rm Proj}[q({\mathbf{x}})]$ denotes the projection operation which
approximates the target distribution $q({\mathbf{x}})$ as Gaussian
distribution with moment matching, i.e., the means and variances of
$q({\mathbf{x}})$ matches with that of the approximated Gaussian distribution.
Let ${\mathbf{h}}_{{\rm A},{\mathcal{M}}}^{\rm post}$ and ${\mathbf{v}}_{{\rm
B},{\mathbf{h}}_{\mathcal{M}}}^{\rm post}$ denote the posterior means and
variances of ${\mathbf{h}}_{\mathcal{M}}$ with respect to $q_{\rm
A}({\mathbf{h}}_{\mathcal{M}})$, then

$\displaystyle{\rm Proj}[q_{\rm
A}({\mathbf{h}}_{\mathcal{M}})]={\mathcal{CN}}({\mathbf{h}}_{\mathcal{M}};{\mathbf{h}}_{{\rm
A},{\mathcal{M}}}^{\rm post},{\rm diag}({\mathbf{v}}_{{\rm
A},{\mathbf{h}}_{\mathcal{M}}}^{\rm post})).$ (38)

From (IV-E), ${\mathbf{h}}_{{\rm A},{\mathcal{M}}}^{\rm ext}$ and
${\mathbf{v}}_{{\rm A},{\mathbf{h}}_{\mathcal{M}}}^{\rm ext}$ can be
calculated as

$\displaystyle{\mathbf{v}}_{{\rm A},{\mathbf{h}}_{\mathcal{M}}}^{\rm
ext}=\left(\frac{1}{{\mathbf{v}}_{{\rm A},{\mathbf{h}}_{\mathcal{M}}}^{\rm
post}}-\frac{1}{{\mathbf{v}}_{{\rm A},{\mathbf{h}}_{\mathcal{M}}}^{\rm
ext}}\right)^{-1},$ (39a) $\displaystyle{\mathbf{h}}_{{\rm
A},{\mathcal{M}}}^{\rm ext}={\mathbf{v}}_{{\rm
A},{\mathbf{h}}_{\mathcal{M}}}^{\rm ext}\odot\left(\frac{{\mathbf{h}}_{{\rm
A},{\mathcal{M}}}^{\rm post}}{{\mathbf{v}}_{{\rm
A},{\mathbf{h}}_{\mathcal{M}}}^{\rm post}}-\frac{{\mathbf{h}}_{{\rm
B},{\mathcal{M}}}^{\rm ext}}{{\mathbf{v}}_{{\rm
B},{\mathbf{h}}_{\mathcal{M}}}^{\rm ext}}\right),$ (39b)

where $\odot$ denotes the Hadarmard product. In the follow, the details of
obtaining the approximated posterior PDFs through VALSE is briefly described
333For further details, please refer to [29, 35]..

The VALSE proceeds as follows: Firstly, the frequencies
$\\{\theta_{i}\\},i\in\hat{S}$ are inferred and the von Mises approximations
$q(\theta_{i}|\widetilde{\mathbf{h}}_{\mathcal{M}}),i\in\hat{S}$ of the PDFs
are obtained. Secondly, the weights $\bm{\beta}$ and support $\mathbf{s}$ are
inferred and the posterior PDFs
$q({\bm{\beta}}|\widetilde{\mathbf{h}}_{\mathcal{M}})$ of the weight and the
PMF $q({\mathbf{s}}|\widetilde{\mathbf{h}}_{\mathcal{M}})$ of the support are
obtained. Finally, the model parameters $\rho$ and $\nu$ are estimated. In the
following, we detail the procedures.

####  IV-E1 Inferring the frequencies

Let $\mathcal{S}$ be the set of indices of the non-zero components of
$\mathbf{s}$, i.e.,

$\displaystyle\mathcal{S}=\\{i|1\leq i\leq L_{\rm max},s_{i}=1\\}.$

Analogously, define $\widehat{\mathcal{S}}$ based on $\widehat{\mathbf{s}}$.
For $i\notin{\mathcal{S}}$, we have
$q(\theta_{i}|\widetilde{\mathbf{h}}_{\mathcal{M}})=p(\theta_{i})$. For
$i\in{\mathcal{S}}$, $q(\theta_{i}|\widetilde{\mathbf{h}}_{\mathcal{M}})$ is

$\displaystyle q(\theta_{i}|\widetilde{\mathbf{h}}_{\mathcal{M}})\propto
p(\theta_{i}){\rm exp}(\Re\\{{\bm{\eta}}_{i}^{\rm
H}{\mathbf{a}}(\theta_{i})\\})\triangleq\exp\\{f(\theta_{i})\\},$ (40)

where the complex vector $\bm{\eta}_{i}$ is given by

$\displaystyle{\bm{\eta}}_{i}=2\left[\left(\widetilde{\mathbf{h}}_{\mathcal{M}}-\sum_{l\in\widehat{\mathcal{S}}\backslash\\{i\\}}{\widehat{\mathbf{a}}}_{{\mathcal{M}},l}{\widehat{\beta}}_{l}\right){\widehat{\beta}}^{*}_{i}-\sum_{l\in\widehat{\mathcal{S}}\backslash\\{i\\}}{\widehat{\mathbf{a}}}_{{\mathcal{M}},l}{\widehat{C}}_{l,i}\right]/\tilde{\bm{\sigma}}_{\mathcal{M}}^{2},$
(41)

where ${\widehat{\mathbf{a}}}_{{\mathcal{M}},l}$ denotes the $l$th column of
${\widehat{\mathbf{A}}}_{{\mathcal{M}}}$, $``\backslash i"$ denotes the
indices $\widehat{\mathcal{S}}$ excluding $i$,
$\hat{\bm{\beta}}_{\widehat{\mathcal{S}}}$ denotes the subvector of
$\hat{\bm{\beta}}$ by choosing the $\widehat{\mathcal{S}}$ rows of
${\bm{\beta}}$. Since it is hard to obtain the analytical results of the
expected value of ${\mathbf{a}}_{\mathcal{M}}(\theta_{i})$ for the PDF (40),
$q(\theta_{i}|\tilde{\mathbf{h}})$ is approximated as a von Mises distribution
via the following two steps 444For further details, please refer to [29,
Algorithm 2: Heurestic 2].: Firstly, search for the most dominant component of
(40) via Heuristic 2 [29, Algorithm 2]. Secondly, use second-order Taylor
approximation to approximate the dominant component as a single von Mises PDF.
It is worth noting that, the number of subcarriers $N$ is usually very large,
for example $N=1024$, thus the peak of the single von Mises PDF should be well
matched to that of (40) in very high accuracy. Let $\bar{\theta}_{i}$ denote
the mean of the dominant component. In contrast with [29] where a single
Newston refinement is adopted, we apply a two step procedure to improve the
robustness. First, apply the Newston step to refine the estimate

$\displaystyle\hat{\mu}_{i}=\bar{\theta}_{i}-\frac{f^{{}^{\prime}}(\bar{\theta}_{i})}{f^{{}^{\prime\prime}}(\bar{\theta}_{i})},$
(42)

where $f(\theta_{i})$ is defined in (40). Since the Newston step may be too
large to across the peak, we set

$\displaystyle\bar{\mu}_{i}=\frac{\hat{\mu}_{i}+\bar{\theta}_{i}}{2}.$ (43)

Then we continue to apply the Newston step to obtain the final estimate as

$\displaystyle\hat{\theta}_{i}=\bar{\mu}_{i}-\frac{f^{{}^{\prime}}(\bar{\mu}_{i})}{f^{{}^{\prime\prime}}(\bar{\mu}_{i})}.$
(44)

Numerical experiments demonstrate that such a two step approach improves the
robustness of the VALSE significantly, compared to the original implementation
[29]. Besides, the concentration parameter $\hat{\kappa}_{i}$ is approximated
as

$\displaystyle\hat{\kappa}_{i}=A^{-1}\left({\rm
e}^{0.5/f^{{}^{\prime\prime}}(\hat{\theta}_{i})}\right),$ (45)

where $A^{-1}(\cdot)$ is the inverse function of
$A(\cdot)=I_{1}(\cdot)/I_{0}(\cdot)$. The expected value of
${\mathbf{a}}_{\mathcal{M}}(\theta_{i})$ is

$\displaystyle\hat{\mathbf{a}}_{\mathcal{M}}(\theta_{i})=\frac{I_{\mathcal{M}}(\hat{\kappa}_{i})}{I_{0}(\hat{\kappa}_{i})}\odot{\mathbf{a}}_{\mathcal{M}}(\hat{\theta}_{i}).$
(46)

It is found that calculating the most dominate component via Heurestic 2 [29,
Algorithm 2] in each iteration is very time consuming. Consequently, we adopt
the warm start strategy where we only implement Heuristic 2 in the
initialization stage. Then we use the previous estimates to perform
refinement, which is very efficiently.

####  IV-E2 Inferring the weights and support

Next by fixing the posterior PDFs
$q(\theta_{i}|\widetilde{\mathbf{h}}_{\mathcal{M}}),i=1,...,L_{\rm max}$, the
posterior PDF $q({\bm{\beta}},{\mathbf{s}}|\tilde{\mathbf{h}}_{\mathcal{M}})$
is obtained. Define the matrices $\mathbf{J}$ and $\mathbf{u}$ as

$\displaystyle{J}_{ij}=\begin{cases}{\mathbf{1}}_{\mathcal{M}}^{\rm
T}\left(\frac{1}{\tilde{\bm{\sigma}}_{\mathcal{M}}^{2}}\right),&i=j\\\
{\widehat{\mathbf{a}}}^{\rm H}_{{\mathcal{M}},i}{\rm
diag}\left(\frac{1}{\tilde{\bm{\sigma}}_{\mathcal{M}}^{2}}\right){\widehat{\mathbf{a}}}_{{\mathcal{M}},j},&i\neq{j}\end{cases},\quad
i,j\in\\{1,2,\cdots,L_{\rm max}\\},$ (47a)
$\displaystyle{\mathbf{u}}=\widehat{\mathbf{A}}_{\mathcal{M}}^{\rm
H}\left(\tilde{\mathbf{h}}_{\mathcal{M}}/\tilde{\bm{\sigma}}_{\mathcal{M}}^{2}\right).$
(47b)

According to (36), the posterior approximation
$q({\bm{\beta}},{\mathbf{s}}|\tilde{\mathbf{h}}_{\mathcal{M}})$ is factored as
the product of $q({\bm{\beta}}|\tilde{\mathbf{h}}_{\mathcal{M}},{\mathbf{s}})$
and $\delta({\mathbf{s}}-{\widehat{\mathbf{s}}})$. For a given
$\widehat{\mathbf{s}}$,
$q({\bm{\beta}}_{\widehat{\mathcal{S}}}|\tilde{\mathbf{h}}_{\mathcal{M}})$ is
a complex Gaussian distribution, and
$q({\bm{\beta}}|\tilde{\mathbf{h}}_{\mathcal{M}};\widehat{\mathbf{s}})$ is

$\displaystyle
q({\bm{\beta}}|\tilde{\mathbf{h}};\widehat{\mathbf{s}})={\mathcal{CN}}({\bm{\beta}}_{\widehat{\mathcal{S}}};\widehat{\bm{\beta}}_{\widehat{\mathcal{S}}},\widehat{\mathbf{C}}_{\widehat{\mathcal{S}}})\prod_{i\not\in\widehat{\mathcal{S}}}\delta(\beta_{i}),$
(48)

where

$\displaystyle\widehat{\mathbf{C}}_{{\mathcal{S}}}=\left({{\mathbf{J}}_{{\mathcal{S}}}}+\frac{{\mathbf{I}}_{|{\mathcal{S}}|}}{\nu}\right)^{-1},$
(49a)
$\displaystyle\widehat{\bm{\beta}}_{{\mathcal{S}}}=\widehat{\mathbf{C}}_{{\mathcal{S}}}{\mathbf{u}}_{{\mathcal{S}}}.$
(49b)

Although (49a) involves a matrix inversion, it can be avoided by utilizing the
rank one update formula.

Then we need to find $\hat{\mathbf{s}}$ which maximizes $\ln Z(\mathbf{s})$,
i.e.,

$\displaystyle\hat{\mathbf{s}}=\underset{\mathbf{s}}{\operatorname{argmax}}\ln
Z(\mathbf{s}),$ (50)

where $\ln Z(\mathbf{s})$ is

$\displaystyle\ln Z(\mathbf{s})\triangleq$
$\displaystyle-\ln\det\left({\mathbf{J}}_{\mathcal{S}}+\frac{1}{\nu}{\mathbf{I}}_{|\mathcal{S}|}\right)+{\mathbf{u}}_{\mathcal{S}}^{\rm
H}\left(\mathbf{J}_{\mathcal{S}}+\frac{1}{{\nu}}\mathbf{I}_{|\mathcal{S}|}\right)^{-1}{\mathbf{u}}_{\mathcal{S}}$
$\displaystyle+||\mathbf{s}||_{0}\ln\frac{\rho}{1-\rho}+||\mathbf{s}||_{0}\ln\frac{1}{{\nu}}+{\rm
const}.$ (51)

Problem (50) is solved via a greedy algorithm, and local optimum is
guaranteed.

####  IV-E3 Inferring the model parameters

After updating the PDFs of the frequencies and weights, the model parameters
$\\{\rho,~{}\nu\\}$ are updated as

$\displaystyle\widehat{\rho}=$
$\displaystyle\frac{||\widehat{\mathbf{s}}||_{0}}{{L_{\rm max}}},$
$\displaystyle\widehat{{\nu}}=$ $\displaystyle\frac{\widehat{\bm{\beta}}^{\rm
H}_{\widehat{\mathcal{S}}}\widehat{\bm{\beta}}_{\widehat{\mathcal{S}}}+{\rm
tr}(\widehat{{\mathbf{C}}}_{\widehat{\mathcal{S}}})}{||\widehat{\mathbf{s}}||_{0}}.$
(52)

Now the the posterior mean and variances of ${\mathbf{h}}$ can be obtained as

$\displaystyle{\mathbf{h}}_{{\rm A},{\mathcal{M}}}^{\rm
post}=\hat{\mathbf{A}}_{{\mathcal{M}},\hat{\mathcal{S}}}\hat{\bm{\beta}}_{\hat{\mathcal{S}}},$
(53a) $\displaystyle{\mathbf{v}}_{{\rm A},{\mathbf{h}}_{\mathcal{M}}}^{\rm
post}={\rm
diag}(\hat{\mathbf{A}}_{{\mathcal{M}},\hat{\mathcal{S}}}\hat{\mathbf{C}}_{\hat{\mathcal{S}}}\hat{\mathbf{A}}_{{\mathcal{M}},\hat{\mathcal{S}}}^{\rm
H})+\left(\hat{\bm{\beta}}_{\hat{\mathcal{S}}}^{\rm
H}\hat{\bm{\beta}}_{\hat{\mathcal{S}}}{\mathbf{1}}_{|{\mathcal{M}}|}-|\hat{\mathbf{A}}_{{\mathcal{M}},\hat{\mathcal{S}}}|^{2}|\hat{\bm{\beta}}_{\hat{\mathcal{S}}}|^{2}\right)$
$\displaystyle+\left[{\rm
tr}(\hat{\mathbf{C}}_{\hat{\mathcal{S}}}){\mathbf{1}}_{|{\mathcal{M}}|}-|\hat{\mathbf{A}}_{{\mathcal{M}},\hat{\mathcal{S}}}|^{2}{\rm
diag}(\hat{\mathbf{C}}_{\hat{\mathcal{S}}})\right],$ (53b)

respectively. Then, the extrinsic message $m_{{\rm A}\rightarrow{\rm
B}}({\mathbf{h}}_{\mathcal{M}})={\mathcal{CN}}({\mathbf{h}}_{\mathcal{M}};{\mathbf{h}}_{{\rm
A},{\mathcal{M}}}^{\rm ext},{\rm diag}({\mathbf{v}}_{{\rm
A},{\mathbf{h}}_{\mathcal{M}}}^{\rm ext}))$ from module A to module B can be
obtained from (39b) and (39a). Remember that we have initialized the extrinsic
message $m_{{\rm A}\rightarrow{\rm B}}({\mathbf{h}}_{\mathcal{M}})$ from
module A to module B at the beginning of Section IV, the algorithm is closed.

###  IV-F Initialization

Performing the joint estimation of the channel, residual CFO and data
detection is a very challenging task, and initialization is a very important
step to have good performance. Since the residual CFO $\omega$ is close to
zero, we initialize it as zero. To obtain a good estimate of the channel, we
run the pilot based VALSE algorithm to obtain the posterior means and
variances of the channel, which are provided as the initialization of the
extrinsic message $m_{{\rm A}\rightarrow{\rm B}}({\mathbf{h}}_{\mathcal{M}})$
from module A to module B. Besides, we also initialize the nuisance parameters
$\sigma^{2}$ and ${\nu}$ via the pilot based VALSE algorithm. The whole
algorithm is summarized as Algorithm 1.

Algorithm 1 JCCD-VALSE algorithm

1: Initialize the extrinsic message $m_{{\rm A}\rightarrow{\rm
B}}({\mathbf{h}}_{\mathcal{M}})$ from module A to module B, the nuisance
parameters $\sigma^{2}$ and $\nu$ via the pilot based VALSE algorithm, as
described in Subsection IV-F; Set the number of iterations $T_{\rm max}$, the
maximum number of paths $L_{\rm max}$;

2: Calculate the message $m_{{\rm B}\rightarrow{\rm
C}}({\mathbf{d}}_{\mathcal{D}})$ (IV-A) and input it to the LDPC decoder
module to obtain the posterior means $d_{{\rm C},i}^{\rm post}$ and variances
$v_{{\rm C},d_{i}}^{\rm post}$ of the data, as described in Subsection IV-B.

3: Calculate the extrinsic message $m_{{\rm B}\rightarrow{\rm
A}}({\mathbf{h}}_{\mathcal{M}})$ through (32) and (IV-D) from module B to
module A.

4: Initialize $q(\theta_{i}|\widetilde{\mathbf{h}}_{\mathcal{M}})$,
$i=1,2,\cdots,L_{\rm max}$; compute
$\hat{\mathbf{a}}_{\mathcal{M}}(\theta_{i})$, $\mathbf{J}$ (47a) and
$\mathbf{u}$ (47b).

5: for $t=1,\cdots,T_{\rm max}$ do

6: Update $\widehat{\mathbf{s}}$ (50),
$\widehat{\bm{\beta}}_{\widehat{\mathcal{S}}}~{}{\rm
and}~{}\widehat{\mathbf{C}}_{\widehat{\mathcal{S}}}$ (49).

7: Update $\widehat{\rho}$, $\widehat{\nu}$ (IV-E3) and Update $\bm{\eta}_{i}$
and $\widehat{\mathbf{a}}_{i}$ for all $i\in\widehat{\mathcal{S}}$.

8: Compute the posterior means and variances of $\mathbf{h}$ as
${\mathbf{h}}_{\rm A}^{\rm post}$ (53a), ${\mathbf{v}}_{{\rm A},{\rm h}}^{\rm
post}$ (53b).

9: Calculate the extrinsic message $m_{{\rm A}\rightarrow{\rm
B}}({\mathbf{h}}_{\mathcal{M}})$ from module B to module A.

10: Calculate the message $m_{{\rm B}\rightarrow{\rm
C}}({\mathbf{d}}_{\mathcal{D}})$ (IV-A) and input it to the LDPC decoder
module to obtain the posterior means $d_{{\rm C},i}^{\rm post}$ and variances
$v_{{\rm C},d_{i}}^{\rm post}$ of the data.

11: Compute the posterior means and variances of $\mathbf{z}$ as
${\mathbf{z}}_{\rm B}^{\rm post}$ (20a) and ${\mathbf{v}}_{{\rm
B},{\mathbf{z}}}^{\rm post}$ (20b), respectively. Implement the Newston step
to update the residual CFO and the EM step to update the noise variance via
(22) and (29), respectively.

12: Calculate the extrinsic message $m_{{\rm B}\rightarrow{\rm
A}}({\mathbf{h}}_{\mathcal{M}})$ through (32) and (IV-D) from module B to
module A.

13: Update $\mathbf{J}$ (47a) and $\mathbf{u}$ (47b).

14: end for

15: Return belief function of data bits $u_{k}$, $k\in{\mathcal{K}}$, channel
estimate $\hat{\mathbf{h}}$.

##  V Numerical Simulation Results

In this section, numerical simulations are conducted to demonstrate the
advantages of the proposed joint processing scheme. A typical underwater
acoustic channel setting is adopted in the simulation: there are $15$ paths,
the arrival time between two adjacent paths follows exponential distribution
with mean value $1$ ms, and the amplitude of each path follows Rayleigh
distribution with mean power decreasing by $20$ dB in $30$ ms delay spread
[7]. The simulation parameters of the OFDM system are listed in Table I. In
addition, a common Doppler scaling factor exists in all paths, which follows
uniform distribution in the range $\left[-2\times 10^{-3},2\times
10^{-3}\right]$. In our simulation, null subcarriers are utilized to estimate
the Doppler scaling factor and the residual CFO [36], based on the results of
which the resampling and CFO compensation are implemented. Performance metrics
of decoding BER and NMSE of channel estimation defined as

$\displaystyle{\rm
NMSE}(\hat{\mathbf{h}})=10\log_{10}\left(\frac{\left\|\hat{\mathbf{h}}-\mathbf{h}_{0}\right\|^{2}}{\left\|\mathbf{h}_{0}\right\|^{2}}\right)$

are compared in the simulation, where $\mathbf{h}_{0}$, $\hat{\mathbf{h}}$ are
the ground truth and estimated value of the frequency-domain channel response
vector, respectively. All the results are averaged over $500$ MC trials.

TABLE I: CP-OFDM parameter settings Bandwidth | $B$ | $4882.8~{}\text{Hz}$
---|---|---
Carrier frequency | $f_{c}$ | $13000~{}\text{Hz}$
Sampling frequency | $f_{s}$ | $39062.5~{}\text{Hz}$
No. subcarriers | $N$ | $1024$
No. data subcarriers | $N_{\mathcal{D}}$ | $672$
No. pilot subcarriers | $N_{\mathcal{P}}$ | $256$
No. null subcarriers | $N_{\mathcal{N}}$ | $96$
Symbol duration | $T$ | $209.72~{}\text{ms}$
Cyclic-prefix length | $T_{\text{cp}}$ | $40.35~{}\text{ms}$

###  V-A Benchmark Algorithms

Several approaches classified as pilot-only based and joint channel and data
decoding based are implemented to make performance comparison. The pilot-only
based are listed as follows:

  * •

The pilot based OMP approach [7]. The dictionary is constructed by $\times 8$
oversampling for time delay search, and the Bayesian information stopping
criterion [37] is adopted, and the algorithm is termed as OMP.

  * •

The pilot based AMP approach [38]. It is numerically found that AMP diverges
with the oversampling factor greater than $2$. We construct the dictionary by
$\times 2$ oversampling. Besides, the Bernoulli Gaussian prior is adopted for
the coefficients and the EM algorithm [38, 39] is incorporated to jointly
estimate the nuisance parameters of the prior distribution and the noise
variance. The algorithm is termed as AMP.

  * •

The pilot based VALSE approach [29]. The observations corresponding to pilots
are directly input to the VALSE algorithm to obtain the path delays and
amplitudes, which are then used to reconstruct the whole channel. The
algorithm is termed as VALSE.

Note that after the whole channel is reconstructed, the estimates (including
means and variances) of the data are input to the single input single output
(SISO) equalizer, i.e, iterative decoding is conducted to obtain the BER.

The GAMP based _j_ oint _c_ hannel estimation and _d_ ata decoding approach
(JCD) named JCD-GAMP is also conducted [14]. Note that the original approach
uses the Markov chain prior and learn the nuisance parameters of the prior via
fitting the simulated channel data. Besides, it uses the true noise variance.
Here we use the Bernoulli Gaussian prior and use the EM to jointly estimate
the nuisance parameters of the prior distribution and the noise variance
instead.

To show the benefits of estimating the residual CFO, we also implement the
VALSE based JCD approach without estimating the CFO, and we name the approach
as JCD-VALSE.

To further show the effectiveness of the proposed algorithm, we also evaluate
the coded BER under the perfect channel state information (PCSI) and the NMSE
of the JCCD-VALSE with data aware, which are termed as PCSI and JCCD-VALSE
(data aware), respectively. Note that the coded BER of PCSI and the JCCD-VALSE
(data aware) will be the lower bounds of the coded BER and the NMSE,
respectively.

###  V-B QPSK Results

Figure 3: Performance metrics evaluated by the various algorithms for QPSK.
(a) coded BER, (b) NMSE.

Fig. 3 presents the performance metrics evaluated by the various algorithms
for QPSK modulation. For the coded BER shown in Fig. 3, PCSI works best,
followed by JCCD-VALSE, JCD-VALSE, JCD-GAMP, VALSE, OMP and AMP. For ${\rm
BER}=10^{-3}$, the SNR needed for the PCSI, JCCD-VALSE, JCD-VALSE, JCD-GAMP,
VALSE, and OMP are $4.21$ dB, $4.36$ dB, $4.47$ dB, $4.81$ dB, $5.11$ dB,
$6.61$ dB. The SNR gap between the JCCD-VALSE and the PCSI method (oracle) is
$4.36-4.21=0.15$ dB, which is very close. Compared to JCD-GAMP, VALSE and AMP
approaches, our JCCD method has SNR gains $4.81-4.36=0.45$ dB,
$5.11-4.36=0.75$ dB, $6.61-4.36=1.75$ dB, respectively. As for the NMSE of the
channel estimation, Fig. 3 shows that JCCD-VALSE (data aware) works best,
followed by JCCD-VALSE, JCD-VALSE (without estimating the residual CFO), JCD-
GAMP, VALSE, OMP and AMP. Besides, JCCD-VALSE asymptotically approaches to
JCCD-VALSE (data aware). To sum up, the proposed approach achieves excellent
data decoding and channel estimation performance for QPSK modulated symbols.

###  V-C $16$ QAM Results

The performance metrics evaluated by the various algorithms for $16$ QAM
modulation are shown in Fig. 4. Note that the results are similar to the QPSK
modulation setting except the following two differences. The first difference
is that because the constellation becomes larger, the SNR needed to achieve
the same coded BER and the NMSE increases (from about $4.36$ dB to $9.7$ dB of
JCCD-VALSE). The second difference is that the gain of the CFO estimation for
the JCCD-VALSE is marginal, compared to JCD-VALSE (without estimating the
residual CFO). The reason is probably that the Doppler scaling factor
estimation accuracy is better in the higher SNR region, which leads to smaller
residual CFO. Still, we conclude that the proposed approach has performance
gains compared to the other approaches for $16$ QAM modulated symbols.

Figure 4: Performance metrics evaluated by the various algorithms for $16$
QAM. (a) coded BER, (b) NMSE.

##  VI Sea Trial Data Decoding Results

To further verify the performance of the proposed JCCD-VALSE scheme in real
underwater acoustic communication channels, sea test data from MACE10 is
adopted for decoding. MACE10 experiment was carried out off the coast of
Marthas Vineyard, Massachusetts, at Jun. 2010. The depth of water was about 80
meters. The source was towed from the ship and towed back at relatively slow
speeds (1-2 m/s). The relative distance between the transmitter and the
receiver changed from 500 m to 7 km. In total 31 transmissions were carried
out during the experiment, each containing 20 OFDM symbols in QPSK and 20
symbols in $16$ QAM, respectively. The transmitted OFDM signal adopts the same
parameters shown in Table I. For data decoding, one transmission data file
recorded during the turn of the source which has quite low SNR is excluded. In
Fig. 5, the estimated typical channel impulse response and the estimated
Doppler speed corresponding to the 30 transmissions in MACE10 are presented.

(a) Estimated channel impulse response

(b) Estimated Doppler speed

Figure 5: Channel condition in MACE10 [40].

###  VI-A QPSK Decoding Results

In total 600 OFDM symbols in QPSK modulation were collected during the sea
experiment. Since all the original received data can be decoded successfully,
additional noise in different level are added in order to show the differences
between different algorithms. Denote the original noise variance in the
received signal to be $\delta^{2}$, the decoding BER results corresponding to
different additional noise levels are shown in Table II. It can be observed
that the proposed JCCD-VALSE algorithm always achieves the lowest BER in all
additional noise levels.

TABLE II: MACE10 decoding BER: QPSK Algorithm | OMP | AMP | VALSE | JCCD-VALSE | JCD-VALSE | JCD-GAMP
---|---|---|---|---|---|---
$0.8\delta^{2}$ | 0.0015 | 0.0021 | 0.0010 | 0 | 0 | 0.0001
$1.6\delta^{2}$ | 0.0085 | 0.0100 | 0.0078 | 0.0056 | 0.0058 | 0.0070
$2.4\delta^{2}$ | 0.0182 | 0.0194 | 0.0153 | 0.0118 | 0.0118 | 0.0125
$3.2\delta^{2}$ | 0.0338 | 0.0382 | 0.0281 | 0.0204 | 0.0214 | 0.0235
$4.0\delta^{2}$ | 0.0511 | 0.0568 | 0.0460 | 0.0358 | 0.0363 | 0.0397
$4.8\delta^{2}$ | 0.0662 | 0.0712 | 0.0616 | 0.0531 | 0.0535 | 0.0578

###  VI-B $16$ QAM decoding results

For 16 QAM modulation, 523 complete OFDM symbols were collected during the
experiment. A direct decoding of the received signals with different
algorithms leads to the results shown in Table III, which clearly shows that
the proposed JCCD-VALSE scheme achieves the best BER performance.

TABLE III: MACE10 decoding BER: 16 QAM Algorithm | OMP | AMP | VALSE | JCCD-VALSE | JCD-VALSE | JCD-GAMP
---|---|---|---|---|---|---
BER | 0.0152 | 0.0144 | 0.0130 | 0.0081 | 0.0084 | 0.0131

##  VII Conclusion

This paper proposes the joint CFO, gridless channel estimation and data
decoding algorithm named as JCCD-VALSE, which automatically estimates the
number of paths, the nuisance parameters of the prior distribution and the
noise variance. The algorithm is decomposed as three modules named the VALSE
module, the minimum mean squared error (MMSE) module and the decoder module,
and by iteratively exchange extrinsic information between modules, the channel
estimation and data decoding performances gradually improve. Numerical
experiments and real data are used to demonstrate the superior performance of
the proposed approach.

## References

  * [1] M. Stojanovic and J. Preisig, “Underwater acoustic communication channels: Propagation models and statistical characterization,” _IEEE Commun. Mag._ , vol. 47, no. 1, pp. 84-89, Jan. 2009.
  * [2] B. S. Sharif, J. Neasham, O. R. Hinton and A. E. Adams, “A computationlly efficient Doppler compensation system for underwater acoustic communications,” _IEEE J. Ocean. Eng._ , vol. 25, no. 1, pp. 52-61, Jan. 2000.
  * [3] M. Stojanovic, “Low complexity OFDM detector for underwater channels,” _in Proc. of MTS/IEEE OCEANS Conference_ , Boston, MA, Sept. 18-21, 2006.
  * [4] B. Li, S. Zhou, M. Stojanovic, L. Freitag and P. Willett, “Multicarrier communication over underwater acoustic channels with nonuniform Doppler shifts,” _IEEE J. Ocean. Eng._ , vol. 33, no. 2, pp. 198-209, April 2009.
  * [5] S. Zhao, S. Yan and J. Xi, “Adaptive turbo equalization for differential OFDM systems in underwater acoustic communications,” _IEEE Trans. Veh. Technol._ , vol. 69, no. 11, pp. 13937-13941, 2020.
  * [6] W. Li and J. C. Preisig, “Estimation of rapidly time-varying sparse channels,” _IEEE J. Ocean. Eng._ , vol. 32, no. 4, pp. 927-939, Oct. 2007.
  * [7] C. R. Berger, S. Zhou, J. C. Preisig, and P. Willett, “Sparse channel estimation for multicarrier underwater acoustic communication: From subspace methods to compressed sensing,” _IEEE Trans. Signal Process._ , vol. 58, no. 3, pp. 1708-1721, 2010.
  * [8] J. A. Tropp and A. C. Gilbert, “Signal recovery from random measurements via orthogonal matching pursuit,” _IEEE Trans. Inform. Theory_ , vol. 53, no. 12, pp. 4655-4666, 2007.
  * [9] B. Peng, P. S. Rossi, H. Dong and K. Kansanen, “Time-domain oversampled OFDM communication in doubly-selective underwater acoustic channels,” _IEEE Commun. Lett._ , vol. 19, no. 6, pp. 1081-1084, June 2015.
  * [10] E. Panayirci, M. T. Altabbaa, M. Uysal, and H. V. Poor, “Sparse channel estimation for OFDM-based underwater acoustic systems in Rician fading with a new OMP-MAP algorithm,” _IEEE Trans. Signal Process._ , vol. 67, no. 6, pp. 1550-1565, 2019.
  * [11] D. L. Donoho, A. Maleki, and A. Montanari, “Message passing algorithms for compressed sensing: I. Motivation and construction,” _in Proc. Inf. Theory Workshop_ , Cairo, Egypt, Jan. 2010, pp. 1-5.
  * [12] C. Wei, H. Liu, Z. Zhang, J. Dang and L. Wu, “Near-optimum sparse channel estimation based on least squares and approximate message passing,” _IEEE Wireless Commun. Lett._ , vol. 6, no. 6, pp. 754-757, Dec. 2017.
  * [13] E. Panayirci, H. Senol and H. V. Poor, “Joint channel Estimation, equalization, and data detection for OFDM systems in the presence of very high mobility,” _IEEE Trans. Signal Process._ , vol. 58, no. 8, pp. 4225-4238, Aug. 2010.
  * [14] P. Schniter, “A message-passing receiver for BICM-OFDM over unknown clustered-sparse channels,” _IEEE J. Sel. Topics Signal Process._ , vol. 5, no. 8, pp. 1462-1474, Dec. 2011.
  * [15] P. Sun, Z. Wang and P. Schniter, “Joint channel-estimation and equalization of single-carrier systems via bilinear AMP,” _IEEE Trans. Signal Process._ , vol. 66, no. 10, pp. 2772-2785, May 2018.
  * [16] Z. Jiang, X. Shen, H. Wang and Z. Ding, “Joint PSK data detection and channel estimation under frequency selective sparse multipath channels,” _IEEE Trans. Commun._ , vol. 68, no. 5, pp. 2726-2739, May 2020.
  * [17] J. F. Sifferlen, H. C. Song, W. S. Hodgkiss, W. A. Kuperman and J. M. Stevenson, “An iterative equalization and decoding approach for underwater acoustic communication,” _IEEE J. Ocean. Eng._ , vol. 33, no. 2, pp. 182-197, April 2008.
  * [18] Z. Yang and Y. R. Zheng, “Iterative channel estimation and Turbo equalization for multiple-input multiple-output underwater acoustic communications,” _IEEE J. Ocean. Eng._ , vol. 41, no. 1, pp. 232-242, Jan. 2016.
  * [19] Z. Chen, J. Wang and Y. R. Zheng, “Frequency-domain Turbo equalization with iterative channel estimation for MIMO underwater acoustic communications,” _IEEE J. Ocean. Eng._ , vol. 42, no. 3, pp. 711-721, July 2017.
  * [20] X. Qin, F. Qu and Y. R. Zheng, “Bayesian iterative channel estimation and Turbo equalization for multiple-input multiple-output underwater acoustic communications,” _IEEE J. Ocean. Eng._ , vol. 46, no. 1, pp. 326-337, Jan. 2021.
  * [21] K. P. Arunkumar and C. R. Murthy, “Iterative sparse channel estimation and data detection for underwater acoustic communications using partial interval demodulation,” _IEEE Trans. Signal Process._ , vol. 66, no. 19, pp. 5041-5055, Oct., 2018.
  * [22] X. Tu, X. Xu and A. Song, “Frequency-domain decision feedback equalization for single-carrier transmissions in fast time-varying underwater acoustic channels,” _IEEE J. Ocean. Eng._ , vol. 46, no. 2, pp. 704-716, April 2021.
  * [23] S. Wang, Z. He, K. Niu, P. Chen and Y. Rong, “New results on joint channel and impulsive noise estimation and tracking in underwater acoustic OFDM systems,” _IEEE Trans. Wireless Commun._ , vol. 19, no. 4, pp. 2601-2612, April 2020.
  * [24] G. Yang, Q. Guo, H. Ding, Q. Yan and D. D. Huang, “Joint message-passing-based bidirectional channel estimation and equalization with superimposed training for underwater acoustic communications,” _IEEE J. Ocean. Eng._ , early access, 2021.
  * [25] Y. Chi, L. L. Scharf, A. Pezeshki and R. Calderbank, “Sensitivity of basis mismatch to compressed sensing,” _IEEE Trans. on Signal Process._ , vol. 59, pp. 2182 - 2195, 2011.
  * [26] T. L. Hansen, B. H. Fleury and B. D. Rao, “Superfast line spectral estimation,” _IEEE Trans. Signal Process._ , vol. 66, no. 10, pp. 2511-2526, 2018.
  * [27] T. L. Hansen, P. B. Jørgensen, M. A. Badiu and B. H. Fleury, “An iterative receiver for OFDM with sparsity-based parametric channel estimation,” _IEEE Trans. Signal Process._ , vol. 66, no. 20, pp. 5454-5469, 2018.
  * [28] X. Meng, S. Wu and J. Zhu, “A unified Bayesian inference framework for generalized linear models,” _IEEE Signal Process. Lett._ , vol. 25, no. 3, pp. 398-402, 2018.
  * [29] M. A. Badiu, T. L. Hansen and B. H. Fleury, “Variational Bayesian inference of line spectra,” _IEEE Trans. Signal Process._ , vol. 65, no. 9, pp. 2247-2261, 2017.
  * [30] K. V. Mardia and P. E. Jupp, _Directional Statistics_. New York, NY, USA: Wiley, 2000.
  * [31] D. J. C. MacKay, “Bayesian interpolation,” _Neural Computation_ , vol. 4, no. 3, pp. 415-447, 1992.
  * [32] J. O. Berger, _Statistical Decision Theory and Bayesian Analysis_. Springer, second edition, 1985.
  * [33] J. Winn and C. M. Bishop, “Variational message passing,” _J. Mach. Learn. Res._ , vol. 6, pp. 661-694, 2005.
  * [34] T. Minka, “A family of algorithms for approximate Bayesian inference,” Ph.D. dissertation, Department of Electrical Engineering and Computer Science, Mass. Inst. Technol., Cambridge, MA, USA, 2001.
  * [35] J. Zhu, Q. Zhang, and X. Meng, “Grid-less variational Bayesian inference of line spectral from quantized samples,” 2019, arXiv:1811.05680, to appear in _China Communications_. [Online]. Available: https://arxiv.org/abs/1811.05680
  * [36] L. Wan, Z.-H. Wang, S. Zhou, T. C. Yang, and Z. Shi, “Performance comparison of Doppler scale estimation methods for underwater acoustic OFDM,” _Journal of Electrical and Computer Engineering_ , vol. 2012, Article ID 703243, 2012.
  * [37] N. Eldarov, G. Tan, and T. Herfet, “Delay-Doppler search for matching pursuit algorithms in time-variant channels,” _in Proc. IEEE Int. Symp. Broadband Multimedia Syst. Broadcast._ , Jun. 2015, pp. 1-5.
  * [38] J. P. Vila and P. Schniter, “Expectation-Maximization Bernoulli-Gaussian approximate message passing,” _Proc. Asilomar Conf. on Signals, Systems, and Computers_ , Pacific Grove, CA, Nov. 2011.
  * [39] J. P. Vila and P. Schniter, “Expectation-Maximization Gaussian-mixture approximate message passing,” _IEEE Trans. Signal Process._ , vol. 61, no. 19, pp. 4658-4672, Oct. 2013.
  * [40] L. Wan, H. Jia, F. Zhou, M. Muzzammil, T. Li and Y. Huang, “Fine Doppler scale estimations for an underwater acoustic CP-OFDM system,” _Signal Process._ , vol. 170, May 2020.