shadows MyMdframed[1][]

# Off-grid Channel Estimation with Sparse Bayesian Learning for OTFS Systems

Zhiqiang Wei, Weijie Yuan, Shuangyang Li, Jinhong Yuan, and Derrick Wing Kwan
Ng Zhiqiang Wei, Weijie Yuan, Shuangyang Li, Jinhong Yuan, and Derrick Wing
Kwan Ng are with the School of Electrical Engineering and Telecommunications,
the University of New South Wales, Australia (email: zhiqiang.wei;
weijie.yuan; shuangyang.li; j.yuan; w.k.ng@unsw.edu.au).

###### Abstract

This paper proposes an off-grid channel estimation scheme for orthogonal time-
frequency space (OTFS) systems adopting the sparse Bayesian learning (SBL)
framework. To avoid channel spreading caused by the fractional delay and
Doppler shifts and to fully exploit the channel sparsity in the delay-Doppler
(DD) domain, we estimate the original DD domain channel response rather than
the effective DD domain channel response as commonly adopted in the
literature. OTFS channel estimation is firstly formulated as a one-dimensional
(1D) off-grid sparse signal recovery (SSR) problem based on a virtual sampling
grid defined in the DD space, where the on-grid and off-grid components of the
delay and Doppler shifts are separated for estimation. In particular, the on-
grid components of the delay and Doppler shifts are jointly determined by the
entry indices with significant values in the recovered sparse vector. Then,
the corresponding off-grid components are modeled as hyper-parameters in the
proposed SBL framework, which can be estimated via the expectation-
maximization method. To strike a balance between channel estimation
performance and computational complexity, we further propose a two-dimensional
(2D) off-grid SSR problem via decoupling the delay and Doppler shift
estimations. In our developed 1D and 2D off-grid SBL-based channel estimation
algorithms, the hyper-parameters are updated alternatively for computing the
conditional posterior distribution of channels, which can be exploited to
reconstruct the effective DD domain channel. Compared with the 1D method, the
proposed 2D method enjoys a much lower computational complexity while only
suffers a slight performance degradation. Simulation results verify the
superior performance of the proposed channel estimation schemes over state-of-
the-art schemes.

##  I Introduction

Future wireless systems are envisioned to enable a wide range of emerging
mobile applications, such as low-earth-orbit satellites (LEOS), autonomous
cars, in-vehicle infotainment, and unmanned aerial vehicles (UAV) [1]. As the
airborne, spaceborne, and vehicles are usually fast moving, the associated
severe Doppler spread experienced by channels in high mobility propagation
environments imposes new challenges on its air interface design. Recently, a
new modulation waveform, namely orthogonal time-frequency space (OTFS), has
been proposed as a promising candidate for realizing high-mobility
communications owing to its intrinsic advantages of handling the problems
caused by Doppler effect[2, 3].

In contrast to the existing one-dimensional (1D) time-domain or frequency-
domain modulation techniques, such as orthogonal frequency-division
multiplexing (OFDM), OTFS is a two-dimensional (2D) modulation scheme, where
the information symbols are carried over 2D localized pulses defined in the
delay-Doppler (DD) domain. In fact, each information symbol multiplexed in the
DD domain is spread across the whole time-frequency (TF) domain, allowing the
potentials of exploiting the full diversity to achieve reliable communications
against various channel impairments. In addition, the channel in the DD domain
exhibits various beneficial properties, such as separability, stability, and
possibly sparsity [1], which can be exploited for efficient channel estimation
and data detection. In particular, processing the received signals in the
Doppler domain of wireless channels allows us to separate the propagation
paths experiencing an identical delay. Besides, OTFS modulation can transform
a time-variant channel in the TF domain into a 2D quasi-time-invariant channel
in the DD domain and thus relieves the impacts caused by channel aging [2, 3].
Additionally, the DD domain effective channel is generally sparse in an open-
space rural propagation environment with a limited number of moving scatters.
As a result, the coupling relationship between data symbols and the channel in
the DD domain is much simpler than the counterpart in the TF domain, which is
crucial for accurate channel acquisition and efficient detection algorithm
design.

Since OTFS multiplexes information symbols in the DD domain, it enables the
possibility of channel probing via designing a new training signaling in the
DD domain[4, 5], which is a challenging but vital requirement for reliable
detection. In particular, channel spreading [6] caused by fractional Doppler
and delay sacrifices the channel sparsity which degrades the channel
estimation performance. Usually, a guard space is inserted between the pilot
and data symbols to mitigate the interference of unknown data symbols for
channel estimation, which incurs a significant signaling overhead[4]. In fact,
the channel estimation problem of OTFS was firstly studied in [4], where a
single impulse was embedded in the DD domain and a threshold-based estimation
method was proposed. Then, its extension to multiple-input multiple-output
OTFS (MIMO-OTFS) systems was investigated in [5] via positioning the pilot
impulses for different transmit antennas with sufficient guard spacing in the
DD domain. However, their channel estimation performances are sensitive to the
availability of guard space [4, 5]. Indeed, to exploit the channel sparsity in
the DD domain, more advanced channel estimation approaches based on compressed
sensing techniques have been proposed in the literature. For instance, the
authors in [7] proposed a three-dimensional (3D) structured orthogonal
matching pursuit (SOMP) algorithm to estimate the delay-Doppler-angle domain
channel via exploiting the 3D structured sparsity. A further extension in [8]
proposed a 3D Newtonized OMP (NOMP) algorithm, which can extract the
fractional component in the Doppler and angle domains via Newton’s method.
Yet, its channel estimation performance is mainly limited by the greedy OMP
approach, which is generally far away from that achieved by the optimal
estimation. Besides, a most recent work in [9] applied the sparse Bayesian
learning approach for the channel estimation in OTFS systems, which
outperforms the OMP methods [7, 8]. Nevertheless, the proposed scheme in [9]
only applies to the integer Doppler case and is essentially an on-grid method,
where all the estimated Doppler shifts are constrained on the coarse DD domain
grid, which is only valid in limited practical scenarios.

In practice, the DD domain channel sparsity may not always hold [6], where
conventional compressed sensing-based channel estimation cannot achieve a
satisfactory performance. In particular, due to the limited bandwidth and the
limited frame duration, the exact Doppler frequency shift or delay shift may
straddle a pair of finite-resolution bins, rather than falling exactly into a
single bin, which is known as off-grid. In this case, the effective DD domain
channel is spread across all the Doppler and delay indices, when off-grid
happens in the Doppler and delay domains, respectively. As a result, although
the original DD domain channel response is sparse, the effective DD domain
channel may not be sparse due to channel spreading. On the other hand,
although increasing the sampling resolution in the DD domain is a viable
approach to suppress the impacts caused by channel spreading, it consumes more
bandwidth and introduces longer latency. As a result, various approaches have
been proposed as a remedy. For instance, our previous work [6] proposed to
apply a Dolph-Chebyshev (DC) window to suppress channel spreading and thus to
enhance the effective DD domain channel sparsity, which substantially improves
the channel estimation performance compared with [4]. Yet, [6] still estimates
the effective DD domain channel, rather than the original DD domain channel
response, and thus remains suffering from the inevitable channel spreading.

In this paper, we formulate an off-grid sparse signal recovery (SSR) problem
and propose to employ a sparse Bayesian learning (SBL) framework to estimate
the original DD domain channel response for OTFS systems. In particular, to
avoid channel spreading and to fully reap the potential sparsity, we aim to
estimate the original DD domain channel response directly rather than the
effective DD domain channel, which is fundamentally different from [7, 9].
Besides, the proposed off-grid channel estimation scheme is essentially a
semi-parametric estimation approach. In particular, the exact Doppler and
delay shifts can be estimated via their nearest grid indices in the recovered
sparse vector while the corresponding off-grid components can be estimated via
some parametric approach, such as the expectation-maximization (EM) algorithm.
There are three advantages of the proposed SBL-based channel estimation
approaches. Firstly, the hyperprior structure of SBL [10] is flexible such
that the channel compactness and sparsity in the DD domain can be easily
exploited. Secondly, SBL can output the posterior distribution of the original
DD domain channel response, where the mean denotes the channel estimate and
the variance determines the channel estimation accuracy. Thirdly, SBL avoids
the use of sparse regularization parameters [11] and it is not sensitive to
the other initialization parameters, such as the number of paths and noise
power. As a result, the proposed channel estimation approach significantly
outperforms the algorithm in [8], which is the only work of compresses
sensing-based channel estimation so far that can handle the fractional Doppler
case. The main contributions of this work are given as follows:

  * •

By characterizing the DD domain input-output relationship, we first
reformulate the channel estimation of OTFS systems as a one-dimensional (1D)
SSR problem, where we directly estimate the original DD domain response rather
than the effective DD domain channel, e.g. [4, 5, 6]. Moreover, we establish a
1D off-grid compressed channel estimation model based on the first-order
linear approximation, which enables us to improve the delay and Doppler shift
estimation accuracy while maintaining a low computational complexity. The
developed model serves as a generalized platform which is applicable to
different pilot and data arrangements, arbitrary windows, as well as pulse
shaping filters adopted at the OTFS transceivers.

  * •

We first propose a 1D off-grid SBL-based channel estimation algorithm, which
can estimate both the on-grid and off-grid components of the delay and Doppler
shifts. In particular, the on-grid component is obtained by the entry indices
with significant values in the recovered sparse vector. The corresponding off-
grid components are treated as hyper-parameters in the SBL framework which can
be estimated via an iterative algorithm based on the expectation-maximization
approach.

  * •

To further reduce the computational complexity, we decouple the delay and
Doppler shift estimation and develop a 2D off-grid compressed channel
estimation model. A two-step algorithm is proposed to estimate the Doppler and
delay shifts, respectively, where the first-step is realized via applying SBL
to multiple snapshot cases and the second step is achieved by solving multiple
1D SSR problems in parallel. The computational complexity order of the
proposed 2D off-grid method is proportional to the summation of the sizes of
delay and Doppler grids, while the counterpart of the proposed 1D off-grid
method is proportional to their product.

  * •

Extensive simulations are conducted to evaluate and compare the performance of
our proposed schemes. We demonstrate that our proposed 1D off-grid SBL-based
scheme can achieve the most accurate channel estimates compared with the
state-of-the-art schemes in the literature [4, 6, 7, 8, 9]. Also, the
performance of our proposed 2D off-grid SBL-based scheme approaches closely to
that of the 1D method, with a significantly reduced computational complexity.

Notations: Boldface capital and lower case letters are reserved for matrices
and vectors, respectively; ${\left(\cdot\right)^{\mathrm{H}}}$ denotes the
Hermitian transpose of a vector or a matrix;
${\left(\cdot\right)^{\mathrm{T}}}$ denotes the transpose of a vector or a
matrix; ${\left(\cdot\right)^{-1}}$ denotes the inverse of a square matrix;
$\mathbb{A}$ denotes the constellation set; $\mathbb{Z}^{+}$ denotes the set
of all non-negative integers; $\mathbb{C}^{M\times N}$ denotes the set of all
$M\times N$ matrices with complex entries; $\lvert\cdot\rvert$ denotes the
absolute value of a complex scalar or the cardinality of a set;
$E_{X}\\{\cdot\\}$ denotes the expectation of the input with respect to
(w.r.t.) the distribution of $X$; $\left(\cdot\right)^{*}$ denotes the
conjugate operation; $\otimes$ is the Kronecker product operator; $\odot$ is
the point-wise product operator; $\left(\cdot\right)_{N}$ denotes the modulus
operation w.r.t. $N$; $\operatorname{\mathrm{diag}}\\{\cdot\\}$ returns a
square diagonal matrix with the elements of input vector on the main diagonal;
$\Re\\{\cdot\\}$ returns the real part of the input complex number;
$\lfloor\cdot\rfloor$ is the floor function which returns the largest integer
smaller than the input value; $\lceil\cdot\rceil$ is the ceiling function
which returns the smallest integer larger than the input value;
$\mathbf{I}_{M}$ denote the identity matrix of size $M\times M$;
$\\{{\mathbf{X}}\\}_{n}$ denotes the $n$-th column of the matrix $\mathbf{X}$
and $\\{{\mathbf{X}}\\}_{m,n}$ denotes the $\left[m,n\right]$-th entry of the
matrix $\mathbf{X}$; $\\{{\mathbf{x}}\\}_{n}$ denotes the $n$-th entry of the
vector $\mathbf{x}$ and $\\{{\mathbf{x}}\\}_{-n}$ denotes the vector without
the $n$-th entry of $\mathbf{x}$. The big-O notation
$\mathcal{O}\left(\cdot\right)$ asymptotically describes the order of
computational complexity. The circularly symmetric complex Gaussian
distribution with mean $\bm{\mu}$ and covariance matrix $\bm{\Sigma}$ is
denoted by ${\cal CN}(\bm{\mu},\bm{\Sigma})$; the uniform distribution between
$a$ and $b$ is denoted by $\mathcal{U}\left[a,b\right]$; $\sim$ stands for
“distributed as”; $\propto$ denotes “proportional to”. The probability density
function (PDF) for a random complex-valued vector ${\bf{x}}\sim{\cal
CN}(\bm{\mu},\bm{\Sigma})$ is defined as
$\mathcal{CN}\left({{\bf{x}}\left|{{\bm{\mu}},{\bf{\Sigma}}}\right.}\right)=\frac{1}{{{\pi^{N}}\left|{\bf{\Sigma}}\right|}}\exp\left\\{{-{{\left({{\bf{x}}-{\bm{\mu}}}\right)}^{\rm{H}}}{{\bf{\Sigma}}^{-1}}\left({{\bf{x}}-{\bm{\mu}}}\right)}\right\\}$.
The PDF for a random variable $x$ follows a Gamma distribution is defined as
$\Gamma\left({x\left|{a,b}\right.}\right)=\frac{{b^{a}}{x^{a-1}}\exp\left\\{{-bx}\right\\}}{{{\Gamma\left(a\right)}}}$,
where constants $a>0$ and $b>0$, and ${\Gamma\left(a\right)}$ is a Gamma
function.

Figure 1: OTFS system model in the DD domain.

##  II System Model

For ease of presentation, as shown in Fig. 1, we directly present the DD
domain input-output relationship of OTFS following the well-known derivations
in the literature[12, 1, 6, 8]:

$y\left[{k,l}\right]=\sum\limits_{k^{\prime}=0}^{N-1}{\sum\limits_{l^{\prime}=0}^{M-1}{x\left[{k^{\prime},l^{\prime}}\right]}}{h_{w}}\left[{k-k^{\prime},l-l^{\prime}}\right]+z\left[k,l\right],\vspace{-2mm}$
(1)

where $x\left[{k^{\prime},l^{\prime}}\right]\in\mathbb{A}$ is the modulated
symbol arranged on the $\left[k^{\prime},l^{\prime}\right]$-th grid point in
the DD domain and $y\left[{k,l}\right]\in\mathbb{C}$ denotes the
$\left[k,l\right]$-th sampled receiving signal in the DD domain. Variable
$z\left[k,l\right]$ denotes the additive noise whose distribution depends on
the imposed receiver window [6]. For example, adopting a rectangular window at
the receiver, we have
$z\left[k,l\right]\sim\mathcal{CN}\left(0,{\sigma}^{2}\right)$, where
${\sigma}^{2}$ denotes the unknown noise power[6]. Integer $N$ denotes the
number of time slots with a slot duration of $T$ for the OTFS system,
resulting in one OTFS frame duration of $NT$. Besides, integer $M$ represents
the number of subcarriers with a subcarrier spacing of $\Delta_{f}$ for the
OTFS system and the total system bandwidth is $M\Delta_{f}$. Interested
readers may refer to our previous work [6] for more derivation details of
(1)111Note that the input-output relationship in (1) was derived for a single-
antenna OTFS transceiver. As we mainly focus on handling the fractional
Doppler and delay estimation issue of OTFS, we would like to adopt a single-
input single-output system in this paper for obtaining important system design
insights. The extension to multiple-input multiple-output systems [13] will be
considered in our future work..

In (1), ${h_{w}}\left[{k-k^{\prime},l-l^{\prime}}\right]$ denotes the
effective DD domain channel, which is obtained by sampling the original DD
domain channel response. In particular, the original DD domain channel
response is given by

${h}\left(\nu,\tau\right)=\sum_{i=1}^{P}{h_{i}}\delta\left(\nu-\nu_{i}\right)\delta\left(\tau-\tau_{i}\right),\vspace{-2mm}$
(2)

where $P$ represents the number of paths in the channel, ${h}_{i}$ is the
channel coefficient of the $i$-th path, and $\delta\left(\cdot\right)$ is the
Dirac delta function. Variables
${\tau_{i}}\in\left(0,\tau_{\mathrm{max}}\right)$ and
${\nu_{i}}\in\left(-\nu_{\mathrm{max}},\nu_{\mathrm{max}}\right)$ denote the
delay and Doppler shifts of the $i$-th path, respectively. In particular,
constants $\nu_{\mathrm{max}}$ and $\tau_{\mathrm{max}}$ denote the maximum
Doppler shift and the maximum delay shift in the propagation environments,
which are related to the maximum relative velocity among the transceivers or
scatters and the maximum propagation distance, respectively. In this paper, we
assume that the number of path $P$ is unknown while $\nu_{\mathrm{max}}$ and
$\tau_{\mathrm{max}}$ are known in advance by our proposed channel estimation
approaches. In fact, the Doppler and delay shift boundary can be obtained
empirically via long-term measurements, i.e., i.e.,
$\nu_{\mathrm{max}}=\max_{i}\\{\left|\nu_{i}\right|\\}$ and
$\tau_{\mathrm{max}}=\max_{i}\\{\tau_{i}\\}$. However, the number of paths for
each channel realization mainly depends on the presenting scatters in the
propagation environment, which is usually unpredictable. As a result, the
effective DD domain channel can be obtained by[12, 1, 6, 8]

${h_{w}}\left[{k-k^{\prime},l-l^{\prime}}\right]=\sum_{i=1}^{P}{h_{i}}w(k-k^{\prime}-k_{\nu_{i}},l-l^{\prime}-l_{\tau_{i}}){e^{-j2\pi{\nu_{i}}l_{\tau_{i}}}}=\sum_{i=1}^{P}{\tilde{h}_{i}}w(k-k^{\prime}-k_{\nu_{i}},l-l^{\prime}-l_{\tau_{i}}),\vspace{-2mm}$
(3)

where $k,k^{\prime}\in\left\\{0,\ldots,N-1\right\\}$,
$l,l^{\prime}\in\left\\{0,\ldots,M-1\right\\}$, and
$k_{\nu_{i}}={\nu_{i}}{NT}\in\mathbb{R}$ and
$l_{\tau_{i}}={\tau_{i}}{M\Delta_{f}}\in\mathbb{R}$ represent the
corresponding normalized version of $\nu_{i}$ and ${\tau_{i}}$, respectively.
To facilitate the channel estimation, the phase
${e^{-j2\pi{\nu_{i}}l_{\tau_{i}}}}$ is absorbed into the channel coefficient
${\tilde{h}_{i}}={{h}_{i}}{e^{-j2\pi{\nu_{i}}l_{\tau_{i}}}}$ since the delay
and Doppler variables are not separable in the exponential function.

In (3), $w\left(k-k^{\prime}-k_{\nu},l-l^{\prime}-l_{\tau}\right)$ is the
sampling function, which captures the joint effects of the Doppler and delay
shift of the original DD domain channel response, the adopted pulse shaping
filters and the TF domain window at the OTFS modulator and demodulator, as
shown in Fig. 1, which can be obtained by[12, 6]

$\displaystyle w\left(k-k^{\prime}-k_{\nu},l-l^{\prime}-l_{\tau}\right)$
$\displaystyle=$
$\displaystyle\frac{1}{{NM}}\sum\limits_{n=0}^{N-1}\sum\limits_{m=0}^{M-1}\sum\limits_{m^{\prime}\neq
m}^{M-1}{A_{{g_{rx}}{g_{tx}}}}\left({-{l_{\tau}}\frac{1}{{M\Delta
f}},\left({m-m^{\prime}}\right)\Delta f-{k_{\nu}}\frac{1}{{NT}}}\right)$
$\displaystyle\times$ $\displaystyle
U\left[{n,m^{\prime}}\right]V\left[{n,m^{\prime}}\right]{e^{-j2\pi\frac{{n\left({k-k^{\prime}-{k_{\nu}}}\right)}}{N}}}{e^{-j2\pi\frac{{\left({ml-m^{\prime}l^{\prime}-m^{\prime}{l_{\tau}}}\right)}}{M}}}$
$\displaystyle+$
$\displaystyle{e^{-j2\pi\frac{{k^{\prime}}}{N}}}\frac{1}{{NM}}\sum\limits_{n=1}^{N-1}{\sum\limits_{m=0}^{M-1}{\sum\limits_{m^{\prime}=0}^{M-1}{{A_{{g_{rx}}{g_{tx}}}}\left({T-{l_{\tau}}\frac{1}{{M\Delta
f}},\left({m-m^{\prime}}\right)\Delta f-{k_{\nu}}\frac{1}{{NT}}}\right)}}}$
$\displaystyle\times$ $\displaystyle
U\left[{n-1,m^{\prime}}\right]V\left[{n-1,m^{\prime}}\right]{e^{-j2\pi\frac{{n\left({k-k^{\prime}-{k_{\nu}}}\right)}}{N}}}{e^{-j2\pi\frac{{\left({ml-m^{\prime}l^{\prime}-m^{\prime}{l_{\tau}}}\right)}}{M}}},$
(4)

where $U\left[{n,m}\right]$ and $V\left[{n,m}\right]$ denote the adopted
transmit window and receiver window in the $\left[{n,m}\right]$-th TF domain
point, respectively. The ambiguity function
${{A_{{g_{rx}}{g_{tx}}}}\left(\tau,\nu\right)}$ for the adopted transmit pulse
shaping filter ${g_{{\rm{tx}}}}\left({t}\right)$ and receiver pulse shaping
filter ${g_{{\rm{rx}}}}\left({t}\right)$ can be defined by

${{A_{{g_{rx}}{g_{tx}}}}\left(\tau,\nu\right)}=\int_{t}{g_{{\rm{tx}}}}\left({t}\right)g_{{\rm{rx}}}^{*}\left({t-\tau}\right){e^{-j2\pi\nu\left({t-\tau}\right)}}dt.\vspace{-2mm}$
(5)

Note that the considered system model in (1)-(II) is a generalized framework,
which are applicable to an arbitrary pulse shaping filter and an arbitrary TF
domain window. In (II), we can observe that the sampling function has a
circular structure with a period of $N$ and $M$ in the Doppler and delay
domains, respectively. As a result, the effective DD domain channel
${h_{w}}\left[{k-k^{\prime},l-l^{\prime}}\right]$ also exhibits a circular
structure with a period of $N$ and $M$ in the Doppler and delay domains,
respectively, i.e.,
${h_{w}}\left[{k-k^{\prime},l-l^{\prime}}\right]={h_{w}}\left[{(k-k^{\prime})_{N},(l-l^{\prime})_{M}}\right]$.

A commonly studied case in the OTFS literature is to apply an ideal pulse
shaping filter[12] to the transceiver, which satisfies the bi-orthogonal
condition[6] as follows

${{A_{{g_{rx}}{g_{tx}}}}\left({\left({n-n^{\prime}}\right)T-{l_{\tau}}\frac{1}{{M\Delta
f}},\left({m-m^{\prime}}\right)\Delta
f-{k_{\nu}}\frac{1}{{NT}}}\right)}=q_{\tau_{\mathrm{max}}}\left(t-nT\right)q_{\nu_{\mathrm{max}}}\left(f-m\Delta
f\right),\vspace{-2mm}$ (6)

with

$q_{\tau_{\mathrm{max}}}\left(t-nT\right)=\left\\{{\begin{array}[]{*{20}{c}}{\delta\left[n\right]}&{\left|{t-nT}\right|\leq\tau_{\mathrm{max}}},\\\\[-2.84526pt]
{q\left(t\right)}&{\mathrm{otherwise}},\end{array}}\right.\;\text{and}\;q_{\nu_{\mathrm{max}}}\left(f-m\Delta
f\right)=\left\\{{\begin{array}[]{*{20}{c}}{\delta\left[m\right]}&{\left|{f-m\Delta
f}\right|\leq\nu_{\mathrm{max}}},\\\\[-2.84526pt]
{q\left(t\right)}&{\mathrm{otherwise}},\end{array}}\right.\vspace{-2mm}$ (7)

where $q\left(x\right)$ is an arbitrary function. Furthermore, with adopting
the rectangular TF domain window[6] at the transceiver, i.e.,
$U\left[{n,m}\right]=V\left[{n,m}\right]=1$, $\forall m,n$, the sampling
function in (II) can be simplified as

$w\left(k-k^{\prime}-k_{\nu},l-l^{\prime}-l_{\tau}\right)=w_{\nu}\left(k-k^{\prime}-k_{\nu}\right)w_{\tau}\left(l-l^{\prime}-l_{\tau}\right),\vspace{-2mm}$
(8)

where the sampling function can be decomposed as the sampling functions in the
Doppler and delay domains, respectively, which are given by

$\displaystyle w_{\nu}\left(k-k^{\prime}-k_{\nu}\right)$
$\displaystyle=\frac{1}{N}\left[{{e^{-j\left({N-1}\right)\pi\frac{{k-k^{\prime}-k_{\nu}}}{N}}}\frac{{\sin\left({\pi\left(k-k^{\prime}-k_{\nu}\right)}\right)}}{{\sin\left({\frac{{\pi\left(k-k^{\prime}-k_{\nu}\right)}}{N}}\right)}}}\right]\;\text{and}$
$\displaystyle w_{\tau}\left(l-l^{\prime}-l_{\tau}\right)$
$\displaystyle=\frac{1}{M}\left[{{e^{-j\left({M-1}\right)\pi\frac{l-l^{\prime}-l_{\tau}}{M}}}\frac{{\sin\left({\pi\left(l-l^{\prime}-l_{\tau}\right)}\right)}}{{\sin\left({\frac{{\pi\left(l-l^{\prime}-l_{\tau}\right)}}{M}}\right)}}}\right],$
(9)

respectively.

For notational simplicity, we focus on the case of ideal pulse shaping and
rectangular window case in this paper. The designed channel estimation
algorithm can be easily extended to the case of arbitrary pulse shaping filter
and arbitrary window with straightforward modifications of the sampling
functions following (II).

To visualize the difference between the original DD domain channel response
and the effective DD domain channel, we plot a snapshot of the corresponding
normalized channel gains in Fig. 3 and Fig. 3, respectively, where we assume
$M=N=16$ and $P=5$. We can observe that the original DD domain channel
response is indeed sparse. However, due to the existence of fractional Doppler
and delay, i.e., non-integer $k_{\nu}$ and $l_{\tau}$, even with ideal pulse
shaping filters, the sampling functions in (II) are non-zero for
$k-k^{\prime}\neq k_{\nu}$ and $l-l^{\prime}\neq l_{\tau}$. It causes channel
spreading in the effective DD domain channel in (3), which is far from sparse.
As a result, estimating the effective DD domain channel using compressed
sensing techniques [7, 9] may not lead to a satisfactory performance.
Therefore, in this paper, we propose to directly estimate the original DD
domain channel response ${h}\left(\nu,\tau\right)$ and reconstruct the
effective DD domain channel according to (3), which is a key for realizing
effective data detection.

Figure 2: Normalized original DD domain channel response.

Figure 3: Normalized effective DD domain channel gain.

##  III 1D Off-grid Compressed Channel Estimation

In this section, with the DD domain input-output relationship in mind, we
formulate the channel estimation problem as a 1D sparse signal recovery (SSR)
problem.

###  III-A Channel Estimation Problem Formulation

To facilitate channel estimation, the system model in (1) can rewritten as

${\bf{y=\Phi}}\left({\mathbf{k}_{{\nu}}},{\mathbf{l}_{{\tau}}}\right){{\tilde{\mathbf{h}}+\mathbf{z}}},\vspace{-2mm}$
(10)

where $\mathbf{y}\in\mathbb{C}^{MN\times 1}$ is the vectorized form of
$y\left[{k,l}\right]$ and vector
$\mathbf{z}\sim\mathcal{CN}\left(\mathbf{0},\sigma^{2}\mathbf{I}_{MN}\right)\in\mathbb{C}^{MN\times
1}$ is the vectorized form of $z\left[{k,l}\right]$. Vectors
${\mathbf{k}_{{\nu}}}=\left[k_{\nu_{1}},\ldots,k_{\nu_{P}}\right]^{\mathrm{T}}\in\mathbb{R}^{P\times
1}$,
${\mathbf{l}_{{\tau}}}=\left[l_{\tau_{1}},\ldots,l_{\tau_{P}}\right]^{\mathrm{T}}\in\mathbb{R}^{P\times
1}$, and
$\tilde{\mathbf{h}}=\left[{\tilde{h}_{1}},\ldots,{\tilde{h}_{P}}\right]^{\mathrm{T}}\in\mathbb{C}^{P\times
1}$ collect the normalized Doppler shifts, the normalized delay shifts, and
the channel coefficients, respectively. Matrix
${\bf{\Phi}}\left({\mathbf{k}_{{\nu}}},{\mathbf{l}_{{\tau}}}\right)\in\mathbb{C}^{MN\times
P}$ is the measurement matrix, which can be formulated as

${\bf{\Phi}}\left({\mathbf{k}_{{\nu}}},{\mathbf{l}_{{\tau}}}\right)=\left[\bm{\phi}\left(k_{\nu_{1}},l_{\tau_{1}}\right),\ldots,\bm{\phi}\left(k_{\nu_{P}},l_{\tau_{P}}\right)\right],\vspace{-2mm}$
(11)

where $\bm{\phi}\left(k_{\nu_{i}},l_{\tau_{i}}\right)\in\mathbb{C}^{MN\times
1}$ and its $\left(kM+l\right)$-th entry can be given by

$\\{\bm{\phi}\left(k_{\nu_{i}},l_{\tau_{i}}\right)\\}_{kM+l}=\sum\limits_{k^{\prime}=0}^{N-1}{\sum\limits_{l^{\prime}=0}^{M-1}{x\left[{k^{\prime},l^{\prime}}\right]w(k-k^{\prime}-k_{\nu_{i}},l-l^{\prime}-l_{\tau_{i}})}}.\vspace{-2mm}$
(12)

The channel estimation problem is to estimate the normalized Doppler shifts
${\mathbf{k}_{{\nu}}}$, the normalized delay shifts ${\mathbf{l}_{{\tau}}}$,
and the channel coefficients $\tilde{\mathbf{h}}$ based on the function
${\bf{\Phi}}\left({\mathbf{k}_{{\nu}}},{\mathbf{l}_{{\tau}}}\right)$ and the
received signal $\bf{y}$. One observation is that when the number of path $P$
is unknown, the length of $\tilde{\mathbf{h}}$ is not fixed and thus
estimating $\tilde{\mathbf{h}}$ is a challenging task. Besides, since the
measurement matrix
${\bf{\Phi}}\left({\mathbf{k}_{{\nu}}},{\mathbf{l}_{{\tau}}}\right)$ depends
on the unknown variables ${\mathbf{k}_{{\nu}}}$ and ${\mathbf{l}_{{\tau}}}$,
conventional methods, such as least square (LS) or minimum mean square error
(MMSE) channel estimation are not applicable to the considered case. Moreover,
the measurement matrix
${\bf{\Phi}}\left({\mathbf{k}_{{\nu}}},{\mathbf{l}_{{\tau}}}\right)$ also
depends on the arrangement of data and pilot symbols in the DD domain
$x\left[{k^{\prime},l^{\prime}}\right]$. For simplicity, we only exploit the
pilot symbols in the measurement matrix
${\bf{\Phi}}\left({\mathbf{k}_{{\nu}}},{\mathbf{l}_{{\tau}}}\right)$ for
channel estimation. In particular, the impacts of data symbols in the
measurement matrix
${\bf{\Phi}}\left({\mathbf{k}_{{\nu}}},{\mathbf{l}_{{\tau}}}\right)$ are
treated as the model error which is assumed to behave as independent Gaussian
noise and can be absorbed into the measurement noise [10]. Thus, the received
signal corresponding to the pilot can be obtained by

$y\left[{k,l}\right]=\sum\limits_{k^{\prime}\in\mathcal{K}_{\mathrm{p}}}{\sum\limits_{l^{\prime}\in\mathcal{L}_{\mathrm{p}}}{x\left[{k^{\prime},l^{\prime}}\right]}}{h_{w}}\left[{k-k^{\prime},l-l^{\prime}}\right]+\tilde{z}\left[k,l\right],\vspace{-2mm}$
(13)

and the measurement matrix can be rewritten as

$\\{\bm{\phi}\left(k_{\nu_{i}},l_{\tau_{i}}\right)\\}_{kM+l}=\sum\limits_{k^{\prime}\in\mathcal{K}_{\mathrm{p}}}{\sum\limits_{l^{\prime}\in\mathcal{L}_{\mathrm{p}}}{x\left[{k^{\prime},l^{\prime}}\right]w(k-k^{\prime}-k_{\nu_{i}},l-l^{\prime}-l_{\tau_{i}})}},\vspace{-2mm}$
(14)

where $\mathcal{K}_{\mathrm{p}}$ and $\mathcal{L}_{\mathrm{p}}$ denote the
index sets of pilot symbols in the Doppler and delay domains, respectively,
and $\tilde{z}\left[k,l\right]$ denotes the new noise variable absorbing the
interference of data symbols. Note that the proposed model does not specify
the pilot arrangement. Specifically, when
$\left|\mathcal{K}_{\mathrm{p}}\right|=\left|\mathcal{L}_{\mathrm{p}}\right|=1$,
the model in (14) degenerates to a single pilot impulse case in [6]. When
$\left|\mathcal{K}_{p}\right|=\left|\mathcal{L}_{p}\right|>1$, it is
equivalent to multiple pulses or the pilot sequence cases in [7]. For
notational simplicity, we consider a single pulse case, i.e.,
$\mathcal{K}_{\mathrm{p}}=\\{k_{\mathrm{p}}\\}$ and
$\mathcal{L}_{\mathrm{p}}=\\{l_{\mathrm{p}}\\}$, where a single pilot impulse
is inserted in the $\left[k_{\mathrm{p}},l_{\mathrm{p}}\right]$-th DD grid
point.

Due to the compactness of the DD domain channel, i.e.,
${\tau_{i}}\in\left(0,\tau_{\mathrm{max}}\right)$ and
${\nu_{i}}\in\left(-\nu_{\mathrm{max}},\nu_{\mathrm{max}}\right)$, only a
fraction of the symbols in $\mathbf{y}$ are affected the single pilot impulse.
To reduce the computational complexity, we consider a truncated
$\mathbf{y}_{\mathrm{T}}$ to estimate the channel, which is actually a
sufficient statistic for channel estimation of OTFS systems. In particular,
given the maximum Doppler shift $\nu_{\mathrm{max}}$ and the maximum delay
shift $\tau_{\mathrm{max}}$, we only need to consider the received signal
$y\left[k,l\right]$ in the range of $k_{\mathrm{p}}-k_{\mathrm{max}}\leq k\leq
k_{\mathrm{p}}+k_{\mathrm{max}}$ and $l_{\mathrm{p}}\leq l\leq
l_{\mathrm{p}}+k_{\mathrm{max}}$, where
$k_{\mathrm{max}}=\left\lceil\nu_{\mathrm{max}}NT\right\rceil$ and
$l_{\mathrm{max}}=\left\lceil\tau_{\mathrm{max}}M\Delta_{f}\right\rceil$. As a
result, we have

$\mathbf{y}_{\mathrm{T}}=\mathbf{\Phi}_{\mathrm{T}}\left({\mathbf{k}_{{\nu}}},{\mathbf{l}_{{\tau}}}\right){{\tilde{\mathbf{h}}+\tilde{\mathbf{z}}_{\mathrm{T}}}},\vspace{-2mm}$
(15)

where
$\mathbf{y}_{\mathrm{T}}\in\mathbb{C}^{M_{\mathrm{T}}N_{\mathrm{T}}\times 1}$
and
$\tilde{\mathbf{z}}_{\mathrm{T}}\in\mathbb{C}^{M_{\mathrm{T}}N_{\mathrm{T}}\times
1}$ are obtained by
$\\{\mathbf{y}_{\mathrm{T}}\\}_{kM_{\mathrm{T}}+l}=y\left[k,l\right]$ and
$\\{\tilde{\mathbf{z}}_{\mathrm{T}}\\}_{kM_{\mathrm{T}}+l}=\tilde{z}\left[k,l\right]$,
respectively, with $N_{\mathrm{T}}=2k_{\mathrm{max}}+1$ and
$M_{\mathrm{T}}=l_{\mathrm{max}}+1$. The truncated measurement matrix is given
by

${\bf{\Phi}}_{\mathrm{T}}\left({\mathbf{k}_{{\nu}}},{\mathbf{l}_{{\tau}}}\right)=\left[\bm{\phi}_{\mathrm{T}}\left(k_{\nu_{1}},l_{\tau_{1}}\right),\ldots,\bm{\phi}_{\mathrm{T}}\left(k_{\nu_{P}},l_{\tau_{P}}\right)\right]\in\mathbb{C}^{M_{\mathrm{T}}N_{\mathrm{T}}\times
P},\vspace{-2mm}$ (16)

where
$\bm{\phi}_{\mathrm{T}}\left(k_{\nu},l_{\tau}\right)\in\mathbb{C}^{M_{\mathrm{T}}N_{\mathrm{T}}\times
1}$ and it is defined by

$\\{\bm{\phi}_{\mathrm{T}}\left(k_{\nu_{i}},l_{\tau_{i}}\right)\\}_{kM_{\mathrm{T}}+l}={{x\left[{k_{\mathrm{p}},l_{\mathrm{p}}}\right]w(k-k_{\mathrm{p}}-k_{\nu_{i}},l-l_{\mathrm{p}}-l_{\tau_{i}})}},\vspace{-2mm}$
(17)

with $k_{\mathrm{p}}-k_{\mathrm{max}}\leq k\leq
k_{\mathrm{p}}+k_{\mathrm{max}}$ and $l_{\mathrm{p}}\leq l\leq
l_{\mathrm{p}}+k_{\mathrm{max}}$. We can observe that the estimation model in
(15) is non-linear due to the coupling of unknown Doppler and delay shifts
with the measurement matrix, which motivates us to propose the first-order
linear approximation in next section.

###  III-B First-Order Linear Approximation

Following the principle of compressed channel sensing[14], we formulate an SSR
problem via a virtual representation based on a virtual sampling grid in the
DD space. Define a virtual sampling grid in the range of
$\left(0,\tau_{\mathrm{max}}\right)$ and
$\left(-\nu_{\mathrm{max}},\nu_{\mathrm{max}}\right)$ in the DD domain, i.e.,
$\overline{\mathbf{k}}_{\nu}=\left[\overline{k}_{\nu_{0}},\ldots,\overline{k}_{\nu_{M_{\tau}N_{\nu}-1}}\right]\in\mathbb{R}^{M_{\tau}N_{\nu}\times
1}$ and
$\overline{\mathbf{l}}_{\tau}=\left[\overline{l}_{\tau_{0}},\ldots,\overline{l}_{\tau_{M_{\tau}N_{\nu}-1}}\right]\in\mathbb{R}^{M_{\tau}N_{\nu}\times
1}$, where $N_{\nu}$ and $M_{\tau}$ denote the virtual sampling grid size in
the Doppler and delay domains, respectively222Note that increasing the virtual
Doppler and delay resolutions can improve the accuracy of the on-grid channel
estimation schemes. However, the computational complexity will be
prohibitively high. Therefore, our proposed off-grid channel estimation can
mitigate the impact of coarse virtual sampling grid and thus can achieve a
satisfactory performance with a much lower computational complexity. The
impacts of virtual sampling resolution for both on-grid and off-grid methods
will be investigated in Section VI through simulations.. Typically considering
an equally spaced sampling grid with a virtual Doppler resolution
$r_{\nu}=\frac{2\nu_{\mathrm{max}}}{N_{\nu}}$ and a virtual delay resolution
$r_{\tau}=\frac{\tau_{\mathrm{max}}}{M_{\tau}}$, we have
$\overline{k}_{\nu_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}}=k^{\prime\prime}r_{\nu}-\nu_{\mathrm{max}}$
and
$\overline{l}_{\tau_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}}=l^{\prime\prime}r_{\tau}$,
$\forall k^{\prime\prime}\in\left\\{0,\ldots,N_{\nu}-1\right\\}$ and $\forall
l^{\prime\prime}\in\left\\{0,\ldots,M_{\tau}-1\right\\}$. To separate the
estimation of the on-grid and off-grid components of Doppler and delay shifts
based on the constructed virtual sampling grid, we can perform a linear
approximation as follows:

$\displaystyle\bm{\phi}_{\mathrm{T}}\left(k_{\nu_{i}},l_{\tau_{i}}\right)$
$\displaystyle\approx\bm{\phi}_{\mathrm{T}}\left(\overline{k}^{i}_{\nu_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}},\overline{l}^{i}_{\tau_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}}\right)+\bm{\phi}^{\prime}_{{\mathrm{T}},\nu}\left(\overline{k}^{i}_{\nu_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}},\overline{l}^{i}_{\tau_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}}\right)\left(k_{\nu_{i}}-\overline{k}^{i}_{\nu_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}}\right)$
$\displaystyle+\bm{\phi}^{\prime}_{{\mathrm{T}},\tau}\left(\overline{k}^{i}_{\nu_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}},\overline{l}^{i}_{\tau_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}}\right)\left(l_{\tau_{i}}-\overline{l}^{i}_{\tau_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}}\right)$
$\displaystyle=\bm{\phi}_{\mathrm{T}}\left(\overline{k}^{i}_{\nu_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}},\overline{l}^{i}_{\tau_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}}\right)+\bm{\phi}^{\prime}_{{\mathrm{T}},\nu}\left(\overline{k}^{i}_{\nu_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}},\overline{l}^{i}_{\tau_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}}\right){\kappa^{i}_{{\nu_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}}}}$
$\displaystyle+\bm{\phi}^{\prime}_{{\mathrm{T}},\tau}\left(\overline{k}^{i}_{\nu_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}},\overline{l}^{i}_{\tau_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}}\right){\iota^{i}_{{\tau_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}}}},$
(18)

where the approximation is due to the exclusion of the higher order
derivatives. Variables
$\overline{k}^{i}_{\nu_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}}$ and
$\overline{l}^{i}_{\tau_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}}$ denote
the unknown nearest grid points of $k_{\nu_{i}}$ and $l_{\tau_{i}}$,
respectively. Vectors
$\bm{\phi}^{\prime}_{{\mathrm{T}},\nu}\left(k_{\nu},l_{\tau}\right)=\frac{{\partial\bm{\phi}_{\mathrm{T}}\left(k_{\nu},l_{\tau}\right)}}{{\partial{k_{{\nu}}}}}\in\mathbb{C}^{M_{\mathrm{T}}N_{\mathrm{T}}\times
1}$ and
$\bm{\phi}^{\prime}_{{\mathrm{T}},\tau}\left(k_{\nu},l_{\tau}\right)=\frac{{\partial\bm{\phi}_{\mathrm{T}}\left(k_{\nu},l_{\tau}\right)}}{{\partial{l_{{\tau}}}}}\in\mathbb{C}^{M_{\mathrm{T}}N_{\mathrm{T}}\times
1}$ denote the first-order gradient of
$\bm{\phi}_{\mathrm{T}}\left(k_{\nu},l_{\tau}\right)$ w.r.t. $k_{\nu}$ and
$l_{\tau}$, respectively. Also, the $\left(kM_{\tau}+l\right)$-th entries of
$\bm{\phi}^{\prime}_{{\mathrm{T}},\nu}\left(k_{\nu},l_{\tau}\right)$ and
$\bm{\phi}^{\prime}_{{\mathrm{T}},\tau}\left(k_{\nu},l_{\tau}\right)$ are
given by

$\displaystyle\\{\bm{\phi}^{\prime}_{{\mathrm{T}},\nu}\left(k_{\nu},l_{\tau}\right)\\}_{kM+l}$
$\displaystyle={{x\left[{k_{\mathrm{p}},l_{\mathrm{p}}}\right]w_{\tau}\left(l-l_{\mathrm{p}}-l_{\tau}\right)w^{\prime}_{\nu}\left(k-k_{\mathrm{p}}-k_{\nu}\right)}}\;\text{and}$
$\displaystyle\\{\bm{\phi}^{\prime}_{{\mathrm{T}},\tau}\left(k_{\nu},l_{\tau}\right)\\}_{kM+l}$
$\displaystyle={{x\left[{k_{\mathrm{p}},l_{\mathrm{p}}}\right]w^{\prime}_{\tau}\left(l-l_{\mathrm{p}}-l_{\tau}\right)w_{\nu}\left(k-k_{\mathrm{p}}-k_{\nu}\right)}},$
(19)

respectively, where

$\displaystyle w^{\prime}_{\nu}\left(k-k_{\mathrm{p}}-{k_{{\nu}}}\right)$
$\displaystyle=\frac{{\partial
w_{\nu}\left(k-k_{\mathrm{p}}-k_{\nu}\right)}}{{\partial{k_{{\nu}}}}}=\frac{1}{N}\sum\limits_{n=0}^{N-1}{\frac{{j2\pi
n}}{N}{e^{-j2\pi
nT\frac{{\left({k-k_{\mathrm{p}}-{k_{\nu}}}\right)}}{{NT}}}}}\;\text{and}$
$\displaystyle w^{\prime}_{\tau}\left(l-l_{\mathrm{p}}-l_{\tau}\right)$
$\displaystyle=\frac{{\partial
w_{\tau}\left(l-l_{\mathrm{p}}-l_{\tau}\right)}}{{\partial{l_{\tau}}}}=-\frac{1}{M}\sum\limits_{m=0}^{M-1}{\frac{{j2\pi
m}}{M}{e^{j2\pi
m\Delta_{f}\frac{{\left({l-l_{\mathrm{p}}-{l_{\tau}}}\right)}}{{M\Delta_{f}}}}}}.$
(20)

Variables ${\kappa^{i}_{{\nu_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}}}}$
and ${\iota^{i}_{{\tau_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}}}}$
represent the off-grid components of the Doppler and delay shifts of the
$i$-th path, respectively, and they are defined as:

${\kappa^{i}_{{\nu_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}}}}=k_{\nu_{i}}-\overline{k}^{i}_{\nu_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}}\in\left[-\frac{1}{2}r_{\nu},\frac{1}{2}r_{\nu}\right]\;\text{and}\;{\iota^{i}_{{\tau_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}}}}=l_{\tau_{i}}-\overline{l}^{i}_{\tau_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}}\in\left[-\frac{1}{2}r_{\tau},\frac{1}{2}r_{\tau}\right],\vspace{-2mm}$
(21)

respectively. Now, substituting (18) into the system model in (15), we have

$\displaystyle\mathbf{y}_{\mathrm{T}}$
$\displaystyle\approx\sum_{i=1}^{P}\left[\bm{\phi}_{\mathrm{T}}\left(\overline{k}^{i}_{\nu_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}},\overline{l}^{i}_{\tau_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}}\right)+\bm{\phi}^{\prime}_{{\mathrm{T}},\nu}\left(\overline{k}^{i}_{\nu_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}},\overline{l}^{i}_{\tau_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}}\right){\kappa^{i}_{{\nu_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}}}}\right.$
$\displaystyle\left.+\bm{\phi}^{\prime}_{{\mathrm{T}},\tau}\left(\overline{k}^{i}_{\nu_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}},\overline{l}^{i}_{\tau_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}}\right){\iota^{i}_{{\tau_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}}}}\right]\tilde{h}_{i}+\tilde{\mathbf{z}}_{\mathrm{T}}.$
(22)

We can observe that the Doppler and delay shifts can be estimated via
estimating the on-grid components, i.e., the nearest grid point
$\left(\overline{k}^{i}_{\nu_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}},\overline{l}^{i}_{\tau_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}}\right)$,
and the off-grid components, i.e.,
$\left({\kappa^{i}_{{\nu_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}}}},{\iota^{i}_{{\tau_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}}}}\right)$,
separately. After obtaining the on-grid and off-grid components, we can
recover the Doppler and delay shift as
$k_{\nu_{i}}=\overline{k}^{i}_{\nu_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}}+{\kappa^{i}_{{\nu_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}}}}$
and
$l_{\tau_{i}}=\overline{l}^{i}_{\tau_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}}+{\iota^{i}_{{\tau_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}}}}$,
respectively. This separation enables the possibility to recast the channel
estimation as a 1D off-grid SSR problem.

###  III-C 1D Off-grid SSR Model

In (III-B), the nearest grid point
$\left(\overline{k}^{i}_{\nu_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}},\overline{l}^{i}_{\tau_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}}\right)$,
the corresponding off-grid component
$\left({\kappa^{i}_{{\nu_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}}}},{\iota^{i}_{{\tau_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}}}}\right)$,
and the channel coefficient $\tilde{h}_{i}$, $\forall i$, are unknown and
remained to be estimated. In other words, the goal is to estimate
$\left\\{\overline{k}^{i}_{\nu_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}},\overline{l}^{i}_{\tau_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}},{\kappa^{i}_{{\nu_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}}}},{\iota^{i}_{{\tau_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}}}},\tilde{h}_{i}\right\\}$
given $\bf{y}$ and the functions
${\bm{\phi}}_{\mathrm{T}}\left({{k}}_{\nu},{{l}}_{{\tau}}\right)$,
$\bm{\phi}^{\prime}_{{\mathrm{T}},\nu}\left(k_{\nu},l_{\tau}\right)$ and
$\bm{\phi}^{\prime}_{{\mathrm{T}},\tau}\left(k_{\nu},l_{\tau}\right)$. As
$P\ll M_{\tau}N_{\nu}$ usually holds, we recast (III-B) as an SSR problem
based on the virtual sampling grid:

$\mathbf{y}_{\mathrm{T}}=\mathbf{\overline{\Phi}}_{\mathrm{T}}\left({\bm{\kappa}_{{\nu}}},{\bm{\iota}_{{\tau}}}\right)\overline{\mathbf{h}}+\overline{\mathbf{z}}_{\mathrm{T}},\vspace{-2mm}$
(23)

where the linear approximation error in (III-B) is also absorbed into the
noise term. The new measurement matrix can be formulated as

${\bf\overline{\Phi}}_{\mathrm{T}}\left({\bm{\kappa}_{{\nu}}},{\bm{\iota}_{{\tau}}}\right)={\bf\Phi}_{\mathrm{T}}+{\bf\Phi}_{\mathrm{T},\nu}\operatorname{\mathrm{diag}}\left({\bm{\kappa}_{{\nu}}}\right)+{\bf\Phi}_{\mathrm{T},\tau}\operatorname{\mathrm{diag}}\left({\bm{\iota}_{{\tau}}}\right)\in\mathbb{C}^{M_{\mathrm{T}}N_{\mathrm{T}}\times
M_{\tau}N_{\nu}},\vspace{-2mm}$ (24)

where ${\bf\Phi}_{\mathrm{T}}\in\mathbb{C}^{M_{\mathrm{T}}N_{\mathrm{T}}\times
M_{\tau}N_{\nu}}$,
${\bf\Phi}_{\mathrm{T},\nu}\in\mathbb{C}^{M_{\mathrm{T}}N_{\mathrm{T}}\times
M_{\tau}N_{\nu}}$,
${\bf\Phi}_{\mathrm{T},\tau}\in\mathbb{C}^{M_{\mathrm{T}}N_{\mathrm{T}}\times
M_{\tau}N_{\nu}}$, $\bm{\kappa}_{{\nu}}\in\mathbb{R}^{M_{\tau}N_{\nu}\times
1}$, and ${\bm{\iota}_{{\tau}}}\in\mathbb{R}^{M_{\tau}N_{\nu}\times 1}$ are
defined based on the virtual sampling grid $\overline{\mathbf{k}}_{\nu}$ and
$\overline{\mathbf{l}}_{\tau}$, and they are given by

$\displaystyle{\bf\Phi}_{\mathrm{T}}$
$\displaystyle=\left[\bm{\phi}_{\mathrm{T}}\left(\overline{k}_{\nu_{0}},\overline{l}_{\tau_{0}}\right),\ldots,\bm{\phi}_{\mathrm{T}}\left(\overline{k}_{\nu_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}},\overline{l}_{\tau_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}}\right),\ldots,\bm{\phi}_{\mathrm{T}}\left(\overline{k}_{\nu_{M_{\tau}N_{\nu}-1}},\overline{l}_{\tau_{M_{\tau}N_{\nu}-1}}\right)\right],$
(25) $\displaystyle{\bf\Phi}_{\mathrm{T},\nu}$
$\displaystyle=\left[\bm{\phi}^{\prime}_{\mathrm{T},\nu}\left(\overline{k}_{\nu_{0}},\overline{l}_{\tau_{0}}\right),\ldots,\bm{\phi}^{\prime}_{\mathrm{T},\nu}\left(\overline{k}_{\nu_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}},\overline{l}_{\tau_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}}\right),\ldots,\bm{\phi}^{\prime}_{\mathrm{T},\nu}\left(\overline{k}_{\nu_{M_{\tau}N_{\nu}-1}},\overline{l}_{\tau_{M_{\tau}N_{\nu}-1}}\right)\right],$
$\displaystyle{\bf\Phi}_{\mathrm{T},\tau}$
$\displaystyle=\left[\bm{\phi}^{\prime}_{\mathrm{T},\tau}\left(\overline{k}_{\nu_{0}},\overline{l}_{\tau_{0}}\right),\ldots,\bm{\phi}^{\prime}_{\mathrm{T},\tau}\left(\overline{k}_{\nu_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}},\overline{l}_{\tau_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}}\right),\ldots,\bm{\phi}^{\prime}_{\mathrm{T},\tau}\left(\overline{k}_{\nu_{M_{\tau}N_{\nu}-1}},\overline{l}_{\tau_{M_{\tau}N_{\nu}-1}}\right)\right],$
$\displaystyle\bm{\kappa}_{{\nu}}$
$\displaystyle=\left[{\kappa_{{\nu_{0}}}},\ldots,{\kappa_{{\nu_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}}}},\ldots,{\kappa_{{\nu_{M_{\tau}N_{\nu}-1}}}}\right]^{\mathrm{T}},\;\text{and}\;{\bm{\iota}_{{\tau}}}=\left[{\iota_{{\tau_{0}}}},\ldots,{\iota_{{\tau_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}}}},\ldots,{\iota_{{\tau_{M_{\tau}N_{\nu}-1}}}}\right]^{\mathrm{T}},$

respectively. The new noise vector is defined as
$\overline{\mathbf{z}}_{\mathrm{T}}\sim\mathcal{CN}\left(\mathbf{0},\overline{\sigma^{2}}\mathbf{I}_{M_{\mathrm{T}}N_{\mathrm{T}}}\right)\in\mathbb{C}^{M_{\mathrm{T}}N_{\mathrm{T}}\times
1}$, where $\overline{\sigma^{2}}$ is the unknown noise power including the
background noise power, the interference power due to the unknown data
symbols, as well as the modeling error induced by our proposed first-order
approximation in (III-B).

In (23), the unknown vectors $\overline{\bf{h}}$, $\bm{\kappa}_{{\nu}}$, and
${\bm{\iota}_{{\tau}}}$ are remained to be estimated. We can observe that when
$k^{\prime\prime}=\overline{k}^{i}_{\nu_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}}$
and
$l^{\prime\prime}=\overline{l}^{i}_{\tau_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}}$,
$\forall i$, these unknown vectors can be defined as
$\\{\overline{\bf{h}}\\}_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}={\tilde{h}_{i}}$,
${\kappa_{{\nu_{{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}}}}}=k_{\nu_{i}}-\overline{k}^{i}_{\nu_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}}$,
and
${\iota_{{\tau_{{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}}}}}=l_{\tau_{i}}-\overline{l}^{i}_{\tau_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}}$.
The remain entries in $\overline{\bf{h}}$, $\bm{\kappa}_{{\nu}}$, and
${\bm{\iota}_{{\tau}}}$ are all zeros. In other words, $\overline{\bf{h}}$,
$\bm{\kappa}_{{\nu}}$, ${\bm{\iota}_{{\tau}}}$ are all $P$-sparse channel
vectors with the same support, i.e., joint sparsity [10]. The on-grid
components are determined by the indices of the non-zero entries in
$\overline{\bf{h}}$, while the off-grid components can be given by the
estimated parameters $\bm{\kappa}_{{\nu}}$ and ${\bm{\iota}_{{\tau}}}$, as
detailed in the next section. We can observe that the proposed model
degenerates to the case of on-grid sparse channel estimation [9] when forcing
the off-grid components as zeros, i.e.,
${\bm{\kappa}_{{\nu}}}={\bm{\iota}_{{\tau}}}=\mathbf{0}$.

##  IV 1D Off-grid SBL-based Channel Estimation Algorithm

In this section, we propose the off-grid SBL algorithm for estimating the
unknown sparse vectors in the 1D SSR problem in (23).

###  IV-A Hierarchical Hyperprior Distribution

To exploit the sparsity of the original DD domain channel response, we
formulate a hierarchical hyperprior distribution following the SBL framework
[10], as illustrated in Fig. 4. In particular, the prior distribution of
unknown channel vector $\overline{\bf{h}}$ is assumed to follow

$p\left({\overline{\bf{h}}\left|{\bm{\alpha}}\right.}\right)=\mathcal{CN}\left({\overline{\bf{h}}\left|{{\bf{0}},{\bf{\Lambda}}}\right.}\right),\vspace{-2mm}$
(26)

where the covariance matrix
${\bf{\Lambda}}=\operatorname{\mathrm{diag}}\left({\bm{\alpha}}\right)$ is a
diagonal matrix and
${\bm{\alpha}}=\left[{{\alpha_{0}},\ldots,{\alpha_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}},\ldots{\alpha_{M_{\tau}N_{\nu}-1}}}\right]^{\mathrm{T}}\in\mathbb{R}^{M_{\tau}N_{\nu}\times
1}$ with ${\alpha_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}}\geq 0$ is the
hyper-parameter to model the sparsity of $\overline{\bf{h}}$. In particular,
when ${\alpha_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}}=0$, the
corresponding channel coefficient
$\\{\overline{\bf{h}}\\}_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}=0$.
Otherwise, the corresponding channel coefficient is non-zero. To avoid the
impact of the initialization of ${\bm{\alpha}}$ on the channel estimation
performance, we further assume that the hyper-parameters ${\bm{\alpha}}$ are
independent and follow the Gamma distribution [10]:

$p\left({{\bm{\alpha}}\left|\rho\right.}\right)=\mathop{\prod}\limits_{k^{\prime\prime}=0}^{N_{\nu}-1}\mathop{\prod}\limits_{l^{\prime\prime}=0}^{M_{\tau}-1}\Gamma\left({{\alpha_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}}\left|{1,\rho}\right.}\right),\vspace{-2mm}$
(27)

where $\rho>0$ is a fixed root parameter for ${\bm{\alpha}}$. The hyperprior
formulation defined in (26) and (27) results in a Laplace prior, see (17) in
[11], and hence a sparse estimation of $\overline{\bf{h}}$ is preferred.

Figure 4: A hierarchical hyperprior distribution for the proposed off-grid SBL
channel estimation.

Moreover, we assume that the noise vector $\overline{\bf{z}}_{\mathrm{T}}$
follows a complex Gaussian distribution

$p\left({\overline{\bf{z}}_{\mathrm{T}}\left|{{\beta_{0}}}\right.}\right)=\mathcal{CN}\left({\overline{\bf{z}}_{\mathrm{T}}\left|{{\bf{0}},\beta_{0}^{-1}{{\bf{I}}_{M_{\mathrm{T}}N_{\mathrm{T}}}}}\right.}\right),\vspace{-2mm}$
(28)

where $\beta_{0}$ denotes the inverse of the unknown noise variance
$\overline{\sigma^{2}}$ and it further follows:

$p\left({{\beta_{0}}\left|{c,d}\right.}\right)=\Gamma\left({{\beta_{0}}\left|{c,d}\right.}\right),\vspace{-2mm}$
(29)

with $c,d>0$ are fixed root parameters. Furthermore, unknown off-grid Doppler
and delay variables follow a uniform distribution within their bounds in (21),
i.e.,

${\kappa_{{\nu_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}}}}\sim\mathcal{U}\left[{-\frac{1}{2}r_{\nu},\frac{1}{2}r_{\nu}}\right]\;\text{and}\;{\iota_{{\tau_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}}}}\sim\mathcal{U}\left[{-\frac{1}{2}r_{\tau},\frac{1}{2}}r_{\tau}\right],\forall
k^{\prime\prime},l^{\prime\prime}.\vspace{-2mm}$ (30)

As shown in Fig. 4, the off-grid Doppler and delay
$\left({\bm{\kappa}_{{\nu}}},{\bm{\iota}_{{\tau}}}\right)$ are modeled as
hyper-parameters, which will be updated concurrently with the posterior
distribution of channel vector $\overline{\bf{h}}$ in our proposed algorithm
in next section.

###  IV-B Off-grid SBL-based Channel Estimation Algorithm

According to the Bayesian principle, the optimal estimation of the unknown
variables in (23) can be obtained from the posterior distribution
$p\left({\overline{\bf{h}},{\bm{\alpha}},{{\bm{\kappa}}_{\nu}},{{\bm{\iota}}_{\tau}},{\beta_{0}}}\left|{\bf{y}}_{\mathrm{T}}\right.\right)$,
which includes the full information of both their prior distribution and the
measurement $\mathbf{y}_{\mathrm{T}}$. However, the posterior distribution
$p\left({\overline{\bf{h}},{\bm{\alpha}},{{\bm{\kappa}}_{\nu}},{{\bm{\iota}}_{\tau}},{\beta_{0}}}\left|{\bf{y}}_{\mathrm{T}}\right.\right)$
cannot be expressed explicitly due to the non-linearly coupled unknown
variables in the measurement matrix
$\overline{\bf{\Phi}}_{\mathrm{T}}\left({{{\bm{\kappa}}_{\nu}},{{\bm{\iota}}_{\tau}}}\right)$
in (23). The unknown hyper-parameters ${\bm{\alpha}}$ and ${{\beta_{0}}}$,
which determine the prior distribution of $\overline{\bf{h}}$ and the
likelihood function
$p\left({{\bf{y}}_{\mathrm{T}}\left|\overline{\bf{h}}\right.}\right)$,
respectively, are also the main obstacles to derive the posterior
distribution. Following the Bayes’ rule, the posterior distribution can be
decomposed as

$p\left({\overline{\bf{h}},{\bm{\alpha}},{{\bm{\kappa}}_{\nu}},{{\bm{\iota}}_{\tau}},{\beta_{0}}}\left|{\bf{y}}_{\mathrm{T}}\right.\right)=p\left(\overline{\bf{h}}\left|{\bf{y}}_{\mathrm{T}},{\bm{\alpha}},{{\bm{\kappa}}_{\nu}},{{\bm{\iota}}_{\tau}},{\beta_{0}}\right.\right)p\left({\bm{\alpha}},{{{\bm{\kappa}}_{\nu}},{{\bm{\iota}}_{\tau}},{\beta_{0}}}\left|{\bf{y}}_{\mathrm{T}}\right.\right).\vspace{-2mm}$
(31)

Inspired by this decomposition, we prefer to derive the conditional posterior
distribution of $\overline{\bf{h}}$ for given hyper-parameters
$\left({\bm{\alpha}},{{\bm{\kappa}}_{\nu}},{{\bm{\iota}}_{\tau}},{\beta_{0}}\right)$,
i.e.,
$p\left({\overline{\bf{h}}\left|{{\bf{y}}_{\mathrm{T}};{\bm{\alpha}},{{\bm{\kappa}}_{\nu}},{{\bm{\iota}}_{\tau}},{\beta_{0}}}\right.}\right)$,
and then update the hyper-parameters
$\left({\bm{\alpha}},{{\bm{\kappa}}_{\nu}},{{\bm{\iota}}_{\tau}},{\beta_{0}}\right)$
to maximize the posterior distribution
$p\left({\bm{\alpha}},{{{\bm{\kappa}}_{\nu}},{{\bm{\iota}}_{\tau}},{\beta_{0}}}\left|{\bf{y}}_{\mathrm{T}}\right.\right)$.
In turn, the updated point estimates of the hyper-parameters can be further
adopted to update the conditional posterior distribution. Such an alternative
inference calculation will be proceeded repeatedly until convergence.

In the $t$-th iteration, given the hyper-parameters
$\left({{{\bm{\alpha}}^{\left(t\right)}},{\bm{\kappa}}_{\nu}^{\left(t\right)},{\bm{\iota}}_{\tau}^{\left(t\right)},\beta_{0}^{\left(t\right)}}\right)$,
the conditional posterior distribution of $\overline{\bf{h}}$ is given by

$\displaystyle
p\left({\overline{\bf{h}}\left|{{\bf{y}}_{\mathrm{T}};{{\bm{\alpha}}^{\left(t\right)}},{{\bm{\kappa}}_{\nu}^{\left(t\right)}},{{\bm{\iota}}_{\tau}^{\left(t\right)}},{\beta_{0}^{\left(t\right)}}}\right.}\right)=\frac{{p\left({{\bf{y}}_{\mathrm{T}},\overline{\bf{h}};{{\bm{\alpha}}^{\left(t\right)}},{{\bm{\kappa}}_{\nu}^{\left(t\right)}},{{\bm{\iota}}_{\tau}^{\left(t\right)}},{\beta_{0}^{\left(t\right)}}}\right)}}{{p\left({{\bf{y}}_{\mathrm{T}};{{\bm{\alpha}}^{\left(t\right)}},{{{\bm{\kappa}}_{\nu}^{\left(t\right)}},{{\bm{\iota}}_{\tau}^{\left(t\right)}},{\beta_{0}^{\left(t\right)}}}}\right)}}$
$\displaystyle=$
$\displaystyle\frac{{p\left({{\bf{y}}_{\mathrm{T}}\left|{\overline{\bf{h}};{{\bm{\kappa}}_{\nu}^{\left(t\right)}},{{\bm{\iota}}_{\tau}^{\left(t\right)}},{\beta_{0}^{\left(t\right)}}}\right.}\right)p\left({\overline{\bf{h}};{\bm{\alpha}^{\left(t\right)}}}\right)}}{{p\left({{\bf{y}}_{\mathrm{T}};{{\bm{\alpha}}^{\left(t\right)}},{{{\bm{\kappa}}_{\nu}^{\left(t\right)}},{{\bm{\iota}}_{\tau}^{\left(t\right)}},{\beta_{0}^{\left(t\right)}}}}\right)}}\propto{p\left({{\bf{y}}_{\mathrm{T}}\left|{\overline{\bf{h}};{{\bm{\kappa}}_{\nu}^{\left(t\right)}},{{\bm{\iota}}_{\tau}^{\left(t\right)}},{\beta_{0}^{\left(t\right)}}}\right.}\right)p\left({\overline{\bf{h}};{\bm{\alpha}^{\left(t\right)}}}\right)},$
(32)

where the prior distribution
$p\left({\overline{\bf{h}};{\bm{\alpha}^{\left(t\right)}}}\right)$ is given by
(26) and
$p\left({{\bf{y}}_{\mathrm{T}}\left|{\overline{\bf{h}};{{\bm{\kappa}}_{\nu}^{\left(t\right)}},{{\bm{\iota}}_{\tau}^{\left(t\right)}},{\beta_{0}^{\left(t\right)}}}\right.}\right)$
denotes the likelihood function for given hyper-parameters
$\left({{{\bm{\kappa}}_{\nu}^{\left(t\right)}},{{\bm{\iota}}_{\tau}^{\left(t\right)}},{\beta_{0}^{\left(t\right)}}}\right)$
and it is given by

$p\left({{\bf{y}}_{\mathrm{T}}\left|{\overline{\bf{h}};{{\bm{\kappa}}_{\nu}^{\left(t\right)}},{{\bm{\iota}}_{\tau}^{\left(t\right)}},{\beta_{0}^{\left(t\right)}}}\right.}\right)=\mathcal{CN}\left({{\bf{y}}_{\mathrm{T}}\left|{\overline{\bf{\Phi}}_{\mathrm{T}}\left({{{\bm{\kappa}}_{\nu}^{\left(t\right)}},{{\bm{\iota}}_{\tau}^{\left(t\right)}}}\right)\overline{\bf{h}},\left(\beta_{0}^{\left(t\right)}\right)^{-1}{{\bf{I}}_{{M_{\mathrm{T}}N_{\mathrm{T}}}}}}\right.}\right).\vspace{-2mm}$
(33)

Note that
${p\left({{\bf{y}}_{\mathrm{T}};{{\bm{\alpha}}^{\left(t\right)}},{{{\bm{\kappa}}_{\nu}^{\left(t\right)}},{{\bm{\iota}}_{\tau}^{\left(t\right)}},{\beta_{0}^{\left(t\right)}}}}\right)}$
in (33) is a normalized term, which is omitted to facilitate the derivation.
In fact,
${p\left({{\bf{y}}_{\mathrm{T}};{{\bm{\alpha}}^{\left(t\right)}},{{{\bm{\kappa}}_{\nu}^{\left(t\right)}},{{\bm{\iota}}_{\tau}^{\left(t\right)}},{\beta_{0}^{\left(t\right)}}}}\right)}$
can be interpreted as the model evidence, which measures the possibility of
the observations ${\bf{y}}_{\mathrm{T}}$ given the hyper-parameters.

Combining (26), (IV-B), and (33), the conditional posterior distribution of
$\overline{\bf{h}}$ can be obtained by

$p\left({\overline{\bf{h}}\left|{{\bf{y}}_{\mathrm{T}};{{\bm{\alpha}}^{\left(t\right)}},{{\bm{\kappa}}_{\nu}^{\left(t\right)}},{{\bm{\iota}}_{\tau}^{\left(t\right)}},{\beta_{0}^{\left(t\right)}}}\right.}\right)=\mathcal{CN}\left(\overline{\bf{h}}\left|\bm{\mu}_{\overline{\bf{h}}}^{\left(t\right)},\bm{\Sigma}_{\overline{\bf{h}}}^{\left(t\right)}\right.\right)\vspace{-2mm}$
(34)

and the conditional posterior mean and covariance matrix in the $t$-th
iteration are given by333The estimation of the posterior mean and covariance
matrix in (35) and (36) are the optimal solution of the inverse problem in
(23) for given
$\left({{{\bm{\alpha}}^{\left(t\right)}},{\bm{\kappa}}_{\nu}^{\left(t\right)},{\bm{\iota}}_{\tau}^{\left(t\right)},\beta_{0}^{\left(t\right)}}\right)$.
Note that some low-complexity approximate inference methods[15], such as
variational Bayesian or message passing algorithm, can achieve a suboptimal
solution, which are out of the scope of this paper.

$\displaystyle\bm{\mu}_{\overline{\bf{h}}}^{\left(t\right)}$
$\displaystyle={\beta_{0}^{\left(t\right)}}{{\bf{\Sigma}}_{\overline{\bf{h}}}^{\left(t\right)}}{\overline{\bf{\Phi}}_{\mathrm{T}}^{\rm{H}}}\left({{{\bm{\kappa}}_{\nu}^{\left(t\right)}},{{\bm{\iota}}_{\tau}^{\left(t\right)}}}\right){\bf{y}}_{\mathrm{T}}\;\text{and}$
(35) $\displaystyle\bm{\Sigma}_{\overline{\bf{h}}}^{\left(t\right)}$
$\displaystyle={\left({{\beta_{0}}{{\overline{\bf{\Phi}}_{\mathrm{T}}\left({{{\bm{\kappa}}_{\nu}^{\left(t\right)}},{{\bm{\iota}}_{\tau}^{\left(t\right)}}}\right)}^{\rm{H}}}\overline{\bf{\Phi}}_{\mathrm{T}}\left({{{\bm{\kappa}}_{\nu}^{\left(t\right)}},{{\bm{\iota}}_{\tau}^{\left(t\right)}}}\right)+{\left({\bf{\Lambda}}^{\left(t\right)}\right)^{-1}}}\right)^{-1}}$
(36)
$\displaystyle={\bm{\Lambda}}^{\left(t\right)}-{\bm{\Lambda}}^{\left(t\right)}\overline{\bm{\Phi}}_{\mathrm{T}}{\left({{{\bm{\kappa}}_{\nu}^{\left(t\right)}},{{\bm{\iota}}_{\tau}^{\left(t\right)}}}\right)^{\rm{H}}}$
(37)
$\displaystyle\times{\left({\left(\beta_{0}^{\left(t\right)}\right)^{-1}{{\bf{I}}_{M_{\mathrm{T}}N_{\mathrm{T}}}}+\overline{\bm{\Phi}}_{\mathrm{T}}\left({{{\bm{\kappa}}_{\nu}^{\left(t\right)}},{{\bm{\iota}}_{\tau}^{\left(t\right)}}}\right){\bm{\Lambda}}^{\left(t\right)}\overline{\bm{\Phi}}_{\mathrm{T}}{{\left({{{\bm{\kappa}}_{\nu}^{\left(t\right)}},{{\bm{\iota}}_{\tau}^{\left(t\right)}}}\right)}^{\rm{H}}}}\right)^{-1}}\overline{\bm{\Phi}}_{\mathrm{T}}\left({{{\bm{\kappa}}_{\nu}^{\left(t\right)}},{{\bm{\iota}}_{\tau}^{\left(t\right)}}}\right){\bm{\Lambda}}^{\left(t\right)},$

respectively. Note that the covariance matrix in (37) is obtained by the
matrix inversion lemma, which can significantly reduce the computational
complexity of the matrix inverse operation in (36) as usually $M_{\tau}\gg
M_{\mathrm{T}}$ and $N_{\nu}\gg N_{\mathrm{T}}$.

Then, given the conditional posterior distribution
$p\left({\overline{\bf{h}}\left|{{\bf{y}}_{\mathrm{T}};{{\bm{\alpha}}^{\left(t\right)}},{{\bm{\kappa}}_{\nu}^{\left(t\right)}},{{\bm{\iota}}_{\tau}^{\left(t\right)}},{\beta_{0}^{\left(t\right)}}}\right.}\right)$
in the $t$-th iteration, the hyper-parameters
$\left({\bm{\alpha}},{{\bm{\kappa}}_{\nu}},{{\bm{\iota}}_{\tau}},{\beta_{0}}\right)$
can be updated by maximizing the posterior distribution
$p\left({\bm{\alpha}},{{{\bm{\kappa}}_{\nu}},{{\bm{\iota}}_{\tau}},{\beta_{0}}}\left|{\bf{y}}_{\mathrm{T}}\right.\right)$.
Note that
$p\left({\bm{\alpha}},{{{\bm{\kappa}}_{\nu}},{{\bm{\iota}}_{\tau}},{\beta_{0}}}\left|{\bf{y}}_{\mathrm{T}}\right.\right)\propto
p\left({{\bf{y}}_{\mathrm{T}},{\bm{\alpha}},{{\bm{\kappa}}_{\nu}},{{\bm{\iota}}_{\tau}},{\beta_{0}}}\right)$
as the observation ${\bf{y}}_{\mathrm{T}}$ is independent of the hyper-
parameters. An expectation-maximization (EM) algorithm[16] is employed for
updating the hyper-parameters in the $\left(t+1\right)$-th iteration. In
particular, the expectation step is w.r.t. the joint distribution function
$p\left({{\bf{y}}_{\mathrm{T}},\overline{\bf{h}},{\bm{\alpha}},{{\bm{\kappa}}_{\nu}},{{\bm{\iota}}_{\tau}},{\beta_{0}}}\right)$
via marginalizing over the conditional posterior distribution of
$\overline{\bf{h}}$. According to the Bayes’ rule, the joint distribution
$p\left({{\bf{y}}_{\mathrm{T}},\overline{\bf{h}},{\bm{\alpha}},{{\bm{\kappa}}_{\nu}},{{\bm{\iota}}_{\tau}},{\beta_{0}}}\right)$
can be decomposed as

$p\left({{\bf{y}}_{\mathrm{T}},\overline{\bf{h}},{\bm{\alpha}},{{\bm{\kappa}}_{\nu}},{{\bm{\iota}}_{\tau}},{\beta_{0}}}\right)=p\left({{\bf{y}}_{\mathrm{T}}\left|{\overline{\bf{h}},{{\bm{\kappa}}_{\nu}},{{\bm{\iota}}_{\tau}},{\beta_{0}}}\right.}\right)p\left({\overline{\bf{h}}\left|{\bm{\alpha}}\right.}\right)p\left({\bm{\alpha}}\right)p\left({{{\bm{\kappa}}_{\nu}}}\right)p\left({{{\bm{\iota}}_{\tau}}}\right)p\left({{\beta_{0}}}\right).\vspace{-2mm}$
(38)

Now, the expectation and maximization steps can be defined as follows:

  * •

Expectation step: In this step, we aim to obtain an objective function of the
hyper-parameters
$\left({\bm{\alpha}},{{\bm{\kappa}}_{\nu}},{{\bm{\iota}}_{\tau}},{\beta_{0}}\right)$
via averaging the joint distribution in (38) over the conditional posterior
distribution of $\overline{\bf{h}}$, i.e.,

$Q\left({{\bm{\alpha}},{{\bm{\kappa}}_{\nu}},{{\bm{\iota}}_{\tau}},{\beta_{0}}\left|{{{\bm{\alpha}}^{\left(t\right)}},{\bm{\kappa}}_{\nu}^{\left(t\right)},{\bm{\iota}}_{\tau}^{\left(t\right)},\beta_{0}^{\left(t\right)}}\right.}\right)={E_{\overline{\bf{h}}\left|{{\bf{y}}_{\mathrm{T}};{{\bm{\alpha}}^{\left(t\right)}},{\bm{\kappa}}_{\nu}^{\left(t\right)},{\bm{\iota}}_{\tau}^{\left(t\right)},\beta_{0}^{\left(t\right)}}\right.}}\left\\{{\ln\left[{p\left({{\bf{y}}_{\mathrm{T}},\overline{\bf{h}},{\bm{\alpha}},{{\bm{\kappa}}_{\nu}},{{\bm{\iota}}_{\tau}},{\beta_{0}}}\right)}\right]}\right\\}.\vspace{-2mm}$
(39)

  * •

Maximization step: In this step, we optimize each hyper-parameter
alternatively to maximize the obtained objective function in the expectation
step , i.e.,

$\displaystyle{{\bm{\alpha}}^{\left(t+1\right)}}$
$\displaystyle=\arg\operatorname{maximize}Q\left({{\bm{\alpha}}\left|{{{\bm{\alpha}}^{\left(t\right)}},{\bm{\kappa}}_{\nu}^{\left(t\right)},{\bm{\iota}}_{\tau}^{\left(t\right)},\beta_{0}^{\left(t\right)}}\right.}\right),$
(40) $\displaystyle{\bm{\kappa}}_{\nu}^{\left(t+1\right)}$
$\displaystyle=\arg\operatorname{maximize}Q\left({{{\bm{\kappa}}_{\nu}}\left|{{{\bm{\alpha}}^{\left(t\right)}},{\bm{\kappa}}_{\nu}^{\left(t\right)},{\bm{\iota}}_{\tau}^{\left(t\right)},\beta_{0}^{\left(t\right)}}\right.}\right),$
(41) $\displaystyle{\bm{\iota}}_{\tau}^{\left(t+1\right)}$
$\displaystyle=\arg\operatorname{maximize}Q\left({{{\bm{\iota}}_{\tau}}\left|{{{\bm{\alpha}}^{\left(t\right)}},{\bm{\kappa}}_{\nu}^{\left(t\right)},{\bm{\iota}}_{\tau}^{\left(t\right)},\beta_{0}^{\left(t\right)}}\right.}\right),\;\text{and}$
(42) $\displaystyle\beta_{0}^{\left(t+1\right)}$
$\displaystyle=\arg\operatorname{maximize}Q\left({{\beta_{0}}\left|{{{\bm{\alpha}}^{\left(t\right)}},{\bm{\kappa}}_{\nu}^{\left(t\right)},{\bm{\iota}}_{\tau}^{\left(t\right)},\beta_{0}^{\left(t\right)}}\right.}\right).$
(43)

Following definitions in (39)-(43), the closed-form updating rule for each
hyper-parameter will be derived as follows. In particular, the updating rule
for the hyper-parameter $\left({\bm{\alpha}},{\beta_{0}}\right)$ are stated in
the following proposition.

###### Proposition 1

The hyper-parameter $\left({\bm{\alpha}},{\beta_{0}}\right)$ can be updated by

$\displaystyle\alpha_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}^{\left({t+1}\right)}$
$\displaystyle=\frac{{\sqrt{1+4\rho\left\\{{{\left|{\\{{{\bm{\mu}}^{\left(t\right)}_{\overline{\bf{h}}}}\\}_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}}\right|}^{2}}+\\{{{\bm{\Sigma}}^{\left(t\right)}_{\overline{\bf{h}}}}\\}_{k^{\prime\prime}M_{\tau}+l^{\prime\prime},k^{\prime\prime}M_{\tau}+l^{\prime\prime}}\right\\}}-1}}{{2\rho}}\;\text{and}$
(44) $\displaystyle{\beta_{0}}^{\left({t+1}\right)}$
$\displaystyle=\frac{{c-1+M_{\mathrm{T}}N_{\mathrm{T}}}}{{d+{A_{{\beta_{0}}}}}},$
(45)

respectively, where

${A_{{\beta_{0}}}}={\left\|{{\bf{y}}_{\mathrm{T}}-\overline{\bm{\Phi}}_{\mathrm{T}}\left({{\bm{\kappa}}_{\nu}^{\left(t\right)},{\bm{\iota}}_{\tau}^{\left(t\right)}}\right){\bm{\mu}}_{\overline{\bf{h}}}^{\left(t\right)}}\right\|^{2}}{\rm{+}}\left(\beta^{\left(t\right)}_{0}\right)^{-1}\mathop{\sum}\limits_{k^{\prime\prime}=0}^{N_{\nu}-1}\mathop{\sum}\limits_{l^{\prime\prime}=0}^{M_{\tau}-1}\left({1-\left(\alpha_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}^{\left(t\right)}\right)^{-1}{{\left\\{{{\bf{\Sigma}}_{\overline{\bf{h}}}^{\left(t\right)}}\right\\}}_{k^{\prime\prime},l^{\prime\prime}}}}\right).\vspace{-2mm}$
(46)

###### proof 1

Please refer to Appendix of Proposition 1.

To update off-grid components of the Doppler shift, i.e.,
${{\bm{\kappa}}_{\nu}}$, the objective function in (41) can be obtained by

$\displaystyle
Q\left({{{\bm{\kappa}}_{\nu}}\left|{{{\bm{\alpha}}^{\left(t\right)}},{\bm{\kappa}}_{\nu}^{\left(t\right)},{\bm{\iota}}_{\tau}^{\left(t\right)},\beta_{0}^{\left(t\right)}}\right.}\right)={E_{\overline{\bf{h}}\left|{{\bf{y}}_{\mathrm{T}};{{\bm{\alpha}}^{\left(t\right)}},{\bm{\kappa}}_{\nu}^{\left(t\right)},{\bm{\iota}}_{\tau}^{\left(t\right)},\beta_{0}^{\left(t\right)}}\right.}}\left\\{{\ln\left[{p\left({{\bf{y}}_{\mathrm{T}}\left|{\overline{\bf{h}},{{\bm{\kappa}}_{\nu}},{{\bm{\iota}}_{\tau}^{\left(t\right)}},{\beta_{0}^{\left(t\right)}}}\right.}\right)p\left({{{\bm{\kappa}}_{\nu}}}\right)}\right]}\right\\}$
$\displaystyle=$
$\displaystyle{E_{\overline{\bf{h}}\left|{{\bf{y}}_{\mathrm{T}};{{\bm{\alpha}}^{\left(t\right)}},{\bm{\kappa}}_{\nu}^{\left(t\right)},{\bm{\iota}}_{\tau}^{\left(t\right)},\beta_{0}^{\left(t\right)}}\right.}}\left\\{{\ln\left[{{\cal
C}{\cal
N}\left({{\bf{y}}_{\mathrm{T}}\left|{\overline{\bm{\Phi}}_{\mathrm{T}}\left({{{\bm{\kappa}}_{\nu}},{{\bm{\iota}}_{\tau}^{\left(t\right)}}}\right)\overline{\bf{h}},\left(\beta_{0}^{\left(t\right)}\right)^{-1}{{\bf{I}}_{M_{\mathrm{T}}N_{\mathrm{T}}}}}\right.}\right)}\right]}\right\\}.$
(47)

It can be seen that maximizing the objective function in (47) is equivalent to

$\mathop{\operatorname{minimize}}\limits_{{\bm{\kappa}}_{\nu}}{E_{\overline{\bf{h}}\left|{{\bf{y}}_{\mathrm{T}};{{\bm{\alpha}}^{\left(t\right)}},{\bm{\kappa}}_{\nu}^{\left(t\right)},{\bm{\iota}}_{\tau}^{\left(t\right)},\beta_{0}^{\left(t\right)}}\right.}}\left\\{{{{\left\|{{\bf{y}}_{\mathrm{T}}-\overline{\bm{\Phi}}_{\mathrm{T}}\left({{{\bm{\kappa}}_{\nu}},{\bm{\iota}}_{\tau}^{\left(t\right)}}\right)\overline{\bf{h}}}\right\|}^{2}}}\right\\}\Leftrightarrow\mathop{\operatorname{minimize}}\limits_{{\bm{\kappa}}_{\nu}}\left({{\bm{\kappa}}_{\nu}^{\rm{T}}{{\bf{A}}_{\nu}}{{\bm{\kappa}}_{\nu}}-2{\bf{b}_{\nu}^{\rm{T}}}{{\bm{\kappa}}_{\nu}}}\right),\vspace{-2mm}$
(48)

where ${{\bf{A}}_{\nu}}\in\mathbb{R}^{M_{\tau}N_{\nu}\times M_{\tau}N_{\nu}}$
and ${{\bf{b}}_{{\nu}}}\in\mathbb{R}^{M_{\tau}N_{\nu}\times 1}$ are given by

$\displaystyle{{\bf{A}}_{\nu}}$
$\displaystyle=\Re\left\\{{\bm{\Phi}}_{{\mathrm{T}},\nu}^{\rm{H}}{{\bm{\Phi}}_{{\mathrm{T}},\nu}}\odot\left({\left({\bm{\mu}}_{\overline{\bf{h}}}^{\left(t\right)}\right)^{*}\left({\bm{\mu}}_{\overline{\bf{h}}}^{\left(t\right)}\right)^{\rm{T}}+\left({\bm{\Sigma}}_{\overline{\bf{h}}}^{\left(t\right)}\right)^{\rm{T}}}\right)\right\\}\;\text{and}$
(49) $\displaystyle{{\bf{b}}_{{\nu}}}$
$\displaystyle=\Re\left\\{\operatorname{\mathrm{diag}}\left({{\bm{\mu}}_{\overline{\bf{h}}}^{\left(t\right)}}\right){\bm{\Phi}}_{{\mathrm{T}},\nu}^{\rm{T}}{{\bf{y}}_{\mathrm{T}}^{*}}-\operatorname{\mathrm{diag}}\left({\left({{\bm{\mu}}_{\overline{\bf{h}}}^{\left(t\right)}\left({\bm{\mu}}_{\overline{\bf{h}}}^{\left(t\right)}\right)^{\rm{H}}+{\bm{\Sigma}}_{\overline{\bf{h}}}^{\left(t\right)}}\right){{\left({{\bm{\Phi}}_{{{\mathrm{T}}}}+{{\bm{\Phi}}_{{\mathrm{T}},\tau}}{\operatorname{\mathrm{diag}}}\left({{{\bm{\iota}}_{\tau}^{\left(t\right)}}}\right)}\right)}^{\rm{H}}}{{\bm{\Phi}}_{{\mathrm{T}},\nu}}}\right)\right\\},$

respectively.

Since ${{\bm{\kappa}}_{\nu}}$ is jointly sparse with $\overline{\bf{h}}$
sharing the same support, to reduce the complexity for optimizing
${{\bm{\kappa}}_{\nu}}$, we can truncate ${{\bf{A}}_{\nu}}$,
${{\bm{\kappa}}_{\nu}}$, and ${{\bf{b}}_{\bm{\nu}}}$ according to
${{\bm{\alpha}}^{\left(t\right)}}$. However, the number of paths $P$, i.e.,
the number of non-zero entries in $\overline{\bf{h}}$ is unknown. As typically
a sparse vector with a length of $M_{\tau}N_{\nu}$ and at most
$\frac{M_{\mathrm{T}}N_{\mathrm{T}}}{\ln\left(M_{\tau}N_{\nu}\right)}$ non-
zero entries can be recovered based on ${M_{\mathrm{T}}N_{\mathrm{T}}}$
measurements by compressed sensing technique[17], we set
$\hat{P}=\left\lfloor{\frac{M_{\mathrm{T}}N_{\mathrm{T}}}{\ln\left(M_{\tau}N_{\nu}\right)}}\right\rfloor$
and truncate ${{\bf{A}}_{\nu}}$, ${{\bm{\kappa}}_{\nu}}$, and
${{\bf{b}}_{\bm{\nu}}}$ accordingly. In particular, if
${\alpha_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}^{\left(t\right)}}$ is one
of the largest $\hat{P}$ entries of ${{\bm{\alpha}}^{\left(t\right)}}$, the
$\left[{k^{\prime\prime}M_{\tau}+l^{\prime\prime}},{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}\right]$-th
entry of ${{\bf{A}}_{\nu}}$ and the
$\left({k^{\prime\prime}M_{\tau}+l^{\prime\prime}}\right)$-th entry of
${{\bf{b}}_{{\nu}}}$ will be selected, resulting in
${{\bf{A}}^{\mathrm{T}}_{\nu}}\in\mathbb{R}^{\hat{P}\times\hat{P}}$ and
${{\bf{b}}^{\mathrm{T}}_{{\nu}}}\in\mathbb{R}^{\hat{P}\times 1}$,
respectively. The index $k^{\prime\prime}M_{\tau}+l^{\prime\prime}$ will be
included in the truncation index set $\mathcal{S}^{\left(t\right)}_{\nu}$ with
$\left|\mathcal{S}^{\left(t\right)}_{\nu}\right|=\hat{P}$. Also, the
$\left({k^{\prime\prime}M_{\tau}+l^{\prime\prime}}\right)$-th entry of
${{\bm{\kappa}}_{\nu}}$ will be selected, resulting a truncated version of
unknown hyper-parameter, i.e.,
${{\bm{\kappa}}^{\mathrm{T}}_{\nu}}\in\mathbb{R}^{\hat{P}\times 1}$. When
${{\bf{A}}^{\mathrm{T}}_{\nu}}$ is invertible,
${{\bm{\kappa}}^{\mathrm{T}}_{\nu}}$ can be updated by

$\left({{\bm{\kappa}}^{\mathrm{T}}_{\nu}}\right)^{\left(t+1\right)}=\left({{\bf{A}}^{\mathrm{T}}_{\nu}}\right)^{-1}{{\bf{b}}^{\mathrm{T}}_{{\nu}}}.\vspace{-3mm}$
(50)

Otherwise, we perform an element-wise update for
${{\bm{\kappa}}^{\mathrm{T}}_{\nu}}$, i.e.,

$\\{\left({{\bm{\kappa}}^{\mathrm{T}}_{\nu}}\right)^{\left(t+1\right)}\\}_{n}=\frac{\\{{\bf{b}}^{\mathrm{T}}_{{\nu}}\\}_{n}-\\{\\{{{\bf{A}}^{\mathrm{T}}_{\nu}}\\}_{n}\\}_{-n}^{\mathrm{T}}\\{{{\bm{\kappa}}^{\mathrm{T}}_{\nu}}\\}_{-n}}{\\{{{\bf{A}}^{\mathrm{T}}_{\nu}}\\}_{n,n}},\forall
n\in\\{1,\ldots,\hat{P}\\}.\vspace{-2mm}$ (51)

Note that
${\kappa_{{\nu_{{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}}}}}\in\left[-\frac{1}{2}r_{\nu},\frac{1}{2}r_{\nu}\right]$
and thus we have

${\\{\overline{\left({{\bm{\kappa}}^{\mathrm{T}}_{\nu}}\right)^{\left(t+1\right)}}\\}_{n}}=\left\\{{\begin{array}[]{*{20}{c}}{{{\\{\left({{\bm{\kappa}}^{\mathrm{T}}_{\nu}}\right)^{\left(t+1\right)}\\}}_{n}}},&{{\rm{if}}\;{{\\{\left({{\bm{\kappa}}^{\mathrm{T}}_{\nu}}\right)^{\left(t+1\right)}\\}}_{n}}\in\left[{-\frac{1}{2}r_{\nu},\frac{1}{2}r_{\nu}}\right]},\\\\[-2.84526pt]
{-\frac{1}{2}r_{\nu}},&{{\rm{if}}\;{{\\{\left({{\bm{\kappa}}^{\mathrm{T}}_{\nu}}\right)^{\left(t+1\right)}\\}}_{n}}<-\frac{1}{2}r_{\nu}},\\\\[-2.84526pt]
{\frac{1}{2}r_{\nu}},&{{\rm{otherwise}}}.\end{array}}\right.\vspace{-2mm}$
(52)

Now, given the estimated
$\overline{\left({{\bm{\kappa}}^{\mathrm{T}}_{\nu}}\right)^{\left(t+1\right)}}$,
we can recover
$\overline{\left({{\bm{\kappa}}_{\nu}}\right)^{\left(t+1\right)}}$ based on
the truncation index set $\mathcal{S}^{\left(t\right)}_{\nu}$. Note that due
to the estimation error of ${{\bm{\alpha}}^{\left(t\right)}}$, the truncation
above may introduce performance loss in estimating the off-grid Doppler shift.
However, when the iterations is close to convergence, the estimation of
${{\bm{\alpha}}^{\left(t\right)}}$ is accurate and thus the performance loss
due to the proposed truncation is marginal.

Similarly, the off-grid components of the delay shift, i.e.,
${{\bm{\iota}}_{\tau}}$, can be updated by

$\mathop{\operatorname{minimize}}\limits_{{\bm{\iota}}_{\tau}}{E_{\overline{\bf{h}}\left|{{\bf{y}}_{\mathrm{T}};{{\bm{\alpha}}^{\left(t\right)}},{\bm{\kappa}}_{\nu}^{\left(t\right)},{\bm{\iota}}_{\tau}^{\left(t\right)},\beta_{0}^{\left(t\right)}}\right.}}\left\\{{{{\left\|{{\bf{y}}_{\mathrm{T}}-\overline{\bm{\Phi}}_{\mathrm{T}}\left({{{\bm{\kappa}}_{\nu}^{\left(t\right)}},{\bm{\iota}}_{\tau}}\right)\overline{\bf{h}}}\right\|}^{2}}}\right\\}\Leftrightarrow\mathop{\operatorname{minimize}}\limits_{{\bm{\iota}}_{\tau}}\left({{\bm{\iota}}_{\tau}^{\rm{T}}{{\bf{A}}_{\tau}}{\bm{\iota}}_{\tau}-2{\bf{b}}_{\tau}^{\rm{T}}{\bm{\iota}}_{\tau}}\right),\vspace{-3mm}$
(53)

where ${{\bf{A}}_{\tau}}\in\mathbb{R}^{M_{\tau}N_{\nu}\times M_{\tau}N_{\nu}}$
and ${{\bf{b}}_{{\tau}}}\in\mathbb{R}^{M_{\tau}N_{\nu}\times 1}$ are given by

$\displaystyle{{\bf{A}}_{\tau}}$
$\displaystyle=\Re\left\\{{\bm{\Phi}}_{{\mathrm{T}},\tau}^{\rm{H}}{{\bm{\Phi}}_{{\mathrm{T}},\tau}}\odot\left({\left({\bm{\mu}}_{\overline{\bf{h}}}^{\left(t\right)}\right)^{*}\left({\bm{\mu}}_{\overline{\bf{h}}}^{\left(t\right)}\right)^{\rm{T}}+\left({\bm{\Sigma}}_{\overline{\bf{h}}}^{\left(t\right)}\right)^{\rm{T}}}\right)\right\\}\;\text{and}$
(54) $\displaystyle{{\bf{b}}_{{\tau}}}$
$\displaystyle=\Re\left\\{\operatorname{\mathrm{diag}}\left({{\bm{\mu}}_{\overline{\bf{h}}}^{\left(t\right)}}\right){\bm{\Phi}}_{{\mathrm{T}},\tau}^{\rm{T}}{{\bf{y}}_{\mathrm{T}}^{*}}-\operatorname{\mathrm{diag}}\left({\left({{\bm{\mu}}_{\overline{\bf{h}}}^{\left(t\right)}\left({\bm{\mu}}_{\overline{\bf{h}}}^{\left(t\right)}\right)^{\rm{H}}+{\bm{\Sigma}}_{\overline{\bf{h}}}^{\left(t\right)}}\right){{\left({{\bm{\Phi}}_{\mathrm{T}}+{{\bm{\Phi}}_{\mathrm{T},\nu}}\operatorname{\mathrm{diag}}\left({{{\bm{\kappa}}_{\nu}^{\left(t\right)}}}\right)}\right)}^{\rm{H}}}{{\bm{\Phi}}_{{\mathrm{T}},\tau}}}\right)\right\\},$

respectively. Similarly, we can truncate ${{\bf{A}}_{\tau}}$,
${{\mathbf{b}}_{{\tau}}}$, and ${{\bm{\iota}}_{\tau}}$ according to
${{\bm{\alpha}}^{\left(t\right)}}$ as
${{\bf{A}}^{\mathrm{T}}_{\tau}}\in\mathbb{R}^{\hat{P}\times\hat{P}}$,
${{\bf{b}}^{\mathrm{T}}_{{\tau}}}\in\mathbb{R}^{\hat{P}\times 1}$, and
${{\bm{\iota}}^{\mathrm{T}}_{\tau}}\in\mathbb{R}^{\hat{P}\times 1}$
respectively, and obtain a truncation index set
$\mathcal{S}^{\left(t\right)}_{\tau}$. When ${{\bf{A}}^{\mathrm{T}}_{\tau}}$
is invertible, ${{\bm{\iota}}^{\mathrm{T}}_{\tau}}$ can be updated by

$\left({{\bm{\iota}}^{\mathrm{T}}_{\tau}}\right)^{\left(t+1\right)}=\left({{\bf{A}}^{\mathrm{T}}_{\tau}}\right)^{-1}{{\bf{b}}^{\mathrm{T}}_{{\tau}}}.\vspace{-2mm}$
(55)

Otherwise, an element-wise update for ${{\bm{\iota}}^{\mathrm{T}}_{\tau}}$ can
be obtained by

$\\{{{\bm{\iota}}^{\mathrm{T}}_{\tau}}\\}_{n}=\frac{\\{{\bf{b}}^{\mathrm{T}}_{{\tau}}\\}_{n}-\\{\\{{\bf{A}}^{\mathrm{T}}_{\tau}\\}_{n}\\}_{-n}^{\mathrm{T}}\\{{{\bm{\kappa}}^{\mathrm{T}}_{\tau}}\\}_{-n}}{\\{{{\bf{A}}^{\mathrm{T}}_{\tau}}\\}_{n,n}},\forall
n\in\\{1,\ldots,\hat{P}\\}.\vspace{-2mm}$ (56)

Since
${\iota_{{\tau_{{k^{\prime\prime}M+l^{\prime\prime}}}}}}\in\left[-\frac{1}{2}r_{\tau},\frac{1}{2}r_{\tau}\right]$,
we have

${\\{\overline{\left({{\bm{\iota}}^{\mathrm{T}}_{\tau}}\right)^{\left(t+1\right)}}\\}_{n}}=\left\\{{\begin{array}[]{*{20}{c}}{{{\\{\left({{\bm{\iota}}^{\mathrm{T}}_{\tau}}\right)^{\left(t+1\right)}\\}}_{n}}},&{{\rm{if}}{{\\{\left({{\bm{\iota}}^{\mathrm{T}}_{\tau}}\right)^{\left(t+1\right)}\\}}_{n}}\in\left[{-\frac{1}{2}r_{\tau},\frac{1}{2}r_{\tau}}\right]},\\\\[-2.84526pt]
{-\frac{1}{2}r_{\tau}},&{{\rm{if}}{{\\{\left({{\bm{\iota}}^{\mathrm{T}}_{\tau}}\right)^{\left(t+1\right)}\\}}_{n}}<-\frac{1}{2}r_{\tau}},\\\\[-2.84526pt]
{\frac{1}{2}r_{\tau}},&{{\rm{otherwise}}}.\end{array}}\right.\vspace{-2mm}$
(57)

Now, we can recover
$\overline{\left({{\bm{\iota}}_{\tau}}\right)^{\left(t+1\right)}}$ based on
the estimated
$\overline{\left({{\bm{\iota}}^{\mathrm{T}}_{\tau}}\right)^{\left(t+1\right)}}$
and the truncation index set $\mathcal{S}^{\left(t\right)}_{\tau}$.

Algorithm 1 1D Off-grid SBL-based Channel Estimation Algorithm

1: InitializationInitialize the convergence tolerance $\epsilon$, the maximum
number of iterations $T_{\mathrm{max}}$, the iteration counter $t=1$, the root
parameters $\rho$, $c$ and $d$, and the noise variance
$\overline{\sigma^{2}}=\frac{\left\|\mathbf{y}_{\mathrm{T}}\right\|^{2}}{100M_{\mathrm{T}}N_{\mathrm{T}}}$.
Initialize the hyper-parameters
${{\bm{\alpha}}^{\left(t\right)}}=\left|\mathbf{\overline{\Phi}}_{\mathrm{T}}^{\mathrm{H}}\left({\bm{\kappa}_{{\nu}}},{\bm{\iota}_{{\tau}}}\right)\mathbf{y}_{\mathrm{T}}\right|$,
${\bm{\kappa}}_{\nu}^{\left(t\right)}={\bm{\iota}}_{\tau}^{\left(t\right)}=\mathbf{0}$,
and $\beta_{0}^{\left(t\right)}=\frac{1}{\overline{\sigma^{2}}}$.

2: repeat

3: Update the conditional posterior mean
$\bm{\mu}_{\overline{\bf{h}}}^{\left(t\right)}$ and covariance matrix
$\bm{\Sigma}_{\overline{\bf{h}}}^{\left(t\right)}$ according to (35) and (36),
respectively

4: Update the hyper-parameters
$\left({{{\bm{\alpha}}^{\left(t+1\right)}},{\bm{\kappa}}_{\nu}^{\left(t+1\right)},{\bm{\iota}}_{\tau}^{\left(t+1\right)},\beta_{0}^{\left(t+1\right)}}\right)$
according to (44), (52), (57), and (45), respectively

5: $t=t+1$

6: until $t=T_{\mathrm{max}}$ or
$\frac{\left\|{{\bm{\alpha}}^{\left(t+1\right)}}-{{\bm{\alpha}}^{\left(t\right)}}\right\|_{2}}{\left\|{{\bm{\alpha}}^{\left(t\right)}}\right\|_{2}}\leq\epsilon$

7: Return the converged solution
$\hat{\overline{\bf{h}}}=\bm{\mu}_{\overline{\bf{h}}}^{\left(t\right)}$,
$\hat{\mathbf{k}}_{{\nu}}=\overline{\mathbf{k}}_{\nu}+{\bm{\kappa}}_{\nu}^{\left(t\right)}$,
and
$\hat{\mathbf{l}}_{{\tau}}=\overline{\mathbf{l}}_{\tau}+{\bm{\iota}}_{\tau}^{\left(t\right)}$.

Now, the proposed 1D off-grid SBL-based channel estimation algorithm for OTFS
can be summarized in Algorithm 1, where the hyper-parameters are updated
alternatively according to (44), (52), (57), and (45), with concurrently
updating the conditional posterior distribution of $\overline{\bf{h}}$
according to (35) and (36). When the maximum number of iterations is reached
or the changes of hyper-parameter ${\bm{\alpha}}$ are smaller than a
predefined convergence tolerance $\epsilon$, the algorithm terminates and
outputs the estimated channel coefficient $\hat{\overline{\bf{h}}}$, the
estimated normalized Doppler shift $\hat{\mathbf{k}}_{{\nu}}$, and the
estimated normalized delay shift $\hat{\mathbf{l}}_{{\tau}}$. By the virtue of
the EM algorithm[16], the proposed algorithm is guaranteed to converge to the
stationary point of
$p\left({\bm{\alpha}},{{{\bm{\kappa}}_{\nu}},{{\bm{\iota}}_{\tau}},{\beta_{0}}}\left|{\bf{y}}_{\mathrm{T}}\right.\right)$
since it is non-decreasing over iterations. The computational complexity of
the proposed algorithm is dominated by computing the conditional posterior
covariance matrix in (36), which has an order of
$\mathcal{O}\left(M^{2}_{\tau}N^{2}_{\nu}M_{\mathrm{T}}N_{\mathrm{T}}\right)$
[10]. We can observe that the 1D off-grid SBL-based channel estimation
algorithm suffers from high computational complexity with a large virtual
sampling grid size, i.e., $M_{\tau}\to\infty$ and $N_{\nu}\to\infty$. This
motivates us to propose the 2D off-grid SBL-based channel estimation as
detailed in the next section, which enjoys a lower computational complexity.

###### Remark 1

The above proposed framework treats the unknown data symbols in the
measurement matrix in (23) as the impact of noise. To mitigate the impacts of
the unknown data symbols, we can further proposed an iterative data-aided
channel estimation scheme. The key idea of data-aided channel estimation [18]
is that the estimated data symbols can be viewed as pilots, which can in turn
further improve the channel estimation performance. In particular, we first
perform the proposed off-grid 1D channel estimation, i.e., Algorithm 1, based
on (23) assuming only pilot symbol available in its measurement matrix. Then,
we reconstruct the effective DD domain channel through (3) and perform data
detection to estimate the data symbols $\hat{x}\left[k,l\right]$ based on (1).
Substituting the estimated symbols $\hat{x}\left[k,l\right]$ into the
measurement matrix in (23), we can perform our proposed Algorithm 1 again to
update the channel estimation. Such an iteration procedure proceeds until the
estimated channel does not change. The performance of the proposed data-aided
channel estimation will be evaluated in the Section VI.

##  V 2D Off-grid SBL-based Channel Estimation

In this section, to reduce the computational complexity of the SBL-based
channel estimation, we propose a 2D off-grid SSR model and a 2D off-grid SBL-
based channel estimation algorithm. Note that all the previous works for
compressed channel sensing of OTFS [7, 9, 8] are based on 1D SSR model, which
are fundamentally different from our proposed 2D model.

###  V-A 2D Off-grid SSR Model

When the pilot sequence is separable in Doppler and delay domains, i.e.,
$\mathbf{X}=\mathbf{x}_{\nu}\otimes\mathbf{x}_{\tau}^{\mathrm{H}}$, where
$\mathbf{X}\in\mathbb{C}^{\left|\mathcal{K}_{\mathrm{p}}\right|\times\left|\mathcal{L}_{\mathrm{p}}\right|}$
collects the pilot symbols $x\left[{k^{\prime},l^{\prime}}\right]$,
$\mathbf{x}_{\nu}\in\mathbb{C}^{\left|\mathcal{K}_{\mathrm{p}}\right|\times
1}$ collects the pilot symbols $x_{\nu}\left[{k^{\prime}}\right]$ in the
Doppler domain, and
$\mathbf{x}_{\tau}\in\mathbb{C}^{\left|\mathcal{L}_{\mathrm{p}}\right|\times
1}$ collects the pilot symbols $x_{\tau}\left[{l^{\prime}}\right]$ in the
delay domain, $\forall k^{\prime}\in\mathcal{K}_{\mathrm{p}}$ and $\forall
l^{\prime}\in\mathcal{L}_{\mathrm{p}}$, the system model in (13) can be
rewritten as

$\displaystyle y\left[{k,l}\right]$
$\displaystyle=\sum\limits_{k^{\prime}\in\mathcal{K}_{\mathrm{p}}}{\sum\limits_{l^{\prime}\in\mathcal{L}_{\mathrm{p}}}x_{\nu}\left[k^{\prime}\right]x_{\tau}\left[l^{\prime}\right]}{h_{w}}\left[{k-k^{\prime},l-l^{\prime}}\right]+z\left[k,l\right]$
$\displaystyle=\sum_{i=1}^{P}{\tilde{h}_{i}}\sum\limits_{k^{\prime}\in\mathcal{K}_{\mathrm{p}}}x_{\nu}\left[k^{\prime}\right]w_{\nu}(k-k^{\prime}-k_{\nu_{i}})\sum\limits_{l^{\prime}\in\mathcal{L}_{\mathrm{p}}}x_{\tau}\left[l^{\prime}\right]w_{\tau}(l-l^{\prime}-l_{\tau_{i}})+z\left[k,l\right]$
$\displaystyle\mathop{\rm{=}}\limits^{\left(a\right)}\sum_{i=1}^{P}{\tilde{h}_{i}}x_{\nu}\left[k_{p}\right]w_{\nu}(k-k_{p}-k_{\nu_{i}})x_{\tau}\left[l_{p}\right]w_{\tau}(l-l_{p}-l_{\tau_{i}})+z\left[k,l\right].$
(58)

Equality $(a)$ is obtained with adopting a single pilot impulse, i.e.,
$\left|\mathcal{K}_{\mathrm{p}}\right|=\left|\mathcal{L}_{\mathrm{p}}\right|=1$.
Note that for the case of a single pilot impulse, the pilot sequence is
separable as
$x\left[{k_{\mathrm{p}},l_{\mathrm{p}}}\right]=x_{\nu}\left[k_{p}\right]x_{\tau}\left[l_{p}\right]$
always holds.

From (V-A), we can observe that it is possible to decouple the estimation of
the normalized Doppler shift $k_{\nu_{i}}$ and the normalized delay shift
$l_{\tau_{i}}$, which is different from the 1D model in Section III.
Benefiting from this decoupling, it is possible to reduce the computational
complexity based on the proposed 2D SSR model as follows. To this end, we
define the virtual sampling grid for Doppler and delay domains individually,
which is different from that in (18). In particular, the Doppler virtual
sampling grid444For notational consistency, in this section, we reuse some
notations of Section III. is defined in the range of
$\left(-\nu_{\mathrm{max}},\nu_{\mathrm{max}}\right)$ with a virtual Doppler
resolution $r_{\nu}=\frac{2\nu_{\mathrm{max}}}{N_{\nu}}$, resulting in
$\overline{\mathbf{k}}_{\nu}=\left[\overline{k}_{\nu_{0}},\ldots,\overline{k}_{\nu_{N_{\nu}-1}}\right]\in\mathbb{R}^{N_{\nu}\times
1}$ and
$\overline{k}_{\nu_{k^{\prime\prime}}}=k^{\prime\prime}r_{\nu}-\nu_{\mathrm{max}}$.
On the other hand, the delay virtual sampling grid is defined in the range of
$\left(0,\tau_{\mathrm{max}}\right)$ with a virtual delay resolution
$r_{\tau}=\frac{\tau_{\mathrm{max}}}{M_{\tau}}$, resulting in
$\overline{\mathbf{l}}_{\tau}=\left[\overline{l}_{\tau_{0}},\ldots,\overline{l}_{\tau_{M_{\tau}-1}}\right]\in\mathbb{R}^{M_{\tau}\times
1}$ and $\overline{l}_{\tau_{l^{\prime\prime}}}=l^{\prime\prime}r_{\tau}$.
Now, similar to (23), we can recast the system model in (V-A) as a 2D SSR
problem based on the defined virtual sampling grids and the first-order
approximation:

$\mathbf{Y}_{\mathrm{T}}=\overline{\mathbf{\Phi}}_{\mathrm{L}}\left(\bm{\kappa}_{\nu}\right)\overline{\mathbf{H}}\overline{\mathbf{\Phi}}_{\mathrm{R}}^{\mathrm{T}}\left(\bm{\iota}_{\tau}\right)+\overline{\mathbf{Z}}_{\mathrm{T}},\vspace{-2mm}$
(59)

where $\mathbf{Y}_{\mathrm{T}}\in\mathbb{C}^{N_{\mathrm{T}}\times
M_{\mathrm{T}}}$ collects the received signal $y\left[k,l\right]$ in the range
of $k_{\mathrm{p}}-k_{\mathrm{max}}\leq k\leq k_{\mathrm{p}}+k_{\mathrm{max}}$
and $l_{\mathrm{p}}\leq l\leq l_{\mathrm{p}}+k_{\mathrm{max}}$ thanks to the
compactness of the DD domain channel. Also, matrix
$\overline{\mathbf{Z}}_{\mathrm{T}}\in\mathbb{C}^{N_{\mathrm{T}}\times
M_{\mathrm{T}}}$ collects the new noise
$\overline{z}\left[k,l\right]\sim\mathcal{CN}\left(\mathbf{0},\overline{\sigma^{2}}\right)$,
where $\overline{\sigma^{2}}$ is the unknown noise power taking into account
the measurement noise, approximation error, as well as the interference from
unknown data symbols. The unknown channel matrix
$\overline{\mathbf{H}}\in\mathbb{C}^{N_{\nu}\times M_{\tau}}$ collects the
channel coefficients $\tilde{h}_{i}$, $\forall i$, i.e.,
$\\{\overline{\mathbf{H}}\\}_{\overline{k}^{i}_{\nu_{k^{\prime\prime}}},\overline{l}^{i}_{\tau_{l^{\prime\prime}}}}=\tilde{h}_{i}$,
where $\overline{k}^{i}_{\nu_{k^{\prime\prime}}}$ and
$\overline{l}^{i}_{\tau_{l^{\prime\prime}}}$ denote the unknown nearest points
of $k_{\nu_{i}}$ and $l_{\tau_{i}}$, respectively. In other words, there are
at most $P$ non-zero entries in the sparse unknown channel matrix
$\overline{\mathbf{H}}$.

In (59), the left measurement matrix
$\overline{\mathbf{\Phi}}_{\mathrm{L}}\left(\bm{\kappa}_{\nu}\right)\in\mathbb{C}^{N_{\mathrm{T}}\times
N_{\nu}}$ and the right measurement matrix
$\overline{\mathbf{\Phi}}_{\mathrm{R}}\left(\bm{\iota}_{\tau}\right)\in\mathbb{C}^{M_{\mathrm{T}}\times
M_{\tau}}$ are defined as

$\overline{\mathbf{\Phi}}_{\mathrm{L}}\left(\bm{\kappa}_{\nu}\right)=\mathbf{\Phi}_{\mathrm{L}}+\mathbf{\Phi}_{\mathrm{L},\nu}\operatorname{\mathrm{diag}}\left(\bm{\kappa}_{\nu}\right)\;\text{and}\;\overline{\mathbf{\Phi}}_{\mathrm{R}}\left(\bm{\iota}_{\tau}\right)=\mathbf{\Phi}_{\mathrm{R}}+\mathbf{\Phi}_{\mathrm{R},\tau}\operatorname{\mathrm{diag}}\left(\bm{\iota}_{\tau}\right),\vspace{-2mm}$
(60)

respectively, where
$\bm{\kappa}_{\nu}=\left[\kappa_{\nu_{0}},\ldots,\kappa_{\nu_{N_{\nu}-1}}\right]\in\mathbb{R}^{N_{\nu}\times
1}$ and
$\bm{\iota}_{\tau}=\left[\iota_{\tau_{0}},\ldots,\iota_{\tau_{M_{\tau}-1}}\right]\in\mathbb{R}^{M_{\tau}\times
1}$. And the matrices
$\mathbf{\Phi}_{\mathrm{L}},\mathbf{\Phi}_{\mathrm{L},\nu}\in\mathbb{C}^{N_{\mathrm{T}}\times
N_{\nu}}$ and
$\mathbf{\Phi}_{\mathrm{R}},\mathbf{\Phi}_{\mathrm{R},\tau}\in\mathbb{C}^{M_{\mathrm{T}}\times
M_{\tau}}$ are defined as
$\\{\mathbf{\Phi}_{\mathrm{L}}\\}_{k,k^{\prime\prime}}=x_{\nu}\left[k_{\mathrm{p}}\right]w_{\nu}(k-k_{\mathrm{p}}-\overline{k}_{\nu_{k^{\prime\prime}}})$,
$\\{\mathbf{\Phi}_{\mathrm{L,\nu}}\\}_{k,k^{\prime\prime}}=x_{\nu}\left[k_{\mathrm{p}}\right]w^{\prime}_{\nu}(k-k_{\mathrm{p}}-\overline{k}_{\nu_{k^{\prime\prime}}})$,
$\\{\mathbf{\Phi}_{\mathrm{R}}\\}_{l,l^{\prime\prime}}=x_{\tau}\left[l_{\mathrm{p}}\right]w_{\tau}(l-l_{\mathrm{p}}-\overline{l}_{\tau_{l^{\prime\prime}}})$,
and
$\\{\mathbf{\Phi}_{\mathrm{R,\tau}}\\}_{l,l^{\prime\prime}}=x_{\tau}\left[l_{\mathrm{p}}\right]w^{\prime}_{\tau}(l-l_{\mathrm{p}}-\overline{l}_{\tau_{l^{\prime\prime}}})$,
respectively, $\forall k^{\prime\prime}\in\left\\{0,\ldots,N_{\nu}-1\right\\}$
and $\forall l^{\prime\prime}\in\left\\{0,\ldots,M_{\tau}-1\right\\}$.

###  V-B 2D Off-grid SBL-based Channel Estimation Algorithm

Based on the built 2D off-grid SSR model in (59), we propose a two-step SBL-
based channel estimation algorithm in this section, where the first step
estimates the composite matrix
$\mathbf{D}=\overline{\mathbf{H}}\overline{\mathbf{\Phi}}_{\mathrm{R}}^{\mathrm{T}}\left(\bm{\iota}_{\tau}\right)$
and the second step aims to recovery $\overline{\mathbf{H}}$ based on the
estimated $\mathbf{D}$.

####  V-B1 First Step

Defining a composite matrix
$\mathbf{D}=\overline{\mathbf{H}}\overline{\mathbf{\Phi}}_{\mathrm{R}}^{\mathrm{T}}\left(\bm{\iota}_{\tau}\right)=\left[\mathbf{d}_{\mathrm{c},0},\ldots,\mathbf{d}_{\mathrm{c},M_{\mathrm{T}}-1}\right]\in\mathbb{C}^{N_{\nu}\times
M_{\mathrm{T}}}$, the system model in (59) can be rewritten as

$\mathbf{Y}_{\mathrm{T}}=\overline{\mathbf{\Phi}}_{\mathrm{L}}\left(\bm{\kappa}_{\nu}\right)\mathbf{D}+\overline{\mathbf{Z}}_{\mathrm{T}},\vspace{-2mm}$
(61)

where $\mathbf{D}$ is a row-sparse channel matrix sharing the same row-
sparsity with $\overline{\mathbf{H}}$, i.e., the sparsity in the Doppler
domain, $\mathbf{d}_{\mathrm{c},l}\in\mathbb{C}^{N_{\nu}\times 1}$ denotes the
$l$-th column of $\mathbf{D}$, and all columns $\mathbf{d}_{\mathrm{c},l}$ are
jointly sparse with the same support, $\forall
l\in\\{0,\ldots,M_{\mathrm{T}}-1\\}$. To recovery $\mathbf{D}$, we can develop
an SBL algorithm for a multiple snapshots case assuming a row sparsity in
$\mathbf{D}$. In particular, the unknown matrix $\mathbf{D}$ and the hyper-
parameters
$\left({\bm{\alpha}}_{\nu},{{\bm{\kappa}}_{\nu}},{\beta_{0}}\right)$ are
updated iteratively, where
${\bm{\alpha}}_{\nu}=\left[{\alpha_{\nu_{0}},\ldots,{\alpha_{\nu_{k^{\prime\prime}}}},\ldots{\alpha_{\nu_{N_{\nu}-1}}}}\right]\in\mathbb{R}^{N_{\nu}\times
1}$ with ${\alpha_{\nu_{k^{\prime\prime}}}}>0$ models the row-sparsity of
$\mathbf{D}$.

Assuming independent columns of $\mathbf{D}$, the prior distribution of the
unknown matrix $\mathbf{D}$ assumes to follow

$p\left({{\bf{D}}\left|{\bm{\alpha}}_{\nu}\right.}\right)=\prod_{l=1}^{M_{\mathrm{T}}}\mathcal{CN}\left({{\mathbf{d}_{\mathrm{c},l}}\left|{{\bf{0}},{\bf{\Lambda}}}\right.}\right),\vspace{-2mm}$
(62)

where
${\bf{\Lambda}}=\operatorname{\mathrm{diag}}\left({\bm{\alpha}}_{\nu}\right)$
and
$p\left({{\bm{\alpha}}_{\nu}\left|\rho\right.}\right)=\prod\limits_{k^{\prime\prime}=0}^{N_{\nu}-1}\mathop{\Gamma}\left({{\alpha_{\nu_{k^{\prime\prime}}}}\left|{1,\rho}\right.}\right)$,
with $\rho>0$. And the unknown off-grid components of the Doppler shift in the
measurement matrix of (61) follows
${\kappa_{{\nu_{k^{\prime\prime}}}}}\sim\mathcal{U}\left[{-\frac{1}{2}r_{\nu},\frac{1}{2}r_{\nu}}\right]$,
$\forall k^{\prime\prime}\in\\{0,\ldots,N_{\nu}-1\\}$. The new noise matrix
follows

$p\left({\overline{\mathbf{Z}}_{\mathrm{T}}\left|{{\beta_{0}}}\right.}\right)=\prod_{l=0}^{M_{\mathrm{T}}-1}\mathcal{CN}\left({{\overline{\mathbf{z}}_{l}}\left|{{\bf{0}},\beta_{0}^{-1}{{\bf{I}}_{N_{\mathrm{T}}}}}\right.}\right),\vspace{-2mm}$
(63)

where
$\overline{\mathbf{Z}}_{\mathrm{T}}=\left[{\overline{\mathbf{z}}_{0}},\ldots,{\overline{\mathbf{z}}_{M_{\mathrm{T}}-1}}\right]$,
$\overline{\mathbf{z}}_{l}\in\mathbb{C}^{N_{\mathrm{T}}\times 1}$, and
$p\left({{\beta_{0}}\left|{c,d}\right.}\right)=\Gamma\left({{\beta_{0}}\left|{c,d}\right.}\right)$,
with $c,d>0$. Now, the likelihood function for given
$\left({{{\bm{\kappa}}_{\nu}},{\beta_{0}}}\right)$ for (61) is given by

$p\left({{\bf{Y}}_{\mathrm{T}}\left|{{\bf{D}};{{\bm{\kappa}}_{\nu}},{\beta_{0}}}\right.}\right)=\prod_{l=0}^{M_{\mathrm{T}}-1}\mathcal{CN}\left({{\mathbf{y}_{l}}\left|{\overline{\mathbf{\Phi}}_{\mathrm{L}}\left(\bm{\kappa}_{\nu}\right)\mathbf{d}_{\mathrm{c},l},\beta_{0}^{-1}{{\bf{I}}_{N_{\mathrm{T}}}}}\right.}\right),\vspace{-2mm}$
(64)

and the joint distribution can be decomposed by

$p\left({{\bf{Y}}_{\mathrm{T}},{\bf{D}},{\bm{\alpha}}_{\nu},{{\bm{\kappa}}_{\nu}},{\beta_{0}}}\right)=p\left({{\bf{Y}}_{\mathrm{T}}\left|{{\bf{D}},{{\bm{\kappa}}_{\nu}},{\beta_{0}}}\right.}\right)p\left({{\bf{D}}\left|{\bm{\alpha}}_{\nu}\right.}\right)p\left({\bm{\alpha}}_{\nu}\right)p\left({{{\bm{\kappa}}_{\nu}}}\right)p\left({{\beta_{0}}}\right),\vspace{-2mm}$
(65)

where ${\mathbf{y}_{l}}\in\mathbb{C}^{N_{\mathrm{T}}\times 1}$ is the $l$-th
column of ${\bf{Y}}_{\mathrm{T}}$.

Now, given the hyper-parameters
$\left({\bm{\alpha}}_{\nu}^{\left(t\right)},{{\bm{\kappa}}^{\left(t\right)}_{\nu}},{\beta^{\left(t\right)}_{0}}\right)$
in the $t$-th iteration, the conditional posterior distribution of ${\bf{D}}$
can be obtained by

$p\left({\bf{D}}\left|{\bf{Y}}_{\mathrm{T}};{{\bm{\kappa}}^{\left(t\right)}_{\nu}},{\beta^{\left(t\right)}_{0}}\right.\right)=\prod_{l=0}^{M_{\mathrm{T}}-1}\mathcal{CN}\left(\mathbf{d}_{\mathrm{c},l}\left|\bm{\mu}^{\left(t\right)}_{\mathbf{d}_{\mathrm{c},l}},\bm{\Sigma}^{\left(t\right)}_{\mathbf{d}_{\mathrm{c}}}\right.\right),\vspace{-2mm}$
(66)

where the conditional posterior mean for $\mathbf{d}_{\mathrm{c},l}$,
$\forall\in\\{0,\ldots,M_{\mathrm{T}}-1\\}$, and their common covariance
matrix are given by

$\displaystyle\bm{\mu}^{\left(t\right)}_{\mathbf{d}_{\mathrm{c},l}}$
$\displaystyle={\beta^{\left(t\right)}_{0}}\bm{\Sigma}^{\left(t\right)}_{\mathbf{d}_{\mathrm{c}}}\overline{\mathbf{\Phi}}_{\mathrm{L}}\left(\bm{\kappa}^{\left(t\right)}_{\nu}\right){\mathbf{y}_{l}}\;\text{and}$
(67) $\displaystyle\bm{\Sigma}^{\left(t\right)}_{\mathbf{d}_{\mathrm{c}}}$
$\displaystyle={\left({{\beta^{\left(t\right)}_{0}}{{\overline{\mathbf{\Phi}}_{\mathrm{L}}^{\rm{H}}\left(\bm{\kappa}^{\left(t\right)}_{\nu}\right)}}\overline{\mathbf{\Phi}}_{\mathrm{L}}\left(\bm{\kappa}^{\left(t\right)}_{\nu}\right)+{\left({\bf{\Lambda}}^{\left(t\right)}\right)^{-1}}}\right)^{-1}}$
$\displaystyle={\bm{\Lambda}}^{\left(t\right)}-{\bm{\Lambda}}^{\left(t\right)}\overline{\mathbf{\Phi}}_{\mathrm{L}}^{\mathrm{H}}\left(\bm{\kappa}^{\left(t\right)}_{\nu}\right){\left({\left(\beta^{\left(t\right)}_{0}\right)^{-1}{{\bf{I}}_{{N_{\mathrm{T}}}}}+\overline{\mathbf{\Phi}}_{\mathrm{L}}\left(\bm{\kappa}^{\left(t\right)}_{\nu}\right){\bm{\Lambda}}\overline{\mathbf{\Phi}}_{\mathrm{L}}^{\mathrm{H}}\left(\bm{\kappa}^{\left(t\right)}_{\nu}\right)}\right)^{-1}}\overline{\mathbf{\Phi}}_{\mathrm{L}}\left(\bm{\kappa}^{\left(t\right)}_{\nu}\right){\bm{\Lambda}}^{\left(t\right)},$
(68)

respectively.

Similar to Proposition 1, following the EM principle [16], the hyper-
parameters $\left({\bm{\alpha}},{\beta_{0}}\right)$ can be updated by the
following rules

$\displaystyle\alpha_{k^{\prime\prime}}^{\left({t+1}\right)}$
$\displaystyle=\frac{{\sqrt{1+4\frac{\rho}{M_{\mathrm{T}}}\frac{1}{M_{\mathrm{T}}}\sum_{l=0}^{M_{\mathrm{T}}-1}\\{{{{\left|\\{\bm{\mu}^{\left(t\right)}_{\mathbf{d}_{\mathrm{c},l}}\\}_{k^{\prime\prime}}\right|}}^{2}}+\\{\bm{\Sigma}^{\left(t\right)}_{\mathbf{d}_{\mathrm{c}}}\\}_{k^{\prime\prime},k^{\prime\prime}}\\}}-1}}{{2\frac{\rho}{M_{\mathrm{T}}}}}\;\text{and}$
(69) $\displaystyle\beta_{0}^{\left({t+1}\right)}$
$\displaystyle=\frac{{c-1+M_{\mathrm{T}}N_{\mathrm{T}}}}{{d+\sum_{l=0}^{M_{\mathrm{T}}-1}{A^{l}_{{\beta_{0}}}}}},$
(70)

respectively, and
${A^{l}_{{\beta_{0}}}}={\left\|{{\bf{y}}_{l}-\overline{\mathbf{\Phi}}_{\mathrm{L}}^{\mathrm{H}}\left(\bm{\kappa}^{\left(t\right)}_{\nu}\right)\bm{\mu}^{\left(t\right)}_{\mathbf{d}_{\mathrm{c},l}}}\right\|^{2}}{\rm{+}}\left(\beta^{\left(t\right)}_{0}\right)^{-1}\mathop{\sum}\limits_{k^{\prime\prime}=0}^{N_{\nu}-1}\left({1-\left(\alpha_{\nu_{k^{\prime\prime}}}\right)^{-1}\\{\bm{\Sigma}^{\left(t\right)}_{\mathbf{d}_{\mathrm{c}}}\\}_{k^{\prime\prime},k^{\prime\prime}}}\right)$.

Also, the hyper-parameter ${{\bm{\kappa}}_{\nu}}$ can be updated by

${{\bm{\kappa}}_{\nu}}=\left({{\bf{A}}_{\nu}}\right)^{-1}{{\bf{b}}_{{\nu}}},\vspace{-2mm}$
(71)

when
${{\bf{A}}_{\nu}}=\frac{1}{{M_{\mathrm{T}}}}\sum_{l=1}^{{{M_{\mathrm{T}}}}}\Re\left\\{\mathbf{\Phi}_{\mathrm{L},\nu}^{\rm{H}}\mathbf{\Phi}_{\mathrm{L},\nu}\odot\left({\left(\bm{\mu}^{\left(t\right)}_{\mathbf{d}_{\mathrm{c},l}}\right)^{*}\left(\bm{\mu}^{\left(t\right)}_{\mathbf{d}_{\mathrm{c},l}}\right)^{\rm{T}}+\left(\bm{\Sigma}^{\left(t\right)}_{\mathbf{d}_{\mathrm{c}}}\right)^{\rm{T}}}\right)\right\\}\in\mathbb{R}^{N_{\nu}\times
N_{\nu}}$ is invertible and ${{\bf{b}}_{{\nu}}}$ is given by

${{\bf{b}}_{{\nu}}}=\frac{1}{{M_{\mathrm{T}}}}\sum_{l=1}^{{{M_{\mathrm{T}}}}}\Re\left\\{{\operatorname{\mathrm{diag}}\left(\bm{\mu}^{\left(t\right)}_{\mathbf{d}_{\mathrm{c},l}}\right)\mathbf{\Phi}_{\mathrm{L},\nu}^{\rm{T}}{{\mathbf{y}}^{*}_{l}}-\operatorname{\mathrm{diag}}\left({\left({\bm{\mu}^{\left(t\right)}_{\mathbf{d}_{\mathrm{c},l}}\left(\bm{\mu}^{\left(t\right)}_{\mathbf{d}_{\mathrm{c},l}}\right)^{\mathrm{H}}+\bm{\Sigma}^{\left(t\right)}_{\mathbf{d}_{\mathrm{c}}}}\right){\mathbf{\Phi}_{\mathrm{L},\nu}}}\right)}\right\\}\in\mathbb{R}^{N_{\nu}\times
1}.\vspace{-2mm}$ (72)

When ${{\bf{A}}_{\nu}}$ is not invertible, we perform an element-wise update
for ${{\bm{\kappa}}_{\nu}}$, i.e.,

$\\{{{\bm{\kappa}}^{\left(t+1\right)}_{\nu}}\\}_{n}=\frac{\\{{\bf{b}}_{{\nu}}\\}_{n}-\\{\\{{\bf{A}}_{\nu}\\}_{n}\\}_{-n}^{\mathrm{T}}\\{{{\bm{\kappa}}_{\nu}}\\}_{-n}}{\\{{{\bf{A}}_{\nu}}\\}_{n,n}}.\vspace{-2mm}$
(73)

Note that
${\kappa_{{\nu_{{k^{\prime\prime}}}}}}\in\left[-\frac{1}{2}r_{\nu},\frac{1}{2}r_{\nu}\right]$
and thus we have

${\\{\overline{{\bm{\kappa}}^{\left(t+1\right)}_{\nu}}\\}_{n}}=\left\\{{\begin{array}[]{*{20}{c}}{{{\\{{\bm{\kappa}}_{\nu}^{\left(t+1\right)}\\}}_{n}}},&{{\rm{if}}\;{{\\{{\bm{\kappa}}_{\nu}^{\left(t+1\right)}\\}}_{n}}\in\left[{-\frac{1}{2}r_{\nu},\frac{1}{2}r_{\nu}}\right]},\\\\[-2.84526pt]
{-\frac{1}{2}r_{\nu}},&{{\rm{if}}\;{{\\{{\bm{\kappa}}_{\nu}^{\left(t+1\right)}\\}}_{n}}<-\frac{1}{2}r_{\nu}},\\\\[-2.84526pt]
{\frac{1}{2}r_{\nu}},&{{\rm{otherwise}}}.\end{array}}\right.\vspace{-2mm}$
(74)

Now, the conditional posterior distribution of the unknown matrix $\mathbf{D}$
and the values of hyper-parameters
$\left({\bm{\alpha}},{{\bm{\kappa}}_{\nu}},{\beta_{0}}\right)$ can be updated
according to (67), (68), (69), (70), and (74) iteratively, as detailed in line
$1-8$ in Algorithm 2 as follows. After convergence, we can obtain
$\hat{\bf{D}}=\left[\bm{\mu}^{\left(t\right)}_{\mathbf{d}_{\mathrm{c},0}},\ldots,\bm{\mu}^{\left(t\right)}_{\mathbf{d}_{\mathrm{c},M_{\mathrm{T}}-1}}\right]$
and the Doppler shift
$\hat{\mathbf{k}}_{{\nu}}=\overline{\mathbf{k}}_{\nu}+{\bm{\kappa}}_{\nu}^{\left(t\right)}$,
where $\hat{\bf{D}}$ will be used to estimate the unknown channel matrix and
the delay shift in the second step.

####  V-B2 Second Step

The estimated composite matrix $\hat{\bf{D}}$ can be rewritten as
$\hat{\bf{D}}=\left[\hat{\mathbf{d}}_{\mathrm{r},0},\ldots,\hat{\mathbf{d}}_{\mathrm{r},N_{\nu}-1}\right]^{\mathrm{T}}$,
where
$\left(\hat{\mathbf{d}}_{\mathrm{r},k}\right)^{\mathrm{T}}\in\mathbb{C}^{1\times{M_{\mathrm{T}}}}$
denotes the $k$-th row vector of $\hat{\bf{D}}$. Besides, based on the
definition of ${\bf{D}}$ in (61), we can assume
$\hat{\bf{D}}=\overline{\mathbf{H}}\overline{\mathbf{\Phi}}_{\mathrm{R}}^{\mathrm{T}}\left(\bm{\iota}_{\tau}\right)$.
Therefore, defining
$\overline{\mathbf{H}}=\left[\overline{\mathbf{h}}_{0},\ldots,\overline{\mathbf{h}}_{N_{\nu}-1}\right]^{\mathrm{T}}$,
where
$\overline{\mathbf{h}}^{\mathrm{T}}_{k}\in\mathbb{C}^{1\times{M_{\tau}}}$
denotes the $k$-th row vector of $\overline{\mathbf{H}}$, we can obtain
${N_{\nu}}$ 1D SSR problems:

$\hat{\mathbf{d}}_{\mathrm{r},k}=\overline{\mathbf{\Phi}}_{\mathrm{R}}\left(\bm{\iota}_{\tau}\right)\overline{\mathbf{h}}_{k},\forall
k\in\\{0,\ldots,{N_{\nu}}-1\\},\vspace{-2mm}$ (75)

which can be solved in parallel using Algorithm 1. The details for recovering
$\overline{\mathbf{h}}_{k}$ are omitted here to concise our presentation.

Algorithm 2 2D Off-grid SBL-based Channel Estimation Algorithm

1: InitializationInitialize the convergence tolerance $\epsilon$, the maximum
number of iterations $T_{\mathrm{max}}$, the iteration counter $t=1$, the root
parameters $\rho$, $c$, and $d$, and the noise variance
$\overline{\sigma^{2}}=\frac{\left\|\mathbf{Y}_{\mathrm{T}}\right\|^{2}}{100M_{\mathrm{T}}N_{\mathrm{T}}}$.
Initialize the hyper-parameters
${\alpha_{\nu_{k^{\prime\prime}}}^{\left(t\right)}}=\frac{1}{M_{\mathrm{T}}}\sum_{l=0}^{M_{\mathrm{T}}-1}\left\\{\left|\\{\overline{\mathbf{\Phi}}^{\mathrm{H}}_{\mathrm{L}}\left(\bm{\kappa}_{\nu}\right)\mathbf{Y}_{\mathrm{T}}\\}_{k^{\prime\prime},l}\right|\right\\}$,
${\bm{\kappa}}_{\nu}^{\left(t\right)}=\mathbf{0}$, and
$\beta_{0}^{\left(t\right)}=\frac{1}{\overline{\sigma^{2}}}$.

2: repeat {First Step}

3: Update the conditional posterior mean
$\bm{\mu}^{\left(t\right)}_{\mathbf{d}_{\mathrm{c},l}}$ and covariance matrix
$\bm{\Sigma}^{\left(t\right)}_{\mathbf{d}_{\mathrm{c}}}$ according to (67) and
(68), respectively

4: Update the hyper-parameters
$\left({{{\bm{\alpha}}_{\nu}^{\left(t+1\right)}},{\bm{\kappa}}_{\nu}^{\left(t+1\right)},\beta_{0}^{\left(t+1\right)}}\right)$
according to (69), (74), and (70), respectively

5: $t=t+1$

6: until $t=T_{\mathrm{max}}$ or
$\frac{\left\|{{\bm{\alpha}}_{\nu}^{\left(t+1\right)}}-{{\bm{\alpha}}_{\nu}^{\left(t\right)}}\right\|_{2}}{\left\|{{\bm{\alpha}}_{\nu}^{\left(t\right)}}\right\|_{2}}\leq\epsilon$

7: Output the converged solution
$\hat{\bf{D}}=\left[\bm{\mu}^{\left(t\right)}_{\mathbf{d}_{\mathrm{c},0}},\ldots,\bm{\mu}^{\left(t\right)}_{\mathbf{d}_{\mathrm{c},M_{\mathrm{T}}-1}}\right]=\left[\hat{\mathbf{d}}_{\mathrm{r},0},\ldots,\hat{\mathbf{d}}_{\mathrm{r},N_{\nu}-1}\right]^{\mathrm{T}}$
and
${\hat{\mathbf{k}}_{{\nu}}}=\overline{\mathbf{k}}_{\nu}+{\bm{\kappa}}_{\nu}^{\left(t\right)}$.

8: Second Step: Recover $\hat{\overline{\bf{h}}}_{k}$ and
$\hat{\mathbf{l}}_{{\tau}}$ according to Algorithm 1, based on the model in
(75), $\forall k$.

Now, the 2D off-grid SBL-based channel estimation algorithm can be summarized
in Algorithm 2. The convergence of the first step can be guaranteed by the
virtue of the EM algorithm[16], where
$p\left({\bm{\alpha}},{{{\bm{\kappa}}_{\nu}},{\beta_{0}}}\left|{\bf{Y}}_{\mathrm{T}}\right.\right)$
will be non-decreasing for each iteration. The computational complexity of the
proposed algorithm is dominated by channel inversion operations in both the
first step and second step, which has an order of
$\mathcal{O}\left(N^{2}_{\nu}N_{\mathrm{T}}+M^{2}_{\tau}M_{\mathrm{T}}\right)$
[10]. We can observe that the computational complexity of Algorithm 2 is much
lower compared with Algorithm 1. In particular, benefiting from the decoupling
of Doppler and delay shift estimation of our proposed 2D method, its
computational complexity order is proportional to the summation of the sizes
of delay and Doppler grids, i.e.,
$\mathcal{O}\left(N^{2}_{\nu}N_{\mathrm{T}}+M^{2}_{\tau}M_{\mathrm{T}}\right)$,
while the counterpart of the proposed 1D off-grid method is proportional to
their product, i.e.,
$\mathcal{O}\left(N^{2}_{\nu}M^{2}_{\tau}N_{\mathrm{T}}M_{\mathrm{T}}\right)$.

##  VI Simulations

In this section, we evaluate the channel estimation performance of our
proposed 1D and 2D off-grid SBL-based channel estimation algorithms. Unless
specified otherwise, the system parameters used in the our simulations are
given as follows. The DD domain grid size is $M=N=32$, the carrier frequency
is $3$ GHz, the subcarrier spacing is $\Delta_{f}=15$ kHz, and the slot
duration is $T=\frac{1}{\Delta_{f}}$. The total number of path is set as
$P=5$, the normalized maximum Doppler shift and maximum delay shift are given
by $k_{\mathrm{max}}=3$ and $l_{\mathrm{max}}=4$, respectively. The Doppler
and delay shifts are randomly and uniformly generated within
$\left[-\nu_{\mathrm{m}ax},\nu_{\mathrm{m}ax}\right]$ and
$\left[0,\tau_{\mathrm{m}ax}\right]$, respectively. Note that the adopted
maximum Doppler and delay shift correspond to the maximum relative velocity of
$506.25$ km/h among the transceivers or scatters and the maximum propagation
distance difference of $2.5$ km between propagating paths, respectively. The
channel coefficients $h_{i}$ are generated according to the distribution
$h_{i}\sim\mathcal{CN}(0,q^{l_{\tau_{i}}})$, where $q^{l_{\tau_{i}}}$ follows
a normalized exponential power delay profile555Here, we assume that the path
loss and shadowing have been compensated following the literature studying
OTFS [12, 19]., i.e.,
$q^{l_{\tau_{i}}}=\frac{\exp(-0.1l_{\tau_{i}})}{\sum_{i}\exp(-0.1l_{\tau_{i}})}$.
A single pilot impulse is inserted at
$\left[k_{p},l_{p}\right]=\left[\frac{M}{2}+1,\frac{N}{2}+1\right]$ in the DD
domain and a normalized constellation is assumed for data symbols, i.e.,
$E\left\\{{{{\left|{x\left[{k,l}\right]}\right|}^{2}}}\right\\}=1$. The power
of the pilot symbol is 30 dB higher than that of the data symbols, i.e.,
$10\log_{10}\frac{\left|x_{p}\right|^{2}}{E\left\\{{{{\left|{x\left[{k,l}\right]}\right|}^{2}}}\right\\}}=30$
dB. The system signal-to-noise ratio (SNR) is defined as
$\mathrm{SNR}=\frac{1}{N_{0}}$ and is set in the range of $\left[0,30\right]$
dB. The virtual Doppler and delay sampling resolutions are set in the range
$r_{\nu},r_{\tau}\in\left[0.2,0.9\right]$. Following [4, 6], the DD domain
grids in the range ${k_{p}}-2{k_{\max}}\leq k\leq{k_{p}}+2{k_{\max}}$ and
${l_{p}}-{l_{\mathrm{max}}}\leq l\leq{l_{p}}+{l_{\mathrm{max}}}$, with
$k\neq{k_{p}}$ and $l\neq{l_{p}}$, remain null serving as guard space to
mitigate the interference from the data symbols to the pilot symbol. In our
proposed Algorithms 1 and 2, we set $\epsilon=10^{-3}$, $\rho=10^{-2}$, and
$c=d=10^{-4}$ [10]. For comparison, we evaluate the performance of the
traditional impulse-based channel estimation method [4], our previously
proposed channel estimation exploiting DC window in the TF domain[6], and
NOMP-based channel estimation in [8], which are denoted as “Traditional
impulse”, “Traditional impulse + DC window”, “Off-grid NOMP”, respectively.
Besides, to demonstrate the performance gain of our proposed off-grid methods,
the channel estimation performance for on-grid schemes with
${\mathbf{k}_{{\nu}}}={\mathbf{l}_{{\tau}}}=\mathbf{0}$ based on OMP[7], 1D
SBL[9], and our proposed 2D SBL are also evaluated, which are denoted as “On-
grid OMP”, “1D on-grid SBL”, “2D on-grid SBL”, respectively. Note that the
additional steps introduced by our proposed off-grid methods is the hyper-
parameter update compared with the on-grid methods, which can be efficiently
computed. However, the computational complexities of both the on-grid and off-
grid methods are dominated by the matrix inversion operation and thus they are
in the same order.

Note that the traditional channel estimation methods of OTFS, such as
“Traditional impulse” and “Traditional impulse + DC window”, are mainly to
estimate the effective DD domain channel rather than the original DD domain
channel response as our proposed scheme in this paper. Therefore, for our
proposed schemes, we reconstruct the effective DD domain channel
${\hat{h}_{w}}\left[{k-k^{\prime},l-l^{\prime}}\right]$ based on the estimated
original DD domain channel response according to (3) and define the normalized
mean square error (NMSE) of channel estimation as

$\mathrm{NMSE}=\frac{\sum_{k-k^{\prime}=0}^{N-1}\sum_{l-l^{\prime}=0}^{M-1}\left|{h_{w}}\left[{k-k^{\prime},l-l^{\prime}}\right]-{\hat{h}_{w}}\left[{k-k^{\prime},l-l^{\prime}}\right]\right|^{2}}{\sum_{k-k^{\prime}=0}^{N-1}\sum_{l-l^{\prime}=0}^{M-1}\left|{h_{w}}\left[{k-k^{\prime},l-l^{\prime}}\right]\right|^{2}},\vspace{-1mm}$
(76)

where the summation is over $k-k^{\prime}\in\left[0,N-1\right]$ and
$l-l^{\prime}\in\left[0,M-1\right]$ as
${h_{w}}\left[{k-k^{\prime},l-l^{\prime}}\right]={h_{w}}\left[{(k-k^{\prime})_{N},(l-l^{\prime})_{M}}\right]$.
All the simulation results are averaged over $1,000$ OTFS frames.

Figure 5: NMSE of effective channel estimation versus SNR with a high virtual
sampling resolution $r_{\nu}=r_{\tau}=0.5$.

Figure 6: NMSE of effective channel estimation versus SNR with a low virtual
sampling resolution $r_{\nu}=r_{\tau}=0.8$.

###  VI-A NMSE versus SNR with Different Virtual Sampling Resolutions

Fig. 6 and Fig. 6 illustrate the performance of our proposed channel
estimation schemes with high and low virtual sampling resolutions,
respectively. In the case of a high virtual sampling resolution with
$r_{\nu}=r_{\tau}=0.5$, we can observe that our proposed 1D off-grid SBL-based
channel estimation achieves the lowest NMSE and the performance gain over the
on-grid method is about 1.5 dB in the high SNR regime. Moreover, the proposed
2D off-grid/on-grid SBL-based schemes achieve a coarser channel estimation
than the proposed 1D methods due to potential estimation error propagation
from the first step to the second step in Algorithm 2. However, owing to the
decoupling of our proposed 2D model, the computational complexity is
significantly reduced. Note that the 2D on-grid SBL-based scheme outperforms
the 2D off-grid SBL-based scheme. This is because the non-linear coupling of
the hyper-parameters $\left({{\bm{\kappa}}_{\nu}},\bm{\iota}_{\tau}\right)$ in
the measurement matrices in our proposed 2D off-grid model in (59) may reduce
the gain offered by estimating the off-grid components for a high virtual
sampling resolution. Additionally, the OMP-based methods, including both “On-
grid OMP” and “Off-grid NOMP”, achieve a better NMSE performance than the
impulse-based methods, including both “Traditional impulse” and “Traditional
impulse + DC window”, since the former methods can exploit the channel
sparsity in the DD domain. However, the OMP-based methods perform worse than
our proposed SBL-based methods as OMP is heuristic and is generally far from
the optimal solution. In fact, our proposed SBL-based methods are derived from
the maximum a posteriori probability estimate in (31) which can effectively
capture the main properties of the channel for estimation. Although the NOMP
method [8] can refine the delay and Doppler shifts estimation through Newton’s
method and achieves a 2 dB NMSE gain compared with the on-grid OMP method [7],
its NMSE is still substantially higher than that of our proposed SBL-based
schemes owing to the superiority of the SBL framework. Note that our previous
work of DC window-based channel estimation [6] can achieve a better
performance than the traditional impulse method in the high SNR regime thanks
to the improved channel sparsity obtained by windowing. However, the work in
[6] still suffers from a high NMSE compared with the compressed sensing-based
channel estimation approaches. In fact, the work in [6] estimates the
effective DD domain channel instead of the original DD domain channel
response, where the sparsity of the former channel is sacrificed compared with
the latter one due to inevitable channel spreading.

For the case of a low virtual sampling resolution with $r_{\nu}=r_{\tau}=0.8$,
we can observe a similar performance trend as Fig. 6 for the compressed
sensing-based channel estimation schemes. The only difference is that for our
proposed 2D model, the off-grid method outperforms the on-grid method since
the former can estimate the off-grid components of delay and Doppler shifts
and thus can achieve a high channel estimation accuracy. We can observe that
the OMP-based methods are even worse than the impulse-based approaches as the
former suffer from the low virtual sampling resolution. However, even with a
low virtual sampling resolution, our proposed 1D off-grid SBL-based channel
estimation can still achieve the best performance.

Figure 7: NMSE of effective channel estimation versus the virtual resolution.

Figure 8: NMSE of effective channel estimation versus SNR without guard space
for $r_{\nu}=r_{\tau}=0.8$.

###  VI-B NMSE versus the Virtual Sampling Resolution

To investigate the impact of the virtual sampling resolution, Fig. 8
illustrates the channel estimation performance versus the virtual sampling
resolution in the Doppler domain with a fixed virtual delay sampling
resolution $r_{\tau}=0.5$. We can observe that the NMSE of the compressed
sensing-based channel estimation increases with decreasing the virtual Doppler
resolution from $r_{\nu}=0.2$ to $r_{\nu}=0.9$. In contrast, the NMSE of the
impulse-based channel estimation methods is invariant since it estimates the
effective DD domain channel which does not depend on the virtual sampling
resolution. However, the computational complexity of the compressed sensing-
based approaches increases significantly with the increasing virtual sampling
resolution. For our proposed 2D model, the on-grid approach outperforms the
off-grid method with a high virtual Doppler resolution while the former scheme
is worse than the latter one with a low virtual Doppler resolution. In fact,
the 2D on-grid SBL-based scheme becomes more accurate for a high virtual
sampling resolution and enjoys a simpler estimation model with fewer unknown
parameters. With reducing the virtual sampling resolution, i.e., increasing
$r_{\nu}$, the NMSE of the proposed 2D on-grid SBL-based scheme becomes higher
since all the estimated delay and Doppler shifts are constrained on the coarse
DD domain grid. In contrast, our proposed off-grid approach can compensate the
loss of the low virtual sampling resolution through estimating the off-grid
components as hyper-parameters and thus enjoys a more accurate channel
estimation.

###  VI-C NMSE versus SNR without Guard Space

In practice, guard spaces are introduced to mitigate the interference from
data symbols to the pilot symbol. However, the insertion causes a
significantly high signaling overhead for OTFS. In particular, nulling the DD
domain space in the range of ${k_{p}}-2{k_{\max}}\leq
k\leq{k_{p}}+2{k_{\max}}$ and ${l_{p}}-{l_{\mathrm{max}}}\leq
l\leq{l_{p}}+{l_{\mathrm{max}}}$ results in
$\left(4{k_{\max}}+1\right)\left(2{l_{\mathrm{max}}}+1\right)$ overhead
symbols out of $MN$ symbols. Therefore, it is interesting to investigate the
robustness of our proposed channel estimation schemes against the interference
via evaluating the channel estimation performance without the use of guard
space. Fig. 8 demonstrates the NMSE of channel estimation without guard space
for $r_{\nu}=r_{\tau}=0.8$. We can observe that our proposed 1D off-grid SBL-
based channel estimation can achieve the lowest NMSE and possesses about 1 dB
gain compared to that of 1D on-grid SBL-based channel estimation, which is
consistent to previous results. Besides, it is worth to note that the on-grid
methods outperform the off-grid methods for both OMP-based and 2D SBL-based
channel estimation schemes. This is because the on-grid methods are based on a
linear SSR model while the off-grid methods are based on a non-linear model,
e.g. (59). As a result, the on-grid methods are more robust against the
interference than that of the off-grid methods. Also, comparing Fig. 6 and
Fig. 8, it can be seen that there is about $1\sim 2$ dB loss in channel
estimation NMSE for all the evaluated schemes due to the higher amount of
interference from data symbols without guard space. This motivates us to
further employ the proposed data-aided channel estimation scheme in next
subsection.

Figure 9: NMSE of effective channel estimation versus SNR with data-aided
channel estimation.

###  VI-D Data-Aided Channel Estimation

Fig. 9 illustrates the NMSE performance for different iterations in our
proposed 1D off-grid SBL-based channel estimation adopting the data-aided
philosophy[18]. We adopt binary phase shift keying (BPSK) modulation for
$x\left[k,l\right]$ and the sum-product algorithm proposed in our previous
work [6] for data detection. To improve the data detection performance, a
half-rate convolutional code and the maximum likelihood sequence estimation
are applied at OTFS transceiver. We can observe that about 1.5 dB decreases in
NMSE can be achieved by the first data-aided iteration, while the performance
improvement is marginal for the remaining iterations. In fact, the first
iteration can retrieve most of the unknown data symbols and provides the most
significant information to reconstruct the measurement matrix in (23). The
following iterations can further improve the accuracy of data detection based
on the refined channel estimates. However, the additional information provided
from the refined channel estimation is marginal compared with the first
iteration. Considering the high computational complexity of data detection and
channel estimation, we only need the first iteration to achieve an acceptable
channel estimation performance in practice.

##  VII Conclusions

In this paper, we proposed an off-grid SBL-based channel estimation scheme,
which estimates the original DD domain channel response instead of the
effective DD domain channel. In particular, the channel estimation problem was
reformulated as 1D and 2D off-grid SSR problems adopting linear approximation,
where the on-grid and off-grid components of both the delay and Doppler shifts
are separated for estimation. Based on the SBL framework, we developed 1D and
2D off-grid SBL-based channel estimation algorithms, where the on-grid
components of the delay and Doppler shifts are estimated based on the
recovered channel vector and the corresponding off-grid components are
estimated by the proposed EM algorithm. Extensive simulations have been
conducted to demonstrate the superior channel estimation performance of our
proposed schemes and revealed some interesting insights. Particularly, (1)
compared with estimating the effective DD domain channel, estimating the
original DD domain channel response can avoid channel spreading due to
fractional delay and Doppler shifts and thus can effectively exploit the DD
domain channel sparsity; (2) compared with the on-grid methods, the proposed
off-grid channel compressed sensing scheme can further improve the channel
estimation accuracy of OTFS systems even with a low virtual sampling
resolution while only introduces a slightly higher computational complexity;
(3) the proposed 1D off-grid SBL-based channel estimation scheme can achieve a
superior performance while the proposed 2D off-grid scheme enjoys a much lower
computational complexity via decoupling the delay and Doppler estimations; (4)
further adopting the data-aided channel estimation in our proposed scheme can
improve its robustness against the interference caused by unknown data symbols
and thus can further improve the channel estimation performance even in the
absence of guard space.

## Appendix: Proof of Proposition 1

To update ${\bm{\alpha}}$, the objective function is obtained by substituting
the items related to ${\bm{\alpha}}$ in (38) into (40), yielding

$\displaystyle
Q\left({{\bm{\alpha}}\left|{{{\bm{\alpha}}^{\left(t\right)}},{\bm{\kappa}}_{\nu}^{\left(t\right)},{\bm{\iota}}_{\tau}^{\left(t\right)},\beta_{0}^{\left(t\right)}}\right.}\right)={E_{\overline{\bf{h}}\left|{{\bf{y}}_{\mathrm{T}};{{\bm{\alpha}}^{\left(t\right)}},{\bm{\kappa}}_{\nu}^{\left(t\right)},{\bm{\iota}}_{\tau}^{\left(t\right)},\beta_{0}^{\left(t\right)}}\right.}}\left\\{{\ln\left[{p\left({\overline{\bf{h}}\left|{\bm{\alpha}}\right.}\right)p\left({\bm{\alpha}}\right)}\right]}\right\\}$
(77) $\displaystyle=$
$\displaystyle{E_{\overline{\bf{h}}\left|{{\bf{y}}_{\mathrm{T}};{{\bm{\alpha}}^{\left(t\right)}},{\bm{\kappa}}_{\nu}^{\left(t\right)},{\bm{\iota}}_{\tau}^{\left(t\right)},\beta_{0}^{\left(t\right)}}\right.}}\left\\{{-\sum\limits_{k^{\prime\prime}=0}^{N_{\nu}-1}{\sum\limits_{l^{\prime\prime}=0}^{M_{\tau}-1}{\ln\left({{\alpha_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}}}\right)}-{{\overline{\bf{h}}}^{\rm{H}}}{{\bm{\Lambda}}^{-1}}\overline{\bf{h}}-\rho\sum\limits_{k^{\prime\prime}=0}^{N_{\nu}-1}{\sum\limits_{l^{\prime\prime}=0}^{M_{\tau}-1}{{\alpha_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}}}}}}\right\\}.$

We can observe that the objective function in (77) is neither convex nor
concave w.r.t. ${\bm{\alpha}}$. In general, there is no systematic and
computational efficient solution to maximize $\eqref{ObjectiveFunctionAlpha}$.
As a compromise, we propose a closed-from updating rule of ${\bm{\alpha}}$,
which is computational efficient and is guaranteed to converge to a stationary
point of (77), as commonly done in the literature[10]. In particular, taking
the derivative of
$Q\left({{\bm{\alpha}}\left|{{{\bm{\alpha}}^{\left(t\right)}},{\bm{\kappa}}_{\nu}^{\left(t\right)},{\bm{\iota}}_{\tau}^{\left(t\right)},\beta_{0}^{\left(t\right)}}\right.}\right)$
w.r.t. ${\alpha_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}}$, we have

$\displaystyle\frac{{\partial
Q\left({{\bm{\alpha}}\left|{{{\bm{\alpha}}^{\left(t\right)}},{\bm{\kappa}}_{\nu}^{\left(t\right)},{\bm{\iota}}_{\tau}^{\left(t\right)},\beta_{0}^{\left(t\right)}}\right.}\right)}}{{\partial{\alpha_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}}}}$
$\displaystyle=-\frac{1}{{{\alpha_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}}}}$
$\displaystyle+\frac{1}{{\alpha_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}^{2}}}\\{{{\left|{\\{{{\bm{\mu}}^{\left(t\right)}_{\overline{\bf{h}}}}\\}_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}}\right|}^{2}}+\\{{{\bm{\Sigma}}^{\left(t\right)}_{\overline{\bf{h}}}}\\}_{k^{\prime\prime}M_{\tau}+l^{\prime\prime},k^{\prime\prime}M_{\tau}+l^{\prime\prime}}\\}-\rho.$
(78)

Imposing $\frac{{\partial
Q\left({{\bm{\alpha}}\left|{{{\bm{\alpha}}^{\left(t\right)}},{\bm{\kappa}}_{\nu}^{\left(t\right)},{\bm{\iota}}_{\tau}^{\left(t\right)},\beta_{0}^{\left(t\right)}}\right.}\right)}}{{\partial{\alpha_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}}}}=0$,
the updating rule of ${\alpha_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}}$
can be obtained by (44).

To update ${\beta_{0}}$, the objective function in (43) can be obtained by

$\displaystyle
Q\left({\beta_{0}\left|{{{\bm{\alpha}}^{\left(t\right)}},{\bm{\kappa}}_{\nu}^{\left(t\right)},{\bm{\iota}}_{\tau}^{\left(t\right)},\beta_{0}^{\left(t\right)}}\right.}\right)={E_{\overline{\bf{h}}\left|{{\bf{y}}_{\mathrm{T}};{{\bm{\alpha}}^{\left(t\right)}},{\bm{\kappa}}_{\nu}^{\left(t\right)},{\bm{\iota}}_{\tau}^{\left(t\right)},\beta_{0}^{\left(t\right)}}\right.}}\left\\{{\ln\left[{p\left({{\bf{y}}_{\mathrm{T}}\left|{\overline{\bf{h}},{{\bm{\kappa}}_{\nu}^{\left(t\right)}},{{\bm{\iota}}_{\tau}^{\left(t\right)}},{\beta_{0}}}\right.}\right)p\left({{\beta_{0}}}\right)}\right]}\right\\}$
$\displaystyle=$
$\displaystyle{E_{\overline{\bf{h}}\left|{{\bf{y}}_{\mathrm{T}};{{\bm{\alpha}}^{\left(t\right)}},{\bm{\kappa}}_{\nu}^{\left(t\right)},{\bm{\iota}}_{\tau}^{\left(t\right)},\beta_{0}^{\left(t\right)}}\right.}}\left\\{\ln\left[{{\cal
C}{\cal
N}\left({{\bf{y}}_{\mathrm{T}}\left|{\overline{\bm{\Phi}}_{\mathrm{T}}\left({{{\bm{\kappa}}_{\nu}^{\left(t\right)}},{{\bf{\iota}}_{\tau}^{\left(t\right)}}}\right)\overline{\bf{h}},\beta_{0}^{-1}{{\bf{I}}_{M_{\mathrm{T}}N_{\mathrm{T}}}}}\right.}\right)}\right]+\ln\left[{\Gamma\left({{\beta_{0}}\left|{c,d}\right.}\right)}\right]\right\\}$
$\displaystyle=$
$\displaystyle\left({c-1+M_{\mathrm{T}}N_{\mathrm{T}}}\right)\ln\left({{\beta_{0}}}\right)-\left({d+{A_{{\beta_{0}}}}}\right){\beta_{0}},$
(79)

where ${A_{{\beta_{0}}}}$ is given by (46). Taking derivative of
$Q\left({\beta_{0}\left|{{{\bm{\alpha}}^{\left(t\right)}},{\bm{\kappa}}_{\nu}^{\left(t\right)},{\bm{\iota}}_{\tau}^{\left(t\right)},\beta_{0}^{\left(t\right)}}\right.}\right)$
w.r.t. $\beta_{0}$ and imposing $\frac{{\partial
Q\left({\beta_{0}\left|{{{\bm{\alpha}}^{\left(t\right)}},{\bm{\kappa}}_{\nu}^{\left(t\right)},{\bm{\iota}}_{\tau}^{\left(t\right)},\beta_{0}^{\left(t\right)}}\right.}\right)}}{{\partial{\alpha_{k^{\prime\prime}M_{\tau}+l^{\prime\prime}}}}}=0$,
we can obtain the updating rule for $\beta_{0}$ as (45).

## References

  * [1] Z. Wei, W. Yuan, S. Li, J. Yuan, G. Bharatula, R. Hadani, and L. Hanzo, “Orthogonal time-frequency space modulation: A full-diversity next generation waveform,” _arXiv preprint arXiv:2010.03344_ , 2020.
  * [2] R. Hadani, S. Rakib, M. Tsatsanis, A. Monk, A. J. Goldsmith, A. F. Molisch, and R. Calderbank, “Orthogonal time frequency space modulation,” in _Proc. IEEE Wireless Commun. and Networking Conf._ , 2017, pp. 1–6.
  * [3] 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 preprint arXiv:2010.13008_ , 2020.
  * [4] P. Raviteja, K. T. Phan, and Y. Hong, “Embedded pilot-aided channel estimation for OTFS in delay-doppler channels,” _IEEE Trans. Veh. Technol._ , vol. 68, no. 5, pp. 4906–4917, May 2019.
  * [5] M. Kollengode Ramachandran and A. Chockalingam, “MIMO-OTFS in high-doppler fading channels: Signal detection and channel estimation,” in _Proc. IEEE Global Commun. Conf._ , Dec. 2018, pp. 206–212.
  * [6] Z. Wei, W. Yuan, S. Li, J. Yuan, and D. W. K. Ng, “Transmitter and receiver window designs for orthogonal time frequency space modulation,” _arXiv preprint arXiv:1709.02505_ , 2020.
  * [7] W. Shen, L. Dai, J. An, P. Fan, and R. W. Heath, “Channel estimation for orthogonal time frequency space (OTFS) massive MIMO,” _IEEE Trans. Signal Process._ , vol. 67, no. 16, pp. 4204–4217, Jul. 2019.
  * [8] M. Li, S. Zhang, F. Gao, P. Fan, and O. A. Dobre, “A new path division multiple access for the massive MIMO-OTFS networks,” _IEEE J. Select. Areas Commun._ , early access, 2020.
  * [9] L. Zhao, W. J. Gao, and W. Guo, “Sparse Bayesian learning of delay-doppler channel for OTFS system,” _IEEE Commun. Lett._ , vol. 24, no. 12, pp. 2766–2769, 2020.
  * [10] Z. Yang, L. Xie, and C. Zhang, “Off-grid direction of arrival estimation using sparse Bayesian inference,” _IEEE Trans. Signal Process._ , vol. 61, no. 1, pp. 38–43, Jan. 2013.
  * [11] S. D. Babacan, R. Molina, and A. K. Katsaggelos, “Bayesian compressive sensing using laplace priors,” _IEEE Trans. Image Process._ , vol. 19, no. 1, pp. 53–63, Jan. 2010.