Progressive Inter-path Interference Cancellation Algorithm for Channel Estimation Using OTFS

1. Introduction

Reliable communication in doubly dispersive channels, where high delay and Doppler spreads due to multipath and mobility are experienced, is becoming an increasingly prominent research topic toward the development of next-generation wireless networks [1,2]. The widely adopted OFDM (Orthogonal Frequency Division Multiplexing) is suitable for static frequency-selective multipath channels but experiences drastic performance degradation in high mobility scenarios due to the loss of orthogonality between the subcarriers as analyzed in [3], inducing Inter-Carrier Interference (ICI). To cater to high-mobility scenarios, several novel waveforms have been introduced in literature. An overview of some of the existing modulation schemes for high mobility wireless channels is reported in [4]. Orthogonal Time Frequency Space (OTFS) modulation was proposed in [5,6] and is garnering a lot of attention due to its robustness against Doppler effect. OTFS is a two-dimensional modulation scheme that multiplexes information symbols in the delay-Doppler domain, where the channel representation becomes sparse and the channel can be considered almost time-invariant [5,6,7]. Conversely, OFDM multiplexes information symbols in the time-frequency (TF) domain where mobility generates a time-varying channel frequency response.

1.1. Problem Context

In this work, we consider a scenario comprising a single antenna transmitter and receiver, focusing on the channel estimation problem with fractional delays and Doppler shifts and no prior information about the number of propagation paths [8,9]. The transmitter sends a pilot symbol in the delay-Doppler (DD) domain, while the receiver processes the received pilot frame to extrapolate the knowledge of the channel parameters.

The channel estimation problem is often considered in the integer DDs scenario, which is facilitated by the assumption of having a channel with delays and Doppler shifts that are multiples of delay and Doppler resolutions. In this case, there is no interference between adjacent DD bins, and a simple threshold method can be used to obtain good performances, as shown in [10]. On the other side, in reality this assumption may not hold, which motivates the study of the the fractional DDs scenario, which is more challenging, since the received replica of the pilot symbol in the DD domain spreads into adjacent bins after 2D-sampling of the DD domain [7,11]. This is the well-known Inter-Path Interference phenomenon. Our focus is on the effects of IPI caused by fractional DDs and IPI cancellation to recover channel state information (CSI) with high accuracy.

1.2. State of the Art

Channel estimation in OTFS systems has been widely investigated in the literature since the detection performance of the OTFS receiver is based on the estimated channel matrix [12,13,14,15]. In [10], an embedded-pilot aided method for channel estimation is proposed, achieving good performances in both integer DDs and fractional Doppler scenarios. Channel estimation in massive multiple input multiple output scenarios (MIMO) has been investigated in [16] where additional sparsity due to a large number of transmit and receive antennas allows the channel estimation problem to become a simpler sparse recovery problem taking advantage of the angular domain. Low-complexity channel estimation schemes were considered in single-antenna systems in [8], where Modified Maximum Likelihood Estimator (M-MLE) and Two-Step Estimator (TSE) algorithms were introduced. The first is based on the asymptotic orthogonality of a delay-Doppler parameter matrix for large OTFS frames, which allows for the simplification of the cost function. On the other hand, the second one has lower complexity and relies on disjoint delay and Doppler estimation. Under the assumption of having good separability of the received pilot replica due to fine delay and Doppler resolutions, in both algorithms the IPI is cancelled and at each iteration a residual received vector is used for estimating channel parameters of the next path. Similarly, IPI cancellation is performed in the low-complexity channel estimation scheme for Discrete Zak Transform based OTFS (DZT-OTFS) proposed in [17]. Sparse Bayesian learning (SBL) algorithms for channel estimation are also considered in [18,19], but M-MLE and TSE algorithms are shown to perform better then SBL in [8] when delay and Doppler resolutions are fine enough. In [20], different channel estimation schemes for the fractional DDs case are proposed with the scope of improving the spectrum efficiency and reduce the peak-to-average-power-ratio (PAPR) by sending pilots in the time-frequency (TF) domain. The problem of high precision channel estimation in IPI scenarios has been addressed in [11], where an iterative delay-Doppler IPI Cancellation (DDIPIC) algorithm was proposed. In this case, the IPI cancellation is performed indirectly through refinement procedures carried out at each iteration. These refinements of the channel parameters take time, increasing the computational burdening of the channel estimation algorithm. In more detail, the computational complexity of the cost function, which is evaluated multiple times during each iteration, increases with the cube of the iteration index due to matrix inversion.

1.3. Contributions

The goal of this work is to proposed a variant of the DDIPIC algorithm from [11], named the Progressive IPI Cancellation (P-IPIC) algorithm, that cancels the IPI by performing residue cost function (RCF) maximization at each step. This algorithm operates in two phases: A search phase, where a certain number of possible paths are detected and a first rough estimate of the channel parameters is provided, A refinement phase, where the found paths are refined, looking for ones declared as false alarms and discarded, until the stopping criterion is met. This procedure allows finding the correct number of paths while the accuracy of the estimates is improved. The contributions of this work can be summarized as follows:

• Progressive IPI cancellation: Instead of discovering possible paths and performing refinements in each iteration maximizing the observation-based cost function as in [11], we perform maximization of a residue-based cost function (RCF) at each step during the discovery phase. This progressive cancellation of the IPI allows providing a first channel estimate aimed to reduce the residue energy as iteration progress achieving high accuracy in channel matrix reconstruction
• Latency reduction due to single global refinement of the estimates: the main drawback of the state-of-the-art DDIPIC algorithm is that it performs refinement procedures at each iteration as discussed in Section 3.5. Progressive IPI cancellation allows to reduce the latency since, as discussed above, the algorithm can operate in two distinct phases and a single global refinement is needed
• Computational complexity reduction in search phase: during the iterative search phase, the proposed algorithm looks for new paths by performing a RCF maximization. This allows us to reduce the computational cost since the computation of the RCF does not require matrix inversion. Moreover, the computational cost does not depend on the iteration index. As a consequence, coarse and fine estimation procedures described in Section 3.4 will take less time
• Resilience towards high-IPI scenarios: in literature it is commonly assumed to have sufficiently fine delay and Doppler resolutions such that all the propagation paths are resolvable in the DD domain. In some cases this may not hold and robust low-complexity channel estimation algorithms must be designed to suite high-IPI scenarios. In this work we consider both low- and high-IPI regimes showing that the proposed approach is suitable for both regimes

The remainder of this paper is organized as follows. In Section 2, the OTFS signal model is presented. In Section 3, OTFS channel estimation is covered, including the baseline DDIPIC method and our proposed P-IPIC method.

2. OTFS Modulation

In this section, the employed models for the OTFS signal, the propagation channel, and receiver-side processing are presented.

2.1. Transmitted Signal Model

In OTFS systems [5,7], the TX arranges $MN$ symbols in the delay-Doppler grid referred to as $I_{DD}$, where M is the number of subcarriers and N is the number of time-slots of the OTFS frame. The DD-domain symbols are then converted to the time-frequency domain through inverse symplectic finite Fourier transform (ISFFT). Denoting with $\mathbf{X}\in {\mathbb{C}}^{M\times N}$ the delay-Doppler symbols multiplexed in the DD domain, the corresponding symbols in the time-frequency grid are obtained as

$$X_{tf}[l,k]=\frac{1}{\sqrt{MN}}\sum _{n=0}^{N-1}\sum _{m=0}^{M-1}X[m,n]e^{j2\pi \left(\frac{nk}{N}-\frac{ml}{M}\right)}.$$

(1)

After that, the time-domain OTFS signal is obtained through the the Heisenberg transform as

$$s\left(t\right)=\sum _{k=0}^{N-1}\sum _{l=0}^{M-1}{X}_{tf}[l,k]g_{tx}(t-kT)e^{j2\pi l\Delta f(t-kT)},$$

(2)

where $g_{tx}\left(t\right)$ is the transmit shaping pulse of duration T. The continuous-time signal in Equation (2) has a total duration of $NT$ and a bandwidth of $M\Delta f$ where T corresponds to the time-slot duration and $\Delta f=1/T$ is the subcarrier spacing. Thus, the delay and Doppler resolutions are, respectively,

$${\tau }_{\mathrm{res}}=\frac{1}{M\Delta f},{\nu }_{\mathrm{res}}=\frac{1}{NT}.$$

(3)

2.2. Channel Model

The transmitted signal passes through the channel that, in the DD-domain, is modeled as

$$h(\tau ,\nu )=\sum _{i=1}^P{g}_i\delta (\tau -{\tau }_i)\delta (\nu -{\nu }_i).$$

(4)

where P is the number of propagation paths and $g_i$, ${\tau }_i$ and ${\nu }_i$ are the complex channel gain, the propagation delay and the Doppler shift of the i-th path, respectively.

2.3. Receiver Processing

The received signal is obtained by computing the Heisenberg transform and adding Additive White Gaussian Noise (AWGN) with monolateral Power Spectral Density (PSD) $N_0$ referred to as $n\left(t\right)$. Thus, the received signal is

$$r\left(t\right)=\int \int h(\tau ,\nu )s(t-\tau )e^{j2\pi \nu (t-\tau )}d\tau d\nu +n\left(t\right),$$

(5)

At the receiver side, inverse operations are performed to recover the received symbols in the DD-domain. Firstly, the receiver computes the Wigner transform by sampling the so-called Cross-Ambiguity Function (CAF) $A_{g_{rx},r}(f,t)$ to convert the received time-domain signal into samples in the time-frequency grid as

$$Y_{tf}[l,k]=A_{g_{rx},r}{(f,t)|}_{f=l\Delta f,t=kT}=\int r\left(t^{\prime }\right)g_{rx}^{*}(t^{\prime }-t)e^{-j2\pi f(t^{\prime }-t)}d{t}^{\prime }{|}_{f=l\Delta f,t=kT}.$$

(6)

Secondly, a SFFT is applied to get the DD-domain received frame

$$Y[m,n]=\frac{1}{\sqrt{MN}}\sum _{k=0}^{N-1}\sum _{l=0}^{M-1}{Y}_{tf}[l,k]e^{-j2\pi \left(\frac{nk}{N}-\frac{ml}{M}\right)}.$$

(7)

From now on we are going to consider OTFS modulation with rectangular shaping filters of duration T and amplitude $1/\sqrt{T}$ at both transmitter and receiver side.

3. OTFS Channel Estimation

The purpose of channel estimation is to determine the parameters of (4), based on pilot signals. In this section, the pilot model and corresponding observation model are presented, followed by the baseline method and the proposed method.

3.1. Pilot Model

In our scenario [8,11], the single antenna transmitter sends a pilot symbol in the DD domain. The corresponding pilot frame is given by

$$X[m,n]=\left\{\begin{array}{cc}\sqrt{E_p}\hfill & m=m_p,n=n_p,\hfill \\ 0\hfill & \mathrm{otherwise},\hfill \end{array}\right.$$

(8)

where $m=m_p$, $n=n_p$ is the delay-Doppler Resource Element (DDRE) in which the pilot is sent and $E_p$ is the energy of the transmitted pilot signal.

At the receiver, after the operations described above and vectorization, the received pilot frame is given by

$$\mathbf{y}=\mathrm{vec}\left(\mathbf{Y}\right)=\mathbf{H}\mathrm{vec}\left(\mathbf{X}\right)+\mathbf{w}=\mathbf{H}\mathbf{x}+\mathbf{w},$$

(9)

where $\mathbf{w}\sim \mathcal{CN}(0,N_0\mathbf{I})\in {\mathbb{C}}^{MN\times 1}$ is the AWGN noise vector and $\mathbf{H}\in {\mathbb{C}}^{MN\times MN}$ is the delay-Doppler domain channel matrix given by

$$\mathbf{H}=\sum _{i=1}^P{g}_i{\mathbf{G}}_i({\tau }_i,{\nu }_i),$$

(10)

where ${\mathbf{G}}_i({\tau }_i,{\nu }_i)\in {\mathbb{C}}^{MN\times MN}$ is a matrix that captures the effect of the propagation delay and the Doppler shift of the i-th path and for $l^{\prime },l^{\prime \prime }=0,1,..,M-1$, $k^{\prime },k^{\prime \prime }=0,1,...,N-1$ is computed as [8,11]

$$G_i[k^{\prime }M+l^{\prime },k^{\prime \prime }M+l^{\prime \prime }]=\frac{e^{-j2\pi {\tau }_i{\nu }_i}}{MN}\sum _{n=0}^{N-1}\sum _{m=0}^{M-1}{r}_{k^{\prime \prime },l^{\prime }}\left(m\right)e^{j2\pi \left(\frac{m}{M}\left(l^{\prime }-l^{\prime \prime }-\frac{M{\tau }_i}{T}\right)-\frac{n}{N}\left(k^{\prime }-k^{\prime \prime }-\frac{N{\nu }_i}{\Delta f}\right)\right)},$$

(11)

where

$$\begin{array}{c}\hfill r_{k^{\prime \prime },l^{\prime }}\left(m\right)=\sum _{s=-m}^{M-1-m}{e}^{j2\pi s\frac{l^{\prime }}{M}}[\left(1-\frac{{\tau }_i}{T}\right)e^{j\pi \left(1+\frac{{\tau }_i}{T}\right)\left(\frac{{\nu }_i}{\Delta f}-s\right)}\mathrm{sinc}\left(\left(1-\frac{{\tau }_i}{T}\right)\left(\frac{{\nu }_i}{\Delta f}-s\right)\right)+\\ \hfill \frac{{\tau }_i}{T}{e}^{-j2\pi \frac{k^{\prime \prime }}{N}}{e}^{j\pi \left(\frac{{\tau }_i}{T}\right)\left(\frac{{\nu }_i}{\Delta f}-s\right)}\mathrm{sinc}\left({\tau }_i\left({\nu }_i-s\Delta f\right)\right)].\end{array}$$

(12)

If the channel coefficients are normalized such that ${\sum }_{i=1}^P{\left|g_i\right|}^2=1$ the received pilot signal power equals the transmitted one. Under this hypothesis, the Pilot Signal-to-Noise Ratio (PSNR) can be computed as in [8]. Hence, the pilot signal power is the total energy $E_p$ divided by the overall time duration $NT$, while the noise power is the product of the power spectral density $N_0$ and the bandwidth $M\Delta f$. Thus,

$$\mathrm{PSNR}=\frac{E_p/NT}{M\Delta f{N}_0}=\frac{E_p}{MN{N}_0}.$$

(13)

During a channel coherence window, the receiver processes the received pilot frame running a channel estimation algorithm and, once the estimated channel matrix $\widehat{\mathbf{H}}$ is obtained, uses the estimate to detect the incoming data frames.

3.2. Observation Model for Channel Estimation

As in [11], channel estimation is performed considering the input-output relation in Equation (9) rewritten as

$$\mathbf{y}=\mathbf{R}\left(\tau ,\nu \right)\mathbf{g}+\mathbf{w},$$

(14)

where $\mathbf{g}={[g_1{g}_2...g_P]}^T\in {\mathbb{C}}^{P\times 1}$ is the channel gains vector and $\mathbf{R}\left(\tau ,\nu \right)=[{\mathbf{r}}_1{\mathbf{r}}_2...{\mathbf{r}}_P]\in {\mathbb{C}}^{MN\times P}$ is a matrix called constituent delay-Doppler parameter matrix (CDDPM). Each column of $\mathbf{R}\left(\tau ,\nu \right)$ is related to a different path and is given by

$${\mathbf{r}}_i({\tau }_i,{\nu }_i)={\mathbf{G}}_i({\tau }_i,{\nu }_i)\mathbf{x}.$$

(15)

Given the input-output relation in Equation (14), the Maximum Likelihood (ML) estimate of the unknowns is given by

$$(\widehat{\mathbf{g}},\widehat{\tau },\widehat{\nu })=arg\underset{\mathbf{g},\tau ,\nu }{min}{\parallel \mathbf{y}-\mathbf{R}\left(\tau ,\nu \right)\mathbf{g}\parallel }^2.$$

(16)

Although the minimization problem depends on three parameters, it can be simplified performing a two-step ML estimate of the unknowns as

$$(\widehat{\tau },\widehat{\nu })=arg\underset{\tau ,\nu }{max}\Phi \left(\mathbf{R}\left(\tau ,\nu \right),\mathbf{y}\right),$$

(17)

and

$$\widehat{\mathbf{g}}={\left({\mathbf{R}}^H\left(\widehat{\tau },\widehat{\nu }\right)\mathbf{R}\left(\widehat{\tau },\widehat{\nu }\right)\right)}^{-1}{\mathbf{R}}^H\left(\widehat{\tau },\widehat{\nu }\right)\mathbf{y}.$$

(18)

Here, the cost function $\Phi \left(\mathbf{R}\left(\tau ,\nu \right),\mathbf{y}\right)$ is computed as

$$\Phi \left(\mathbf{R}\left(\tau ,\nu \right),\mathbf{y}\right)={\mathbf{y}}^H\mathbf{R}\left(\tau ,\nu \right){\left({\mathbf{R}}^H\left(\tau ,\nu \right)\mathbf{R}\left(\tau ,\nu \right)\right)}^{-1}{\mathbf{R}}^H\left(\tau ,\nu \right)\mathbf{y}.$$

(19)

The DDIPIC algorithm proposed in [11] performs an iterative maximization of the cost function in Equation (17) in order to jointly estimate delays and Dopplers of the paths. The operation of the algorithm is summarized below.

3.3. DDIPIC Algorithm

The algorithm proceeds iteratively for a maximum of $P_{max}$ iterations and the CDDPM matrix is initialized as $\mathbf{R}\left(\tau ,\nu \right)=\mathbf{0}$. The details of the DDIPIC algorithm are provided in Algorithm 1. At each iteration, the parameters of one path are estimated as follows:

• Coarse estimation (line 2): the cost function in (19) is maximized over the DD grid of integer multiples of delay and Doppler resolutions providing a coarse estimate of the delay and Doppler shift. The search area is limited to $[0,L{\tau }_{\mathrm{res}}]\times [-K{\nu }_{\mathrm{res}},K{\nu }_{\mathrm{res}}]$ where $L=\lceil \frac{{\tau }_{max}}{{\tau }_{\mathrm{res}}}\rceil$ and $K=\lceil \frac{{\nu }_{max}}{{\nu }_{\mathrm{res}}}\rceil$
• Fine estimation (lines 5–11): an iterative procedure runs for a maximum of $s_{max}$ iterations. In each iteration, the cost function in (19) is maximized over a search area comprising fractional delays and Dopplers around, for $s=1$, the coarse estimate obtained in step 1 or, for $s>1$, the fine estimate obtained in the previous iteration. The search space is narrowed down as the iteration progresses and the fine estimation procedure ends when the stopping criteria for $\tau$ and $\nu$ are met. The search space for the fine estimation procedure is called $F_{DD}^{\left(s\right)}$ and we define $\Gamma =\{-\lfloor \frac{m_{\tau }}{2}\rfloor ,...,0,...,\lfloor \frac{m_{\tau }}{2}\rfloor \}$ and $\Lambda =\{-\lfloor \frac{n_{\nu }}{2}\rfloor ,...,0,...,\lfloor \frac{n_{\nu }}{2}\rfloor \}$ where $m_{\tau }$ and $n_{\nu }$ are design parameters. Notice that after the Cartesian product (denoted ×) of the two sets, only positive delays need to be considered. The parameters ${ϵ}_{\tau }$ and ${ϵ}_{\nu }$ are chosen to achieve the desired accuracy compatibly with $s_{max}$.
• Check stopping criterion (line 15-19): a residue vector is computed as ${\mathcal{E}}^{\left(i\right)}=\mathbf{y}-\mathbf{R}(\widehat{\tau },\widehat{\nu })\widehat{\mathbf{g}}$ and if $\parallel {\mathcal{E}}^{\left(i\right)}\parallel <ϵ$ is satisfied the algorithm stops and a number of $\widehat{P}$ paths are detected, otherwise, before starting the search for a new path, all previous estimates are refined as described below. The convergence tolerance parameter $ϵ$ for the stopping criterion has to be chosen in order not to underestimate the number of propagation paths, not to miss the ones with smaller amplitudes, and not to overestimate the number of propagation paths due to residual IPI or noise. A reasonable value is $ϵ=3\sqrt{MN{N}_0}$.
• Refinement procedure (line 20-21): all coarse and fine estimates are performed again for all the previously detected paths until the current one. In particular, at the i-th iteration lines 2–13 are run again iteratively for $i_{\mathrm{ref}}=1,...,i$.

3.4. Proposed P-IPIC Algorithm

The proposed variant of the DDIPIC algorithm aims to reduce the complexity of DDIPIC and is based on a RCF maximization and performs an iterative search as in the original algorithm until the maximum number of iterations is reached or the stopping criterion is met. The pseudocode of the proposed P-IPIC algorithm is provided in Algorithm 2. The main parameters are as for the DDIPIC algorithm.

The proposed algorithm features several key differences compared to Algorithm 1:

• Coarse and fine estimation based on the residuals (lines 3 and 10): The key idea is that, once a path is detected and the residue vector is obtained, during the next search, a new path is discovered looking at the residual pilot signal of the previous iteration. Thus, at the i-th iteration, the algorithm looks for a new path maximizing an RCF given by

  $$\begin{array}{cc}\hfill \Phi \left(\mathbf{r}(\tau ,\nu ),{\mathcal{E}}^{(i-1)}\right)& ={\left({\mathcal{E}}^{(i-1)}\right)}^H\mathbf{r}(\tau ,\nu )\left({\mathbf{r}}^H(\tau ,\nu )\mathbf{r}(\tau ,\nu )\right){)}^{-1}{\mathbf{r}}^H(\tau ,\nu ){\mathcal{E}}^{(i-1)}\hfill \end{array}$$

  (20)

  $$\begin{array}{cc}& =\frac{{\left({\mathbf{r}}^H(\tau ,\nu ){\mathcal{E}}^{(i-1)}\right)}^{*}{\mathbf{r}}^H(\tau ,\nu ){\mathcal{E}}^{(i-1)}}{{\parallel \mathbf{r}(\tau ,\nu )\parallel }^2}=\frac{|{\mathbf{r}}^H(\tau ,\nu ){\mathcal{E}}^{(i-1)}{|}^2}{{\parallel \mathbf{r}(\tau ,\nu )\parallel }^2}.\hfill \end{array}$$

  (21)
  In this way, at each iteration, the IPI between the paths is directly canceled facilitating the new search.
• Global refinement (line 21): It is important to notice that this solution is prone to error propagation across all estimates. This may lead to residue vectors with residual peaks that are interpreted as small-amplitude paths by the algorithm. In order to avoid false alarms or multidetections due to residual uncancelled IPI, the P-IPIC algorithm needs a final refinement step in which all the estimates are refined while the residue in monitored. If the residue satisfies the stopping criterion, all remaining unrefined estimates are declared false alarms and, therefore, discarded. Specifically, during the search phase a number of paths equal to $\widehat{P}$ is detected. After that, during the refinement phase the paths are refined as in the DDIPIC algorithm with the only difference in using the cost function in Equation (19) with a Tikhonov regularization, also known as ridge-regression, and in performing channel gain refinement as in line 15. As discussed above, during this phase some paths may be discarded. Therefore, the refinement procedure provides a number of estimated paths equal to $\widehat{\widehat{P}}\le \widehat{P}$.
• Regularization (line 15): Another difference from the original version of the IPI cancellation algorithm is that the Ordinary Least Squares (OLS) estimate of channel gains in Equation (18) is replaced by a Regularized Least Squares (RLS) estimate to make the matrix inversion more stable. This is due to the fact that, due to residual IPI, some paths are detected multiple times and as a result the CDDPM matrix may have some almost collinear columns. Therefore, the RLS estimate of the channel gains vector is given by

  $$\widehat{\mathbf{g}}={\left({\mathbf{R}}^H(\widehat{\tau },\widehat{\nu })\mathbf{R}(\widehat{\tau },\widehat{\nu })+N_0\mathbf{I}\right)}^{-1}{\mathbf{R}}^H(\widehat{\tau },\widehat{\nu })\mathbf{y}$$

  (22)

  where $\mathbf{I}$ is the identity matrix.

3.5. Latency/Complexity Comparison

The key advantage of the proposed approach is that it requires a single final refinement procedure compared to the $\widehat{P}-1$ refinements performed by the DDIPIC algorithm. A rough latency comparison of the two approaches is shown in Table 1 where we assumed that the P-IPIC algorithm does not overestimate the number of paths significantly during the first iterative procedure, thus $\widehat{\widehat{P}}\approx \widehat{P}$. Under this hypothesis, and denoting by $T_{\mathrm{coarse}}$ the average time needed to perform a coarse estimate and with $T_{\mathrm{fine}}$ the average time needed to perform a fine estimate, we can make the following observations.

• DDIPIC: it performs a coarse and fine estimation for every path in a time equal to $\widehat{P}(T_{\mathrm{coarse}}+T_{\mathrm{fine}})$. Moreover, from the second iteration to the end it performs $\widehat{P}-1$ refinement procedures of the estimates. Therefore, at each refinement procedure, it performs i coarse and fine estimates.
• P-IPIC: it performs a coarse and fine estimation for every path in a time equal to $\widehat{P}(T_{\mathrm{coarse}}+T_{\mathrm{fine}})$ during the search phase. Moreover, in the refinement phase performs a second coarse and fine estimate for each path.

During the search phase, the proposed P-IPIC algorithm runs coarse and fine estimation procedures of channel delays and Dopplers performing maximization of the cost function in Equation (21), while during the refinement phase maximizes the regularized version of Equation (19). The complexity comparison of the cost function in Equation (19) and the RCF is shown in Table 2. It can be noticed that the first advantage in computing the RCF is that it doesn’t require matrix inversion but only multiplications and additions, instead the computation of the cost function in Equation (19) does. In order to compare the two cost functions we can make the following observations.

• DDIPIC cost function in (19): during the i-th iteration the number of required operations is

  (a)

  ${\mathbf{y}}^H\mathbf{R}\left(\tau ,\nu \right)$, $\mathbf{y}\in {\mathbb{C}}^{MN\times 1}$, $\mathbf{R}\left(\tau ,\nu \right)\in {\mathbb{C}}^{MN\times i}$ requires $MNi$ multiplications and $(MN-1)i$ additions. Then, ${\mathbf{R}}^H\left(\tau ,\nu \right)\mathbf{y}$ is simply obtained applying the Hermitian operator

  (b)

  ${\mathbf{R}}^H\left(\tau ,\nu \right)\mathbf{R}\left(\tau ,\nu \right)\in {\mathbb{C}}^{i\times i}$ requires $MN{i}^2$ multiplications and $(MN-1)i^2$ additions. After that, matrix inversion has a complexity in the order of $\mathcal{O}\left(i^3\right)$

  (c)

  ${\mathbf{y}}^H\mathbf{R}\left(\tau ,\nu \right){\left({\mathbf{R}}^H\left(\tau ,\nu \right)\mathbf{R}\left(\tau ,\nu \right)\right)}^{-1}$ requires $i^2$ multiplications and $i(i-1)$ additions

  (d)

  The matrix multiplication of the $1\times i$ matrix obtained in the previous step by ${\left({\mathbf{y}}^H\mathbf{R}\left(\tau ,\nu \right)\right)}^H$ requires i multiplications and $i-1$ additions

• P-IPIC RCF: during the i-th iteration the number of required operations is

  (a)

  ${\mathbf{r}}^H(\tau ,\nu ){\mathcal{E}}^{(i-1)}$,$\mathbf{r}\in {\mathbb{C}}^{MN\times 1}$,${\mathcal{E}}^{(i-1)}\in {\mathbb{C}}^{MN\times 1}$ requires $MN$ multiplications and $MN-1$ additions. After that, an additional multiplication is needed to compute the squared module

  (b)

  ${\parallel \mathbf{r}(\tau ,\nu )\parallel }^2$ requires $MN$ multiplications and $MN-1$ additions

  (c)

  One final division is required to compute the ratio between numerator and denominator

The computation of the RCF is independent on the iteration index i, while the computational complexity of the cost function in Equation (19) increases as $\mathcal{O}\left(i^3\right)$ due to the inversion of a $i\times i$ matrix.

4. Simulation Results

In this section, the proposed method is evaluated and compared against the baseline DDIPIC method.

4.1. Simulation Parameters

To investigate the performance of the proposed approach, we consider OTFS modulation with $M=16$ and $N=16$ with subcarrier spacing of $\Delta f=25\mathrm{kHz}$ and a carrier frequency $f_c=5.1\mathrm{GHz}$. Thus, delay and Doppler resolutions are ${\tau }_{\mathrm{res}}=2.5\mu \mathrm{s},$ and ${\nu }_{\mathrm{res}}=1.5625\mathrm{kHz}$. As in [11], the parameters for running channel estimation algorithms are $P_{max}=15,s_{max}=10,{ϵ}_{\tau }={10}^{-10},{ϵ}_{\nu }={10}^{-2},m_{\tau }=10,n_{\nu }=10$ . Finally, as discussed above, the convergence tolerance parameter for the stopping criterion is set according to the noise variance and OTFS frame parameters.

4.2. Simulation Scenarios

The high mobility multipath channel is assumed to have $P=4$ Rayleigh faded paths with exponential Power Delay Profile (PDP) with decaying time constant of $10\mu \mathrm{s}$ and a maximum delay of $10\mu \mathrm{s}$. We consider two different case studies corresponding to different scenarios:

• Low-IPI Regime: propagation delays and Doppler shifts assume fractional values but the paths are far enough to guarantee low IPI. Thus, all paths are separable in the DD domain.
  The delays and Doppler shifts of the paths are $\tau =\left[00.86{\tau }_{\mathrm{res}}2.11{\tau }_{\mathrm{res}}3.2{\tau }_{\mathrm{res}}\right]\mathrm{s}$, $\nu =[-1.87{\nu }_{\mathrm{res}}2.24{\nu }_{\mathrm{res}}-3{\nu }_{\mathrm{res}}1.3{\nu }_{\mathrm{res}}]\mathrm{Hz}$, respectively.
• High-IPI Regime: propagation delays and Doppler shifts assume fractional values, and the paths are close enough to create high IPI.
  The propagation delays of the paths are $\tau =\left[02.33.157.7\right]\mu \mathrm{s}$. The maximum User Equipment (UE) speed is $v_{UE}=190\frac{\mathrm{m}}{\mathrm{s}}\left(684\frac{\mathrm{km}}{\mathrm{h}}\right)$ resulting into a maximum Doppler spread of ${\nu }_{max}=\frac{v_{UE}}{c}{f}_c=3.23\mathrm{kHz}$, where $c=3·{10}^8\frac{\mathrm{m}}{\mathrm{s}}$ is the speed of light.
  The Doppler shifts of all paths are generated assuming Jakes’ Doppler Power Spectrum (DPS) using ${\nu }_i={\nu }_{max}cos\left({\theta }_i\right)$ where ${\theta }_i\sim \mathcal{U}[0,2\pi ]$.

An example of a pilot frame received in both regimes is provided in Figure 1.

4.3. Performance Metrics

We will consider the following performance metrics:

• Channel estimation: The accuracy of the channel estimate is measured through the Normalized Mean Square Error (NMSE) computed as

  $$\mathrm{NMSE}=\mathbb{E}\left[\frac{\parallel \mathbf{H}-\widehat{\mathbf{H}}{\parallel }_F^2}{{\parallel \mathbf{H}\parallel }_F^2}\right].$$

  (23)
• Channel parameter estimation: The accuracy of the channel parameter estimates is also considered evaluating the Normalized Root Mean Square Error (NRMSE) as

  $$\mathrm{NRMSE}\left(\vartheta \right)=\sqrt{\frac{\mathbb{E}\left[{\sum }_{i=1}^P{|{\vartheta }_i-\widehat{{\vartheta }_i}|}^2\right]}{{\sum }_{i=1}^P{\left|{\vartheta }_i\right|}^2}},$$

  (24)

  where $\vartheta$ is a variable that represents delay, Doppler or channel gain.
• Model order performance: The estimation of P will be evaluated by the average number of estimated paths and the Probability Mass Function (PMF) of the number of detected paths.
• Communication: Communication performances of the proposed approach are evaluated by measuring the Bit Error Rate (BER) as a function of the SNR/bit given by $E_b/N_0$ where $E_b$ is the average energy per bit. We consider a Quadrature Amplitude Modulation (QAM) with $Q=16$ constellation points and $E_b=\frac{2(Q-1)}{3lo{g}_{2}Q}$. The input-output relation for data frames is the one in Equation 9. At the receiver, a Linear Minimum Mean Square Error (LMMSE) equalizer provides the estimated symbols as

  $$\widehat{\mathbf{x}}={\left({\widehat{\mathbf{H}}}^H\widehat{\mathbf{H}}+\frac{1}{\mathrm{SNR}}\mathbf{I}\right)}^{-1}{\widehat{\mathbf{H}}}^H\mathbf{y},$$

  (25)

  where $\mathrm{SNR}=\frac{MN{E}_s/NT}{M\Delta f{N}_0}=\frac{E_s}{N_0}$ is the communication Signal-to-Noise Ratio and $E_s=E_{b}lo{g}_{2}Q$ is the average energy per symbol. After that, a symbol-by-symbol Maximum Likelihood (ML) detector takes the decision.

4.4. Results and Discussion

4.4.1. Channel Estimation Performance

Figure 2 shows the NMSE as a function of the pilot SNR in both low- and high-IPI regimes. In the Low-IPI case at low PSNR the performances of the proposed variant and the DDIPIC algorithm are comparable, while at higher values of PSNR the proposed variant outperforms the DDIPIC algorithm, achieving much lower NMSE. On the other hand, in the high-IPI regime, the P-IPIC algorithm achieves a NMSE 1-2 dB lower than the DDIPIC, except for the case of PSNR=35 dB where the NSME performances become comparable. Performances of the threshold method in [10] are also presented in both cases as a baseline for comparison and it can be noticed they are drastically degraded due to the presence of IPI since it assumes integer DDs.

4.4.2. Channel Parameter Estimation Performance

Figure 3 shows the NRMSE of channel gains, delays and Doppler shift as a function of pilot SNR in a channel randomly generated according to channel models 1 and 2. In both low- and high-IPI regimes, the P-IPIC is capable of estimating channel parameters with higher accuracy. It can be noticed that, in the High-IPI case, at low PSRN the NRMSE is high due to the fact that small amplitude paths are missed. On the other side, at high PSNR the NRMSE for the P-IPIC increases because small errors in the parameter estimates during the search phase lead to higher residual IPI that is interpreted as small-amplitude paths, reducing the overall accuracy of the estimates. However, the P-IPIC still outperforms the DDIPIC algorithm for high PSNR.

4.4.3. Number of Detected Paths

Figure 4 shows the average number of detected paths as a function of PSNR. In the Low-IPI case, performances of both algorithms are comparable. At low PSNR, the average number of detected paths is below the true value ($P=4$) because small-amplitude paths are missed. At higher PSNR values, both algorithms are able to detect the right number of paths. On the other hand, in the High-IPI case detection capabilities are still comparable but at high PSNR values both algorithms overestimate a bit the number of propagation paths due to high residual interference effects. Figure 5 shows the probability distribution among all the possible values for P in the case of PSNR=35 dB. In the Low-IPI regime, both algorithms detect the right number of paths with a probability close to $\sim 0.98$. Thus, the probability of miss detection and the probability of false alarm are close to zero (probability of miss dtection $<0.02$). In the high-IPI regime, the probability that the algorithm finds the correct number of paths is reduced to $\sim 0.73$ for the DDIPIC and $\sim 0.65$ for the P-IPIC. The probability of miss detection is close to zero in both cases ($\sim 0.015$ for the P-IPIC). In the end, the probability of false alarm is $\sim 0.27$ for DDIPIC and $\sim 0.34$ for P-IPIC.

4.4.4. Communication Performance

BER performance with the estimated channel (PSNR=30 dB) and with perfect Channel State Information at the Receiver (CSIR) are provided in Figure 5, showing the BER performance as a function of the SNR/bit in both low- and high-IPI regimes. As expected, the proposed P-IPIC achieves lower bit error rates in both regimes. Further, BER curves are closer to the ideal curves obtained with perfect channel state information at receiver.

5. Conclusions and Future Works

In this work, we addressed the challenge of low-complexity fractional delay-Doppler (DD) channel estimation in Orthogonal Time Frequency Space (OTFS) systems due to Inter-Path Interference (IPI). A variant of the Delay-Doppler Inter-Path Interference Cancellation (DDIPIC) algorithm is proposed, which progressively cancels IPI as estimates are obtained, requiring only a final refinement procedure. This approach reduces the latency of the algorithm, with the time difference in latency being almost proportional to the square of the number of estimated paths. Numerical results demonstrate that the proposed algorithm achieves both lower Normalized Mean Square Error (NMSE) and Bit Error Rate (BER) compared to other schemes with comparable detection capabilities. The high-IPI scenario is also taken into account, in which the delay and Doppler resolutions are not sufficiently fine, and IPI cancellation becomes tricky.

In future work, the usage of Deep Learning (DL) approaches, as in [21], can be considered to further reduce computational complexity while maintaining good performance.