# Orthogonal Time Frequency Space Modulation – Part I: Fundamentals and
Challenges Ahead

Zhiqiang Wei, Shuangyang Li, Weijie Yuan, Robert Schober, and Giuseppe Caire,
(_Invited Paper_) Z. Wei is with the School of Mathematics and Statistics,
Xi’an Jiaotong University, Xi’an 710049, China (e-mail:
zhiqiang.wei@xjtu.edu.au).S. Li was with the School of Electrical Engineering
and Telecommunications, University of New South Wales, Sydney, NSW 2052,
Australia, when this letter was submitted. He is now with the Department of
Electrical, Electronic, and Computer Engineering, University of Western
Australia, Perth, WA 6009, Australia (e-mail: shuangyang.li@uwa.edu.au). W.
Yuan is with the Department of Electrical and Electronic Engineering, Southern
University of Science and Technology, Shenzhen 518055, China (e-mail:
yuanwj@sustech.edu.cn).R. Schober is with the Institute for Digital
Communications (IDC), the Friedrich-Alexander University Erlangen-Nuremberg,
Erlangen 91054, Germany (e-mail: robert.schober@fau.de).G. Caire is with the
Electrical Engineering and Computer Science Department, the Technische
Universität Berlin, Berlin 10587, Germany (e-mail: caire@tu-berlin.de).

###### Abstract

This letter is the first part of a three-part tutorial on orthogonal time
frequency space (OTFS) modulation, which is a promising candidate waveform for
future wireless networks. This letter introduces and compares two popular
implementations of OTFS modulation, namely the symplectic finite Fourier
transform (SFFT)- and discrete Zak transform (DZT)-based architectures. Based
on these transceiver architectures, fundamental concepts of OTFS modulation,
including the delay-Doppler (DD) domain, DD domain information multiplexing,
and its potential benefits, are discussed. Finally, the challenges ahead for
OTFS modulation are highlighted. Parts II and III of this tutorial on OTFS
modulation focus on transceiver designs and integrated sensing and
communication (ISAC), respectively.

###### Index Terms:

OTFS, delay-Doppler domain, SFFT, DZT.

##  I Introduction

Next-generation wireless networks are required to support new applications and
challenging deployment scenarios, including mobile communications on board
aircraft (MCA), low-earth-orbit (LEO) satellites[1], high speed trains,
unmanned aerial vehicles (UAVs)[2], and vehicle-to-vehicle (V2V) networks[3].
In this context, high-mobility communications have become an essential
technology for future wireless networks[4]. Wireless communications in high-
mobility propagation environments experience fast time-varying channels, which
are typically caused by the Doppler effect[5]. If the Doppler effect was not
considered for system design, wireless communications over high-mobility
channels would suffer from significant performance degradation. The widely-
adopted orthogonal frequency division multiplexing (OFDM) modulation may not
be able to support efficient and reliable communications in high-mobility
environments due to the severe inter-carrier interference (ICI) induced by
Doppler spread. Besides, the limited coherence time of the time-frequency
domain channel will incur a large channel estimation overhead for OFDM
systems. For example, for an OFDM system operating at a carrier frequency of
$f_{c}=3.5$ GHz, employing a subcarrier spacing of $\Delta f=15$ kHz, and
supporting a relative velocity of $v=300$ km/h, the maximum Doppler shift is
$\nu_{\mathrm{max}}=972.22$ Hz while the OFDM symbol duration including a 20%
cyclic prefix (CP) is 80 ms. The channel’s coherence time is around
$\frac{1}{4\nu_{\mathrm{max}}}=257.14$ ms. Thus, one channel coherence
interval can only accommodate at most 3 OFDM symbols, i.e., a pilot symbol
needs to be inserted once every three OFDM symbols, causing a large overhead.

Recently, a new waveform, namely orthogonal time frequency space (OTFS)
modulation, was proposed to cater to high-mobility applications and has
attracted significant attention since its introduction[6, 4]. In particular,
for OTFS modulation, the information symbols are placed in the delay-Doppler
(DD) domain, wherein the time-varying channels are quasi-static and sparse[7].
Besides, DD domain information multiplexing simplifies the coupling between
the information symbols and the channels, which facilitates efficient and
reliable transceiver design. More importantly, owing to the two-dimensional
(2D) spreading inherent to OTFS modulation, each DD domain symbol experiences
the entire time-frequency (TF) domain channel response across one OTFS frame,
which enables OTFS to exploit the full degrees of freedom (DoFs) of high-
mobility channels. Moreover, delay and Doppler are typical sensing parameters
and the DD domain channel directly mirrors the geometry of the scatterers in
high-mobility environments[3, 8]. Therefore, OTFS modulation is regarded as a
promising waveform for integrated sensing and communication (ISAC) systems[9].

TABLE I: Notations for Main System Parameters. Notation | Physical meaning | Notation | Physical meaning
---|---|---|---
${h}_{\mathrm{DD}}\left(\tau,\nu\right)$ | Continuous DD domain channel response | ${H}_{\mathrm{DD}}\left[l,k,l^{\prime},k^{\prime}\right]$ | Discrete DD domain channel response
$w\left({l,k,l^{\prime},k^{\prime},l_{i},k_{i}}\right)$ | Sampling function | $\alpha\left[l,k,l_{i},k_{i}\right]$ | Phase rotation term
$\mathcal{DZ}_{x}\left[l,k\right]$ |  Discrete Zak transform of $x$ | ${{{A_{{g_{\rm tx}}{g_{\rm rx}}}}\left(\tau,\nu\right)}}$ | Ambiguity function

The research on OTFS modulation, including the derivation of its input-output
characteristics[10, 11, 12, 13], performance analysis[14, 15, 16, 17],
transceiver design[18, 19, 20, 21, 22], and its interplay with other promising
techniques[23, 24], is quickly evolving and attracts an increasing number of
researchers. However, it might be challenging for new researchers entering
this field to grasp the relevant mathematical tools and abstract concepts. The
authors of [4] have provided an overview of the basic concepts, transceiver
architectures, and waveforms of OTFS modulation, but did not provide any
mathematical details. This motivates this three-part tutorial, which aims to
provide an accessible introduction for new researchers interested in OTFS
modulation. In particular, Part I introduces two popular implementations of
OTFS modulation, reviews some related fundamental concepts, and discusses
future research challenges. Part II provides an in-depth discussion of OTFS
transceiver design, which is a key research direction in the OTFS literature.
Part III is dedicated to OTFS-based ISAC systems, which are regarded as an
enabling technology for the sixth-generation (6G) wireless networks, and
highlights a range of general research directions for OTFS modulation.

We note that there are alternative waveforms [25, 26], which are closely
related to OTFS modulation and can also facilitate high-mobility
communications. However, in this tutorial, we focus on OTFS modulation due to
the rapidly growing interest in this modulation format.

Notations: $\mathbb{C}^{M\times N}$, $\mathbb{R}^{M\times N}$, and
$\mathbb{Z}^{M\times N}$ denote the sets of $M\times N$ matrices with complex,
real, and integer entries, respectively; $\left(\cdot\right)^{*}$ denotes the
conjugate operation; $\left(\cdot\right)_{N}$ denotes the modulus operation
with respect to (w.r.t.) $N$; $\delta\left(\cdot\right)$ and
$\delta\left[\cdot\right]$ are the continuous- and discrete-time delta
functions, respectively. For clarity, we summarize the notations for the main
system parameters adopted in this paper in Table I.

##  II Two Implementations of OTFS Modulation

The two most commonly used implementations of OTFS modulation are symplectic
finite Fourier transform (SFFT)- [10, 11] and discrete Zak transform
(DZT)-based OTFS [12, 13]. In this section, we review the details of both
transceiver architectures and discuss their similarities and differences.

###  II-A SFFT-based OTFS

Figure 1: SFFT-based OTFS transceiver architecture.

The SFFT-based OTFS transceiver architecture is shown in Fig. 1[10, 11]. We
assume a bandwidth of $B_{\mathrm{OTFS}}$ and a time duration of
$T_{\mathrm{OTFS}}$, defined in the TF domain, are available for accommodating
one OTFS frame. Bandwidth $B_{\mathrm{OTFS}}$ is divided into $M$ subcarriers
with subcarrier spacing $\Delta f=\frac{B_{\mathrm{OTFS}}}{M}$ and time
duration $T_{\mathrm{OTFS}}$ is divided into $N$ time slots with slot duration
$T=\frac{T_{\mathrm{OTFS}}}{N}$. At the transmitter (TX) side, the information
symbol for the $\left(l,k\right)$-th DD resource element (DDRE) is denoted as
$X_{\mathrm{DD}}\left[{l,k}\right]$, where $l\in\\{0,\ldots,M-1\\}$ is the
delay index and $k\in\\{0,\ldots,N-1\\}$ is the Doppler index. The DD domain
symbols $X_{\mathrm{DD}}\left[{l,k}\right]$ are transformed into TF domain
signal $X_{\mathrm{TF}}\left[{m,n}\right]$ via the inverse symplectic finite
Fourier transform (ISFFT) [19], i.e.,

$X_{\mathrm{TF}}\left[{m,n}\right]=\frac{1}{{\sqrt{NM}}}\sum\limits_{k=0}^{N-1}{\sum\limits_{l=0}^{M-1}{X_{\mathrm{DD}}\left[{l,k}\right]{e^{j2\pi\left({\frac{{nk}}{N}-\frac{{ml}}{M}}\right)}}}},\vspace{-2mm}$
(1)

where $m\in\\{0,\ldots,M-1\\}$ is the subcarrier index and
$n\in\\{0,\ldots,N-1\\}$ is the slot index. The OTFS transmit signal
${s}\left(t\right)$ is obtained by performing multicarrier modulation on the
TF domain signal ${X}_{\mathrm{TF}}\left[{m,n}\right]$, also known as
Heisenberg transform [10], i.e.,

$s\left(t\right)=\sum\limits_{n=0}^{N-1}{\sum\limits_{m=0}^{M-1}{X_{\mathrm{TF}}\left[{m,n}\right]{{g_{{\rm{tx}}}}\left({t-nT}\right){e^{j2\pi
m\Delta f\left({t-nT}\right)}}}}},\vspace{-2mm}$ (2)

where ${{g_{{\rm{tx}}}}\left({t}\right)}$ is the TX pulse shaping filter. In
practice, the TF domain signal $X_{\mathrm{TF}}\left[{m,n}\right]$ is first
transformed to a time domain sequence via the inverse finite Fourier transform
(IFFT), is then pulse shaped in the discrete time domain, and is finally
transformed to the analog TX signal $s\left(t\right)$ by digital-to-analog
conversion (DAC).

A deterministic channel model in the DD domain has been widely adopted for
linear time-varying (LTV) channels [5, 27] and the continuous DD domain
channel response is modeled as

${h}_{\mathrm{DD}}\left(\tau,\nu\right)=\sum\nolimits_{i=1}^{P}{h_{i}}\delta\left(\tau-\tau_{i}\right)\delta\left(\nu-\nu_{i}\right),\vspace{-2mm}$
(3)

where $P\in\mathbb{Z}$ is the number of resolvable paths and
${h}_{i}\in\mathbb{C}$ is the channel coefficient of the $i$-th path.
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, where
$\tau_{\mathrm{max}}$ and $\nu_{\mathrm{max}}$ denote the maximum delay and
Doppler shifts in the considered high-mobility propagation environment,
respectively. Note that DD domain channel responses fluctuate much slower than
time-delay domain or TF domain channel responses, since only a drastic change
of the propagation path lengths or the moving speeds would cause DD domain
channel response variations.

After passing through the LTV channel, the received signal at the OTFS
receiver (RX) is given by

$\displaystyle r\left(t\right)$
$\displaystyle=\int_{-\infty}^{\infty}{\int_{-\infty}^{\infty}{h_{\mathrm{DD}}\left({\tau,\nu}\right)}}{{e^{j2\pi\nu\left({t-\tau}\right)}}}s\left({t-\tau}\right)\mathrm{d}\tau\mathrm{d}\nu+z\left(t\right)$
$\displaystyle=\sum\nolimits_{i=1}^{P}{h_{i}}s\left({t-\tau_{i}}\right){{e^{j2\pi\nu_{i}\left({t-\tau_{i}}\right)}}}+z\left(t\right),$
(4)

where $z\left(t\right)$ denotes the additive white Gaussian noise (AWGN) in
the time domain. At the OTFS RX, after receive filtering, we first perform
multicarrier demodulation of $r\left(t\right)$ to obtain the TF domain
received signal. This operation is also known as Wigner transform[10] and
given by

${Y}_{\mathrm{TF}}\left[m,n\right]=\int_{-\infty}^{\infty}{r\left(t\right)g_{{\rm{rx}}}^{*}\left({t-nT}\right){e^{-j2\pi
m\Delta f\left({t-nT}\right)}}\mathrm{d}t},\vspace{-3mm}$ (5)

where ${{g_{{\rm{rx}}}}\left({t}\right)}$ is the RX filter. Then, the TF
domain received signal is transformed into the DD domain by the SFFT. The
resulting DD domain received signal is given by

$Y_{\mathrm{DD}}\left[{l,k}\right]=\frac{1}{{\sqrt{NM}}}\sum\limits_{n=0}^{N-1}{\sum\limits_{m=0}^{M-1}{{Y}_{\mathrm{TF}}\left[{m,n}\right]{e^{-j2\pi\left({\frac{{kn}}{N}-\frac{{lm}}{M}}\right)}}}}.\vspace{-3mm}$
(6)

Note that (2) and (5) degenerate to conventional OFDM modulation and
demodulation, respectively, when $N=1$ and rectangular TX and RX filters are
adopted. Moreover, when $N=1$ and no restrictions on the form of the TX and RX
filters are imposed, (2) and (5) degenerate to pulse-shaped OFDM modulation
and demodulation, respectively [28].

$\displaystyle w\left({l,k,l^{\prime},k^{\prime},l_{i},k_{i}}\right)$
$\displaystyle=\frac{1}{{NM}}\sum\nolimits_{n=0}^{N-1}\sum\nolimits_{m=0}^{M-1}\left[\sum\nolimits_{m^{\prime}\neq
m}^{M-1}{A_{{g_{\rm tx}}{g_{\rm rx}}}}\left({-{l_{i}}\frac{1}{{M\Delta
f}},\left({m-m^{\prime}}\right)\Delta f-{k_{i}}\frac{1}{{NT}}}\right)\right.$
$\displaystyle\left.+{{e^{-j2\pi\frac{{k^{\prime}}}{N}}}\sum\nolimits_{m^{\prime}=0}^{M-1}{{A_{{g_{\rm
tx}}{g_{\rm rx}}}}\left({T-{l_{i}}\frac{1}{{M\Delta
f}},\left({m-m^{\prime}}\right)\Delta
f-{k_{i}}\frac{1}{{NT}}}\right)}}\right]{e^{-j2\pi\frac{{n\left({k-k^{\prime}-{k_{i}}}\right)}}{N}}}{e^{j2\pi\frac{{\left({ml-m^{\prime}l^{\prime}-m^{\prime}{l_{i}}}\right)}}{M}}}$
(9)
$\alpha\left[l,k,l_{i},k_{i}\right]=\left\\{{\begin{array}[]{*{20}{c}}{e^{j2\pi\frac{\left(l-l_{i}\right)_{M}k_{i}}{MN}}}&{l\in\left\\{l_{i},\ldots,M-1\right\\}}\\\
{e^{j2\pi\frac{\left(l-l_{i}\right)_{M}k_{i}-k_{i}M-\left(k-k_{i}\right)_{N}M}{MN}}}&{l\in\left\\{0,\ldots,l_{i}-1\right\\}}\end{array}}\right..\vspace{-2mm}$
(12)

Combining (1)-(6), the DD domain input-output relationship for OTFS modulation
is given by

$Y_{\mathrm{DD}}\left[{l,k}\right]=\sum\limits_{k^{\prime}=0}^{N-1}{\sum\limits_{l^{\prime}=0}^{M-1}{X_{\mathrm{DD}}\left[{l^{\prime},k^{\prime}}\right]}}{H}_{\mathrm{DD}}\left[{l,k,l^{\prime},k^{\prime}}\right]+Z_{\mathrm{DD}}\left[l,k\right],\vspace{-3mm}$
(7)

where $Z_{\mathrm{DD}}\left[l,k\right]$ denotes the effective DD domain noise.
The discrete DD domain channel
${H}_{\mathrm{DD}}\left[{l,k,l^{\prime},k^{\prime}}\right]$ is given by

${H}_{\mathrm{DD}}\left[{l,k,l^{\prime},k^{\prime}}\right]=\sum\nolimits_{i=1}^{P}{h_{i}}w\left({l,k,l^{\prime},k^{\prime},l_{i},k_{i}}\right){e^{-j2\pi{\nu_{i}}{\tau_{i}}}},\vspace{-2mm}$
(8)

where $k_{i}={\nu_{i}}NT\in\mathbb{R}$ and $l_{i}={\tau_{i}}M\Delta
f\in\mathbb{R}$. We distinguish between integer and fractional Doppler and
delay scenarios based on whether $k_{i}$ and $l_{i}$ are integer or fractional
numbers, respectively[10]. In (8),
$w\left({l,k,l^{\prime},k^{\prime},l_{i},k_{i}}\right)$ is the sampling
function for analog DD domain channel response, which captures joint effects
of delay and Doppler shifts, TX pulse shaping filter, and RX filter, and is
given by (II-A) at the top of next page, where the ambiguity function
${{A_{{g_{\rm tx}}{g_{\rm rx}}}}\left(\tau,\nu\right)}$ is defined as

${{{A_{{g_{\rm tx}}{g_{\rm
rx}}}}\left(\tau,\nu\right)}}=\int_{-\infty}^{\infty}{g_{{\rm{tx}}}}\left({t}\right)g_{{\rm{rx}}}^{*}\left({t-\tau}\right){e^{-j2\pi\nu\left({t-\tau}\right)}}\mathrm{d}t.\vspace{-3mm}$
(10)

For the widely adopted rectangular TX and RX filters and assuming integer
$k_{i}$ and $l_{i}$[10], (7) simplifies to the well-known 2D circular
convolution input-output relationship as follows:

$\displaystyle Y_{\mathrm{DD}}\left[{l,k}\right]$
$\displaystyle=\sum\nolimits_{i=1}^{P}h_{i}X_{\mathrm{DD}}\left[\left(l-l_{i}\right)_{M},\left(k-k_{i}\right)_{N}\right]\alpha\left[l,k,l_{i},k_{i}\right]$
$\displaystyle+Z_{\mathrm{DD}}\left[l,k\right],$ (11)

where $\alpha\left[l,k,l_{i},k_{i}\right]$ is a phase rotation term defined in
(12) at the top of next page.

###  II-B DZT-based OTFS

Figure 2: DZT-based OTFS transceiver architecture.

SFFT-based OTFS modulation maps the DD domain information symbols first to the
TF domain and then to the time domain[10, 11]. In contrast, DZT-based OTFS
modulation directly transforms the DD domain information symbols into the time
domain, as shown in Fig. 2[12, 13]. The DZT and inverse DZT (IDZT) define a
mapping between a one-dimensional (1D) sequence and a 2D sequence. Their
definitions are provided in the following.

######  Definition 1 (Discrete Zak Transform)

Let $x$ be a periodic sequence with period $MN$. The DZT of $x$ is defined as

$\mathcal{DZ}_{x}\left[l,k\right]\buildrel\Delta\over{=}\frac{1}{\sqrt{N}}\sum\nolimits_{n=0}^{N-1}x\left[l+nM\right]e^{-j2\pi\frac{n}{N}k},\vspace{-2mm}$
(13)

for $l\in\\{0,\ldots,M-1\\}$ and $k\in\\{0,\ldots,N-1\\}$.

######  Definition 2 (Inverse Discrete Zak Transform)

Given the $MN$-periodic sequence $x$, whose DZT is given by
$\mathcal{DZ}_{x}\left[l,k\right]$, the IDZT is defined as

$x\left[l+nM\right]\buildrel\Delta\over{=}\frac{1}{\sqrt{N}}\sum\nolimits_{k=0}^{N-1}\mathcal{DZ}_{x}\left[l,k\right]e^{j2\pi\frac{n}{N}k},\vspace{-2mm}$
(14)

for $l\in\\{0,\ldots,M-1\\}$ and $n\in\\{0,\ldots,N-1\\}$.

As shown in [13], a time-delay (TD) domain sequence can be obtained by
applying the IDZT to the DD domain symbols at the OTFS TX, i.e.,

${x}_{\mathrm{TD}}\left[{l+nM}\right]=\sqrt{M}\sum\nolimits_{k=0}^{K-1}{X_{\mathrm{DD}}\left[{l,k}\right]\mathcal{DZ}_{g_{\mathrm{tx}}}\left[l,k\right]{e^{j2\pi\frac{n}{N}k}}},\vspace{-1mm}$
(15)

where $\mathcal{DZ}_{g_{\mathrm{tx}}}\left[l,k\right]$ is the DZT of the
periodically extended sampled TX pulse shape
$g_{\mathrm{tx}}\left(\frac{[k]_{MN}}{M}T\right)$, $k\in\mathbb{Z}$. The OTFS
transmit signal $s(t)$ can be obtained by DAC of
${x}_{\mathrm{TD}}\left[{l+nM}\right]$. Interpolating with a rectangular
kernel, the OTFS transmit signal is given by

$s(t)={x}_{\mathrm{TD}}\left[{l+nM}\right],t\in\left[\frac{l-0.5}{M}T+nT,\frac{l+0.5}{M}T+nT\right].\vspace{-1.5mm}$
(16)

At the OTFS RX, the channel output $r(t)$ in (II-A) is subjected to an analog-
to-digital conversion (ADC) and the resulting TD domain received signal is
given by

$y_{\mathrm{TD}}\left[{l+nM}\right]=r(t)|_{t=\frac{l}{M}T+nT}.\vspace{-2mm}$
(17)

Then, the DD domain received signal is obtained by performing DZT on the TD
domain received signal, i.e.,

$Y_{\mathrm{DD}}\left[{l,k}\right]=\sqrt{MN}\mathcal{DZ}_{y_{\mathrm{TD}}}\left[l,k\right]\mathcal{DZ}^{*}_{g_{\mathrm{rx}}}\left[l,k\right],\vspace{-1.5mm}$
(18)

where $\mathcal{DZ}_{y_{\mathrm{TD}}}\left[l,k\right]$ is the DZT of
$y_{\mathrm{TD}}\left[{l+nM}\right]$ and
$\mathcal{DZ}_{g_{\mathrm{rx}}}\left[l,k\right]$ is the DZT of the
periodically extended sampled RX filter
$g_{\mathrm{rx}}\left(\frac{[k]_{MN}}{M}T\right)$, $k\in\mathbb{Z}$.

A common choice for the TX and RX filters are rectangular filters. In this
case, we have
$\mathcal{DZ}_{g_{\mathrm{tx}}}\left[l,k\right]=\mathcal{DZ}_{g_{\mathrm{rx}}}\left[l,k\right]=\frac{1}{\sqrt{MN}}$,
$\forall l,k$. Moreover, assuming both integer delays and integer Doppler
shifts, the DD domain input-output relationship can be obtained as [13]

$\displaystyle Y_{\mathrm{DD}}\left[{l,k}\right]$
$\displaystyle=\sum\nolimits_{i=1}^{P}h_{i}X_{\mathrm{DD}}\left[\left(l-l_{i}\right)_{M},\left(k-k_{i}\right)_{N}\right]\alpha\left[k,l,k_{i},l_{i}\right]$
$\displaystyle+Z_{\mathrm{DD}}\left[l,k\right],$ (19)

which is identical to (II-A). We refer to [13] for a more detailed discussion
of the input-output relationship for non-rectangular TX and RX filters,
fractional delays, and fractional Doppler shifts.

###  II-C Their Similarities and Differences

The SFFT- and DZT-based OTFS implementations are both suitable for high-
mobility communications. Besides, for rectangular TX and RX filters, integer
delays, and integer Doppler shifts, SFFT- and DZT-based OTFS lead to the same
input-output relationship, see (II-A) and (II-B). Unveiling their connections
for the case of non-rectangular filters, fractional delays, and fractional
Doppler shifts is still an open problem in the literature.

SFFT-based OTFS is compatible with conventional OFDM by concatenating the
ISFFT/SFFT module and an OFDM modulator/demodulator, while DZT-based OTFS
enjoys a lower computational complexity as only $N$-point inverse finite
Fourier transform (IFFT) and finite Fourier transform (FFT) are needed with a
symbol-wise interleaver/de-interleaver at the OTFS transceivers, respectively,
see (15) and (18). Moreover, the SFFT-based OTFS architecture facilitates TF
domain signal processing, such as window function design[19] and design of
multiple access schemes[29]. On the other hand, for DZT-based OTFS, the TX and
RX filters can be implemented digitally in the DD domain based on their DZTs.
This not only reduces the implementation complexity, but also facilitates a
systematic pulse shape design in the DD domain. A detailed discussion on pulse
shape design can be found in Part II of this tutorial.

##  III Discussions on OTFS Modulation

Based on the transceiver architectures introduced in Section II, we shed light
on some fundamental concepts of OTFS modulation in this section.

###  III-A What is the Delay-Doppler Domain

In the context of signal processing, we define and describe a signal in a
certain domain. Describing a signal in the time domain characterizes the
fluctuation of its amplitude and phase w.r.t. the time variable, while
representing a signal in the frequency domain characterizes the corresponding
fluctuation rate. Furthermore, the TF domain is expanded by the time duration
and frequency bandwidth occupied by a signal. The characterization of a signal
in the TF domain reveals the temporal evolution of its spectral content
(within a certain time duration). On the other hand, the delay and Doppler
parameters are usually used to characterize a time-varying system rather than
a signal. For example, $h\left(\tau,\nu\right)$ in (3) characterizes the
response of a time-varying system that introduces delays $\tau_{i}$, $\forall
i$, and Doppler shifts $\nu_{i}$, $\forall i$, to the input signal. OTFS
modulation is based on signal representation theory[7], cf. (1), (2), and
(15), which allows to characterize a signal directly in the DD domain. In
particular, a DD domain impulse $\delta\left[l-l_{0},k-k_{0}\right]$,
$l,l_{0}\in\\{0,\ldots,M-1\\}$, $k,k_{0}\in\\{0,\ldots,N-1\\}$ corresponds to
a time-varying pulsetone in the time domain[7], since it can be viewed as the
output of a virtual time-varying system that introduces delay
$l_{0}\frac{1}{M\Delta f}$ and Doppler shift $k_{0}\frac{1}{NT}$ w.r.t. an
input DD domain impulse $\delta\left[l,k\right]$. At the OTFS RX, when
transforming a signal from the time or TF domain to the DD domain, cf. (5),
(6), and (18), we are inverting the channel operations the TX signal has
experienced. According to the uncertainty principle of the Fourier transform,
the longer the TX signal duration, the higher the resolution for
distinguishing different Doppler shifts in the channel, while the larger the
TX signal frequency bandwidth, the higher the resolution for distinguishing
different delay shifts in the channel [12].

###  III-B How to Embed Information into the DD Domain

As shown in (15), the OTFS TX modulates information symbols
$X_{\mathrm{DD}}\left[{l,k}\right]$ onto the DD domain pulse samples
$\mathcal{DZ}_{g_{\mathrm{tx}}}\left[l,k\right]$, $l\in\\{0,\ldots,M-1\\}$,
$k\in\\{0,\ldots,N-1\\}$, which can be interpreted as the generation of a set
of time-varying carrier signals in the time domain to carry the corresponding
information symbol. A set of 2D impulses with zero bandwidth in the DD domain,
i.e., $\delta\left[l-l_{0},k-k_{0}\right]$, $l,l_{0}\in\\{0,\ldots,M-1\\}$,
$k,k_{0}\in\\{0,\ldots,N-1\\}$, would be ideal for carrying information since
they are perfectly localized[12]. Unfortunately, this is not possible in
practice since such impulses would require infinite time duration and infinite
frequency bandwidth. As shown in [12], 2D sinc pulses (their samples are given
by $\mathcal{DZ}_{g_{\mathrm{tx}}}\left[l,k\right]$ in (15),
$l\in\\{0,\ldots,M-1\\}$, $k\in\\{0,\ldots,N-1\\}$) with finite time duration
and frequency bandwidth, which are approximately orthogonal in the DD domain,
constitute an efficient set of DD domain information-carrying impulses.

###  III-C Benefits of DD Domain Modulation

Due to the Heisenberg uncertainty principle, it is impossible to find a 2D
signal perfectly localized in the TF domain, while in the DD domain, we can
obtain a set of 2D localized signals (their samples are given by
$\mathcal{DZ}_{g_{\mathrm{tx}}}\left[l,k\right]$ in (15)) that are
approximately orthogonal [12]. Moreover, as mentioned before, the OTFS
modulator employs time-varying carrier signals for carrying information, which
caters to the time-varying features of high-mobility channels. In particular,
the high-mobility channel, corresponding to DD domain channel response,
$h_{\mathrm{DD}}\left(\tau,\nu\right)$ in (3) can be modeled as a linear time-
varying system concatenated with the virtual time-varying system formed by the
OTFS modulator. As a result, the high-mobility channel additively changes the
time-varying features of the carrier signals, i.e., the received carrier
signal experiences a shorter or a longer delay ($l\pm l_{i}$) and a lower or a
higher Doppler shift ($k\pm k_{i}$), compared to the transmitted carrier
signal in the $\left(l,k\right)$-th DDRE. As a result, the additive channel
operations on the carrier signals, i.e., the delay and Doppler shifts, can be
simply described by a shift operation and a phase rotation in the DD domain,
i.e., (II-A) and (II-B), which is beneficial for efficient and reliable
channel equalization.

##  IV Challenges Ahead

So far, the research on OTFS modulation has focused on the characterization of
the input-output relationship, transceiver design, and potential synergies
with other communication and signal processing techniques, as elaborated in
Parts I, II, and III, respectively. However, some aspects of OTFS modulation
are not well understood or have not been addressed by the research community,
yet. In this section, we discuss important open problems, which might be the
main challenges ahead for OTFS research.

###  IV-A Latency Concern

Compared to OFDM, OTFS modulation spreads information over a longer block to
combat channel fluctuations, which may increase the overall latency. In
particular, the duration of each OTFS signal is $N$ times of that of an OFDM
symbol, with the same number of subcarriers and the same subcarrier
spacing[10]. Fortunately, multiplexing information in the DD domain turns a
time-varying channel into a quasi-static channel in the DD domain, which might
significantly reduce the amount of signaling overhead required for channel
estimation. A comprehensive comparative study of the latency of OTFS and OFDM
is still missing. Besides, the latency of OTFS can be made identical to that
of OFDM by increasing the subcarrier spacing or reducing $N$ at the expense of
a reduced Doppler resolution. A detailed performance comparison of OTFS and
OFDM for a given latency is an interesting research direction.

###  IV-B Fractional Doppler

To limit the latency of OTFS modulation, the signal duration $NT$ cannot be
chosen exceedingly large, which gives rise to the well-known fractional
Doppler issue[10]. Existing works in the literature have addressed the DD
domain channel estimation problem in the case of fractional Doppler by various
means, e.g., using windowing[19], oversampling[30], and off-grid
estimation[20]. However, detection performance may significantly degrade in
the case of fractional Doppler. In fact, the commonly used message passing-
based detection algorithms [10] are prone to diverge, due to the loops in the
factor graph caused by fractional Doppler shift. How to improve detection
performance and robustness in the case of fractional Doppler is a critical
issue. We note that the authors in [21, 31] have proposed a cross domain
message passing detector and a variational Bayes detector, respectively, which
can improve the convergence performance in the case of fractional Doppler.

###  IV-C Multiple Access

Serving multiple users in high-mobility environments is a challenging task.
For OTFS modulation, it is not clear yet in which domain the user should be
multiplexed. Conventional TF domain multiple access schemes may suffer from
severe multi-user interference (MUI) caused by doubly selective fading
channels. On the other hand, DD domain multiple access requires careful
allocation of the DDREs and the guard space between users to minimize MUI.
Moreover, the 2D circular convolution input-output relationship of OTFS leads
to special MUI patterns that might be exploited for user scheduling and
resource allocation design. Furthermore, compared to orthogonal multiple
access, power domain and/or code domain non-orthogonal multiple access schemes
[32] can accommodate more users while imposing a critical challenge for
interference management.

###  IV-D MU MIMO-OTFS

Multiple-input multiple-output (MIMO) systems can provide additional spatial
DoFs for OTFS systems, which may be exploited to improve communication
reliability and the achievable rate. Moreover, these additional DoFs can be
used to accommodate multiple users in the spatial domain, which would allow to
avoid the guard resource elements needed in existing DD and TF domain multiple
access schemes[18]. However, incorporating OTFS modulation in MU MIMO systems
is a challenging task since the signals have to be optimized simultaneously in
the delay, Doppler, and spatial domains. The authors of [33] investigated the
precoding design for massive fully-digital MIMO-OTFS systems, which may incur
exceedingly high hardware complexity and power consumption. In fact, a general
and concise MU MIMO-OTFS model, based on which channel estimation, precoding,
and multiple access schemes can be designed, has not been reported so far.

###  IV-E Cross-layer Design

High-mobility communication networks, such as UAV and V2V networks, may suffer
not only from time-varying channels but also from a highly dynamic network
topology. However, existing works on OTFS modulation mainly focus on the
physical layer design to improve communication performance. Cross-layer design
spanning both the physical and network layers for OTFS systems is a critical
but challenging task. In fact, cross-layer design facilitates coordination,
interaction, and joint optimization of communication protocols across
different layers, which is beneficial for mobility management and energy-
efficient and reliable communications. In high-mobility environments, it is
critical to sense the fast time-varying network state information (NSI), such
as the location and velocity of the network elements exploiting the OTFS
modulation waveform, and to design the OTFS transceivers accordingly.
Furthermore, the network topology may be actively adjusted depending on the
state of the high-mobility environment to improve the communication
performance.

##  V Conclusions

As the first part of this three-part tutorial, this letter presented and
compared two common implementations of OTFS, namely SFFT- and DZT-based OTFS.
Based on the introduced OTFS transceiver architectures, some fundamental
concepts of OTFS modulation were discussed. Subsequently, several important
research challenges regarding OTFS modulation were highlighted. Part II and
Part III will discuss OTFS transceiver design and the benefits of OTFS for
ISAC and other applications, respectively.

## References

  * [1] N. Cheng, J. He, Z. Yin, C. Zhou, H. Wu, F. Lyu, H. Zhou, and X. Shen, “6G service-oriented space-air-ground integrated network: A survey,” __to appear in_ Chinese Journal of Aeronautics_, Dec. 2021.
  * [2] Y. Cai, Z. Wei, R. Li, D. W. K. Ng, and J. Yuan, “Joint trajectory and resource allocation design for energy-efficient secure UAV communication systems,” _IEEE Trans. Commun._ , vol. 68, no. 7, pp. 4536–4553, Mar. 2020.
  * [3] W. Yuan, Z. Wei, S. Li, J. Yuan, and D. W. K. Ng, “Integrated sensing and communication-assisted orthogonal time frequency space transmission for vehicular networks,” _IEEE J. Select. Topics Signal Process._ , vol. 15, no. 6, pp. 1515–1528, Nov. 2021.
  * [4] Z. Wei, W. Yuan, S. Li, J. Yuan, G. Bharatula, R. Hadani, and L. Hanzo, “Orthogonal time-frequency space modulation: A promising next-generation waveform,” _IEEE Wireless Commun._ , vol. 28, no. 4, pp. 136–144, Aug. 2021\.
  * [5] F. Hlawatsch and G. Matz, _Wireless Communications Over Rapidly Time-varying Channels_. Academic Press, 2011.
  * [6] 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.
  * [7] R. Hadani and A. Monk, “OTFS: A new generation of modulation addressing the challenges of 5G,” _arXiv preprint arXiv:1802.02623_ , 2018.
  * [8] L. Gaudio, M. Kobayashi, G. Caire, and G. Colavolpe, “On the effectiveness of OTFS for joint radar parameter estimation and communication,” _IEEE Trans. Wireless Commun._ , vol. 19, no. 9, pp. 5951–5965, Sep. 2020.
  * [9] S. Li, W. Yuan, C. Liu, Z. Wei, J. Yuan, B. Bai, and D. W. K. Ng, “A novel ISAC transmission framework based on spatially-spread orthogonal time frequency space modulation,” __to appear in_ IEEE J. Select. Areas Commun._, 2022.
  * [10] P. Raviteja, K. T. Phan, Y. Hong, and E. Viterbo, “Interference cancellation and iterative detection for orthogonal time frequency space modulation,” _IEEE Trans. Wireless Commun._ , vol. 17, no. 10, pp. 6501–6515, Oct. 2018.
  * [11] A. Farhang, A. RezazadehReyhani, L. E. Doyle, and B. Farhang-Boroujeny, “Low complexity modem structure for OFDM-based orthogonal time frequency space modulation,” _IEEE Commun. Lett._ , vol. 7, no. 3, pp. 344–347, Jun. 2018.
  * [12] S. K. Mohammed, “Derivation of OTFS modulation from first principles,” _IEEE Trans. Veh. Technol._ , vol. 70, no. 8, pp. 7619–7636, Aug. 2021.
  * [13] F. Lampel, A. Alvarado, and F. M. Willems, “Orthogonal time frequency space modulation: A discrete Zak transform approach,” _arXiv preprint arXiv:2106.12828_ , 2021.
  * [14] G. D. Surabhi, R. M. Augustine, and A. Chockalingam, “On the diversity of uncoded OTFS modulation in doubly-dispersive channels,” _IEEE Trans. Wireless Commun._ , vol. 18, no. 6, pp. 3049–3063, Jun. 2019.
  * [15] 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,” _IEEE Trans. Wireless Commun._ , vol. 20, no. 9, pp. 6033–6048, Sep. 2021.
  * [16] A. RezazadehReyhani, A. Farhang, M. Ji, R. R. Chen, and B. Farhang-Boroujeny, “Analysis of discrete-time MIMO OFDM-Based orthogonal time frequency space modulation,” in _Proc. IEEE Intern. Commun. Conf._ , May 2018, pp. 1–6.
  * [17] R. Chong, S. Li, J. Yuan, and D. W. K. Ng, “Achievable rate upper-bounds of uplink multiuser OTFS transmissions,” __to appear in_ IEEE Wireless Commun. Lett._, 2022.
  * [18] 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.
  * [19] Z. Wei, W. Yuan, S. Li, J. Yuan, and D. W. K. Ng, “Transmitter and receiver window designs for orthogonal time-frequency space modulation,” _IEEE Trans. Commun._ , vol. 69, no. 4, pp. 2207–2223, Apr. 2021.
  * [20] ——, “Off-grid channel estimation with sparse Bayesian learning for OTFS systems,” __to appear in_ IEEE Trans. Wireless Commun._, 2022\.
  * [21] S. Li, W. Yuan, Z. Wei, and J. Yuan, “Cross domain iterative detection for orthogonal time frequency space modulation,” _IEEE Trans. Wireless Commun._ , vol. 21, no. 4, pp. 2227–2242, Apr. 2022.
  * [22] Z. Yuan, F. Liu, W. Yuan, Q. Guo, Z. Wang, and J. Yuan, “Iterative detection for orthogonal time frequency space modulation with unitary approximate message passing,” _IEEE Trans. Wireless Commun._ , vol. 21, no. 2, pp. 714–725, Jul. 2022.
  * [23] Z. Ding, R. Schober, P. Fan, and H. Vincent Poor, “OTFS-NOMA: An efficient approach for exploiting heterogenous user mobility profiles,” _IEEE Trans. Commun._ , vol. 67, no. 11, pp. 7950–7965, Nov. 2019.
  * [24] L. Gaudio, M. Kobayashi, G. Caire, and G. Colavolpe, “Hybrid digital-analog beamforming and MIMO radar with OTFS modulation,” _arXiv preprint arXiv:2009.08785_ , 2020.
  * [25] T. Thaj, E. Viterbo, and Y. Hong, “Orthogonal time sequency multiplexing modulation: Analysis and low-complexity receiver design,” _IEEE Trans. Wireless Commun._ , vol. 20, no. 12, pp. 7842–7855, Dec. 2021.
  * [26] T. Zemen, M. Hofer, D. Löschenbrand, and C. Pacher, “Iterative detection for orthogonal precoding in doubly selective channels,” in _Proc. IEEE Personal, Indoor and Mobile Radio Commun. Sympos._ , Sep. 2018, pp. 1–7.
  * [27] Y. Hong, T. Thaj, and E. Viterbo, _Delay-Doppler Communications: Principles and Applications_. Academic Press, 2022.
  * [28] T. Strohmer and S. Beaver, “Optimal OFDM design for time-frequency dispersive channels,” _IEEE Trans. Commun._ , vol. 51, no. 7, pp. 1111–1122, Jul. 2003.
  * [29] V. Khammammetti and S. K. Mohammed, “OTFS-based multiple-access in high Doppler and delay spread wireless channels,” _IEEE Wireless Commun. Lett._ , vol. 8, no. 2, pp. 528–531, Apr. 2019.
  * [30] Y. Ge, Q. Deng, P. C. Ching, and Z. Ding, “Receiver design for OTFS with a fractionally spaced sampling approach,” _IEEE Trans. Wireless Commun._ , vol. 20, no. 7, pp. 4072–4086, Jul. 2021.
  * [31] W. Yuan, Z. Wei, J. Yuan, and D. W. K. Ng, “A simple variational Bayes detector for orthogonal time frequency space (OTFS) modulation,” _IEEE Trans. Veh. Technol._ , vol. 69, no. 7, pp. 7976–7980, Jul. 2020.
  * [32] Z. Wei, L. Zhao, J. Guo, D. W. K. Ng, and J. Yuan, “Multi-beam NOMA for hybrid mmwave systems,” _IEEE Trans. Commun._ , vol. 67, no. 2, pp. 1705–1719, Feb. 2019.
  * [33] B. C. Pandey, S. K. Mohammed, P. Raviteja, Y. Hong, and E. Viterbo, “Low complexity precoding and detection in multi-user massive MIMO OTFS downlink,” _IEEE Trans. Veh. Technol._ , vol. 70, no. 5, pp. 4389–4405, May 2021.


# Machine Learning Based Forward Solver: An Automatic Framework in gprMax

Utsav Akhaury1, Iraklis Giannakis2, Craig Warren3, and Antonios Giannopoulos4

<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>1Laboratoire
d’astrophysique, École Polytechnique Fédérale de Lausanne (EPFL), Switzerland
2School of Geosciences, University of Aberdeen, Aberdeen, United Kingdom
3Department of Mechanical and Construction Engineering, Northumbria
University, Newcastle, United Kingdom  4School of Engineering, The University
of Edinburgh, Edinburgh, United Kingdom

###### Abstract

General full-wave electromagnetic solvers, such as those utilizing the finite-
difference time-domain (FDTD) method, are computationally demanding for
simulating practical GPR problems. We explore the performance of a near-real-
time, forward modeling approach for GPR that is based on a machine learning
(ML) architecture. To ease the process, we have developed a framework that is
capable of generating these ML-based forward solvers automatically. The
framework uses an innovative training method that combines a predictive
dimensionality reduction technique and a large data set of modeled GPR
responses from our FDTD simulation software, gprMax. The forward solver is
parameterized for a specific GPR application, but the framework can be
extended in a straightforward manner to different electromagnetic problems.

###### Index Terms:

Full-Waveform Inversion (FWI), Machine Learning (ML), Principle Component
Analysis (PCA), Singular Value Decomposition (SVD), Random Forest, XGBoost
(Extreme Gradient Boosting)

##  I Introduction

Recent progress in the field of GPR simulations [1] enabled research on
developing advanced modelling tools based on big data and machine learning
(ML). In this paper, we explore the concept of a near-real-time, forward
modeling approach for GPR that is based on ML architectures. The initial idea
for our ML approach was developed in [2] and further enhanced in [3]. Our
approach greatly eases the ML-based framework generation by introducing the
Random Parameter Generation feature in gprMax [4], our bespoke FDTD simulation
software. Using this feature, we generate a large number of GPR models with
randomly varying parameters. Subsequently, we reduce the dimensionality of
these data similar to [3], and finally, we compare the performance of
different ML regressors that take as inputs these random parameters and the
compressed A-Scans of the models for training. The ultimate goal is to predict
the output responses of GPR models depending on the selected parameters used
as inputs by the user. The whole process is automatic, and although the
training is time-consuming and computationally demanding, the final output is
a near real-time forward solver capable of predicting the resulting A-Scan
subject to the given inputs.

A fast numerical solver can greatly accelerate full-waveform inversion (FWI),
since the most computationally intensive module of the latter is the emended
numerical modelling in each iteration. This would make FWI commercially
appealing and attainable using minimum computational resources without the
need for high performance computing.

##  II Dataset Generation

###  II-A Random Parameter Generation - New Module in gprMax

This new feature facilitates the generation of random parameters for a
specific gprMax model and uses the _numpy.random_ module from NumPy [5] as
backend. It allows the user to specify the Probability Distribution Function
(PDF) from which these random parameters are drawn. Subsequently, for every
numerical parameter defining the GPR scenario, the user must enter two values
(in pairs) that define the specified PDF’s parameters. The following
convention is followed to activate the random parameter generation mode:

  * •

_#command_name: distr parameter_1.1 parameter_1.2 parameter_2.1 parameter_2.2
parameter_3.1 parameter_3.2 …_

_distr_ specifies the PDF from which random numbers are drawn. For example, to
specify a uniform distribution, the user must enter _u_ , and the two
subsequent values (_parameter_x.1 & parameter_x.2_) entered would define the
lower & upper bounds for the $x^{th}$ parameter. Note that if _parameter_x.1 =
parameter_x.2_ , then the random number generation is skipped and the constant
value is used.

As a quick demonstration, the following command creates a material called
_my_sand_ which has a relative permittivity $\epsilon_{r}$ drawn from a
uniform distribution within the range $[2,5]$, a conductivity of
$\sigma=0.01S/m$, $\mu_{r}=1$ and $\sigma_{*}=0$.

  * •

_#material: u 2 5 0.01 0.01 1 1 0 0 my_sand_

In case the generated random parameter exceeds the model domain bounds, it is
automatically constrained to fit inside the domain, which ensures that the
FDTD execution is not stopped midway. Finally, all the randomly generated
parameters for every GPR model generated are saved to a pickle file. We
automatically compress the pickle file by removing redundant features (i.e.
those features that have a constant value for every GPR model generated) in
order to improve the training efficiency.

###  II-B Metal Cylinder buried in a Dielectric Medium

Using the Random Parameter Generation feature detailed above, we generated a
dataset of 6250 2D GPR models for a metal cylinder buried in a dielectric
half-space. We then randomly split the dataset into train-test subsets - 5000
models for training and 1250 for testing. The geometry of the simple 2D
scenario is straight-forward and an image from the geometry view is shown in
Fig. 1. The transparent region around the boundary of the domain depicts the
PML (Perfectly matched layer) region. The line source used for the model’s
excitation is polarized along the $z$ direction, and is excited using a Ricker
waveform with amplitude of one and centre frequency of 1.5 GHz. Table I
summarizes the parameters that remain constant for every model.

In each of these models, we vary the following defining parameters, which are
drawn from a Uniform Distribution.

  * •

The cylinder radius

  * •

The depth at which the cylinder is placed in the dielectric medium

  * •

The electromagnetic permittivity ($\epsilon_{r}$) of the dielectric half-space

These are summarized in Table II.

The electromagnetic response (output A-Scans) of all these models are
simulated using the gprMax software and saved along with all the randomly
generated parameters. These data-pairs would eventually be fed to the ML
scheme. Note that since the line source is aligned along the $z$ direction,
only the $E_{z}$ component of the electric field is non-zero. And hence, we
only use the $E_{z}$ component for training.

TABLE I: Constant Model Parameters Parameter | Values
---|---
Domain Bounds $(x,y,z)$ | 0-240 mm, 0-210 mm, 0-2 mm
Dielectric half-space bounds $(x,y,z)$ | 0-240 mm, 0-170 mm, 0-2 mm
Spatial Discretization $(dx,dy,dz)$ | 2 mm, 2 mm, 2 mm
Time Window | 3 ns
Source Location $(x,y,z)$ | 100 mm, 170 mm, 0 mm
Receiver Location $(x,y,z)$ | 140 mm, 170 mm, 0 mm
$\sigma$ \- dielectric half-space | 0
$\sigma_{*}$ \- dielectric half-space | 0
$\mu_{r}$ \- dielectric half-space | 1
TABLE II: Randomly Varying Parameters Parameter | Lower Bound | Upper Bound
---|---|---
Cylinder Radius (r) | 5 mm | 20 mm
Cylinder Depth | 20 mm | 140 mm
$\epsilon_{r}$ \- dielectric half-space | 4 | 8

Figure 1: Geometry of a sample GPR model generated.

##  III Dimensionality Reduction

To increase the efficiency of our ML-based forward solver, we pre-process our
data in a suitable manner. To that end, we reduce the dimensionality of the
A-Scans before feeding them to the ML model. It is observed that the A-Scans
can effectively be represented using much lesser number of components without
any significant loss in signal quality. We compared the performance of two
popular algorithms, namely Principle Component Analysis (PCA) and Singular
Value Decomposition (SVD), which give us very similar results in terms of
degree of compression.

In either case, we start by using 10 components and computing the Normalized
Mean Squared Error (NMSE) between the original and reconstructed signals. We
then iteratively increase the number of components used for representing the
signal until the difference between successive NMSEs falls below a certain
threshold ($10^{-12}$). We noticed that most of the signal information is
captured by the first 10 to 15 components. Obviously, adding more components
would further minimize the NMSE. However, even with around 30 components, we
get an NMSE $\approx 10^{-12}$ on 1250 test samples, which is sufficiently low
for our purpose.

###  III-A Principle Component Analysis (PCA)

Principal component analysis (PCA) is among the oldest methods for reducing
the dimensionality of datasets [6], and it does so by increasing
interpretability while minimizing information loss. It performs this operation
by creating new uncorrelated variables that consecutively maximize variance.
These new variables are called the principal components. Finding these
principle components reduces to solving an eigenvalue/eigenvector problem, and
the new variables are defined by the dataset at hand, not a priori, thus
rendering PCA an adaptive data analysis technique. PCA has already been
successfully applied to GPR data for clutter reduction in [7] \- [8], and more
precisely for dimensionality reduction in [9].

Fig. 2 shows a few comparison plots between the A-Scans reconstructed from
their compressed representations (by applying inverse PCA transform) and the
original A-Scans.

Figure 2: Dimensionality Reduction using PCA. The dimension is reduced from
637 to 30, and yet the reconstruction is in excellent agreement with the
original signal. Notice that most of the signal information is captured by the
first 10 components, which are much larger in magnitude compared to the ones
that follow.

###  III-B Singular Value Decomposition (SVD)

Singular value decomposition (SVD) is a factorization method of a real or
complex matrix that extends the concept of the eigen-decomposition of a square
normal matrix with an orthonormal eigenbasis to any $m\times n$ matrix [10].
Truncated SVD is a quicker and more compact type of SVD in which the data
matrix is truncated and the rest of it is discarded. This can be much quicker
and more economical than simple SVD.

For our work, we use scikit-learn’s TruncatedSVD transformer. As opposed to
PCA, this estimator does not center the data before performing SVD, which
implies that it can efficiently handle sparse matrices. Fig. 3 shows a few
comparison plots between the A-Scans reconstructed from their compressed
representations (by applying inverse SVD transform) and the original A-Scans.

Figure 3: Dimensionality Reduction using TruncatedSVD. The dimension of the
A-Scans is reduced from 637 to 30.

##  IV Machine Learning based Forward Solver

Once we have the compressed representations of the output responses of our GPR
models, we feed them to our ML-based model for training. Supervised machine
learning will unravel the physical causal relationship between the inputs and
the resulting trace, and effectively approximate the underlying physical laws
of the given problem. To that end, we compare the performance of a few popular
ML regressors, such as Random Forest [11] and XGBoost [12], that we extend to
perform multi-variate multi-output regression. We had also evaluated the
performance of other regressors such as Support Vector Machines (SVM) [13],
etc. However, they do not even come close to the accuracy achieved by Random
Forest or XGBoost, and hence, we only focus on these two algorithms in our
work.

The ML solver outputs the compressed form of the output A-Scan, which can be
reconstructed by inverting the same dimensionality reduction algorithm that
was used to compress the training data.

###  IV-A Random Forest

Random Forest [11] is a supervised machine learning algorithm that is
constructed using decision trees [14] and invokes the concept of ensemble
learning [15], a technique that combines predictions from multiple ML
algorithms to make a more accurate prediction compared to a single model. As
an improvement to the decision tree algorithm, it reduces overfitting to
datasets and increases precision. It is a bagging-based algorithm in which
only a subset of features is selected at random. The maximum tree depth, which
is a tunable input parameter, controls the tendancy of the model to overfit to
the input data. In case of regression, Random Forest returns the mean or
average prediction of the individual trees. It can handle large datasets with
ease but does not perform well on very sparse data [16]. The idea of using
Random Forest regression for GPR-based predictions has been introduced in
[17].

We evaluate the performance of our trained Random Forest model on our test
dataset. A few output predictions using Random Forest are shown in Fig. 4.

Figure 4: Output predictions for Random Forest.

###  IV-B XGBoost

XGBoost [12], or Extreme Gradient Boosting, is also a decision-tree-based
ensemble Machine Learning algorithm that uses a gradient boosting framework.
It builds upon the concept of gradient boosting algorithm and provides a more
efficient implementation by invoking parallelization, tree pruning, and
hardware optimization. For most regression problems, XGBoost provably has the
best combination of prediction performance and processing time compared to
other algorithms. A few output predictions using XGBoost are shown in Fig. 5.

Figure 5: Output predictions for XGBoost.

###  IV-C Chain Regression

A Chain Regressor is a multi-label model that arranges regressions into a
chain. Every model makes a prediction within the order determined by the chain
using all of the available features provided to the model and the predictions
of models that precede in the chain. Chain Regressors are suitable where
multi-output regression is required and can improve the training performance,
as shown in [18]. However, they also take significantly longer to train.

We test the performance of Random Forest and XGBoost arranged in regression
chains. However, it is observed that there is no visible improvement when the
two cases are applied to our test dataset. On the contrary, it leads to a
degradation in performance. A few output predictions using Chain Regression
for XGBoost are shown in Fig. 6.

Figure 6: Output predictions for XGBoost with Chain Regression.

##  V Results and Comparison

Table III shows a comparison of the performance of the various ML schemes we
have tested. We applied all the methods discussed above to our test dataset of
1250 GPR models. All algorithms were tested on a 2.3 GHz Quad-Core Intel Core
i7 processor. For a quantitative comparison, we computed the Normalized Mean
Squared Error (NMSE) between the predicted and the simulated ground-truth
output responses (that were obtained using gprMax). Note that a lower value of
NMSE implies better performance. Alongside, we also noted the time taken by
the algorithm for training.

Random Forest achieves the least NMSE on our test dataset, followed by
XGBoost. As mentioned earlier, introducing Chain Regression does not show any
improvement compared to the original multi-output regressors. Furthermore, it
significantly consumes more time during training.

TABLE III: Performance Comparison Method | NMSE | Training Time
---|---|---
Random Forest | 0.0182 | 2.8s
Random Forest + Chain Regression | 0.0809 | 179.5s
XGBoost | 0.0285 | 8.8s
XGBoost + Chain Regression | 0.1011 | 17.4s

##  VI Conclusion

We successfully implemented a framework for generating ML-based forward
solvers in gprMax, our FDTD simulation software. This was greatly eased by
introducing the Random Parameter Generation feature. We observed that for our
test dataset, Random Forest gives us the best performance, taking into
consideration both, the NMSE between the ground truth & output predictions,
and the training time. Ultimately, gprMax gives the user the flexibility to
train & test the performance of any other suitable ML model to better fit a
specific dataset. Furthermore, as shown in [3], the performance of neural
networks [19] could also be investigated on our dataset in detail, which
remains the future scope of this project. The ML solver can also easily handle
3D geometries, in which case, there would only be a greater number of
parameters. The underlying principle remains the same.

## Acknowledgment

The project was funded via the Google Summer of Code (GSoC) 2021 programme.
GSoC initiative is a global program focused on bringing student developers
into open source software development. The source code for this project can be
found at https://github.com/gprMax/gprMax/pull/294.

## References

  * [1] C. Warren et al., “A CUDA-based GPU engine for gprMax: Open source FDTD electromagnetic simulation software,” Comput. Phys. Commun., to be published. [Online]. Available: https://www.sciencedirect. com/science/article/pii/S0010465518303990.
  * [2] I. Giannakis, A. Giannopoulos, and C. Warren, “A machine learning approach for simulating ground penetrating radar,” in Proc. 17th Int. Conf. Ground Penetrating Radar, Rapperswil-Jona, Switzerland, 2018, pp. 1–4.
  * [3] I. Giannakis, A. Giannopoulos, and C. Warren, ”A Machine Learning Based Fast Forward Solver for Ground Penetrating Radar With Application to Full Waveform Inversion,” IEEE Transactions on Geoscience and Remote Sensing. 2019 ; Vol. 57, No. 7. pp. 4417-4426.
  * [4] C. Warren, A. Giannopoulos, and I. Giannakis, “gprMax: Open source software to simulate electromagnetic wave propagation for ground penetrating radar,” Comput. Phys. Commun., vol. 209, pp. 163–170, Dec. 2016.
  * [5] C.R. Harris, K.J. Millman, S.J. van der Walt, et al., ”Array programming with NumPy”, Nature 585, 357–362, 2020.
  * [6] Y. Bengio, A. Courville, and P. Vincent, “Representation learning: A review and new perspectives,” IEEE Trans. Pattern Anal. Mach. Intell., vol. 35, no. 8, pp. 1798–1828, Aug. 2013.
  * [7] G. B. Kaplan, O. Içoglu, A. B. Yoldemir and M. Sezgin, “Real-time object detection using dynamic principal component analysis,” in Proc. 13th Int. Conf. Ground Penetrating Radar, Lecce, Italy, Jun. 2010, pp. 1–6.
  * [8] E. Tebchrany, F. Sagnard, V. Baltazart, J. Tarel and X. Dérobert, “Assessment of statistical-based clutter reduction techniques on ground-coupled GPR data for the detection of buried objects in soils,” in Proc. 15th Int. Conf. Ground Penetrating Radar, Brussels, Belgium, Jun./Jul. 2014, pp. 604–609.
  * [9] P. Kaczmarek and J. Pietrasinski, “Principal component analysis in interpretation of A-Scan measurements in GPR system,” in Proc. 15th Int. Radar Symp. (IRS), Gdansk, Poland, Aug. 2014, pp. 1–5.
  * [10] V. Klema and A. Laub, ”The singular value decomposition: Its computation and some applications,” in IEEE Transactions on Automatic Control, vol. 25, no. 2, pp. 164-176, April 1980.
  * [11] L. Breiman, ”Random Forests”, Machine Learning 45, 5–32, 2001.
  * [12] T. Chen and C. Guestrin, ”XGBoost: A scalable tree boosting system”, Proc. 22nd ACM SIGKDD Int. Conf. Knowl. Discovery Data Mining, pp. 785-794, 2016.
  * [13] M. A. Hearst, S. T. Dumais, E. Osuna, J. Platt and B. Scholkopf, ”Support vector machines,” in IEEE Intelligent Systems and their Applications, vol. 13, no. 4, pp. 18-28, July-Aug. 1998.
  * [14] A. Navada, A. N. Ansari, S. Patil and B. A. Sonkamble, ”Overview of use of decision tree algorithms in machine learning,” 2011 IEEE Control and System Graduate Research Colloquium, 2011, pp. 37-42.
  * [15] F. Huang, G. Xie and R. Xiao, ”Research on Ensemble Learning,” 2009 International Conference on Artificial Intelligence and Computational Intelligence, 2009, pp. 249-252.
  * [16] C. Tang, D. Garreau, and U. V. Luxburg, ”When do random forests fail?”, In Proceedings of the 32nd International Conference on Neural Information Processing Systems (NIPS’18). Curran Associates Inc., Red Hook, NY, USA, 2987–2997.
  * [17] I. Giannakis, A. Giannopoulos and C. Warren, ”A Machine Learning Scheme for Estimating the Diameter of Reinforcing Bars Using Ground Penetrating Radar,” in IEEE Geoscience and Remote Sensing Letters, vol. 18, no. 3, pp. 461-465, March 2021.
  * [18] H. Borchani, G. Varando, C. Bielza, P. Larrañaga, ”A survey on multi-output regression”, WIREs Data Mining Knowl Discov, 5: 216-233, 2015.
  * [19] C. M. Bishop and G. Hinton, Neural Networks for Pattern Recognition. London, U.K.: Oxford Univ. Press, 1996.