Iterative Detection for Orthogonal Time Frequency Space Modulation with Unitary Approximate Message Passing

Recently the orthogonal time frequency space (OTFS) modulation has attracted much attention due to its capability of achieving reliable communications for high mobility applications [1, 2, 3, 4, 5]. OTFS offers both time and frequency diversity as each symbol is spread over the time and frequency domains through the two dimensional (2D) inverse symplectic finite Fourier transform (ISFFT) [1, 2]. Compared with orthogonal frequency division multiplexing (OFDM), OTFS can achieve significant performance gains in high mobility scenarios [6]. In addition, when the number of channel paths is small, the effective channel in the delay-Doppler (DD) domain is sparse, which allows efficient channel estimation and data detection using message passing techniques [2].

A practical and powerful detector is essential to harvest the full time and frequency diversity promised by OTFS. The optimal maximum a posteriori (MAP) detector is impractical due to its complexity growing exponentially with the length of the OTFS block. In [7], an effective channel matrix in the DD domain was derived, based on which a low-complexity two-stage detector was proposed. The authors in [8] proposed a detection scheme, which uses minimum mean squared error (MMSE) equalization in the first iteration, followed by parallel interference cancellation with a soft-output sphere decoder in subsequent iterations. A low complexity iterative rake decision feedback detector was proposed in [9], [10], which extracts and coherently combines the multiple copies of the symbols (due to multipath propagation) in the DD grid using maximal ratio combining (MRC). Multiple modified MRC detectors with enhanced performance were also developed in [10]. They performs similarly or even better than the message passing based detector [2] while with a linear complexity. The design of low-complexity detectors were also investigated based on the message passing techniques [2], [11], [12]. A variational Bayes (VB) based detector was proposed in [13] to achieve better convergence compared with the existing message passing based detectors. The detectors in [2], [13], and [14] take advantage of the sparsity of the channel matrix in the DD domain, and their complexity depends on the number of nonzero elements in each row of the channel matrix, which is denoted by $S$. Without considering fractional Doppler shifts, $S$ is equal to the number of channel paths. In general, a wideband system is able to provide sufficient delay resolution, however, the Doppler resolution depends on the time duration of the OTFS block. To fulfill the low latency requirement in future wireless communications, the time duration should be relatively small, where it is necessary to consider fractional Doppler shifts [2, 11]. In this case, the value of $S$ is proportional to the number of paths and can be significantly larger than the number of channel paths. This leads to two issues: the complexity of the existing detectors is a concern especially in the case of rich scattering environments, and the message passing based detectors may suffer from performance loss due to the presence of short loops in the corresponding system graph model.

This work aims to address the above issues by designing efficient detectors based on the approximate message passing (AMP) algorithm [15], [16]. The use of AMP in this work is due to two facts: AMP was developed based on loopy belief propagation in dense factor graphs, i.e., short loops can be handled by AMP, and AMP has a low complexity. It is known that AMP works well for i.i.d. (sub-)Gaussian system transfer matrix, but it may suffer from performance loss or even diverge for a general system transfer matrix [17]. Inspired by [18], the work in [19] showed that AMP can still work well for a general system transfer matrix when a unitary transform of the original model is used, where the unitary matrix for transformation can be the conjugate transpose of the left singular matrix of the general system transfer matrix [19, 20] obtained through singular value decomposition (SVD). This variant to AMP is called unitary AMP (UAMP), which is formerly known as UTAMP. UAMP has been used for low complexity robust sparse Bayesian learning [21], bilinear recovery [20], inverse synthetic aperture radar (ISAR) imaging with high Doppler resolution [22], direction of arrival (DOA) estimation [23], etc. It will be shown in this paper that UAMP is well suitable for OTFS because the channel matrix in the DD domain is a block circulant matrix with circulant block (BCCB), and the 2D discrete Fourier transform can be used for the unitary transformation, thereby allowing an efficient implementation using the 2D fast Fourier transform (FFT) algorithm. This leads to an attractive UAMP-based detector in two aspects. First, the complexity of the UAMP-based detector is in the order of the logarithm of the OTFS block length per symbol per iteration, which is independent of $S$. Second, thanks to the robustness of UAMP, the UAMP-based detector delivers significantly better performance. In addition, as the noise variance is normally unknown in practice, the noise variance estimation is incorporated in the design of the UAMP-based detector in this work, so that an extra noise variance estimator is not required (existing OTFS detectors except the one in [9] often assume a perfect knowledge of the noise variance). It is shown that the UAMP-based detector can significantly outperform the AMP-based detector and other state-of-the-art detectors [2], [9, 10], [13]. In this paper, we also extend our investigations to coded OTFS systems, and integrate (U)AMP-based detectors into a powerful turbo receiver to achieve joint iterative detection and decoding.

The remainder of the paper is organized as follows. In Section II, the OTFS system model is introduced. We design (U)AMP-based detectors with bi-orthogonal waveform in Section III, and rectangular waveform in Section IV. The discussion is extended to iterative joint detection and decoding in a coded OFTS system and the performance prediction of the OTFS system with the proposed detectors is discussed in Section V. Simulations results are provided in Section VI, followed by conclusions in Section VII.

Notations- Boldface lower-case and upper-case letters denote vectors and matrices, respectively. The superscripts $(\cdot)^{H}$ and $(\cdot)^{T}$ represent conjugate transpose and transpose operations, respectively. The superscript $*$ is used to denote the conjugate operation. We define $[\cdot]_{M}$ as the mod $M$ operation. The probability density function of a complex Gaussian variable with mean $\hat{x}$ and variance $\nu_{x}$ is represented by $\mathcal{N}(x|{\hat{x}},\nu_{x})$. The notation $\left\langle f(\mathbf{x})\right\rangle_{q(\mathbf{x})}$ denotes the expectation of the function $f(\mathbf{x})$ with respect to probability density $q(\mathbf{x})$. The relation $f(x)=cg(x)$ for some positive constant $c$ is written as $f(x)\propto g(x)$. Notation $\otimes$ represents the Kronecker product, and $\bm{a}\cdot\bm{b}$ and $\bm{a}\cdot/\bm{b}$ represent the component-wise product and division between vectors $\bm{a}$ and $\bm{b}$, respectively. The notation $\bm{X}=\text{reshape}_{M}\left(\bm{x}\right)$ represents that the vector $\bm{x}$ is reshaped as an $M\times N$ matrix $\bm{X}$ column by column, where the length of vector $\bm{x}$ is $MN$. We use $\bm{x}=\text{vec}\left(\bm{X}\right)$ to represent vectorization of matrix $\bm{X}$ column by column. The notation $\text{diag}(\bm{a})$ represents a diagonal matrix with the elements of $\bm{a}$ as its diagonal. We use $|\bm{A}|^{2}$ to denote the element-wise magnitude squared operation for matrix $\bm{A}$. We use 1 and 0 to denote an all-ones vector and an all-zeros vector with a proper length, respectively. We use $q_{j}$ to denoted the $j$th entry of $\bm{q}$. The superscript $t$ of $\textbf{s}^{t}$ denotes the iteration index in an iterative algorithm.

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

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

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

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

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

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

where $l_{i}$ and $k_{i}$ are the delay index and Doppler index of the $i$th path, and $\kappa_{{i}}$ $\in[-\frac{1}{2},\frac{1}{2}]$ is a fractional Doppler associated with the $i$th path. Here, $M\Delta f$ is the total bandwidth of the system and $NT$ is the duration of an OTFS block.

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

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

The rectangular waveform has been used in the literature as a more practical waveform. When a length-$T$ rectangular waveform is used as the transmitter and receiver filters, as in [29], [30, 31], we append a cyclic prefix (CP) to each length-$M$ sub-block of the signal $\bm{s}$ in the time domain before transmitting over the time-varying channel.

After removing CPs at the receiver side, the received signal in the time domain can be expressed as [29, 30, 31]

where $\bm{\omega}$ is an additive Gaussian white noise vector, and $\bm{H}_{T}$ is an $NM\times NM$ block diagonal matrix

where $\bm{H}_{n}\in\mathbb{C}^{M\times M}$ is the channel matrix corresponding to the $n$-th sub-block of the transmitted signal. Then the received signal in the DD domain is obtained by performing block-wise DFTs, followed by the SFFT operation, i.e.,

Finally, the channel input-output relationship in the DD domain can be expressed as

where $\bm{H}$ represents the effective channel matrix in the DD domain.

Based on model (35), we can directly apply Algorithm 1 to detect $\bm{x}$, leading to an AMP-based OTFS detector. However, Algorithm 2 cannot be directly applied. Next, we design a detector based on the UAMP algorithm.

Since the channel matrix $\bm{H}_{T}$ is a block diagonal matrix, we can perform SVD block-by-block, i.e.,

where

and each $\bm{\Lambda}_{n}$ is an $M\times M$ diagonal matrix (so does $\bm{\Lambda}$ itself). By applying the unitary transformation with the unitary matrix $\bm{U}^{H}(\bm{F}_{N}^{H}\otimes\bm{I}_{M})$ to model (35), we have

where $\bm{r}=\bm{U}^{H}(\bm{F}_{N}^{H}\otimes\bm{I}_{M})\bm{y}\in\mathbb{C}^{MN\times 1}$, $\bm{\Phi}=\bm{V}(\bm{F}_{N}^{H}\otimes\bm{I}_{M})$ and $\bar{\bm{\omega}}=\bm{U}^{H}(\bm{F}_{N}^{H}\otimes\bm{I}_{M})\bm{\omega}$ is still a zero-mean Gaussian noise vector with the same covariance matrix as $\bm{\omega}$ because $\bm{U}^{H}(\bm{F}_{N}^{H}\otimes\bm{I}_{M})$ is a unitary matrix. So, UAMP can be applied to model (40), leading to an algorithm similar to Algorithm 2 with a few equations modified, as detailed in the following:

• a.

  Vector $\bm{d}$ consists of the diagonal elements of matrix $\bm{\Lambda}$.

• b.

  The computation of $\bm{p}$ in Line 2 of Algorithm 2 is replaced by

  	$\displaystyle\bm{p}$	$\displaystyle=$	$\displaystyle\bm{\Lambda}\bm{V}(\bm{F}_{N}^{H}\otimes\bm{I}_{M})\hat{\bm{x}}-\bm{\nu}_{p}\cdot\bm{s}^{t-1}$		(41)
  		$\displaystyle=$	$\displaystyle\bm{d}\cdot\bm{V}\text{vec}(\hat{\bm{X}}\bm{F}_{N}^{H})-\bm{\nu}_{p}\cdot\bm{s}^{t-1},$

  where $\hat{\bm{X}}\triangleq\text{reshape}_{M}(\hat{\bm{x}})$, and $\hat{\bm{X}}\bm{F}_{N}^{H}$ can be computed using the FFT algorithm.

• c.

  The computation of $\bm{q}$ in Line 9 of Algorithm 2 is replaced by

  	$\displaystyle\bm{q}$	$\displaystyle=$	$\displaystyle\hat{\bm{x}}^{t}+\nu_{q}(\bm{\Lambda}\bm{\Phi})^{H}\bm{s}$		(42)
  		$\displaystyle=$	$\displaystyle\hat{\bm{x}}^{t}+\nu_{q}(\bm{F}_{N}\otimes\bm{I}_{M})\bm{V}^{H}\bm{\Lambda}^{H}\bm{s}$
  		$\displaystyle=$	$\displaystyle\hat{\bm{x}}^{t}+\nu_{q}\text{vec}(\bm{S}\bm{F}_{N}),$

  where $\bm{S}\triangleq\text{reshape}_{M}(\bm{V}^{H}(\bm{d}^{*}\cdot\bm{s}))$, and $\bm{S}\bm{F}_{N}$ can be computed using the FFT algorithm.

The UAMP-based detector with rectangular waveform only involves element-wise product (division) except (41) and (42). Note that $\bm{V}$ is a block diagonal matrix, so the computational complexity of the detector is $\mathcal{O}\left(M^{2}N\right)+\mathcal{O}(MN|\mathcal{A}|)$ per OTFS block. which is independent of $S$. Compared to the MRC detector in [9, 10], the UAMP-based detector has a higher complexity, but as shown in Section VI, it can achieve significantly better performance, especially in the case of high order modulations.

In this section, we consider a coded OFTS system and investigate the design of a turbo receiver for joint detection and decoding.

The turbo system is shown in Fig. 3, where the information bits are firstly encoded and then interleaved before mapping. Each symbol $x_{j}\in\mathcal{A}=\{\alpha_{1},...,\alpha_{|\mathcal{A}|}\}$ in the DD domain is mapped from a sub-sequence of the coded bit sequence, which is denoted by $\bm{c}_{j}=[c_{j}^{1},...,c_{j}^{log|\mathcal{A}|}]$. Each $\alpha_{a}$ corresponds to a length-$\mathrm{log}|\mathcal{A}|$ binary sequence denoted by $\{\alpha_{a}^{1},...,\alpha_{a}^{log|\mathcal{A}|}\}$. The turbo receiver consists of a soft-in-soft-out (SISO) detector and a SISO decoder, which exchange information in an iterative manner. Based on the log-likelihood ratios (LLRs) provided the SISO decoder and the output of the OTFS demodulator as shown in Fig. 3, the task of the SISO detector is to compute the extrinsic LLR for each coded bit, i.e.,

where $L^{a}(c_{j}^{q})$ is the ouput extrinsic LLR of the decoder in last iteration. The extrinsic LLR $L^{e}(c_{j}^{q})$ is passed to the decoder. The derivation for $L^{e}(c_{j}^{q})$ in terms of extrinsic mean and variance can be find in [32], and $L^{e}(c_{j}^{q})$ can be expressed as

where $m_{j}^{e}$ and $v_{j}^{e}$ are the extrinsic mean and variance of $x_{j}$, and $\mathcal{A}_{q}^{0}$ and $\mathcal{A}_{q}^{0}$ represent the subsets of all $\alpha_{a}$ corresponding to $c_{j}^{q}=0$ and $c_{j}^{q}=1$, respectively. The extrinsic variance and mean are defined as [32]

where $m_{j}$ and $v_{j}$ are the a priori mean and variance of $x_{j}$ calculated based on the output LLRs of the SISO decoder [33], [34], [35] and $m_{j}^{p}$ and $v_{j}^{p}$ are the a posteriori mean and variance of $x_{j}$. By examining the derivation of the (U)AMP algorithm, we can find that $\bm{q}$ and $\bm{\nu}_{q}$ consists of the extrinsic means and extrinsic variances of the symbols in vector $\bm{x}$ as they are the messages passed from the observation side and do not contain the immediate a priori information about $\bm{x}$. Therefore, we have

in Algorithm 1 and

in Algorithm 2. Then (44) can be readily used to compute the extrinsic LLRs of the coded bits.

It is noted that, with the LLRs provided by the SISO decoder, we can compute the probability $p(x_{j}=\alpha_{a})$ for each $x_{j}$, which is no longer $1/|\mathcal{A}|$ in Algorithm 1 and Algorithm 2. Therefore, $\xi_{j,a}$ in Line 7 of Algorithm 1 is amended to

and Line 10 of Algorithm 2 is changed to