Spectrally Efficient Time-Frequency Space Modulation: A Delay-Doppler Domain FTN Waveform for ISAC Systems

Spectrally Efficient Time-Frequency Space Modulation: A Delay-Doppler Domain FTN Waveform for ISAC Systems Michele Mirabella∗†, Pasquale Di Viesti∗†, and Giorgio Matteo Vitetta∗† ∗Dept. of Engineering “Enzo Ferrari”, University of Modena and Reggio Emilia, Modena, Italy †Consorzio Nazionale Interuniversitario per le Telecomunicazioni (CNIT) Email: {michele.mirabella, pasquale.diviesti, giorgio.vitetta}@unimore.it Abstract—This paper introduces the spectrally efficient timefrequency space (SETFS) modulation, a two-dimensional (2D) modulation extending the concepts of spectrally efficient frequency division multiplexing (SEFDM) and orthogonal time–frequency space (OTFS) modulation by increasing the spectral overlap at two levels, i.e., between subcarriers and sub-bands. This double compression improves the use of the available spectrum. Moreover, in SETFS, the periodicity in both time and frequency is maintained through the use of a double cyclic prefix (DCP) in order to simplify channel equalization and avoid 2D intersymbol interference (ISI). However, the presence of a spectral overlap induces 2D inter-carrier interference (ICI), that makes symbol detection more complicated. To mitigate channel estimation errors and also support high-resolution sensing, a novel channel estimator, dubbed Newton-based estimation and spectral cancellation algorithm (NESCA) is proposed. Numerical results show that SETFS achieves up to a 19% gain in the achievable communication rate with respect to OTFS at the cost of higher receiver complexity. Moreover, SETFS combined with NESCA provides near optimum performance in delay and Doppler estimation compared to other two techniques, thus offering a favorable trade-off between spectral efficiency, computational complexity, and sensing accuracy. Index Terms—Channel Estimation, Detection, Faster-thanNyquist, Orthogonal Time-Frequency Space, Spectral Efficiency I. INTRODUCTION With the advent of the sixth-generation (6G) of wireless communications, integrated sensing and communication (ISAC) has emerged as a key paradigm to enhance spectral efficiency (SE) while simultaneously enabling a wide range of sensing-oriented applications. A crucial challenge in ISAC is the design of waveforms able to offer good performance in both communication and sensing under specific spectral and complexity constraints [1]. Additionally, as carrier frequencies continue to increase, signal propagation may be substantially affected by Doppler shifts and time-variant channel effects; this can severely degrade the performance of traditional mulThis work has been supported by the European Union under the Italian National Recovery and Resilience Plan (PNRR). Specifically, it has been carried out within the ”Telecommunications of the Future” partnership (PE00000001 - program ”RESTART”), CUP E63C22002040007 - D.D. n.1549 of 11/10/2022, and within the initiatives of Mission 4 Component 2, Investment 1.4 (D.D. 1033 17/06/2022, CN00000023). Furthermore, the project has involved the MOST – Sustainable Mobility National Research Center. ticarrier modulations such as orthogonal frequency division multiplexing (OFDM) [2]. Recently, to overcome the last limitation, orthogonal time–frequency space (OTFS) modulation has been proposed as a promising alternative to OFDM, thanks to its inherent robustness against Doppler effects and its ability to provide full time–frequency (TF) diversity [3], [4]. OTFS maps information symbols in the delay–Doppler (DD) domain and relies on the discrete symplectic Fourier transform (DSFT) to spread each symbol across all TF resources. When implemented with a double cyclic prefix (DCP) [5], the channel filtering results in atwo-dimensional (2D) cyclic convolution in the DD domain, equivalent to an element-wise multiplication in the TF domain. This property ensures the absence of inter-symbol interference (ISI) between adjacent OTFS symbols and enables simple single-tap equalization. Moreover, channel estimation can be efficiently performed through the transmission of a single OTFS pilot symbol at the beginning of each frame [6]. In this work, we extend the concept of spectrally efficient multicarrier modulation, originally introduced in onedimensional (1D) frequency-domain (FD) modulations such as spectrally efficient frequency division multiplexing (SEFDM) [7]–[9], to the 2D OTFS framework. More specifically, we introduce novel modulation format, dubbed spectrally efficient time–frequency space (SETFS), which generalizes SEFDM to the DD domain. Unlike conventional OTFS, SETFS intentionally increases the spectral overlap at two levels, namely between subcarriers and sub-bands, thus breaking the orthogonality in both time and frequency while achieving higher SE. To the best of our knowledge, no prior work has investigated such a non-orthogonal 2D modulation. The closest related studies concern multiple-input multiple-output (MIMO) fasterthan-Nyquist (FTN) signaling for OTFS systems, where the lack of orthogonality is introduced across antennas, while the underlying modulation remains orthogonal [10], [11]. In contrast, in this work we introduce SETFS modulation as a generalization of the SEFDM concept to the 2D OTFS framework, including the use of both time-domain (TD) and FD cyclic prefixes (CPs) in its generation (i.e., a double cyclic prefix, DCP). The resulting modulation format ensures cyclicity in both time and frequency, and, similarly as OTFSDCP, allows for simple element-wise equalization, so avoiding the need for higher-complexity decision-feedback equalization (DFE) techniques (e.g. see [12] for zero-padded, ZP OTFS). However, unlike OTFS-DCP, the spectral compression introduces intentional 2D ICI, which breaks the full orthogonality of the modulation. As a result, additional detection complexity is required at the receiver to properly recover the transmitted symbols while preserving the SE gain [8], [9]. The main contribution of this manuscript is threefold and can be summarized as follows: 1) We propose a novel modulation scheme, named SETFS, which extends SEFDM through the use of fractional discrete symplectic Fourier transform (FrDSFT) and of proper pulse shaping (root-raised cosine, RRC, pulses). 2) We develop a low-complexity SETFS baseband transceiver architecture, including channel estimation via a single SETFS pilot symbol, single-tap minimum mean square error (MMSE) equalization and message-passing (MP) detection adapted from SEFDM systems [9]. 3) We develop a novel channel refinement algorithm, named Newton-based estimation and spectral cancellation algorithm (NESCA), which estimates the channel parameters (namely, the delays, Doppler shifts, amplitudes, and overall number of paths) to reconstruct a denoised channel representation with the aim of improving the subsequent equalization accuracy. The rest of this manuscript is organized as follows: in Section II, the main steps for the generation of the SETFS signal are sketched and the baseband models for the received (RX) signal in the presence of a doubly selective channel are illustrated. Section III is devoted to the development of the NESCA, whereas its performance is compared with the conventional and widely employed periodogram method (also employed in the OFDM and OTFS literature; e.g., see [1], [3]) in Section IV. Data detection performance for SETFS are compared with that achieved by OTFS-DCP in terms of achievable capacity when the parallel-branch message passing algorithm (PBMPA) developed in [9] is employed for SETFS after channel equalization to recover the intentional ICI introduced by compression. Finally, some conclusions are offered in Section V. Notation: In this paper, the following notation is adopted: 1) (·)∗and (·)Hdenote complex conjugate and complex conjugate transpose, respectively; 2) the symbols ⊙and ⊘ represent the Hadamard and product and division operators, respectively; 3) the symbols ⊗and ⊛represent the Kronecker and Khatri-Rao product operators, respectively; 4) modN[·] indicates the modulo Noperator; 5) X≜[xm,n]defines a matrix Xof proper size and xm,n denotes the element appearing on its mth row and nth column; 6) y= vec(Y) defines an (MN)-dimensional column vector resulting from the ordered concatenation of the columns of the M×Nmatrix Y;7)ℜ{x}indicates the real part of the complex variable x;8)0Ndenotes the zero vector of size N;9)ΞNindicates the unitary discrete Fourier transform (DFT) matrix of size N, whose (p, q)element is exp(−j2πpq/N)√N. II. SIGNAL AND SYSTEM MODELS In this section, we focus on the derivation of an analytical model for the complex envelope of the transmitted SETFS signal and of the corresponding received signal in the presence of a doubly selective fading channel. The baseband architecture of the proposed SETFS-based communication system is illustrated in Fig. 1. In the following, we assume that the SETFS signal conveys the M×Nmatrix C≜[cm,n](with m= 0,1, ..., M −1and n= 0,1, ..., N −1) containing Nuchannel symbols belonging to an Mcary constellation and Nsc ≜MN −Nuzeroes (associated with to the suppressed carriers, SCs). Our derivation of the SETFS signal model parallels that provided in [5, Sec. II-C] for an OTFS-DCP signal; the main difference is represented by the fact that the matrix C undergoes an order (M, N)fractional inverse discrete symplectic Fourier transform (FrIDSFT) of (positive) parameters (β(FD), β(TD))both less than or equal to one (note that, if (β(FD), β(TD)) = (1,1), SETFS coincides with OTFS); this produces X≜[xm,n]=FM,β(FD) C FH N,β(TD) , (1) where FX,β is a square matrix of order X, whose (p, q)th element is equal to exp(−j2πpqβ/X)/√X. The matrix X, in (1), is cyclically extended along both dimensions, with period Nin the index n(i.e., in the TD) and period Min the index m(i.e., in the FD); the new elements are generated as xk,l =xmodM[k],modN[l]for k /∈ {0,1, ..., M −1}and l /∈ {0,1, ..., N −1}. This results in the introduction of a TD CP, a FD CP and a FD cyclic postfix (CPO) having sizes N(TD) cp ,N(FD) cp and N(FD) cpo , respectively. The complex envelope of the SETFS signal conveying the cyclically extended symbol can be expressed in a similar way as that of OTFS-DCP, i.e., as (see [5, Eq. (62)]) s(t;C) = M/2+N(FD) cpo −1 X k=−(M/2+N(FD) cp ) s(TD) k(t;C), (2) where s(TD) k(t;C)≜ N−1 X l=−N(TD) cp xk,l pt−lTs/β(TD) ·expj2πkBsβ(FD)t−lTs/β(TD). (3) It is important to point out that s(TD) k(t;C)(3) can be interpreted as the frequency-shifted version of a SEFDM signal (e.g., see [9, Eq. (2)]), characterized by a CP of size N(TD) cp , TX pulse p(t), symbol interval Tsand a frequency shift equal to kBsβ(FD), where Bsrepresents, up to the compression factor β(FD), the frequency spacing between adjacent SEFDM signals. This observation suggests to adopt the same pulse and the same rules for subcarrier suppression as in SEFDM (see [9]). This leads to: 1) selecting p(t)as a pulse whose spectrum P(f)is the RRC with roll-off factor α∈[0,1] [5, Eq. (17)]; 2) Symbol Mapping Transmitter SC Insertion FrIDSFT DCP Insertion Receiver ( ; )stC   ()pt ( , )ht  Filter Bank ( ; )rtC Sampling & S/P TD CP Removal Channel Estimator Equalizer Data Detection Data out Data in ˆd    m,nc FrDSFT ˆd  Fig. 1. Baseband model of the considered SETFS-based communication system. meeting the condition Bsβ(FD) =β(TD)/Ts; 3) setting to zero multiple (i.e., Nsc) elements of the matrix C. Based on the last choices and on its cyclic structure, the signal s(TD) k(t;C) (3) can be represented through its Fourier series (FS) as s(TD) k(t;C) = +∞ X q=−∞ S(TD) k,q (C) exp(j2πfqt)(4) in the interval [0, T ′], with T′≜NTs/β(TD); here, fq≜q/T′ the qth SEFDM subcarrier frequency. Moreover, in (4), S(TD) k,q (C) = β(TD) √NTs Xk,q Pq−kN T′(5) is the qth Fourier coefficient of s(TD) k(t;C)(3) and Xk,q ≜1 √N N−1 X n=0 xk,n exp−j2πq n N. (6) is the qth coefficient of the order NDFT of the sequence {xk,n}with respect to the index n. Let us assume now the signal s(t;C)(2) is sent over a doubly selective fading channel, having channel impulse response (CIR) h(t, τ) = L−1 X l=0 alexp(j2πνlt)δ(τ−τl), (7) characterized by the Ltriplets {(al, τl, νl); l= 0,1, ..., L − 1}; here, al,τland νlrepresent the gain, the delay1and the Doppler shift associated with the lth path, respectively and Lis the overall number of paths. The baseband equivalent of the received signal r(t;C)feeds a filter bank, consisting of M distinct matched filters (see [5, eq. (99)]). The response of the ˜ kth matched filter (with ˜ k=−M/2,−M/2+1, ..., M/2−1), having impulse response ϕ˜ k(t)=p∗(−t) exp(j2π˜ ktβ(TD)/Ts), (8) is sampled at the instant t˜n=τL−1+ ˜nTs, with ˜n= 0,1, ..., N −1; this produces the sample r[˜ k, ˜n]. The MN samples acquired at the output of the filter bank are collected in the M×Nmatrix R≜[r[˜ k, ˜n]]. In [5, Sec. II-C.4], it is proved that the property of double periodicity, ensured by the 1We assume the CIR components are arranged in ascending order of delays, so that τ0and τL−1represent the minimum and maximum delays, respectively. introduction of a DCP, and the specific choice of the pulse p(t)allows to express Ras R≜H⊙T(C)+W; (9) here, H≜[h[˜ k, ˜n]] and W≜[w[˜ k, ˜n]] are M×Nmatrices that represent the channel state information (CSI) and the complex noise matrix affecting R(an additive white Gaussian noise, AWGN, model is assumed for the 2D sequence {w[˜ k, ˜n]}; the variance of noise samples is denoted σ2 w), respectively. Additionally, h[˜ k, ˜n]≜ L−1 X l=0 Alexp(−j2π˜ kFτl) exp(j2π˜nFνl), (10) with ˜ k= 0,1, ..., M −1and ˜n= 0,1, ..., N −1,Al≜ alexp(j2πνl(τl+τL−1)),Fτl≜(τl−τL−1)β(TD)/Tsand Fνl≜νl/(Bsβ(FD))represent the complex gain, the normalized delay and normalized Doppler shift associated with the lth channel path, respectively and (see (9)) T(C) = FM,β(FD) (C⊙¯ G)FH N,β(TD) , (11) is an M×Nmatrix accounting for the channel symbols through Cand the pulse shape through ¯ G≜[¯gp,q]. As shown in [5, eq. (82)], ¯gp,q exhibits a weak dependence on the normalized frequencies Fτland Fνl. For this reason, the coefficients {¯gp,q}corresponding to the Nuunsuppressed channel symbols are characterized by similar amplitudes (in particular, by a unitary amplitude). This property significantly simplifies channel estimation, equalization and subsequent data detection. These steps are described in the next section. III. RECEIVER OPERATIONS This section is organized into two parts. The first part revisits the channel estimation, equalization, and data detection stages, which constitute the core of the receiver processing chain. The second part focuses on the derivation of a novel Newton-based algorithm for the estimation of channel parameters, which is able to provide a refined reconstruction of the channel response. While the first part is essential to ensure reliable data recovery, the second provides a more detailed channel characterization. Note that, such parameter estimation not only improves communication performance at the receiver, but can also be exploited in co-located ISAC systems to enable radar sensing capabilities. A. Channel Estimation, Equalization & Detection We first consider the transmission of a SETFS frame consisting of DSETFS symbols, the first serving as a pilot. In the following the SETFS pilot symbol2is denoted C0, whereas the remaining SETFS symbols are {Cd≜[c(d) m,n]; d= 1,2, ..., D− 1}; moreover, the RX signal matrix R(9) associated with the dth SETFS symbol is denoted Rd. In addition, we assume that: 1) all the elements of the matrix Rd(9) are affected by AWGN with variance σ2 wfor any d; 2) C0contains a single non-null pilot channel symbol, belonging to the unsuppressed region of C(the remaining elements of the SETFS pilot symbol are equal to zero); 3) the channel variations over each SETFS frame can be deemed negligible; 4) thanks to the last assumption, CSI can be estimated on the basis of C0and, then, be exploited for channel equalization over the other (D−1) SETFS symbol intervals of the same frame. The estimate ˆ Hof His evaluated as ˆ H≜ˆ h[˜ k, ˜n]=R0⊘T(C(0)), (12) where T(C(0))is known at receive side (see (11) with C= C(0)). Note that the CSI estimate achieved through (12) can also be directly exploited to equalize the remaining (D−1) SETFS symbols in the TF domain. In particular, equalization can be accomplished in the TF domain to compensate for the distortions introduced by the channel (see (9)). The dth MMSE equalized matrix is computed as ˆ Td≜ˆ T(C(d))=Rd⊙Qmmse, (13) where Qmmse ≜[q[˜ k, ˜n]] is an M×Nmatrix, whose (˜ k, ˜n)th element is q[˜ k, ˜n]] ≜ˆ h∗[˜ k, ˜n]/(ˆ h∗[˜ k, ˜n]ˆ h[˜ k, ˜n]+ ˆσ2 w), with ˆσ2 w being an estimate of the noise variance in (9). To ease our description of symbol detection, we rewrite the expression (11) in vector form for the dth SETFS data symbol as γd≜vec(T(C(d)))=ΨH β(FD),β(TD) ¯ cd, (14) with d= 1,2, ..., D −1; here, Ψβ(FD),β(TD) ≜FN,β(TD) ⊗ FH M,β(FD) is an (MN)×(MN)matrix and ¯ cd≜vec(C(d)⊙ ¯ G)is an (MN)-dimensional vector. In fact, detection can be accomplished on a symbol-by-symbol basis by first applying the order (M, N)FrDSFT of parameters (β(FD), β(TD))to ˆ γd≜vec( ˆ Td)(see (13)); this produces ˆ cd=Γ¯ cd+˘ w, (15) where Γ≜Ψβ(FD),β(TD) ΨH β(FD),β(TD) is the (MN)×(MN) matrix accounting the 2D intentional ICI introduced by spectral compression. Detecting the channel symbols in (15) requires compensating for the intentional 2D ICI represented by Γ. In this work, the detection stage is not explicitly addressed; instead, we rely on the PBMPA detector proposed in [9] for a SEFDM-based communication system, since the signal model and the structure of the ICI matrix Γin (15) are similar (see [9, Eq. (7)]). 2The considered pilot rate is thus 1/D. B. Channel parameter estimation algorithm In this section, we focus on the problem of estimating the parameters of the communication channel, namely a) the set of parameters {(Al, Fτl, Fνl); l= 0,1, ..., L −1}and b) the number of paths L, which completely describe the matrix H (see (10)). The availability of an estimate ˆ Hof the channel H(see (10)) is assumed. To begin, we rewrite the model (12) in vector form as ˆ h≜vec( ˆ H)=h+w, (16) where wis the (MN)-dimensional column vector of noise samples affecting ˆ h, whereas h≜vec(H) = B(fτ,fν)a; (17) here, a≜[A0, A1, ..., AL−1]T,fτ≜[Fτ0, Fτ1, ..., FτL−1]T and fν≜[Fν0, Fν1, ..., FνL−1]Tcollect the Lcomplex gains, normalized delays and normalized Doppler frequencies, respectively. Moreover, in (17), B(fτ,fν)≜C∗ N(fν)⊛CM(fτ)(18) is an (MN)×Lmatrix, Cp(f)≜[¯ b0(f),¯ b1(f), ..., ¯ bp−1(f)]T is a p×Lmatrix for any positive integer p,¯ bx(f)≜ [¯ bx(f0),¯ bx(f1), ...,¯ bx(fL−1)]Tis an L-dimensional column vector and ¯ bx(fl)≜exp(−j2πxfl)for any l. Given ˆ h(16), our goal is estimating a,fτ,fνand L. The algorithm we propose is dubbed NESCA and consists of four steps: 1) an initialization step based on the computation of the 2D spectrum of ˆ h; 2) a residual spectrum maximization step; 3) a Newton-based frequency refinement and 4) a spectral cancellation procedure. Each of the four steps is briefly described in the following Step 1 - The initialization step consists in: 1) evaluating the 2D spectral matrix (through a discrete symplectic Fourier transform, DSFT) Y≜[y[ ˜m, ˜n]] = LrLcΞH M0ˆ HZP ΞN0(19) where ˆ HZP ≜ˆ H,(0M0T N0−N) (0M0−M0T N) (0M0−M0T N0−N). (20) is the ZP version of ˆ H; the last matrix has size M0×N0, where M0≜LrMand N0≜LcN, with Lrand Lcbeing the oversampling factors, for the rows of ˆ Hand its columns, respectively (12). 2) Initializing the residual spectrum Yres ≜[yres[˜ k, ˜n]] = Yand the parameter ¯ L= 0. Step 2 - The residual spectrum energy ∥Yres∥2is evaluated; if this quantity is greater than or equal to3a given detection threshold ϵ,¯ Lis increased by one and the parameters associated with a new (i.e., the (¯ L−1)th) 2D tone are estimated as ˆ Fτ¯ L−1= ˆm/M0, (21) 3If this condition is not met, the algorithm stops and the ¯ L-dimensional estimate vectors of the channel parameters are provided as output. ˆ Fν¯ L−1= ˆn/N0−0.5(22) and ˆ A¯ L−1=yres[ ˆm, ˆn], (23) where (see (19)) ( ˆm, ˆn) = arg max ˜m,˜nyres[ ˜m, ˜n]. (24) Step 3 - The third step is fed by ¯ L-dimensional vectors ˆ a≜[ˆ A0,ˆ A1, ..., ˆ A¯ L−1]T,ˆ fτ≜[ˆ Fτ0,ˆ Fτ1, ..., ˆ Fτ¯ L−1]Tand ˆ fν≜[ˆ Fν0,ˆ Fν1, ..., ˆ Fν¯ L−1]T, collecting the estimates of the complex gains, normalized delays and normalized Doppler frequencies, respectively, with ¯ L≤L. It refines these estimates through an iterative procedure, which maximizes, in an approximate fashion, the maximum likelihood (ML) cost function p(ˆ h|˜ a,˜ fτ,˜ fν) = 1 (πσ2 w)MN exp−L(ˆ h|˜ a,˜ fτ,˜ fν)/σ2 w. (25) In the last expression, the quantity L(ˆ h|˜ a,˜ fτ,˜ fν)≜∥ˆ h−B(˜ fτ,˜ fν)˜ a∥2(26) is the log-likelihood (LL) function of the considered optimization problem. The refinement procedure is initialized by setting ˆ a(0) =ˆ a,ˆ f(0) τ=ˆ fτand ˆ f(0) ν=ˆ fν,ˆ f(0) = [(ˆ f(0) τ)T,(ˆ f(0) ν)T]T and the iteration index iset to 1. Then, an iterative procedure, is started; it consists in: 1) A Complex amplitude update – The new complex amplitude vector ˆ a(i)=B†ˆ f(i−1) τ,ˆ f(i−1) νˆ h(27) is computed; here, ˆ f(i−1) τand ˆ f(i−1) νdenote the estimates of ˆ fτand ˆ fνat the end of the (i−1)th iteration. 2) A Frequency update – The new estimate ˆ f(i)=ˆ f(i−1) −µ¨ H−1 L(ˆ a(i),ˆ f(i−1))∇L(ˆ a(i),ˆ f(i−1))(28) is evaluated; here, ˆ f(i)= [(ˆ f(i) τ)T,(ˆ f(i) ν)T]T,µis a (positive) step size (in our simulations, this parameter has been set to one), and ∇L(˜ a,˜ f)and ¨ HL(˜ a,˜ f)represent, for a given couple (˜ a,˜ f)the gradient vector and Hessian matrix of the LL function L(26). The gradient and the Hessian of Lcan be expressed as4 ∇L(˜ a,˜ f)=−2hℜ˜ a∗⊙˙ BH ˜ fτ(ˆ h−B˜ a)T, ℜ˜ a∗⊙˙ BH ˜ fν(ˆ h−B˜ a)TiT , (29) and as ¨ HL(˜ a,˜ f)≜T˜ fτ,˜ fτU˜ fτ,˜ fν U˜ fτ,˜ fνV˜ fν,˜ fν, (30) 4The dependence of the matrix B(18) on the variables (˜ fτ,˜ fν)is not explicitly shown in the following to ease reading. respectively; here5, T˜ fτ,˜ fτ≜"∂2L ∂˜ fτl∂˜ fτl′#= 2ℜn(˜ a˜ aH)⊙˙ BH ˜ fτ ˙ B∗ ˜ fτ −˜ a∗ˆ h−B˜ aH¨ B˜ fτ,˜ fτo, (31) U˜ fτ,˜ fν≜"∂2L ∂˜ fτl∂˜ fνl′#= 2ℜn(˜ a˜ aH)⊙˙ BH ˜ fτ ˙ B∗ ˜ fν −˜ a∗ˆ h−B˜ aH¨ B˜ fτ,˜ fνo (32) and V˜ fν,˜ fν≜"∂2L ∂˜ fνl∂˜ fνl′#= 2ℜn(˜ a˜ aH)⊙˙ BH ˜ fν ˙ B∗ ˜ fν −˜ a∗ˆ h−B˜ aH¨ B˜ fν,˜ fνo (33) are ¯ Lׯ Lmatrices. Moreover, ˙ B˜ fxis an (MN)ׯ Lmatrix, whose lth column contains the partial derivative of the lth column of B(see (18) with (fτ,fν) = (˜ fτ,˜ fν)), evaluated with respect to ˜ fx. Similarly, ¨ B˜ fx,˜ fyisa(MN)ׯ Lmatrix, whose lth column contains the second-order derivative of the lth column of B(18) evaluated with respect to ˜ fxand ˜ fy. Note that, thanks to the linearity property of the derivative, the equality ¨ B˜ fx,˜ fy=¨ B˜ fy,˜ fxholds. A proof of the results illustrated above is sketched in Appendix A. Step 4 - Once the refinement step for the ¯ L2D tones is completed, a spectral cancellation procedure is carried out; this step aims at computing the residual 2D M0×N0spectrum Yres ≜yres[ ˜m, ˜n]=Y−Ycˆ a,ˆ fτ,ˆ fν, (34) where Yis defined by (19). Moreover, in (34), Yc(ˆ a,ˆ fτ,ˆ fν)≜ [yc[ ˜m, ˜n]] is an M0×N0matrix, whose ( ˜m, ˜n)th element is computed as yc[ ˜m, ˜n] = 1 MN ¯ L−1 X l=0 ˆ Alϕ˜m(ˆ Fτl, M0)ϕ∗ ˜n(ˆ Fνl, N0)(35) where ϕx(ˆ F, X0)≜1−expj2πx−F X0 1−expj2πx X0−F (36) with x= 0,1, ..., X0−1for any positive integer X0and normalized frequency F. This completes our description of last step of the NESCA. Note that: a) after spectral cancellation, the residual spectrum energy is evaluated again in Step 2; b) this step aims at improving estimation accuracy and may a have a significant impact in the presence of multiple closelyspaced channel paths. The NESCA algorithm is summarized in Algorithm 1. Its overall complexity, is NNESCA ∼ =N2D-FFT +NitL2MN, where N2D-FFT =M0N0log2(M0N0)is the complexity of the initialization (due to the evaluation of the 2D spectrum Y (19)). 5Note that the dependence of the matrices T˜ fτ,˜ fτ,U˜ fτ,˜ fνand V˜ fν,˜ fνon the variable ˜ ais not explicitly shown to ease notation. Algorithm 1: Newton-based estimation and spectral cancellation algorithm (NESCA) Input: The M×Nmatrix ˆ H(12), the detection threshold ϵ, the oversampling factors (Lr, Lc), the step size µand the overall number of refinement iterations Nit. 1Initialization: Compute Yfrom ˆ Hthrough (19), set ˆ a,ˆ fτland ˆ fνlto empty vectors and ¯ L= 0. Initialize Yres =Y. 2Energy criterion:If∥Yres∥2≥ϵ, then increase ¯ Lby one and go to step 3, otherwise go to Output. 3Residual spectrum maximization: Evaluate the new quantities (ˆ A¯ L−1,ˆ Fτ¯ L−1,ˆ Fν¯ L−1),through (23), (21) and (22), using Yres in the spectral maximization (24). Then, append ˆ Fτ¯ L−1,ˆ Fν¯ L−1and ˆ A¯ L−1=yres[ ˆm, ˆn]to ˆ fτ,ˆ fνand ˆ a, respectively. 4Refinement: Initialize ˆ a(0) =ˆ aand ˆ f(0) = [ˆ fT τ,ˆ fT ν]T. for i= 1 to Nit do 3aAmplitude update: compute ˆ a(i)through (27). 3bFrequency update: Compute ˆ f(i)= [(ˆ f(i) τ)T(ˆ f(i) ν)T]Tusing (28). end Set ˆ a=ˆ a(Nit),ˆ fτ=ˆ f(Nit) τ,ˆ fν=ˆ f(Nit) νand ˆ L=¯ L. 5Spectral cancellation: Update Yres through (34) using ˆ a,ˆ fτand ˆ fν; then, go to step 2. Output: The estimates ˆ a,ˆ fτ,ˆ fνand ˆ Lof a,fτ,fνand L, respectively. TABLE I COMPUTATIONAL COMPLEXITY ORDER OF THE CONSIDERED CHANNEL ESTIMATORS Algorithm O(·) 2D-FFT N2D−FFT 2D-FFTi N2D−FFT +L(IrIc+¯ M¯ N) NESCA N2D−FFT +Nit L2MN IV. NUMERICAL RESULTS Computer simulations have been run to assess the performance of the NESCA and to compare it with a) the 2D periodogram method and b) with its variant including a refinement procedure based on interpolation of the 2D periodogram [13, Sec. IV-A]; these are dubbed 2D-FFT and 2D-FFTi in the following. The computational complexity orders (O(·)) of the 2D-FFT algorithm, the 2D-FFTi algorithm and the NESCA are listed in Table I. Note that: 1) The 2D-FFT algorithm represents the baseline technique for the other two algorithms. 2) The 2D-FFTi refines each 2D frequency by interpolating the 2D spectral components around each peak detected by the 2D-FFT algorithm. The interpolation orders for the rows and columns of the 2D spectrum Y(19) are denoted Irand Ic, respectively. The interpolation procedure generates Lfiner 2D grids, each of size ¯ Mׯ Nand for which spectral maximization is carried out to refine the L2D frequencies and amplitudes. 3) The NESCA operates in an off-grid fashion and its complexity scales linearly with the variables Mand N, and grows quadratically with the number of channel paths. However, since Lis generally a small quantity, the major contribution to the complexity of the NESCA is given by its initialization step thorugh the DSFT used to evaluate the 2D spectrum Y (19). In our simulations, Gaussian quadrature rules (GQR) [14] have been employed to generate a channel model approximating a doubly selective channel characterized by a truncated exponential power delay profile (maximum delay τmax = 0.1µs) and Jake’s power spectrum for the Doppler effect (maximum Doppler shift νmax = 178 kHz). The resulted channel model has L= 3 paths, whose complex amplitudes are represented as the superposition of three complex exponentials. The following parameters have been selected for the SETFS frame and modulation format: 1) SETFS frame size D= 8;2)M=N= 32 in the generation of each data matrix C; 2) subcarrier spacing ∆f= 60 kHz; 3) TD CP size N(TD) cp =N/4; 4) FD CP and CPO sizes N(FD) cp =N(FD) cpo =M/16; 5) carrier frequency fc= 24 GHz; 6) roll-off factor for the RRC pulse α= 0.15; 7) cardinality of the phase-shift keying (PSK) constellation Mc= 4; 8) compression parameters β(TD) =β(FD) = 0.9. These choices result in: 1) channel symbol time Ts∼ =0.5µs; 2) frequency spacing between adjacent SEFDM symbols Bs= 1.92 MHz; 3) number of SCs for each SEFDM symbol Nsc = 5 (see comments above [5, Sec. II-A, Eq. (17)]); 4) useful bandwidth occupied by each SETFS symbol B≜MBs= 61.4MHz; 5) delay and Doppler resolutions, τbin ≜Ts/(Mβ(FD))∼ =16 ns and νbin ≜Bsβ(TD)/N ∼ =60 kHz. Note also that, from a communication perspective, the values selected for the parameters (β(FD), β(TD))have allowed us to achieve a good trade-off between SE and data detection complexity; in fact, reducing these two parameters improves SE, but increases the complexity of data detection because of the stronger ICI [8]. From a sensing perspective, instead, the values of (β(FD), β(TD))influence the 2D resolution grid; in fact, selecting smaller values for (β(FD), β(TD))entails a reduction in the spacing between normalized delays and an increase in that between normalized Doppler frequencies. This directly affects the achievable parameter estimation accuracy, as it modifies the related Cram´ er-Rao lower bound (CRLB). The performance of the considered methods for channel parameter estimation has been compared in terms of root mean square error (RMSE) achieved for delay (RMSEτ) and Doppler shift (RMSEν) estimation. The outcome {(ˆ Al,ˆ Fτl,ˆ Fνl); l= 0,1, ..., ˆ L−1}of these algorithms is employed to reconstruct the channel estimates through (10), which are then used by the MMSE equalizer. After that, the PBMPA detector is used to recover the intentional 2D ICI due to compression. The communication performance is evaluated in terms of effective capacity, defined as Ceff ≜ MI/(β(TD) β(FD)), where the mutual information (MI) is computed on the basis of the soft-output information generated by the considered detector. Note that [4, Eq. (54)] has been employed for the evaluation of the MI. The performance of the proposed SETFS modulation has also been compared with that of a reference OTFS-DCP system, adopting the same SETFS parameters except for β(TD) and β(FD), which are both set to one. Additionally, a conventional log-likelihood ratio (LLR) detector operating after MMSE equalization has been employed for OTFS-DCP, since it allows symbol-by-symbol detection, providing softoutput information with low computational complexity. This setup provides a fair benchmark to assess the communication and sensing trade-offs of SETFS with respect to OTFS-DCP. In generating our numerical results, the following choices have been made for all the considered algorithms: 1) the oversampling factors Lr=Lc= 4 have been used for 2DFFT, 2D-FFTi and NESCA; 2) the overall number of paths L and the noise variance σ2 whave been assumed to be known at the receiver; 3) Nit = 5 and µ= 1 have been selected for the NESCA. Additionally, the detection threshold ϵhas been always set to the 60% of the initial spectrum energy, i.e., p∥Y∥2(see (19)); 4) in 2D-FFTi method, Ir=Ic= 7 have been selected and the 2D spline interpolation has been employed to generate finer spectral samples over ¯ Mׯ N grid, with ¯ M=¯ N= 256; 5) the PBMPA detector has been configured to perform 4iterations. Some results about RMSEτand RMSEνare shown in Fig. 2, for an SNR ∈[−10,30] dB; the CRLB is also provided as a reference. The Ceff performance achievable in the same SNR range using all the considered methods is shown in Fig. 3. The results shown in Fig. 2 lead to the following conclusions: 1) An SNR threshold can be easily identified in all the RMSE curves; this threshold is approximately equal to 12 dB for all the considered algorithms. 2) Above the threshold, all the considered techniques achieve lower RMSEτwhen using SETFS with respect to OTFS-DCP. However, their RMSEνvalues are generally higher, suggesting less accurate Doppler estimation when employing SETFS. This trend is also reflected in the corresponding CRLBs; in particular, when OTFS-DCP is employed, the CRLB for delay estimation is higher, while that for Doppler estimation is lower than those observed with SETFS. 3) Both 2D-FFT and 2D-FFTi algorithms exhibit a RMSE floor at high SNRs (say, for an SNR higher than 20 dB). Keep in mind that, on the one hand, the 2D-FFT restricts the peak search to the 2D spectrum Y(19); thus, its accuracy is limited by the Fourier transform orders. On the other hand, the 2DFFTi employs a refinement step based on spectral interpolation and its accuracy is limited by the grid size. 4) Compared to the other methods, the NESCA achieves superior estimation accuracy (close to the CRLB) in terms of both RMSEτand RMSEνwhile maintaining manageable complexity; this is mainly due to the fact that the NESCA algorithm operates in an off-grid fashion. It is also worth noting that its RMSEτtends to deviate from the respective CRLBs at approximately 30 dB, while its RMSEνcontinues to follow the CRLB trend. This is due to the limited number of iterations carried out by the algorithm. Further simulations have evidenced that increasing Nit improves the NESCA RMSE performance at the price, however, of a higher computational complexity. 5) The computational complexity of the 2D-FFTi algorithm and NESCA is respectively 1.85 and 1.27 times higher than that of the 2D-FFT algorithm (see Table I); this means that the NESCA computational requirements are 27% higher than 2D-FFT and 46% lower than 2D-FFTi. Finally, the results illustrated in Fig. 3 show that: 1) The LLR detector in OTFS-DCP is less sensitive to channel reconstruction errors and achieves similar performance when combined with 2D-FFT, 2D-FFTi and NESCA algorithm. Its complexity is approximately O(Mc). 2) SETFS modulation saturates approximately to a 19% higher Ceff compared to OTFS-DCP thanks to spectral compression. However, this gain comes at the cost of a notably higher computational complexity6of the PBMPA detector compared to the LLR detector; note that the last detector cannot be used for SETFS as it does not compensate for the 2D intentional ICI. 3) The NESCA combined with the PBMPA detector outperforms the other techniques for SETFS within the SNR range [10,20] dB. V. CONCLUSIONS In this manuscript, a novel modulation format, dubbed SETFS and generalizing the SEFDM paradigm to the 2D DD domain of OTFS, and a channel estimator, called NESCA, to be used in a SETFS-based communication system operating over a doubly selective fading channel, have been proposed. The SETFS modulation achieves significant SE gains, but requires an higher detection complexity than OTFS because of the presence of 2D ICI. Our simulation results show that NESCA achieves accurate channel parameter estimation, outperforming both 2D-FFT and 2D-FFTi (with a complexity 27% higher than 2D-FFT and 46% lower than 2D-FFTi), and closely approaching the CRLB in delay and Doppler estimation. It also outperforms the other two techniques, in terms of achievable communication rate, when combined with the PBMPA detector. In conclusion, SETFS modulation offers a promising trade-off between SE, complexity, and sensing performance, making it an potential candidate for future 6G ISAC systems. Our ongoing research activities in this field focus on the application of SETFS to more complex ISAC scenarios. APPENDIX A. Gradient and Hessian of the LL function In this appendix, the derivation of the gradient (29) and of the Hessian matrix (30) for the cost function L(26) is sketched. Before doing so, some basic properties of the 6Complexity increases approximately according to a factor NMMc. −10 −5 0 5 10 15 20 25 30 35 10−4 10−3 10−2 10−1 100 SNR (dB) RMSEτ(µs) −10 −5 0 5 10 15 20 25 30 35 100 101 102 103 SNR (dB) RMSEν(kHz) 2D-FFT 2D-FFTi NESCA CRLB Fig. 2. Root mean square error performance achieved by 2D-FFT, the 2D-FFTi and NESCA in a) delay and b) Doppler shift estimation in the cases of SETFS with (β(FD), β(TD)) = (0.9,0.9) (solid lines) and OTFS-DCP (dashed lines), for a variable SNR ∈[−10,30] dB; the CLRBs are also shown. −10 −5 0 5 10 15 20 25 30 35 0.5 1 1.5 2 2.5 SNR (dB) Ceff , bit/s/Hz 2D-FFT + PBMPA (SETFS) 2D-FFTi + PBMPA (SETFS) NESCA + PBMPA (SETFS) 2D-FFT + LLR (OTFS) 2D-FFTi + LLR (OTFS) NESCA + LLR (OTFS) Shannon limit Fig. 3. Effective capacity, Ceff , achieved by the PBMPA detector for SETFS with (β(FD), β(TD)) = (0.9,0.9) (solid lines) when the MMSE equalizer employs the reconstructed channel generated through the estimates of the 2F-FFT, 2D-FFTi and NESCA. The performance offered by a LLR detector when OTFS-DCP is employed (i.e., β(TD) =β(FD) = 1) is also shown as a reference (dashed lines); moreover, Shannon capacity is included as it provides an upper bound. The SNR ∈[−10,30] dB is considered. Dirac vectors as well as of the stacking operations are illustrated. To this aim, we define: a) the L-dimensional vectors x≜[x0, x1, ..., xL−1]T,y≜[y0, y1, ..., yL−1]Tand z≜[z0, z1, ..., zL−1]T(the lth element is denoted xl,yland zl, respectively); b) the L×Lmatrix A≜[al,l′](its (l, l′)th element is denoted al,l′); c) the L-dimensional vector δL,l, called Dirac vector, having its lth element equal to unity and its remaining (L−1) elements equal to zero. Given these definitions, the following properties can be proved (here, Pi indicates the ith property): •P1:xTδL,l =xl •P2:δL,l AδL,l′=al,l′ •P3:ifzl=xlylfor any l, then z=x⊙y •P4:ifal,l′=xlyl′, for any land l′, then A=xyT Our derivations start with substituting the expression of h(16) in the RHS of (26) and evaluating the partial derivative of the resulting expression with respect to the scalar frequency ˜ Fτl (with l= 0,1, ..., L −1); this yields ∂L ∂˜ fτl =−2ℜn˜ aHδL,l δT L,l ˙ BH ˜ fτˆ h−B˜ ao, (37) where ˙ BH ˜ fτis an (MN)×Lmatrix, whose lth column contains the partial derivative of the lth column of B(see (18) with (fτ,fν)=(˜ fτ,˜ fν)), with respect to ˜ fτl. The partial derivative of L(26) with respect to the scalar frequency ˜ fνl(with l= 0,1, ..., L −1) can be computed in a similar way; this results in ∂L ∂˜ fνl =−2ℜn˜ aHδL,l δT L,l ˙ BH ˜ fνˆ h−B˜ ao, (38) where ˙ B˜ fνis an (MN)×Lmatrix, whose lth column contains the partial derivative of the lth column of B(see (18) with (fτ,fν) = (˜ fτ,˜ fν)) with respect to ˜ fνl. Given the partial derivatives in (37) and (38), the gradient vector ∇L(˜ a,˜ f) (having length equal to 2L) can be obtained, since ∇L(˜ a,˜ f)≜"∂L ∂˜ fτ0 ,∂L ∂˜ fτ1 , ..., ∂L ∂˜ fτL−1 , ∂L ∂˜ fν0 ,∂L ∂˜ fν1 , ..., ∂L ∂˜ fνL−1#T . (39) Exploiting the property P1 in both (37) and (38) and the property P3 in (39) easily leads to (29). The evaluation of the (2L)×(2L)Hessian matrix ¨ HL(˜ a,˜ f) requires computing the partial derivatives of L(26) with respect to the couple of variables (˜ fτl,˜ fτl′),(˜ fτl,˜ fνl′), (˜ fνl,˜ fτl′)and (˜ fνl,˜ fνl′). Following the same line of reasoning as the one illustrated for the evaluation of the gradient of L, it can be easily shown that ∂2L ∂˜ fτl∂˜ fτl′ =−2ℜn˜ aHδL,l δT L,l′¨ BH ˜ fτ,˜ fτ(ˆ h−B˜ a)o + 2ℜn˜ aHδL,l δT L,l ˙ BH ˜ fτ ˙ B˜ fτδL,l′δT L,l′˜ ao, (40) ∂2L ∂˜ fτl∂˜ fνl′ =−2ℜn˜ aHδL,l δT L,l′¨ BH ˜ fτ,˜ fν(ˆ h−B˜ a)o + 2ℜn˜ aHδL,l δT L,l ˙ BH ˜ fτ ˙ B˜ fνδL,l′δT L,l′˜ ao, (41) ∂2L ∂˜ fνl∂˜ fτl′ =−2ℜn˜ aHδL,l δT L,l′¨ BH ˜ fν,˜ fτ(ˆ h−B˜ a)o + 2ℜn˜ aHδL,l δT L,l ˙ BH ˜ fν ˙ B˜ fτδL,l′δT L,l′˜ ao(42) and ∂2L ∂˜ fνl∂˜ fνl′ =−2ℜn˜ aHδL,l δT L,l′¨ BH ˜ fν,˜ fν(ˆ h−B˜ a)o + 2ℜn˜ aHδL,l δT L,l ˙ BH ˜ fν ˙ B˜ fνδL,l′δT L,l′˜ ao. (43) Given equations (40)-(43), the (2L)×(2L)Hessian matrix ¨ HL(˜ a,˜ f)can be easily put in the form (30), by exploiting: 1) P1 for all the terms ˜ aHδL,l and δT L,l′˜ ain (40)-(43); 2) P1 for the terms δT L,l′¨ BH ˜ fτ,˜ fτ y,δT L,l′¨ BH ˜ fτ,˜ fν y, δT L,l′¨ BH ˜ fν,˜ fτ yand δT L,l′¨ BH ˜ fν,˜ fν y, in (40), (41), (42) and (43), respectively (here, y≜ˆ h−B˜ a); 3) P2 for the terms δT L,l ˙ BH ˜ fτ ˙ B˜ fτδL,l′,δT L,l ˙ BH ˜ fτ ˙ B˜ fνδL,l′, δT L,l ˙ BH ˜ fν ˙ B˜ fτδL,l′and δT L,l ˙ BH ˜ fν ˙ B˜ fνδL,l′, in (40), (41), (42) and (43), respectively. 4) P3 and P4 in the resulting expressions. REFERENCES [1] Z. Wei et al., “Integrated Sensing and Communication Signals Toward 5G-A and 6G: A Survey,” IEEE Internet Things J., vol. 10, no. 13, pp. 11 068–11 092, Jul. 2023. [2] F.-M. Han and X.-D. Zhang, “Wireless multicarrier digital transmission via weyl-heisenberg frames over time-frequency dispersive channels,” IEEE Trans. Commun., vol. 57, no. 6, pp. 1721–1733, Jun. 2009. [3] 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. [4] 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. [5] M. Mirabella, P. Di Viesti, and G. M. Vitetta, “On the Use of a TwoDimensional Cyclic Prefix in OTFS Modulation and Its Implications,” IEEE Open J. Commun. Soc., vol. 5, pp. 3340–3367, May 2024. [6] ——, “On the Use of a Double Cyclic Prefix in Orthogonal TimeFrequency Space Modulation for Communication and Sensing,” in 2024 IEEE 25th International Workshop on Signal Processing Advances in Wireless Communications (SPAWC), 2024, pp. 726–730. [7] I. Darwazeh and M. Rodrigues, “A Spectrally Efficient Frequency Division Multiplexing Based Communications System,” in Proceedings of the 8th International OFDM-Workshop (InOWo’03), Sep. 2003. [8] Y. Ma et al., “A Low-Complexity Receiver for Multicarrier Faster-ThanNyquist Signaling Over Frequency Selective Channels,” IEEE Commun. Lett., vol. 24, no. 1, pp. 81–85, Jan. 2020. [9] M. Mirabella, P. Di Viesti, and G. M. Vitetta, “A Novel Message Passing Algorithm for Soft-Output Detection in Faster-than-Nyquist Multicarrier Systems,” in 2025 IEEE 26th Int. Workshop on Signal Process. and Artificial Intell. for Wireless Commun. (SPAWC), 2025, pp. 1–5. [10] X. Wang et al., “Multi-Carrier Faster-Than-Nyquist Signaling for OTFS Systems,” IEEE Trans. Veh. Technol., pp. 1–15, Jul. 2025. [11] Z. Hong and S. Sugiura, “Reduced-Complexity MIMO Faster-ThanNyquist Signaling Transceiver for OTFS Modulation,” IEEE Open J. Commun. Soc., vol. 6, pp. 6126–6141, 2025. [12] T. Thaj and E. Viterbo, “Low Complexity Iterative Rake Decision Feedback Equalizer for Zero-Padded OTFS Systems,” IEEE Trans. Veh. Technol., vol. 69, no. 12, pp. 15 606–15 622, Dec. 2020. [13] M. Mirabella, P. Di Viesti, and G. M. Vitetta, “Deterministic Algorithms for Four-Dimensional Imaging in Colocated MIMO OFDM-Based Radar Systems,” IEEE Open J. Commun. Soc., vol. 4, pp. 1516–1543, Jul. 2023. [14] E. Chiavaccini and G. M. Vitetta, “GQR models for multipath Rayleigh fading channels,” IEEE J. Sel. Areas Commun., vol. 19, no. 6, pp. 1009– 1018, Jun. 2001.