# Message Passing Based Structured Sparse Signal Recovery for Estimation of
OTFS Channels with Fractional Doppler Shifts

Fei Liu, Zhengdao Yuan, Qinghua Guo, , Zhongyong Wang and Peng Sun The work of
F. Liu, Z. Yuan and Z. Wang was supported by National Natural Science
Foundation of China (61901417,U1804152) and Science and Technology Research
Project of Henan Province (212102210556, 202102210313, 212102210566).
Corresponding authors: Qinghua Guo and Zhongyong Wang.F. Liu, Z. Wang and P.
Sun are with the School of Information Engineering, Zhengzhou University,
Zhengzhou 450002, China, (e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>iepengsun@zzu.edu.cn)Z. Yuan is with the Artificial Intelligence Technology
Engineering Research Center, Open University of Henan, and also with the
School of Information Engineering, Zhengzhou University, Zhengzhou 450002,
China (e-mail: yuan_zhengdao@163.com).Q. Guo is with the School of Electrical,
Computer and Telecommunications Engineering, University of Wollongong,
Wollongong, NSW 2522, Australia (e-mail: qguo@uow.edu.au).

###### Abstract

The orthogonal time frequency space (OTFS) modulation has emerged as a
promising modulation scheme for high mobility wireless communications. To
enable efficient OTFS detection in the delay-Doppler (DD) domain, the DD
domain channels need to be acquired accurately. To achieve the low latency
requirement in future wireless communications, the time duration of the OTFS
block should be small, therefore fractional Doppler shifts have to be
considered to avoid significant modelling errors due to the assumption of
integer Doppler shifts. However there lack investigations on the estimation of
OTFS channels with fractional Doppler shifts in the literature. In this work,
we develop a channel estimator for OTFS with particular attention to
fractional Doppler shifts, and both bi-orthogonal waveform and rectangular
waveform are considered. Instead of estimating the DD domain channel directly,
we estimate the channel gains and (fractional) Doppler shifts that
parameterize the DD domain channel. The estimation is formulated as a
structured sparse signal recovery problem with a Bayesian treatment. Based on
a factor graph representation of the problem, an efficient message passing
algorithm is developed to recover the structured sparse signal (thereby the
OTFS channel). The Cramer-Rao Lower Bound (CRLB) for the estimation is
developed and the effectiveness of the algorithm is demonstrated through
simulations.

###### Index Terms:

Orthogonal time frequency space modulation, message passing, channel
estimation, fractional Doppler shifts.

##  I Introduction

Rencently, the orthogonal time frequency space (OTFS) modulation has been
proposed to achieve reliable communications in high-mobility scenarios
[1],[2],[3]. OTFS provides both time and frequency diversity because each
symbol is spread over the time and frequency domains through the two
dimensional inverse symplectic finite Fourier transform (ISFFT) [1],[2]. When
the number of channel paths is small, the effective channel in the delay-
Doppler (DD) domain is sparse, which allows efficient data detection using the
message passing techniques [2]. A variety of OTFS detection methods have been
proposed in the literature to harvest the time and frequency diversity
promised by OTFS [4, 5, 6, 7, 8, 9, 10, 11]. However, all the detection
methods assume perfect channel state information, which has to be estimated in
practice.

In OTFS, the delay shifts and Doppler shifts are discretized in the DD domain.
In general, a wideband system is able to provide sufficient delay resolution,
so that fractional delay shifts do not need to be considered [12]. However,
the Doppler resolution depends on the time duration of the OTFS block. To
fulfill the low latency requirement in future wireless communications, the
time duration of the OTFS block should be relatively small, hence fractional
Doppler shifts have to be considered to avoid significant modelling errors due
to the assumption of integer Doppler shifts [2],[13]. Therefore factional
Doppler shifts have to be considered in OTFS channel estimation.

A number of channel estimation methods have been proposed in the literature.
In [14], an entire OTFS block is used to accommodate pilot symbols for channel
estimation, and the estimated channel is used for data detection in the
subsequent OTFS block. This results in a considerable loss in spectrum
efficiency and the detection performance may deteriorate due to the channel
variation between two OTFS blocks. To solve this problem, pilot and data
symbols are placed in the same OTFS block in [13], where guard interval is
used to avoid interference between pilot and data symbols. With this scheme,
channel estimation and data detection can be performed for the same OTFS
block. In [15], channel estimation in massive multiple input and multiple
output (MIMO)-OTFS systems is studied, where the downlink time-varying massive
MIMO channels are transformed to the delay-Doppler-angle domain, and the
channel estimation is formulated as a sparse signal recovery problem. In [16],
the uplink-aided downlink massive MIMO-OTFS channel estimation over the delay-
Doppler-angle domain is studied. However, in these works, only integer Doppler
shifts are considered, except [13]. In [13], a single pilot symbol is used to
facilitate the design of a low complexity thresholding method to acquire the
channel in the DD domain. Although factional Doppler shifts are considered for
OTFS with the bi-orthogonal waveform in [13], it is not clear how to estimate
the channel with fractional Doppler shifts for OTFS with a more practical
waveform, e.g., the rectangular waveform. In addition, the method estimates
the DD domain channel directly, and does not provide the estimates of the
channel gains and fractional Doppler shifts.

In this work, we address the issue of OTFS channel estimation with particular
attention to fractional Doppler shifts. As the OTFS channel is parameterized
by the channel gains and fractional Doppler shifts, we estimate these
parameters rather than the DD domain channel directly, mainly due to two
reasons. Firstly, this can lead to much better performance as the number of
variables to be estimated is significantly smaller. Secondly, the estimates of
Doppler shifts can be useful, e.g., for estimating the velocity of mobile
users. In addition, both the bi-orthogonal waveform and the rectangular
waveform are considered. We formulate the OTFS channel estimation as a
structured signal recovery problem, which is solved using Bayesian inference.
With a factor graph representation of the problem, a message passing algorithm
is developed to estimate the channel gains and fractional Doppler shifts
jointly. In contrast to [13], we show that using a few pilot symbols rather
than a single pilot symbol is more preferable as the peak to average power
ratio (PAPR) of OTFS signals can be significantly reduced. Our proposed
algorithm can work with either single pilot symbol or multiple pilot symbols,
and significantly outperforms the state-of-the-art method. The Cramer-Rao
lower bound (CRLB) for the channel estimation is derived to verify the
performance of the proposed algorithm. The bit error rates (BER) of the OTFS
system with perfect channel and estimated channel are also compared to
demonstrate the effectiveness of the proposed algorithm.

The main contributions of this work are summarized as follows:

  * •

The OTFS DD domain channel is acquired by estimating the relevant parameters,
i.e., the channel gains and fractional Doppler shifts. The estimation of the
parameters is formulated as a novel structured sparse signal recovery problem
with a Bayesian treatment.

  * •

Both the bi-orthogonal waveform and the rectangular waveform are considered in
this work. To the best of our knowledge, the investigations on the fractional
Doppler shift estimation are very limited in the literature. The work in [13],
where the DD domain channel is estimated directly, does not provide fractional
Doppler estimate explicitly, and it is not applicable to the case of the
rectangular waveform.

  * •

Due to the uniqueness of the formulated structured sparse signal recovery
problem, a dedicated message passing algorithm is proposed to efficiently
solve the problem. The algorithm is able to work with either a single pilot
symbol or multiple pilot symbols, which can significantly reduce the PAPR of
the OTFS signal.

  * •

To evaluate the performance of the estimator, the CRLB is derived. Comparison
results with existing methods are provided to demonstrate the advantages of
the proposed method.

For simplicity, this work is focused on a plain OTFS system, but the proposed
algorithm can be extended for a more complex system, such as MIMO-OTFS.

The remainder of the paper is organized as follows. In Section II, we
introduce OTFS modulation and demodulation, and OTFS input-output relation in
the DD domain. Then we formulate the channel estimation as a structured sparse
signal recovery problem in Section III, and develop the message passing based
algorithm in Section IV. The CLRB for channel estimation is derived in Section
V. Simulation results are provided in Section VI, followed by conclusions in
Section VII.

Notations\- Boldface lower-case and upper-case letters denote column vectors
and matrices, respectively. The superscripts $(\cdot)^{T}$ and $(\cdot)^{*}$
represent transpose and conjugate operations, respectively. We use
$[\cdot]_{M}$ to denote the modulo-$M$ operation. The probability density
function of a complex Gaussian variable with mean $\hat{x}$ and variance
$\nu_{x}$ is represented by $\mathcal{CN}(x;{\hat{x}},\nu_{x})$. The Gamma
distribution with scale $\epsilon$ and rate $\eta$ is represented as
$Ga(x;\epsilon,\eta)$. The uniform distribution over the range $[a,b]$ is
represented by $U[a,b]$. The relation $f(x)=cg(x)$ for some positive constant
$c$ is written as $f(x)\propto g(x)$. The notation $\otimes$ represents the
Kronecker product, and $\bm{a}\cdot\bm{b}$ and $\bm{a}\cdot/\bm{b}$ represent
the element-wise product and division between vectors $\bm{a}$ and $\bm{b}$,
respectively. We use $Tr(\bm{A})$ and $\bm{I}_{N}$ to denote the trace of
$\bm{A}$ and an $N\times N$ identity matrix, respectively. We use $|x|^{2}$ to
denote the magnitude squared operation for $x$ and use $||\bm{x}||_{2}^{2}$ to
denote the squared norm of vector $\bm{x}$. We use $q_{j}$ to denoted the
$j$th element of $\bm{q}$. The superscript $t$ of $\textbf{s}^{t}$ denotes the
iteration index in an iterative algorithm. The notation $\mathbb{E}(x)$ is
used to denote the expectation of the random variable $x$. The notations
$\mathcal{R}[a]$ and $\mathcal{I}[a]$ denote the real and imaginary parts of
the complex variable $a$, respectively.

ISFFT Heisenberg transform Channel $h(\tau,\nu)$ Wigner transform SFFTTime-
Frequency DomainDelay-Doppler
Domain$x[k,l]$$X[n,m]$$s(t)$$r(t)$$Y[n,m]$$y[k,l]$ Figure 1: OTFS modulation
and demodulation [2].

##  II OTFS System Model

###  II-A System Model in the DD Domain

As shown in Fig. 1, the OTFS modulation and demodulation are implemented with
2D inverse SFFT (ISFFT) and SFFT at the transmitter and receiver, respectively
[1][17]. A (coded) bit sequence is mapped to symbols
$\\{x[k,l],k=0,\cdots,N-1,l=0,\cdots M-1\\}$ in the DD domain, where
$x[k,l]\in\mathcal{A}=\\{\alpha_{1},....\alpha_{|\mathcal{A}|}\\}$ with
$|\mathcal{A}|$ being the cardinality of $\mathcal{A}$, and $l$ and $k$ denote
the indices of the delay and Doppler shifts, respectively. As shown in Fig. 1,
ISFFT is performed to convert the symbols to signals in the time-frequency
(TF) domain, i.e.,

$\displaystyle
X_{tf}[n,m]=\frac{1}{\sqrt{MN}}\sum_{k=0}^{N-1}\sum_{l=0}^{M-1}x[k,l]e^{j2\pi(\frac{nk}{N}-\frac{ml}{M})}.$
(1)

Then $\\{X_{tf}[m,n]\\}$ are converted to a continuous-time waveform $s(t)$
using the Heisenberg transform with a transmit waveform $g_{tx}(t)$, i.e.,

$\displaystyle
s(t)=\sum_{n=0}^{N-1}\sum_{m=0}^{M-1}X_{tf}[n,m]g_{tx}(t-nT)e^{j2\pi m\Delta
f(t-nT)},$ (2)

where $\Delta f$ is the subcarrier spacing and $T=1/\Delta f$.

The signal $s(t)$ is then transmitted through a time-varying channel and the
received signal in the time domain can be expressed as

$\displaystyle r(t)=\int\int h(\tau,\nu)s(t-\tau)e^{j2\pi\nu(t-\tau)}d\tau
d\nu,$ (3)

where $h(\tau,\nu)$ is the channel impulse response in the (continuous) DD
domain. The channel impulse response can be expressed as

$\displaystyle
h(\tau,\nu)=\sum_{i=1}^{P}h_{i}\delta(\tau-\tau_{i})\delta(\nu-\nu_{i}),$ (4)

where $\delta(\cdot)$ is the Dirac delta function, $P$ is the number of
resolvable propagation paths, and $h_{i}$, $\tau_{i}$ and $\nu_{i}$ represent
the gain, delay shift and Doppler shift associated with the $i$th path,
respectively. The delay and Doppler shift taps for the $i$th path are

$\displaystyle\tau_{i}=\frac{l_{i}}{M\Delta f},$ (5)
$\displaystyle\nu_{i}=\frac{k_{i}+\kappa_{i}}{NT},$ (6)

where $0\leq l_{i}\leq l_{max}$ and $-k_{max}\leq k_{i}\leq k_{max}$ are the
delay index and Doppler index of the $i$th path, $l_{max}$ and $k_{max}$
represent the largest indices of the delay taps and Doppler taps,
respectively, and $\kappa_{i}$ $\in[-0.5,0.5]$ is the fractional Doppler shift
associated with the $i$th path.

At the receiver side, a receive waveform $g_{rx}(t)$ is used to transform the
received signal $r(t)$ to the TF domain i.e.,

$\displaystyle Y(t,f)=\int g_{rx}^{*}(t^{\prime}-t)r(t^{\prime})e^{-j2\pi
f(t^{\prime}-t)}dt^{\prime},$ (7)

which is then sampled at $t=nT$ and $f=m\Delta f$, yielding $Y[n,m]$. Finally,
an SFFT is applied to $\\{Y[n,m]\\}$ to obtain the signal $y[k,l]$ in the DD
domain, i.e.,

$\displaystyle
y[k,l]=\frac{1}{\sqrt{MN}}\sum_{n=0}^{N-1}\sum_{m=0}^{M-1}Y[n,m]e^{-j2\pi(\frac{nk}{N}-\frac{ml}{M})}.$
(8)

If the transmit waveform $g_{tx}(t)$ and receive waveform $g_{rx}(t)$ satisfy
the bi-orthogonal property [1], the channel input-output relationship in the
DD domain can be expressed as [2, 9]

$\displaystyle
y[k,l]=\sum_{i=1}^{P}\sum_{q=-N_{i}}^{N_{i}}h_{i}\frac{1-e^{-j2\pi(-q-\kappa_{i})}}{N-Ne^{-j\frac{2\pi}{N}(-q-\kappa_{i})}}e^{-j2\pi\frac{l_{i}(k_{i}+\kappa_{i})}{MN}}$
$\displaystyle\times
x\big{[}[k-k_{i}+q]_{N},[l-l_{i}]_{M}\big{]}+\omega[k,l],$ (9)

where $N_{i}\ll N$ is an integer, and $\omega[k,l]$ denotes the Gaussian noise
in the DD domain with mean 0 and variance $\gamma^{-1}$ (or precision
$\gamma$). We can see that for each path, the transmitted signal is circularly
shifted, and scaled by a channel gain. We stack $\\{x[k,l]\\}$ to form a
vector $\bm{x}\in\mathbb{C}^{MN\times 1}$, where the $j$th element $x_{j}$ is
$x[k,l]$ with $j=kM+l$. Similarly, a vector $\bm{y}\in\mathbb{C}^{MN\times 1}$
can also be constructed from $\\{y[k,l]\\}$. Then (9) can be rewritten in a
vector form as

$\displaystyle\bm{y}=\bm{H}_{bi}\bm{x}+\bm{\omega},$ (10)

where $\bm{\omega}$ is the corresponding noise vector, and
$\bm{H}_{bi}\in\mathbb{C}^{MN\times MN}$ represents the effective channel in
the DD domain, which can be expressed as

$\displaystyle\bm{H}_{bi}=\sum_{i=1}^{P}\sum_{q=-N_{i}}^{N_{i}}\bm{I}_{N}(-[q-k_{i}]_{N})\otimes\Big{[}\bm{I}_{M}(l_{i})h_{i}$
$\displaystyle\times\left(\frac{1-e^{-j2\pi(-q-\kappa_{i})}}{N-Ne^{-j\frac{2\pi}{N}(-q-\kappa_{i})}}\right)e^{-j2\pi\frac{l_{i}(k_{i}+\kappa_{i})}{MN}}\Big{]},$
(11)

where $\bm{I}_{N}(-[q-k_{i}]_{N})$ denotes an $N\times N$ matrix obtained by
circularly shifting the rows of an identity matrix by $-[q-k_{i}]_{N}$, and
$\bm{I}_{M}(l_{i})$ is obtained similarly.

When the rectangular waveform is used for both $g_{tx}(t)$ and $g_{rx}(t)$,
the received signal $y[k,l]$ in the DD domain can be expressed as [2, 9]

$\displaystyle
y[k,l]=\sum_{i=1}^{P}\sum_{q=-N_{i}}^{N_{i}}h_{i}e^{j2\pi\left(\frac{l-l_{i}}{M}\right)\left(\frac{k_{i}+\kappa_{i}}{N}\right)}\alpha_{i}(k,l,q)$
$\displaystyle\times x[[k-k_{i}+q]_{N},[l-l_{i}]_{M}]+\omega[k,l],$ (12)

where

$\alpha_{i}(k,l,q)=\begin{cases}\frac{1}{N}\beta_{i}(q)&l_{i}\leq l<M\\\
\frac{1}{N}(\beta_{i}(q)-1)e^{-j2\pi\frac{[k-k_{i}+q]_{N}}{N}}&0\leq
l<l_{i}\end{cases},$ (13)

and

$\displaystyle\beta_{i}(q)=\frac{e^{-j2\pi(-q-\kappa_{i})}-1}{e^{-j\frac{2\pi}{N}(-q-\kappa_{i})}-1}.$
(14)

Similarly, we can also rewrite (12) in a vector form as

$\displaystyle\bm{y}=\bm{H}_{rect}\bm{x}+\bm{\omega},$ (15)

and the channel matrix $\bm{H}_{rect}\in\mathbb{C}^{MN\times MN}$ can be
written as

$\displaystyle\bm{H}_{rect}=$
$\displaystyle\sum_{i=1}^{P}\sum_{q=-N_{i}}^{N_{i}}\bm{I}_{N}(-[q-k_{i}]_{N})\otimes$
$\displaystyle\left(\bm{\Lambda}\bm{I}_{M}(l_{i})h_{i}e^{-j2\pi\frac{l_{i}(k_{i}+\kappa_{i})}{MN}}\right)\cdot\bm{\Delta}^{-[q-k_{i}]_{N}},$
(16)

where $\bm{\Lambda}\in\mathbb{C}^{M\times M}$ is a diagonal matrix and the
$l$th diagonal element $\bm{\Lambda}_{ll}$ can be expressed as

Figure 2: An OTFS block with pilot symbols inserted in the DD domain [15].

$\displaystyle\bm{\Lambda}_{ll}=\begin{cases}e^{j2\pi\frac{l(k_{i}+\kappa_{i})}{MN}}\beta_{i}(q)/N&l_{i}\leq
l<M\\\
e^{j2\pi\frac{l(k_{i}+\kappa_{i})}{MN}}\left(\beta_{i}(q)-1\right)/N&0\leq
l<l_{i}\\\ \end{cases}.$ (17)

Note that $\bm{\Delta}^{-[q-k_{i}]_{N}}$ is an $MN\times MN$ block matrix,
which is obtained by circularly shifting the blocks in a matrix $\bm{\Delta}$
by $-[q-k_{i}]_{N}$. The matrix $\bm{\Delta}$ is a block diagonal matrix

$\displaystyle\bm{\Delta}=\left(\begin{array}[]{cccccc}\bm{\Delta}_{0}&\bm{0}&\cdots&\bm{0}\\\
\bm{0}&\bm{\Delta}_{1}&\cdots&\bm{0}\\\ \vdots&\vdots&\ddots&\vdots\\\
\bm{0}&\bm{0}&\cdots&\bm{\Delta}_{N-1}\\\ \end{array}\right),$ (22)

where $\bm{\Delta}_{n}=\bm{\Psi I}_{M}(l_{i})$ and
$\bm{\Psi}\in\mathbb{C}^{M\times M}$ is a diagonal matrix with the $m$th
diagonal element ${\Psi}_{mm}$ given as

$\displaystyle{\Psi}_{mm}=\begin{cases}1&l_{i}\leq m<M\\\ e^{-j2\pi m/N}&0\leq
m<l_{i}\\\ \end{cases}.$ (23)

Figure 3: Received OTFS signal block where the signals in the orange area are
used for channel estimation [15].

###  II-B Data and Pilot Symbols Arrangement

In this paper, we use the scheme in [13] and [15], where pilot symbols are
inserted to the DD plane with guard interval to avoid the interference between
data and pilot symbols, as shown in Fig. 2. This allows that channel
estimation and data detection are carried out for the same OTFS block.

As shown in Fig. 2, the size of the pilot symbol block is $M_{p}\times N_{p}$.
To avoid the interference between pilot and data symbols, the size of guard
interval should be at least $2(k_{max}+\widehat{N})$ along the Doppler
dimension and $l_{max}$ along the delay dimension, where we assume
$\\{N_{i},i=1,2,\cdots,P\\}$ take a same value $\widehat{N}$. With this
scheme, the overhead due to pilot symbols and guard interval is
$(2l_{max}+M_{p})(4k_{max}+4\widehat{N}+N_{p})/MN$. The received signal in the
DD domain is shown in Fig. 3, where the signals in the orange area correspond
to the pilot symbols, and they are used for channel estimation.

In order to minimize the spectrum overhead due to pilot symbols, the use of a
small number of pilot symbols is highly desirable. To guarantee the channel
estimation quality, the power of the pilot symbols should be sufficiently
high, which can be much larger than the power of data symbols. Fortunately,
the operation of ISFFT at the transmitter can spread the power of pilot
symbols, therefore greatly alleviating the PAPR problem [13]. In [13], a
single pilot symbol with a very high power is used to facilitate the design of
a low complexity channel estimation algorithm based on thresholding. However,
the use of a few pilot symbols can be more preferable because this can
significantly reduce the PAPR at the cost of a small loss in spectrum
efficiency. In Table I, we compare the average PAPR of OTFS signals in the
time domain, where we assume that the rectangle waveform is used. The pilot
symbols are arranged in a single column and placed along the delay dimension,
i.e., $N_{p}=1$. We assume that the signal to noise ratio (SNR) of data
symbols is 14 dB, and compare the average PAPR of the OTFS signals versus the
number of pilot symbols and the SNR of pilot symbols (denoted by SNRp). It can
be seen that more pilot symbols are used, a lower PAPR can be achieved.
Moreover, with a fixed power budget for pilot, the use of multiple pilot
symbols leads to considerable PAPR reduction. Take an example: the use of a
single pilot symbol with SNRp = 50dB has the same power budget as the use of
10 pilot symbols with SNRp = 40dB. From Table I, a significant PAPR reduction
of about 8dB can be achieved. It is worth mentioning that, with the proposed
method in this paper, almost the same estimation performance can be achieved
for both cases (as shown later), making the use of multiple pilot symbols very
attractive.

TABLE I: PAPR of OTFS signals in the time domain, where $M=128$, $N=32$ and the SNR of data symbol is 14 dB. SNRp |  # of pilot symbols |  PAPR
---|---|---
40dB |  1 |  12.5527 dB
|  10 |  10.4095 dB
50dB |  1 |  18.7619 dB
|  10 |  11.3900 dB

In this work, we design a flexible channel estimation method, which is able to
work with a single or multiple pilot symbols. Moreover, the proposed method is
able to estimate the channel gains and the fractional Doppler shifts, which
are crucial to improving the estimation performance as demonstrated later.

##  III Structured Sparse Signal Recovery For OTFS Channel Estimation

The OTFS channel can be estimated using the pilot symbols and the received
signals in the orange area shown in Fig. 3. It can be seen from (9) and (29)
that the OTFS channel is parameterized by the nonzero channel gains
$\\{h_{i}\\}$, the fractional Doppler shifts $\\{\kappa_{i}\\}$ and their
indices $\\{l_{i},k_{i}\\}$, $i=1,\cdots,P$. These parameters will be
estimated, based on which the OTFS channel in the DD domain can be
constructed.

###  III-A Bi-Orthogonal Waveform

Model (9) involves a number of unknown variables, including $h_{i}$,
$\kappa_{i}$, $l_{i}$ and $k_{i}$ for $i=1,\cdots,P$. In addition, $P$ is
unknown as well. These make the estimation very difficult. To overcome this,
we define two variables $t$ and $d$ to represent the indices of delay shifts
and Doppler shifts, respectively. As only the pilot symbols need to be
considered, we have $t\in[0,l_{max}]$ and $d\in[-k_{max},k_{max}]$. Then (9)
can be rewritten as

$\displaystyle
y[k,l]=\sum_{t=0}^{l_{max}}\sum_{d=-k_{max}}^{k_{max}}h_{t,d}e^{-j2\pi\frac{t(d+\kappa_{d})}{MN}}\sum_{q=-\widehat{N}}^{\widehat{N}}f(q,\kappa_{d})$
$\displaystyle x\big{[}[k-d+q]_{N},[l-t]_{M}\big{]}+\omega[k,l],$ (24)

where

$\displaystyle
f(q,\kappa_{d})=\frac{1-e^{-j2\pi(-q-\kappa_{d})}}{N-Ne^{-j\frac{2\pi}{N}(-q-\kappa_{d})}}.$
(25)

The received signals $\\{y[k,l]\\}$ in (24) corresponding to the pilot symbols
are in the orange area in Fig. 3. Our aim is to find the non-zero elements in
$\\{h_{t,d},\kappa_{d},t\in[0,l_{max}],d\in[-k_{max},k_{max}]\\}$ based on
$\\{y[k,l]\\}$. Therefore, the estimation can be formulated as a sparse signal
recovery problem.

To facilitate the design of the estimation algorithm, we rewrite (24) in a
matrix form as

$\displaystyle\bm{y}={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\bm{X}_{bi}}\bm{c}+\bm{\omega},$
(26)

where $\bm{y}\in\mathbb{C}^{Z\times 1}$,
$Z=(l_{max}+M_{p})(N_{p}+2k_{max}+2\widehat{N})$ is formed by stacking
$\\{y[k,l]\\}$ as a vector,
${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\bm{X}_{bi}}\in\mathbb{C}^{Z\times(l_{max}+1)(2k_{max}+1)(2\widehat{N}+1)}$
is constructed based on the pilot symbols, $\bm{\omega}\in\mathbb{C}^{Z\times
1}$ is obtained by stacking $\\{\omega[k,l]\\}$, and
$\bm{c}\in\mathbb{C}^{(l_{max}+1)(2k_{max}+1)B\times 1}$ with
$B=2\widehat{N}+1$ can be expressed as

$\displaystyle\bm{c}=[\bm{c}^{T}_{0,-k_{max}},~{}\bm{c}^{T}_{1,-k_{max}},~{}\cdots,~{}\bm{c}^{T}_{t,d},\cdots,\bm{c}^{T}_{l_{max},k_{max}}]^{T},$
(27)

where $\bm{c}_{t,d}\in\mathbb{C}^{B\times 1}$ is given as

$\displaystyle\bm{c}_{t,d}=h_{t,d}e^{-j2\pi\frac{(d+\kappa_{d})t}{MN}}\cdot$
$\displaystyle\ \
[f(-\widehat{N},\kappa_{d}),f(-\widehat{N}+1,\kappa_{d}),\cdots,f(\widehat{N},\kappa_{d})]^{T}.$
(28)

We can see that the vector $\bm{c}$ has a structure, i.e., $\bm{c}$ is block
sparse, and for a non-zero block (subvector) $\bm{c}_{t,d}$, it is
parameterized by only two parameters $h_{t,d}$ and $\kappa_{d}$. We will
recover the structured sparse vector $\bm{c}$ and obtain the estimates of the
non-zero channel gains and fractional Doppler shifts.

###  III-B Rectangular Waveform

To facilitate channel estimation, we place the pilot symbols to ensure that
their delay index $l\geq l_{max}$, so that (12) is reduced to

$\displaystyle y[k,l]=\\!\\!\\!\\!\\!\\!\\!\\!$
$\displaystyle\sum_{i=1}^{P}\sum_{q=-N_{i}}^{N_{i}}h_{i}\frac{1-e^{-j2\pi(-q-\kappa_{i})}}{N-Ne^{-j\frac{2\pi}{N}(-q-\kappa_{i})}}e^{j2\pi\frac{(l-l_{i})(k_{i}+\kappa_{i})}{MN}}$
(29) $\displaystyle~{}~{}\times
x\big{[}[k-k_{i}+q]_{N},[l-l_{i}]_{M}\big{]}+\omega[k,l].$

Comparing (29) to (9), we can find that their difference lies in that each
transmitted symbol in (29) is rotated by an additional phase
$e^{j2\pi\frac{l(k_{i}+\kappa_{i})}{MN}}$.

Similar to the case of bi-orthogonal waveform, model (29) can be rewritten in
a matrix form

$\displaystyle\bm{y}={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\bm{X}_{rect}}\bm{c}+\bm{\omega},$
(30)

where $\bm{c}$ is the structured sparse vector to be recovered, and the
$(z,n)$th element of matrix $\bm{X}_{rect}$ is given as

$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{X}_{rect}^{z,n}=e^{j2\pi\frac{l_{z}(d_{n}+\kappa_{n})}{MN}}{X}_{bi}^{z,n}},$
(31)

where $l_{z}$ is the delay index corresponding to the $z$th element of
$\bm{y}$, $d_{n}$ and $\kappa_{n}$ are the parameters $d$ and $\kappa_{d}$ of
the $n$th element of vector $\bm{c}$ (refer to the definition of vector
$\bm{c}$ in (27), and ${X}_{bi}^{z,n}$ denotes the the $(z,n)$th element of
matrix $\bm{X}_{bi}$. It is noted that when $z$ and $n$ are given, $l_{z}$ and
$d_{n}$ are known, but $\kappa_{n}$ is unknown. So it is different from the
case of the bi-orthogonal waveform in (26) that the matrix $\bm{X}_{rect}$
depends on the fractional Doppler shifts, which are unknown. Later we will
show that this can be solved by using an iterative estimation strategy.

As discussed above, the channel estimation problem can be formulated as
structured sparse signal recovery. Next, we will develop a Bayesian method to
recover the block sparse vector $\bm{c}$ and obtain the estimates of
$\\{h_{t,d}\\}$ and $\\{\kappa_{d}\\}$. In particular, the factor graph
techniques [18] are used, based on which an efficient message passing
algorithm is developed.

##  IV Bayesian Approach and Message Passing Algorithm

In this section, we first focus on OTFS with the bi-orthogonal waveform, and
then extend our discussion to the case of the rectangular waveform.

For the convenience of notation, we define $j=(l_{max}+1)(k_{max}+d)+t+1$. As
$t\in[0,l_{max}]$ and $d\in[-k_{max},k_{max}]$, there is one-to-one map
between the index $j$ and the index pair $(t,d)$. Then we define

$\displaystyle\bm{c}_{j}=\bm{c}_{t,d},~{}h_{j}=h_{t,d},~{}\kappa_{j}=\kappa_{d},$
(32)

which will be used hereafter.

Inspired by the sparse Bayesian learning [19], we assume a Gaussian prior with
mean 0 and precision $\lambda_{j}$ for $h_{j}$ to promote sparsity, i.e.,

$\displaystyle p(h_{j}|\lambda_{j})=\mathcal{CN}(h_{j};0,\lambda_{j}^{-1}),$
(33)

where the hyperparameter $\lambda_{j}$ is Gamma distributed, i.e.,
$p(\lambda_{j})=Ga(\lambda_{j};\epsilon_{j},\eta_{j})$. The hyper-parameters
$\epsilon_{j}$ and $\eta_{j}$ are respectively set to be 1 and 0, leading to a
uniform or noninformative distribution for $\lambda_{j}$ [20]. As the noise
precision $\gamma$ is normally unknown, it will also be estimated with an
improper prior $p(\gamma)$[19].

To facilitate the algorithm design, we introduce an auxiliary vector

$\displaystyle\bm{g}_{j}\triangleq$
$\displaystyle[f(-\widehat{N},\kappa_{j}),f(-\widehat{N}+1,\kappa_{j}),\cdots,f(\widehat{N},\kappa_{j})]^{T}e^{-j2\pi\frac{t(d+\kappa_{j})}{MN}}$
$\displaystyle\triangleq$
$\displaystyle[\Phi(-\widehat{N},\kappa_{j}),\cdots,\Phi(q,\kappa_{j}),\cdots,\Phi(\widehat{N},\kappa_{j})]^{T},$
(34)

where

$\Phi(q,\kappa_{j})=f(q,\kappa_{j})e^{-j2\pi\frac{t(d+\kappa_{j})}{MN}}.$ (35)

Then we have

$\displaystyle\bm{c}_{j}=h_{j}\bm{g}_{j}.$ (36)

Define $\bm{h}=[h_{1},\cdots h_{J}]^{T}$,
$\bm{g}=[\bm{g}_{1}^{T},\cdots,\bm{g}_{J}^{T}]^{T}$,
$\bm{\kappa}=[\kappa_{1},\cdots,\kappa_{J}]^{T}$ and
$\bm{\lambda}=[\lambda_{1},\cdots,\lambda_{J}]^{T}$, where
$J=(l_{max}+1)(2k_{max}+1)$. Then the joint conditional distribution of the
unknown variables can be factorized as

$\displaystyle p$
$\displaystyle(\bm{c},\bm{h},\bm{g},\bm{\kappa},\bm{\lambda},\gamma|\bm{y})$
$\displaystyle\propto$
$\displaystyle~{}p(\bm{y}|\bm{c},\gamma)p(\bm{c}|\bm{h},\bm{g})p(\bm{h}|\bm{\lambda})p(\bm{\lambda})p(\bm{g}|\bm{\kappa})p(\bm{\kappa})p(\gamma)$
$\displaystyle=$
$\displaystyle~{}p(\gamma)p(\bm{y}|\bm{c},\gamma){\prod}_{j,b}p(c_{jb}|{h}_{j},g_{jb})p(h_{j}|\lambda_{j})p(\lambda_{j})$
$\displaystyle\qquad\qquad\times p(g_{jb}|\kappa_{j})p(\kappa_{j})$
$\displaystyle\triangleq$
$\displaystyle~{}f_{\gamma}(\gamma)f_{\bm{y}}(\bm{c},\gamma){\prod}_{j,b}f_{c_{jb}}(c_{jb},h_{j},g_{jb})f_{h_{j}}(h_{j},\lambda_{j})f_{\lambda_{j}}(\lambda_{j})$
$\displaystyle\qquad\qquad\times
f_{g_{jb}}(g_{jb},\kappa_{j})f_{\kappa_{j}}(\kappa_{j}),$ (37)

where $1\leq j\leq J$, $1\leq b\leq B$, and $c_{jb}$ and $g_{jb}$ represent
the $b$th element of $\bm{c}_{j}$ and $\bm{g}_{j}$, respectively. The
correspondence between the factors and distributions (and their detailed
forms) are listed in Table II. We aim to obtain the approximate marginals of
the parameters (thereby their estimates) by performing approximate inference.
In particular, we use the factor graph techniques, and the variational message
passing (VMP) [21, 22] and belief propagation (BP) message passing [18] are
combined to achieve efficient approximate inference. We note that variational
inference has been used for channel estimation, e.g., the works in [23] and
[24].

$c_{j1}$$c_{jb}$$c_{jB}$$\bm{\vdots}$$\bm{\vdots}$$c_{J1}$$c_{Jb}$$c_{JB}$$f_{c_{j1}}$$f_{c_{jb}}$$f_{c_{jB}}$$f_{c_{J1}}$$f_{c_{Jb}}$$f_{c_{JB}}$$g_{j1}$$g_{jb}$$g_{jB}$$h_{j}$$g_{J1}$$g_{Jb}$$g_{JB}$$h_{J}$$f_{g_{j1}}$$f_{g_{jb}}$$f_{g_{jB}}$$f_{h_{j}}$$f_{g_{J1}}$$f_{g_{Jb}}$$f_{g_{JB}}$$f_{h_{J}}$$\kappa_{j}$$\lambda_{j}$$\kappa_{J}$$\lambda_{J}$$f_{\kappa_{j}}$$f_{\lambda_{j}}$$f_{\kappa_{J}}$$f_{\lambda_{J}}$$f_{\bm{y}}$$\gamma$$f_{\gamma}$$II$$I$ Figure 4:  Factor graph representation of (37). TABLE II: Correspondence between the factors and distributions in (37). Factor |  Distribution |  Function
---|---|---
$f_{\gamma}$ |  $p(\gamma)$ |  $\gamma^{-1}$
$f_{\bm{y}}$ |  $p(\bm{y}|\bm{c},\gamma)$ |  $\mathcal{CN}(\bm{y};\bm{X}_{bi}\bm{c},\gamma^{-1}\bm{I}_{Z})$
$f_{c_{jb}}$ |  $p(c_{jb}|h_{j},g_{jb})$ |  $\delta(c_{jb}-h_{j}g_{jb})$
$f_{h_{j}}$ |  $p(h_{j}|\lambda_{j})$ |  $\mathcal{CN}(h_{j};0,\lambda_{j}^{-1})$
$f_{\lambda_{j}}$ |  $p(\lambda_{j})$ |  $Ga(\lambda_{j};\epsilon_{j},\eta_{j})$
$f_{g_{jb}}$ |  $p(g_{jb}|\kappa_{j})$ |  $\delta(g_{jb}-\Phi(-\widehat{N}+b-1,\kappa_{j}))$
$f_{\kappa_{j}}$ |  $p(\kappa_{j})$ |  $U[-0.5,0.5]$

The factorization of (37) can be visualized by the factor graph shown in Fig.
4, where we partition the factor graph into two parts: Part I and Part II. To
keep the graph clear, we only show the function nodes and variables nodes
associated with the $j$th and $J$th blocks of $\bm{c}$. In the following, we
derive the message computations for the forward (from left to right) and
backward (from right to left) passing in both parts. We use $m_{A\rightarrow
B}(x)$ to denote a message passed from a function node $A$ to a variable node
$B$, which is a function of $x$, and use $n_{B\rightarrow A}(x)$ to denote a
message passed from a variable node $B$ to a function node $A$ , which is also
a function of $x$. Meanwhile, the arrows above the mean and the variance of a
Gaussian message indicate the direction of the message passing. In addition,
we use $b(x)$ to denote the belief of a variable $x$. Note that, if a forward
computation requires backward messages, the relevant messages in the previous
iteration is used by default.

###  IV-A Message Computations in Part I

Assume that the belief of $\bm{c}$ is known, which turns out to be Gaussian,
i.e., $b(\bm{c})\propto\mathcal{CN}(\bm{c};\bm{c}^{p},\bm{V}_{c}^{p})$, as
given later in (43) and (44). The VMP or mean field (MF) is used at the
function node $f_{\bm{y}}(\bm{c},\gamma)$, and the message from
$f_{\bm{y}}(\bm{c},\gamma)$ to the variable node $\gamma$ can be computed as

$\displaystyle m_{f_{\bm{y}}\rightarrow\gamma}(\gamma)$
$\displaystyle\propto\exp\int\ln f_{\bm{y}}(\bm{c},\gamma)b(\bm{c})d\bm{c}$
$\displaystyle\propto{\gamma^{Z}}\exp\\{-\gamma[(\bm{y}-\bm{X}_{bi}\bm{c}^{p})^{H}(\bm{y}-\bm{X}_{bi}\bm{c}^{p})$
$\displaystyle\qquad\qquad+Tr\\{\bm{X}_{bi}\bm{V_{c}}^{p}\bm{X}_{bi}^{H}\\}]\\}.$
(38)

With the prior $f_{\gamma}(\gamma)\propto 1/\gamma$, the belief $b(\gamma)$
can be expressed as

$\displaystyle
b(\gamma)\propto\gamma^{-1}m_{f_{\bm{y}}\rightarrow\gamma}(\gamma).$ (39)

Then, with the MF rule and the belief $b(\gamma)$, we can compute the outgoing
message $m_{f_{\bm{y}}\rightarrow{c}_{jb}}(c_{jb})$ as

$\displaystyle m_{f_{\bm{y}}\rightarrow{c}_{jb}}(c_{jb})$
$\displaystyle\propto\exp\int\ln
f_{\bm{y}}(\bm{c},\gamma)b(\gamma)b(\bm{c}_{\sim c_{jb}})d\gamma d\bm{c}_{\sim
c_{jb}},$ (40)

where $\bm{c}_{\sim c_{jb}}$ denotes a vector obtained by removing $c_{jb}$
from $\bm{c}$. However, the computation of the message involves high
complexity. An efficient way is to first compute the belief of $c_{jb}$, then
compute the outgoing extrinsic message. According to the MF rule, the message
$m_{f_{\bm{y}}\rightarrow\bm{c}}(\bm{c})$ can be expressed as

$\displaystyle m_{f_{\bm{y}}\rightarrow\bm{c}}(\bm{c})$
$\displaystyle\propto\exp\int\ln f_{\bm{y}}(\bm{c},\gamma)b(\gamma)d\gamma$
$\displaystyle\propto~{}\mathcal{CN}\left(\bm{c};(\bm{X}_{bi}^{H}\bm{X}_{bi})^{-1}\bm{X}_{bi}^{H}\bm{y},\hat{\gamma}^{-1}(\bm{X}_{bi}^{H}\bm{X}_{bi})^{-1}\right),$
(41)

where

$\displaystyle\hat{\gamma}$ $\displaystyle=\int\gamma b(\gamma)d\gamma$
$\displaystyle=\frac{Z}{(\bm{y}-\bm{X}_{bi}\bm{c}^{p})^{H}(\bm{y}-\bm{X}_{bi}{c}^{p})+Tr\\{\bm{X}_{bi}{V_{c}}^{p}\bm{X}_{bi}^{H}\\}}.$
(42)

As the incoming message $m_{f_{c_{jb}}\rightarrow
c_{jb}}(c_{jb})\propto\mathcal{CN}(c_{jb};\overset{\scriptscriptstyle\leftarrow}{c_{jb}},\overset{\scriptscriptstyle\leftarrow}{\nu}_{c_{jb}})$,
the belief of $\bm{c}$ is Gaussian with covariance matrix $\bm{V}_{c}^{p}$ and
mean vector $\bm{c}^{p}$, which can be computed as

$\displaystyle\bm{V}_{c}^{p}=$
$\displaystyle(\bm{V}_{c}^{-1}+\hat{\gamma}\bm{X}_{bi}^{H}\bm{X}_{bi})^{-1},$
(43) $\displaystyle\bm{c}^{p}=$
$\displaystyle\bm{V}^{p}_{c}(\bm{V}_{c}^{-1}\overset{\scriptscriptstyle\leftarrow}{\bm{c}}+\hat{\gamma}\bm{X}_{bi}^{H}\bm{y}),$
(44)

where $\bm{V}_{\bm{c}}$ is a diagonal matrix with the diagonal elements given
by $\\{\overset{\scriptscriptstyle\leftarrow}{\nu}_{c_{jb}}\\}$ and
$\overset{\scriptscriptstyle\leftarrow}{\bm{c}}$ is a column vector that
consists of $\\{\overset{\scriptscriptstyle\leftarrow}{c}_{jb}\\}$. The
computations of the mean $\overset{\scriptscriptstyle\leftarrow}{c}_{jb}$ and
the variance $\overset{\scriptscriptstyle\leftarrow}{\nu}_{c_{jb}}$ are
delayed to (90) and (91). Hence the outgoing message
$m_{f_{\bm{y}}\rightarrow{c}_{jb}}(c_{jb})$ is also Gaussian, and can be
expressed as [25]

$\displaystyle
m_{f_{\bm{y}}\rightarrow{c}_{jb}}(c_{jb})=\mathcal{CN}(c_{jb};\overset{\scriptscriptstyle\rightarrow}{c}_{jb},\overset{\scriptscriptstyle\rightarrow}{\nu}_{c_{jb}})$
(45)

with

$\displaystyle\overset{\scriptscriptstyle\rightarrow}{\nu}_{c_{jb}}=$
$\displaystyle\left(\frac{1}{\nu_{c_{jb}}^{p}}-\frac{1}{\overset{\scriptscriptstyle\leftarrow}{\nu}_{c_{jb}}}\right)^{-1},$
(46) $\displaystyle\overset{\scriptscriptstyle\rightarrow}{c}_{jb}=$
$\displaystyle\overset{\scriptscriptstyle\rightarrow}{\nu}_{c_{jb}}\left(\frac{{c}_{jb}^{p}}{\nu_{c_{jb}}^{p}}-\frac{\overset{\scriptscriptstyle\leftarrow}{c}_{jb}}{\overset{\scriptscriptstyle\leftarrow}{\nu}_{c_{jb}}}\right),$
(47)

where ${c}_{jb}^{p}$ is the $b$th element of the $j$th block of
${\bm{c}}^{p}$, and $\nu_{c_{jb}}^{p}$ is the $b$th element of the $j$th block
of the vector that consists of the diagonal elements of $\bm{V}_{\bm{c}}^{p}$.

###  IV-B Message Computations in Part II

####  IV-B1 Forward Message Passing

With the incoming message from Part I $n_{c_{jb}\rightarrow
f_{c_{jb}}}(c_{jb})=m_{f_{\bm{y}}\rightarrow{c}_{jb}}(c_{jb})=\mathcal{CN}(c_{jb};\overset{\scriptscriptstyle\rightarrow}{c}_{jb},\overset{\scriptscriptstyle\rightarrow}{\nu}_{c_{jb}})$
and the factor $f_{c_{jb}}(c_{jb},h_{j},g_{jb})=\delta(c_{jb}-h_{j}g_{jb})$,
we can obtain an intermediate function node $\tilde{f}_{c_{jb}}(h_{j},g_{jb})$
according to BP [18], i.e.,

$\displaystyle\tilde{f}_{c_{jb}}(h_{j},g_{jb})=$ $\displaystyle\int
f_{c_{jb}}(c_{jb},h_{j},g_{jb})n_{c_{jb}\rightarrow
f_{c_{jb}}}(c_{jb})d{c_{jb}}$ $\displaystyle=$
$\displaystyle\mathcal{CN}(h_{j}g_{jb};\overset{\scriptscriptstyle\rightarrow}{c}_{jb},\overset{\scriptscriptstyle\rightarrow}{\nu}_{c_{jb}}).$
(48)

With the local function $\tilde{f}_{c_{jb}}(h_{j},g_{jb})$, we can compute the
message from $f_{c_{jb}}$ to $h_{j}$ by treating $g_{jb}$ as a constant, i.e.,

$\displaystyle m_{f_{c_{jb}}\rightarrow
h_{j}}(h_{j})=\mathcal{CN}(h_{j};\overset{\scriptscriptstyle\rightarrow}{h}_{jb},\overset{\scriptscriptstyle\rightarrow}{\nu}_{h_{jb}}),$
(49)

where

$\displaystyle\overset{\scriptscriptstyle\rightarrow}{h}_{jb}=$
$\displaystyle\frac{\overset{\scriptscriptstyle\rightarrow}{c}_{jb}}{\hat{g}_{jb}},$
(50) $\displaystyle\overset{\scriptscriptstyle\rightarrow}{\nu}_{h_{jb}}=$
$\displaystyle\frac{\overset{\scriptscriptstyle\rightarrow}{\nu}_{c_{jb}}}{|\hat{g}_{jb}|^{2}}$
(51)

with $\hat{g}_{jb}$ being the mean of the Gaussian belief of $g_{jb}$, which
is computed in (89). The product of the Gaussian messages
$\\{m_{f_{c_{jb}}\rightarrow h_{j}}(h_{j}),\forall b\\}$ is still Gaussian
[25], i.e.,

$\displaystyle q_{j}(h_{j})=\prod_{b}m_{f_{c_{jb}}\rightarrow
h_{j}}(h_{j})\propto\mathcal{CN}(h_{j};\hat{q}_{j},\nu_{q_{j}}),$ (52)

where

$\displaystyle{\nu}_{q_{j}}=$
$\displaystyle\left(\sum_{b}\frac{1}{\overset{\scriptscriptstyle\rightarrow}{\nu}_{h_{jb}}}\right)^{-1},$
(53) $\displaystyle\hat{q}_{j}=$
$\displaystyle{\nu}_{q_{j}}\sum_{b}\frac{\overset{\scriptscriptstyle\rightarrow}{h}_{jb}}{\overset{\scriptscriptstyle\rightarrow}{\nu}_{h_{jb}}}.$
(54)

With the message $m_{f_{h_{j}}\rightarrow
h_{j}}(h_{j})\propto\mathcal{CN}(h_{j};0,\widehat{\lambda}_{j}^{-1})$, which
is given later in (77), the belief $b(h_{j})$ of $h_{j}$ is obtained as

$\displaystyle b(h_{j})=q_{j}(h_{j})m_{f_{h_{j}}\rightarrow
h_{j}}(h_{j})\propto\mathcal{CN}(h_{j};\hat{h}_{j},\nu_{h_{j}}),$ (55)

where

$\displaystyle\hat{h}_{j}=$
$\displaystyle\frac{\hat{q}_{j}}{1+{\nu}_{q_{j}}\widehat{\lambda}_{j}},$ (56)
$\displaystyle\nu_{h_{j}}=$
$\displaystyle\left(\frac{1}{{\nu}_{q_{j}}}+\widehat{\lambda}_{j}\right)^{-1}.$
(57)

The message $m_{f_{h_{j}}\rightarrow\lambda_{j}}(\lambda_{j})$ is then
computed by using the MF rule, i.e.,

$\displaystyle m_{f_{h_{j}}\rightarrow\lambda_{j}}(\lambda_{j})\propto$
$\displaystyle exp\left\\{\int\ln
f_{h_{j}}(h_{j},\lambda_{j})b(h_{j})d{h_{j}}\right\\}$ $\displaystyle\propto$
$\displaystyle\lambda_{j}exp\\{-\lambda_{j}(|\hat{h}_{j}|^{2}+\nu_{h_{j}})\\},$
(58)

so the belief $b(\lambda_{j})$ of the hyperparameter $\lambda_{j}$ is given as

$\displaystyle b(\lambda_{j})=$ $\displaystyle
m_{f_{h_{j}}\rightarrow\lambda_{j}}(\lambda_{j})f_{\lambda_{j}}(\lambda_{j})$
$\displaystyle\propto$
$\displaystyle\lambda_{j}^{\epsilon_{j}}exp\left\\{-\lambda_{j}(\eta_{j}+|\hat{h}_{j}|^{2}+\nu_{h_{j}})\right\\}.$
(59)

The message from $f_{c_{jb}}$ to $g_{jb}$ is computed with the MF rule, i.e.,

$\displaystyle m_{f_{c_{jb}}\rightarrow g_{jb}}(g_{jb})\propto$ $\displaystyle
exp\left\\{\int\ln\tilde{f}_{c_{jb}}(h_{j},g_{jb})b(h_{j})d{h_{j}}\right\\}$
$\displaystyle=$
$\displaystyle\mathcal{CN}(g_{jb};\overset{\scriptscriptstyle\rightarrow}{g}_{jb},\overset{\scriptscriptstyle\rightarrow}{\nu}_{g_{jb}}),$
(60)

where

$\displaystyle\overset{\scriptscriptstyle\rightarrow}{g}_{jb}=$
$\displaystyle\frac{\overset{\scriptscriptstyle\rightarrow}{c}_{jb}\hat{h}_{j}^{*}}{|\hat{h}_{j}|^{2}+\nu_{h_{j}}},$
(61) $\displaystyle\overset{\scriptscriptstyle\rightarrow}{\nu}_{g_{jb}}=$
$\displaystyle\frac{\overset{\scriptscriptstyle\rightarrow}{\nu}_{c_{jb}}}{|\hat{h}_{j}|^{2}+\nu_{h_{j}}}.$
(62)

It is noted that the local function node $f_{g_{j}b}$ includes a nonlinear
function $\Phi(q,\kappa_{j})$, which makes the message computation about
$\kappa_{j}$ intractable. Inspired by the extended Kalman filter [26] , to
solve this problem, $\Phi(q,\kappa_{j})$ is linearized by using the first
order Taylor expansion with the estimate of $\kappa_{j}$ in last iteration,
i.e.,

$\displaystyle\Phi(q,\kappa_{j})\approx\Phi(q,\hat{\kappa}_{j}^{{}^{\prime}})+\Phi^{\prime}(q,\hat{\kappa}_{j}^{{}^{\prime}})(\kappa_{j}-\hat{\kappa}_{j}^{{}^{\prime}}),$
(63)

with

$\displaystyle\Phi^{\prime}(q,\hat{\kappa}_{j}^{{}^{\prime}})=$
$\displaystyle\left(\frac{-j2\pi
t}{MN}\right)\Phi(q,\hat{\kappa}_{j}^{{}^{\prime}})+$ $\displaystyle
e^{-j2\pi\frac{t(d+\hat{\kappa}_{j}^{{}^{\prime}})}{MN}}\frac{1}{N}\sum_{n=1}^{N-1}j\frac{2n\pi}{N}e^{j\frac{2n\pi}{N}(q+\hat{\kappa}_{j}^{{}^{\prime}})},$
(64)

where $\hat{\kappa}_{j}^{{}^{\prime}}$ denotes the estimates of $\kappa_{j}$
in last iteration. Then, with the BP rule and the approximation in (63), the
message $m_{f_{g_{jb}}\rightarrow\kappa_{j}}(\kappa_{j})$ can be expressed as

$\displaystyle m_{f_{g_{jb}}\rightarrow\kappa_{j}}(\kappa_{j})=$
$\displaystyle\int f_{g_{jb}}(g_{jb},\kappa_{j})m_{f_{c_{jb}}\rightarrow
g_{jb}}(g_{jb})d{g_{jb}}.$ (65)

Note that $\kappa_{j}$ is a real valued variable. To ensure that the message
$m_{f_{g_{jb}}\rightarrow\kappa_{j}}(\kappa_{j})$ is real, we rewrite the
function $f_{g_{jb}}(g_{jb},\kappa_{j})$ as

$\displaystyle f_{g_{jb}}(g_{jb},\kappa_{j})$
$\displaystyle=\delta\left(\mathcal{R}[g_{jb}]-\mathcal{R}[\Phi(Q,\kappa_{j})]\right)\delta\left(\mathcal{I}[g_{jb}]-\mathcal{I}[\Phi(Q,\kappa_{j})]\right),$
(66)

where $Q=-\widehat{N}+b-1$. Then the message
$m_{f_{g_{jb}}\rightarrow\kappa_{j}}(\kappa_{j})$ in (65) can be obtained as

$\displaystyle m_{f_{g_{jb}}\rightarrow\kappa_{j}}(\kappa_{j})\propto$
$\displaystyle\mathcal{N}(\kappa_{j};\overset{\scriptscriptstyle\rightarrow}{\kappa}_{jb}^{\mathcal{R}},\overset{\scriptscriptstyle\rightarrow}{\nu}_{\kappa_{jb}}^{\mathcal{R}})\mathcal{N}(\kappa_{j};\overset{\scriptscriptstyle\rightarrow}{\kappa}_{jb}^{\mathcal{I}},\overset{\scriptscriptstyle\rightarrow}{\nu}_{\kappa_{jb}}^{\mathcal{I}})$
$\displaystyle\propto$
$\displaystyle\mathcal{N}(\kappa_{j};\overset{\scriptscriptstyle\rightarrow}{\kappa}_{jb},\overset{\scriptscriptstyle\rightarrow}{\nu}_{\kappa_{jb}}),$
(67)

where

$\displaystyle\overset{\scriptscriptstyle\rightarrow}{\kappa}_{jb}^{\mathcal{R}}=$
$\displaystyle\frac{\mathcal{R}[\overset{\scriptscriptstyle\rightarrow}{g}_{jb}]-\mathcal{R}[\Phi(Q,\hat{\kappa}_{j}^{{}^{\prime}})]+\mathcal{R}[\Phi^{\prime}(Q,\hat{\kappa}_{j}^{{}^{\prime}})]\hat{\kappa}_{j}^{{}^{\prime}}}{\mathcal{R}[\Phi^{\prime}(Q,\hat{\kappa}_{j}^{{}^{\prime}})]},$
(68)
$\displaystyle\overset{\scriptscriptstyle\rightarrow}{\nu}_{\kappa_{jb}}^{\mathcal{R}}=$
$\displaystyle\frac{\overset{\scriptscriptstyle\rightarrow}{\nu}_{g_{jb}}}{2|\mathcal{R}[\Phi^{\prime}(Q,\hat{\kappa}_{j}^{{}^{\prime}})]|^{2}},$
(69)
$\displaystyle\overset{\scriptscriptstyle\rightarrow}{\kappa}_{jb}^{\mathcal{I}}=$
$\displaystyle\frac{\mathcal{I}[\overset{\scriptscriptstyle\rightarrow}{g}_{jb}]-\mathcal{I}[\Phi(Q,\hat{\kappa}_{j}^{{}^{\prime}})]+\mathcal{I}[\Phi^{\prime}(Q,\hat{\kappa}_{j}^{{}^{\prime}})]\hat{\kappa}_{j}^{{}^{\prime}}}{\mathcal{I}[\Phi^{\prime}(Q,\hat{\kappa}_{j}^{{}^{\prime}})]},$
(70)
$\displaystyle\overset{\scriptscriptstyle\rightarrow}{\nu}_{\kappa_{jb}}^{\mathcal{I}}=$
$\displaystyle\frac{\overset{\scriptscriptstyle\rightarrow}{\nu}_{g_{jb}}}{2|\mathcal{I}[\Phi^{\prime}(Q,\hat{\kappa}_{j}^{{}^{\prime}})]|^{2}},$
(71)

and

$\displaystyle\overset{\scriptscriptstyle\rightarrow}{\nu}_{\kappa_{jb}}=$
$\displaystyle\left(\frac{1}{\overset{\scriptscriptstyle\rightarrow}{\nu}_{\kappa_{jb}}^{\mathcal{R}}}+\frac{1}{\overset{\scriptscriptstyle\rightarrow}{\nu}_{\kappa_{jb}}^{\mathcal{I}}}\right)^{-1},$
(72) $\displaystyle\overset{\scriptscriptstyle\rightarrow}{\kappa}_{jb}=$
$\displaystyle\overset{\scriptscriptstyle\rightarrow}{\nu}_{\kappa_{jb}}\left(\frac{\overset{\scriptscriptstyle\rightarrow}{\kappa}_{jb}^{\mathcal{R}}}{\overset{\scriptscriptstyle\rightarrow}{\nu}_{\kappa_{jb}}^{\mathcal{R}}}+\frac{\overset{\scriptscriptstyle\rightarrow}{\kappa}_{jb}^{\mathcal{I}}}{\overset{\scriptscriptstyle\rightarrow}{\nu}_{\kappa_{jb}}^{\mathcal{I}}}\right).$
(73)

The prior of $\kappa_{j}$ is a uniform distribution over $[-0.5,0.5]$, i.e.,
$-0.5\leq\kappa_{j}\leq 0.5$. To simplify the computation of the belief
$b(\kappa_{j})$, we simply carry out the clipping operation, i.e.,

$b(\kappa_{j})\propto\begin{cases}\mathcal{N}(\kappa_{j};-0.5,0);&\hat{\kappa}_{j}\leq-0.5\\\
\mathcal{N}(\kappa_{j};0.5,0);&\hat{\kappa}_{j}\geq 0.5\\\
\mathcal{N}(\kappa_{j};\hat{\kappa}_{j},\nu_{\kappa_{j}});&otherwise\\\
\end{cases}$ (74)

where

$\nu_{\kappa_{j}}=\left(\sum_{b}\frac{1}{\overset{\scriptscriptstyle\rightarrow}{\nu}_{\kappa_{jb}}}\right)^{-1},$
(75)
$\hat{\kappa}_{j}=\nu_{\kappa_{j}}\sum_{b}\frac{\overset{\scriptscriptstyle\rightarrow}{\kappa}_{jb}}{\overset{\scriptscriptstyle\rightarrow}{\nu}_{\kappa_{jb}}}.$
(76)

####  IV-B2 Backward Message Passing

We firstly compute the message $m_{f_{h_{j}}\rightarrow h_{j}}(h_{j})$ from
the function node $f_{h_{j}}$ to $h_{j}$ by using the MF rule as follows

$\displaystyle m_{f_{h_{j}}\rightarrow h_{j}}(h_{j})\propto$ $\displaystyle
exp\left\\{\int\ln
f_{h_{j}}(h_{j},\lambda_{j})b(\lambda_{j})d{\lambda_{j}}\right\\}$
$\displaystyle\propto$
$\displaystyle\mathcal{CN}(h_{j};0,\hat{\lambda}_{j}^{-1}),$ (77)

where

$\displaystyle\hat{\lambda}_{j}=\int\lambda_{j}b(\lambda_{j})d\lambda_{j}=\frac{\epsilon_{j}+1}{\eta_{j}+|\hat{h}_{j}|^{2}+\nu_{h_{j}}}.$
(78)

Then the message $n_{h_{j}\rightarrow f_{c_{jb}}}(h_{j})$ from variable node
$h_{j}$ to $f_{c_{jb}}$ is updated by the BP rule, i.e.,

$\displaystyle n_{h_{j}\rightarrow
f_{c_{jb}}}(h_{j})=\frac{b(h_{j})}{m_{f_{c_{jb}}\rightarrow
h_{j}}}\propto\mathcal{CN}(h_{j};\overset{\scriptscriptstyle\leftarrow}{h}_{jb},\overset{\scriptscriptstyle\leftarrow}{\nu}_{h_{jb}}),$
(79)

where

$\displaystyle\overset{\scriptscriptstyle\leftarrow}{\nu}_{h_{jb}}=$
$\displaystyle\left({\frac{1}{\nu_{h_{j}}}}-\frac{1}{\overset{\scriptscriptstyle\rightarrow}{\nu}_{h_{jb}}}\right)^{-1},$
(80) $\displaystyle\overset{\scriptscriptstyle\leftarrow}{h}_{jb}=$
$\displaystyle\overset{\scriptscriptstyle\leftarrow}{\nu}_{h_{jb}}\left(\frac{\hat{h}_{j}}{\nu_{h_{j}}}-\frac{\overset{\scriptscriptstyle\rightarrow}{h}_{jb}}{\overset{\scriptscriptstyle\rightarrow}{\nu}_{h_{jb}}}\right).$
(81)

The backward message $n_{\kappa_{j}\rightarrow f_{g_{jb}}}(\kappa_{j})$ is
Gaussian with mean $\overset{\scriptscriptstyle\leftarrow}{\kappa}_{jb}$ and
variance $\overset{\scriptscriptstyle\leftarrow}{\nu}_{\kappa_{jb}}$ and can
be calculated as

$\displaystyle\overset{\scriptscriptstyle\leftarrow}{\nu}_{\kappa_{jb}}=$
$\displaystyle(1/\nu_{\kappa_{j}}-1/\overset{\scriptscriptstyle\rightarrow}{\nu}_{\kappa_{jb}})^{-1},$
(82) $\displaystyle\overset{\scriptscriptstyle\leftarrow}{\kappa}_{jb}=$
$\displaystyle\overset{\scriptscriptstyle\leftarrow}{\nu}_{\kappa_{jb}}(\hat{\kappa}_{j}/\nu_{\kappa_{j}}-\overset{\scriptscriptstyle\rightarrow}{\kappa}_{jb}/\overset{\scriptscriptstyle\rightarrow}{\nu}_{\kappa_{jb}}).$
(83)

For the cases that $\hat{\kappa}_{j}=0.5$ or $-0.5$ and $\nu_{\kappa_{j}}=0$
in (74), we set

$\displaystyle\overset{\scriptscriptstyle\leftarrow}{\nu}_{\kappa_{jb}}=\nu_{\kappa_{j}},\quad\overset{\scriptscriptstyle\leftarrow}{\kappa}_{jb}=\hat{\kappa}_{j}.$
(84)

Then the message $m_{f_{g_{jb}}\rightarrow g_{jb}}(g_{jb})$ is calculated by
using the BP rule, i.e.,

$\displaystyle m_{f_{g_{jb}}\rightarrow g_{jb}}(g_{jb})=$
$\displaystyle\int{f_{g_{jb}}(g_{jb},\kappa_{j})n_{\kappa_{j}\rightarrow
f_{g_{jb}}}(\kappa_{j})}d{\kappa_{j}}$ $\displaystyle\propto$
$\displaystyle\mathcal{CN}(g_{jb};\overset{\scriptscriptstyle\leftarrow}{g}_{jb},\overset{\scriptscriptstyle\leftarrow}{\nu}_{g_{jb}}),$
(85)

where

$\displaystyle\overset{\scriptscriptstyle\leftarrow}{g}_{jb}=$
$\displaystyle\Phi(Q,\hat{\kappa}_{j}^{{}^{\prime}})+\Phi^{\prime}(Q,\hat{\kappa}_{j}^{{}^{\prime}})(\overset{\scriptscriptstyle\leftarrow}{\kappa}_{jb}-\hat{\kappa}_{j}^{{}^{\prime}}),$
(86) $\displaystyle\overset{\scriptscriptstyle\leftarrow}{\nu}_{g_{jb}}=$
$\displaystyle\overset{\scriptscriptstyle\leftarrow}{\nu}_{\kappa_{jb}}|\Phi^{\prime}(Q,\hat{\kappa}_{j}^{{}^{\prime}})|^{2}.$
(87)

Hence we can obtain the belief
$b(g_{jb})\propto\mathcal{CN}(g_{jb};\hat{g}_{jb},\nu_{g_{jb}})$ with

$\nu_{g_{jb}}=\left(\frac{1}{\overset{\scriptscriptstyle\leftarrow}{\nu}_{g_{jb}}}+\frac{1}{\overset{\scriptscriptstyle\rightarrow}{\nu}_{g_{jb}}}\right)^{-1},\\\
$ (88)
$\hat{g}_{jb}=\nu_{g_{jb}}\left(\frac{\overset{\scriptscriptstyle\leftarrow}{g}_{jb}}{\overset{\scriptscriptstyle\leftarrow}{\nu}_{g_{jb}}}+\frac{\overset{\scriptscriptstyle\rightarrow}{g}_{jb}}{\overset{\scriptscriptstyle\rightarrow}{\nu}_{g_{jb}}}\right).$
(89)

Finally, by combining the incoming message $n_{g_{jb}\rightarrow
f_{c_{jb}}}(g_{jb})=m_{f_{g_{jb}}\rightarrow g_{jb}}(g_{jb})$ and
$n_{h_{j}\rightarrow
f_{c_{jb}}}(h_{j})\propto\mathcal{CN}(h_{j};\overset{\scriptscriptstyle\leftarrow}{h}_{jb},\overset{\scriptscriptstyle\leftarrow}{\nu}_{h_{jb}})$,
the message $m_{f_{c_{jb}}\rightarrow c_{jb}}(c_{jb})$ is Gaussian with mean
$\overset{\scriptscriptstyle\leftarrow}{c}_{jb}$ and variance
$\overset{\scriptscriptstyle\leftarrow}{\nu}_{c_{jb}}$, which are computed as
[27]

$\displaystyle\overset{\scriptscriptstyle\leftarrow}{c}_{jb}=$
$\displaystyle\overset{\scriptscriptstyle\leftarrow}{h}_{jb}\overset{\scriptscriptstyle\leftarrow}{g}_{jb},$
(90) $\displaystyle\overset{\scriptscriptstyle\leftarrow}{\nu}_{c_{jb}}=$
$\displaystyle|\overset{\scriptscriptstyle\leftarrow}{h}_{jb}|^{2}+|\overset{\scriptscriptstyle\leftarrow}{g}_{jb}|^{2}\overset{\scriptscriptstyle\leftarrow}{\nu}_{h_{jb}}+\overset{\scriptscriptstyle\leftarrow}{\nu}_{h_{jb}}.$
(91)

The message passing algorithm is summarized in Algorithm 1. The algorithm can
be terminated when it reaches a maximum number of iteration or the difference
between the estimates of parameters of two consecutive iterations is less than
a threshold.

###  IV-C Extension to Rectangular Waveform

The derivation of Algorithm 1 in the above is for OTFS with the bi-orthogonal
waveform. Thanks to the iterative process of the algorithm, it can be readily
extended for OTFS with the rectangular waveform. As shown in (31), the
difference between $\bm{X}_{bi}$ and $\bm{X}_{rect}$ is that there is a
Doppler shift-dependent factor for each of the elements in $\bm{X}_{rect}$,
and the parameter $\kappa_{n}$ is unknown. Here we note that when $z$ and $n$
are given, the values of $l_{z}$ and $d_{n}$ are known.

This problem can be solved by taking advantage of the iterative estimation
strategy. We can start the iterative process with initialization
$\kappa_{n}=0$ and treat ${\bm{X}}_{rect}$ as a known matrix with the
estimated Doppler shifts in last iteration plugged in (31) (here we abuse the
use of the notation ${\bm{X}}_{rect}$, which is actually an estimate of the
true ${\bm{X}}_{rect}$ with the estimate of $\kappa_{n}$). Then the algorithm
developed for the OTFS with the bi-orthogonal waveform can be used to recover
$\bm{c}$, i.e., $\kappa_{d}$ is estimated, so that $\kappa_{n}$ can be found
based on the corresponding $\kappa_{d}$. Then with the estimated $\kappa_{n}$,
${\bm{X}}_{rect}$ is updated for the next round iteration. Hence, we only need
to add an extra step after step 11 of Algorithm 1 to update ${\bm{X}}_{rect}$
for OTFS with the rectangular waveform.

Algorithm 1 Message Passing Algorithm for OTFS Channel Estimation

Initialize $\overset{\scriptscriptstyle\leftarrow}{c}_{jb}=0$,
$\overset{\scriptscriptstyle\leftarrow}{\nu}_{c_{jb}}=1$,
$\hat{\lambda}_{j}=1$ for $\forall j,b$; $\hat{\gamma}=1$, and $t=1$.
Repeat

1: update $\bm{V}_{c}^{p}$ and $\bm{c}^{p}$ with (43) and (44);

2: update noise precision $\hat{\gamma}$ with (42);

3: update the $\overset{\scriptscriptstyle\rightarrow}{c}_{jb}$ and
$\overset{\scriptscriptstyle\rightarrow}{\nu}_{c_{jb}}$ with (47) and (46);

4: update $\overset{\scriptscriptstyle\rightarrow}{h}_{jb}$ and
$\overset{\scriptscriptstyle\rightarrow}{\nu}_{h_{jb}}$ with (50) and (51);

5: update the belief $b(h_{j})$ of $h_{j}$ with (56) and (57);

6: update $\hat{\lambda}_{j}$ with (78);

7: update the belief $b(h_{j})$ of $h_{j}$ with (56) and (57);

8: update $\overset{\scriptscriptstyle\leftarrow}{h}_{jb}$ and
$\overset{\scriptscriptstyle\leftarrow}{\nu}_{h_{jb}}$ with (81) and (80);

9: update $\overset{\scriptscriptstyle\rightarrow}{g}_{jb}$ and
$\overset{\scriptscriptstyle\rightarrow}{\nu}_{g_{jb}}$ with (61) and (62);

10: update $\overset{\scriptscriptstyle\rightarrow}{\kappa}_{jb}$ and
$\overset{\scriptscriptstyle\rightarrow}{\nu}_{\kappa_{jb}}$ with (72) and
(73);

11: update the belief $b(\kappa_{j})$ of $\kappa_{j}$ with (74); update
${\bm{X}}_{rect}$ with estimated Doppler shifts plugged in (31) in the case of
rectangular waveform;

12: update $\overset{\scriptscriptstyle\leftarrow}{\nu}_{\kappa_{jb}}$ and
$\overset{\scriptscriptstyle\leftarrow}{\kappa}_{jb}$ with (82) and (83);

13: update $\overset{\scriptscriptstyle\leftarrow}{g}_{jb}$ and
$\overset{\scriptscriptstyle\leftarrow}{\nu}_{g_{jb}}$ with (87) and (86);

14: update the belief $b(g_{jb})$ of $g_{jb}$ with (88) and (89);

15: update $\overset{\scriptscriptstyle\leftarrow}{c}_{jb}$ and
$\overset{\scriptscriptstyle\leftarrow}{\nu}_{c_{jb}}$ with (90) and (91);

16: $t=t+1$.

Until terminated

##  V Cramer-Rao Lower Bound

To evaluate the performance of the proposed algorithm, we derive the CRLB for
the channel gain and Doppler shift estimation. We assume that the locations of
the nonzero elements in $\bm{c}$ are known to facilitate the derivation of the
CRLB. Hence, the CRLB derived in this section is a loose performance bound for
the estimation.

We firstly rewrite (9) and (29) as

$\displaystyle y[k,l]=u[k,l]+\omega[k,l],$ (92)

where

$\displaystyle
u[k,l]=\sum_{i=1}^{P}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\sum_{q=-\hat{N}}^{\hat{N}}}$
$\displaystyle
h_{i}f(q,\kappa_{i})e^{-j2\pi\frac{l_{i}(k_{i}+\kappa_{i})}{MN}}$
$\displaystyle\times x([k-k_{i}+q]_{N},[l-l_{i}]_{M})$ (93)

for the bi-orthogonal waveform, and

$\displaystyle
u[k,l]=\sum_{i=1}^{P}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\sum_{q=-\hat{N}}^{\hat{N}}}$
$\displaystyle
h_{i}f(q,\kappa_{i})e^{j2\pi\frac{(l-l_{i})(k_{i}+\kappa_{i})}{MN}}$
$\displaystyle\times x([k-k_{i}+q]_{N},[l-l_{i}]_{M})$ (94)

for the rectangular waveform. In order to derive the CRLB, we define
$\bm{\theta}=[h_{1},h_{p},\cdots,h_{P},\kappa_{1},\kappa_{p},\cdots,\kappa_{P}]^{T}$.
According to [28], the CRLB of the $j$th element in $\bm{\theta}$ is the $j$th
diagonal element of the inverse of the Fisher information matrix, i.e.,

$\displaystyle\theta_{j}^{CRLB}=[\bm{I}^{-1}(\bm{\theta})]_{jj},$ (95)

where the Fisher information matrix $\bm{I}(\bm{\theta})$ has a size of
$2P\times 2P$ with the $(i,j)$th element given as

$\displaystyle[\bm{I}(\bm{\theta})]_{ij}=-\mathbb{E}\left[\frac{\partial^{2}\ln
p(\bm{y};\bm{\theta})}{\partial{\theta}_{i}\partial{\theta}_{j}}\right]$ (96)

for $i=1,2,\cdots,2P$, $j=1,2,\cdots,2P$, and the expectation is taken with
respect to $p(\bm{y};\bm{\theta})$. The logarithm of the likelihood function
$\ln p(\bm{y};\bm{\theta})$ can be expressed as

$\displaystyle\ln
p(\bm{y};\bm{\theta})=-Z\ln(\pi\gamma^{-1})-\gamma\sum_{z=0}^{Z-1}|y_{z}-u_{z}|^{2},$
(97)

where $y_{z}$ denotes the $z$th element in $\bm{y}$ and $u_{z}$ denotes the
corresponding $u[k,l]$. Then, (96) can be rewritten as

$\displaystyle[\bm{I}(\bm{\theta})]_{ij}=\gamma\sum_{z=0}^{Z-1}\left[-\frac{\partial
u_{z}}{\partial{\theta}_{i}}\frac{\partial
u_{z}^{*}}{\partial{\theta}_{j}}-\frac{\partial
u_{z}^{*}}{\partial{\theta}_{i}}\frac{\partial
u_{z}}{\partial{\theta}_{j}}\right].$ (98)

For simplicity, we define $k^{\prime}\triangleq[k-k_{p}+q]_{N}$ and
$l^{\prime}\triangleq[l-l_{p}]_{M}$. When $1\leq p\leq P$ ($p$ represents $i$
or $j$), the derivative $\frac{\partial{u}_{z}}{\partial{\theta}_{p}}$ is
about the channel gain $h_{p}$. For the bi-orthogonal waveform, we have

$\displaystyle\frac{\partial{u}_{z}}{\partial
h_{p}}=\sum_{q=-\widehat{N}}^{\widehat{N}}$ $\displaystyle
f(q,\kappa_{p})e^{-j2\pi\frac{l_{p}(k_{p}+\kappa_{p})}{MN}}x[k^{\prime},l^{\prime}],$
(99)

while for the rectangular waveform we have

$\displaystyle\frac{\partial{u}_{z}}{\partial
h_{p}}=\sum_{q=-\widehat{N}}^{\widehat{N}}$ $\displaystyle
f(q,\kappa_{p})e^{j2\pi\frac{(l-l_{p})(k_{p}+\kappa_{p})}{MN}}x[k^{\prime},l^{\prime}].$
(100)

When $P<p^{\prime}\leq 2P$ ($p^{\prime}$ represents $i$ or $j$), the
derivative $\frac{\partial{u}_{z}}{\partial{\theta}_{p^{\prime}}}$ is about
$\kappa_{p}$, where $p=p^{\prime}-P$. For the bi-orthogonal waveform, we have

$\displaystyle\frac{\partial{u}_{z}}{\partial\kappa_{p}}=$
$\displaystyle\sum_{q=-\widehat{N}}^{\widehat{N}}h_{p}\biggl{[}\left(\frac{1}{N}\sum_{n=1}^{N-1}j\frac{2n\pi}{N}e^{j\frac{2n\pi}{N}(q+\kappa_{p})}\right)e^{-j2\pi\frac{l_{p}(k_{p}+\kappa_{p})}{MN}}$
$\displaystyle+$ $\displaystyle
f(q,\kappa_{p})\left(e^{-j2\pi\frac{l_{p}(k_{i}+\kappa_{p})}{MN}}\frac{-j2\pi
l_{p}}{MN}\right)\biggl{]}x[k^{\prime},l^{\prime}],$ (101)

while for the rectangular waveform

$\displaystyle\frac{\partial{u}_{z}}{\partial\kappa_{p}}=$
$\displaystyle\sum_{q=-\widehat{N}}^{\widehat{N}}h_{p}\biggl{[}\left(\frac{1}{N}\sum_{n=1}^{N-1}j\frac{2n\pi}{N}e^{j\frac{2n\pi}{N}(q+\kappa_{p})}\right)e^{j2\pi\frac{(l-l_{p})(k_{p}+\kappa_{p})}{MN}}$
$\displaystyle+$ $\displaystyle
f(q,\kappa_{p})\left(e^{j2\pi\frac{(l-l_{p})(k_{p}+\kappa_{p})}{MN}}\frac{j2\pi(l-l_{p})}{MN}\right)\biggl{]}x[k^{\prime},l^{\prime}].$
(102)

To evaluate the performance of the proposed algorithm in terms of the
normalized mean squared error (NMSE), we use the average normalized CRLB for
the estimation of channel gains and fractional Doppler shifts, which are
defined as

$\displaystyle\hat{\bm{h}}_{bound}=$
$\displaystyle\frac{\sum_{j=1}^{P}{\theta}_{j}^{CRLB}}{||\bm{h}_{\bm{\theta}}||_{2}^{2}},$
(103) $\displaystyle\hat{\bm{\kappa}}_{bound}=$
$\displaystyle\frac{\sum_{j=P+1}^{2P}{\theta}_{j}^{CRLB}}{||\bm{\kappa}_{\bm{\theta}}||_{2}^{2}},$
(104)

where $\bm{h}_{\bm{\theta}}=(h_{1},h_{2},\cdots,h_{P})^{T}$ and
$\bm{\kappa}_{\bm{\theta}}=(\kappa_{1},\kappa_{2},\cdots,\kappa_{P})^{T}$.

##  VI Simulation Results

In this section, we evaluate the performance of the proposed message passing
algorithm in terms of the NMSE of the estimated channel parameters
$\hat{\bm{h}}$ and $\hat{\bm{\kappa}}$, and the reconstructed channel matrices
$\widehat{\bm{H}}_{bi}$ and $\widehat{\bm{H}}_{rect}$. The NMSE is defined as

$\displaystyle
NMSE(\bm{x})=\frac{1/L\sum_{l=1}^{L}||\hat{\bm{x}}_{l}-\bm{x}||_{2}^{2}}{||\bm{x}||_{2}^{2}},$
(105)

where $\hat{\bm{x}}_{l}$ denotes the estimate of a variable $\bm{x}$ in the
$l$th trial, and $L$ is the number of trials. In addition, the bit error rate
(BER) of data detection based on the reconstructed channel matrix is also
evaluated. We set $M=128$ and $N=32$, i.e., there are $32$ time slots and
$128$ subcarriers in the TF domain. The carrier frequency is 3 GHz, the
subcarrier spacing is 2 kHz, and the quadrature phase shift keying (QPSK) is
adopted for modulation. The velocity of the mobile user is set to be
$v=120km/h$, leading to a maximum Doppler frequency shift index $k_{max}=4$.
We assume that the maximum delay index $l_{{max}}$ = 10, the Doppler index of
the $i$th path is uniformly distributed over $[-k_{max},k_{max}]$, and the
delay index is uniformly distributed over $[1,l_{max}]$, excluding the first
path ($\l_{i}=0$). As mentioned before, the fractional Doppler $\kappa_{i}$
has a uniform distribution over $[-0.5,0.5]$. The channel path gains
$\\{h_{j}\\}$ are independently drawn from the complex Gaussian distribution
$\mathcal{N}(0,1/P)$. Furthermore, we define the SNR of pilot and data as

$\displaystyle{\rm SNRp}=$ $\displaystyle 10\log
10\left(\frac{\mathcal{P}_{p}}{\gamma^{-1}}\right),$ (106) $\displaystyle{\rm
SNRd}=$ $\displaystyle 10\log
10\left(\frac{\mathbb{E}|x_{d}|^{2}}{\gamma^{-1}}\right),$ (107)

respectively, where $\mathcal{P}_{p}$ denotes the average power of pilot
symbols.

###  VI-A NMSE performance comparison

We examine the NMSE performance for the estimation of ${\bm{h}}$,
${\bm{\kappa}}$ and ${\bm{H}}$ for both the bi-orthogonal and the rectangular
waveforms. For comparison, we also include the threshold based channel
estimation method in [13] for the case of the bi-orthogonal waveform (as it is
not clear how to perform the channel estimation in the case of the rectangular
waveform with fractional Doppler shifts). The CRLB derived in the previous
section is also shown for reference.

Figure 5: NMSE of $\widehat{\bm{H}}$ of the proposed algorithm and the threshold based method in [13]. TABLE III: NMSE and PAPR comparisons: 1 pilot symbol versus 10 pilot symbols with the same power budget. SNRp | NMSE | PAPR
---|---|---
$\hat{\bm{h}}$ | $\hat{\bm{\kappa}}$ | $\widehat{\bm{H}}$
50dB |  |  |  |
1 pilot | -37.95 dB | -34.48 dB | -40.17 dB | 18.76 dB
40dB |  |  |  |
10 pilots | -36.95 dB | -32.10 dB | -39.19 dB | 10.41 dB
Figure 6: NMSE comparison for the bi-orthogonal waveform with different $P$:
(a) NMSE of $\hat{\bm{h}}$; (b) NMSE of $\hat{\bm{\kappa}}$. Figure 7: NMSE
comparison for the rectangular waveform with different $P$: (a) NMSE of
$\hat{\bm{h}}$; (b) NMSE of $\hat{\bm{\kappa}}$. Figure 8: NMSE comparison for
the bi-orthogonal waveform with different number of pilot symbols: (a) NMSE of
$\hat{\bm{h}}$; (b) NMSE of $\hat{\bm{\kappa}}$. Figure 9: NMSE comparison for
the rectangular waveform with different number of pilot symbols: (a) NMSE of
$\hat{\bm{h}}$; (b) NMSE of $\hat{\bm{\kappa}}$.

In Fig. 5, we compare the NMSE performance of the reconstructed channel matrix
$\widehat{\bm{H}}$ between our proposed algorithm and the threshold based
method in [13]. Only a single pilot symbol is used and the number of paths
$P=6$ and $10$ are considered. Moreover, with the recommendation in [13], the
threshold is set to be $3\sigma_{p}$ for the threshold based method, where
$\sigma_{p}=(10^{{\rm SNRp}/10})^{-0.5}$. Note that the NMSE of $\hat{\bm{h}}$
and $\hat{\bm{\kappa}}$ cannot be compared as the threshold based method do
not provide these estimates. It can be seen that the proposed algorithm
significantly outperforms the threshold based method in all the cases. The
reason is that the threshold based method does not estimate the channel gains
and fractional Doppler shifts, and it is equivalent to the estimate of the
vector $\bm{c}$ without considering its structure. The results indicate that
exploiting the structure of the vector $\bm{c}$ is crucial to improving the
estimation performance.

In Tab. III, we compare the PAPR and NMSE performance of $\hat{\bm{h}}$,
$\hat{\bm{\kappa}}$ and $\widehat{\bm{H}}$ by using our proposed algorithm in
two cases: a single pilot symbol and 10 pilot symbols. We assume a same power
budget, so a single pilot symbol with SNRp = 50dB corresponds to 10 pilot
symbols with SNRp = 40dB. From the table we can see that, the NMSE performance
of using 10 pilot symbols is very close to that of using a single pilot
symbol, but the PAPR with 10 pilot symbols is only 10.41 dB, in contrast to
18.76 dB for the case of using a single pilot symbol.

The NMSE performances of the algorithm for different SNRp and number of paths
$P$ are shown in Fig. 6 and Fig. 7, where the number of pilot symbols is set
to be 10. The performance of the threshold based method is not included as it
only works with a single pilot symbol. It can be seen that a smaller $P$ leads
to a better performance as the vector $\bm{c}$ is more sparse and the number
of parameters to be estimated is smaller. Comparing them to the CRLB, we can
see that the NMSE performance of the proposed algorithm is fairly good as the
CRLB is a loose lower bound, which is derived with the assumption that the
locations and the number of the paths are known. In Fig. 8 and Fig. 9, we show
the performance of the proposed method with different numbers of pilot
symbols. As expected, with more pilot symbols, the proposed algorithm delivers
better performance.

###  VI-B BER performance comparison

(a) $P=6$; SNRp = 35dB; 1 pilot symbol

(b) $P=6$; SNRp = 40dB; 1 pilot symbol

(c) $P=10$; SNRp = 40dB; 1 pilot symbol

(d) $P=10$; SNRp = 40dB; 10 pilot symbols

Figure 10: BER performance comparison versus SNRp.

We compare the BER performance of the OTFS system with the UTAMP based
detector proposed in [4]. Assume that the bi-orthgonal waveform is used so
that we can compare the performance of the system with the threshold based
channel estimation method [13]. As a benchmark, we also show the BER
performance of the OTFS system with perfect channel matrix. In addition, in
order to demonstrate the significance of considering fractional Doppler
shifts, we include the performance of the OTFS system with the assumption of
integer Doppler shifts (although fractional Doppler shifts exist in the
generation of the OTFS channels). The results are shown in Fig. 10. It can be
seen that, when the fractional Doppler shifts are ignored, the system
performance is degraded severely due to the modelling errors. In addition, in
the case of single pilot symbol, the system with the proposed channel
estimation algorithm delivers much better BER performance than that with the
threshold based channel estimation method. Meanwhile, the performance of the
system with our proposed channel estimation algorithm can be very close to
that of the system with perfect channel matrix. By comparing Fig. 10(b) with
Fig. 10(c), we can see that a larger $P$ leads to a better BER performance
because more diversity gain can be achieved.

##  VII Conclusions

In this paper, we have addressed the issue of the estimation of OTFS channels
with fractional Doppler shifts, where both the bi-orthogonal waveform and the
rectangular waveform are considered. The estimation is formulated as a sparse
structured signal recovery problem and a Bayesian treatment is investigated.
With a factor graph representation of the problem, we have derived an message
passing based algorithm to estimate the channel gains and fractional Doppler
shifts. The CRLB has also been derived to evaluate the estimation performance
of the proposed algorithm. It has been shown that the proposed algorithm
significantly outperforms the existing algorithm, and it is able to work with
multiple pilot symbols to achieve significant PAPR reduction.

## Acknowledgment

The authors would like to thank Prof. Jinhong Yuan at the University of New
South Wales for his valuable comments and suggestions.

## References

  * [1] R. Hadani, S. Rakib, M. Tsatsanis, A. Monk, and R. Calderbank, “Orthogonal Time Frequency Space Modulation,” in _2017 IEEE Wireless Communications and Networking Conference (WCNC)_ , 2017.
  * [2] P. Raviteja, P. K. T., H. Yi, and V. Emanuele, “Interference Cancellation and Iterative Detection for Orthogonal Time Frequency Space Modulation,” _IEEE Transactions on Wireless Communications_ , vol. 17, no. 10, pp. 6501–6515, 2018.
  * [3] R. Hadani and A. Monk, “OTFS: A New Generation of Modulation Addressing the Challenges of 5G,” _ArXiv_ , vol. abs/1802.02623, 2018.
  * [4] Z. Yuan, F. Liu, W. Yuan, Q. Guo, Z. Wang, and J. Yuan, “Iterative Detection for Orthogonal Time Frequency Space Modulation with Unitary pproximate Message Passing ,” _ArXiv_ , vol. abs/2008.06688, 2020.
  * [5] P. Raviteja, E. Viterbo, and Y. Hong, “OTFS Performance on Static Multipath Channels,” _IEEE Wireless Communications Letters_ , vol. 8, no. 3, pp. 745–748, June 2019.
  * [6] K. R. Murali and A. Chockalingam, “On OTFS Modulation for High-Doppler Fading Channels,” in _2018 Information Theory and Applications Workshop (ITA)_ , 2018, pp. 1–10.
  * [7] S. Li, J. Yuan, W. Yuan, Z. Wei, B. Bai, and D. W. K. Ng, “Performance Analysis of Coded OTFS Systems over High-Mobility Channels,” _arXiv e-prints_ , p. arXiv:2010.13008, Oct. 2020.
  * [8] L. Li, H. Wei, Y. Huang, Y. Yao, W. Ling, G. Chen, P. Li, and Y. Cai, “A Simple Two-stage Equalizer With Simplified Orthogonal Time Frequency Space Modulation Over Rapidly Time-varying Channels,” _CoRR_ , vol. abs/1709.02505, 2017. [Online]. Available: http://arxiv.org/abs/1709.02505
  * [9] P. Raviteja, Y. Hong, E. Viterbo, and E. Biglieri, “Practical pulse-shaping waveforms for reduced-cyclic-prefix OTFS,” _IEEE Transactions on Vehicular Technology_ , vol. 68, no. 1, pp. 957–961, 2019.
  * [10] T. Zemen, M. Hofer, and D. Loeschenbrand, “Low-Complexity Equalization for Orthogonal Time and Frequency Signaling (OTFS),” 2017. [Online]. Available: http://arxiv.org/abs/1710.09916
  * [11] S. Tiwari, S. S. Das, and V. Rangamgari, “Low complexity LMMSE Receiver for OTFS,” _IEEE Communications Letters_ , vol. 23, no. 12, pp. 2205–2209, 2019.
  * [12] D. Tse and P. Viswanath, _Fundamentals of Wireless Communication_. USA: Cambridge University Press, 2005.
  * [13] P. Raviteja, K. T. Phan, and Y. Hong, “Embedded Pilot-Aided Channel Estimation for OTFS in Delay–Doppler Channels,” _IEEE Transactions on Vehicular Technology_ , vol. 68, no. 5, pp. 4906–4917, 2019.
  * [14] M. Kollengode Ramachandran and A. Chockalingam, “MIMO-OTFS in High-Doppler Fading Channels: Signal Detection and Channel Estimation,” in _2018 IEEE Global Communications Conference (GLOBECOM)_ , 2018, pp. 206–212.
  * [15] W. Shen, L. Dai, J. An, P. Fan, and R. W. Heath, “Channel Estimation for Orthogonal Time Frequency Space (OTFS) Massive MIMO,” _IEEE Transactions on Signal Processing_ , vol. 67, no. 16, pp. 4204–4217, 2019.
  * [16] Y. Liu, S. Zhang, F. Gao, J. Ma, and X. Wang, “Uplink-Aided High Mobility Downlink Channel Estimation Over Massive MIMO-OTFS System,” _IEEE Journal on Selected Areas in Communications_ , vol. 38, no. 9, pp. 1994–2009, 2020.
  * [17] A. Monk, R. Hadani, M. Tsatsanis, and S. Rakib, “OTFS - Orthogonal Time Frequency Space,” _CoRR_ , vol. abs/1608.02993, 2016. [Online]. Available: http://arxiv.org/abs/1608.02993
  * [18] F.R. Kschischang and B.J. Frey and H.-A. Loeliger, “Factor graphs and the sum-product algorithm,” _IEEE Trans. Inform. Theory_ , vol. 47, no. 2, pp. 498–519, Feb. 2001.
  * [19] M. E. Tipping, “Sparse Bayesian Learning and the Relevance Vector Machine,” _Journal of Machine Learning Research_ , vol. 1, no. 3, pp. 211–244, 2001\.
  * [20] S. D. Babacan, R. Molina, and A. K. Katsaggelos, “Bayesian Compressive Sensing Using Laplace Priors,” _IEEE Transactions on Image Processing_ , vol. 19, no. 1, pp. 53–63, 2010.
  * [21] Winn, John, Bishop, M. Christopher, Jaakkola, and Tommi, “Variational Message Passing,” _Journal of Machine Learning Research_ , vol. 6, pp. 661–694, 2005.
  * [22] E. P. Xing, M. I. Jordan, and S. Russell, “A Generalized Mean Field Algorithm for Variational Inference in Exponential Families,” in _Nineteenth Conference on Uncertainty in Artificial Intelligence_ , 2002.
  * [23] L. He, Y. Wu, S. Ma, T. Ng, and H. V. Poor, “Superimposed Training-Based Channel Estimation and Data Detection for OFDM Amplify-and-Forward Cooperative Systems Under High Mobility,” _IEEE Transactions on Signal Processing_ , vol. 60, no. 1, pp. 274–284, 2012.
  * [24] K. Zhong, Y. Wu, and S. Li, “Signal Detection for OFDM-Based Virtual MIMO Systems under Unknown Doubly Selective Channels, Multiple Interferences and Phase Noises,” _IEEE Transactions on Wireless Communications_ , vol. 12, no. 10, pp. 5309–5321, 2013.
  * [25] Q. Guo and D. D. Huang, “A Concise Representation for the Soft-in Soft-out LMMSE Detector,” _IEEE Communications Letters_ , vol. 15, no. 5, pp. 566–568, 2011.
  * [26] S. J. Julier and J. K. Uhlmann, “Unscented filtering and nonlinear estimation,” _Proceedings of the IEEE_ , vol. 92, no. 3, pp. 401–422, 2004\.
  * [27] Z. Yuan, Q. Guo, and M. Luo, “Approximate message passing with unitary transformation for robust bilinear recovery,” _IEEE Transactions on Signal Processing_ , vol. 69, pp. 617–630, 2021.
  * [28] S. K. Sengijpta, “Fundamentals of Statistical Signal Processing: Estimation Theory,” _Technometrics_ , vol. 37, no. 4, pp. 465–466, 1993.