1
Real-Time and Scalable Zak-OTFS Receiver Processing on GPUs
arXiv:2604.02266v1 [eess.SP] 2 Apr 2026
Junyao Zheng∗ , Chung-Hsuan Tung∗ , Yuncheng Yao, Nishant Mehrotra, Sandesh Mattu, Zhenzhou Qi, Danyang Zhuo, Robert Calderbank, and Tingjun Chen Abstract—Orthogonal time frequency space (OTFS) modulation offers superior robustness to high-mobility channels compared to conventional orthogonal frequency-division multiplexing (OFDM) waveforms. However, its explicit delay-Doppler (DD) domain representation incurs substantial signal processing complexity, especially with increased DD domain grid sizes. To address this challenge, we present a scalable, real-time Zak-OTFS receiver architecture on GPUs through hardware–algorithm co-design that exploits DD-domain channel sparsity. Our design leverages compact matrix operations for key processing stages, a branchless iterative equalizer, and a structured sparse channel matrix of the DD domain channel matrix to significantly reduce computational and memory overhead. These optimizations enable low-latency processing that consistently meets the 99.9-th percentile realtime processing deadline. The proposed system achieves up to 906.52 Mbps throughput with a DD grid size of (16384,32) using 16QAM modulation over 245.76 MHz bandwidth. Extensive evaluations under a Vehicular-A channel model demonstrate strong scalability and robust performance across CPU (Intel Xeon) and multiple GPU platforms (NVIDIA Jetson Orin, RTX 6000 Ada, A100, and H200), highlighting the effectiveness of compute-aware Zak-OTFS receiver design for next-generation (NextG) high-mobility communication systems. Index Terms—Zak-OTFS, real-time signal processing, delayDoppler communication, GPUs.
I. I NTRODUCTION RTHOGONAL time frequency space (OTFS) modulation has received significant attention in NextG networks due to its robustness on channel Doppler shift compared to traditional orthogonal frequency-division multiplexing (OFDM) modulation. OTFS was first proposed as multi-carrier OTFS (MC-OTFS) [1] by adding an inverse symplectic finite Fourier transform (ISFFT) precoder to existing OFDM architectures for compatibility, but it suffers from spectral efficiency loss due to the cyclic prefix (CP). In contrast, Zak-OTFS [2], [3] directly uses the Zak transform to map the delay-Doppler (DD) domain symbols to time-domain pulse-train waveforms, i.e., pulsones, that span the entire frame duration, enhancing channel predictability and enabling more structured processing. Since OTFS modulates information symbols as a 2D grid in the DD domain, its baseband digital signal processing (DSP) involves matrix operations that incur high computational complexity, especially with increased DD grid sizes. Hence, recent research efforts have focused on developing low-complexity algorithms and practical OTFS systems, as summarized in Table I. Several studies have investigated
O
This work was supported in part by NSF under grants CNS-2211944, CCF2342690, ECCS-2434131, CNS-2443137, CNS-2450567, and OAC-2503010, the AFOSR under grants FA 8750-20-2-0504 and FA 9550-23-1-0249, and an NVIDIA Academic Grant. ∗ indicates equal contribution. J. Zheng, C.-H. Tung, N. Mehrotra, S. Mattu, Z. Qi, R. Calderbank, and T. Chen are with the Department of Electrical and Computer Engineering, Duke University, NC 27708, USA (email: {junyao.zheng, chunghsuan.tung, tingjun.chen}@duke.edu). Y. Yao and D. Zhuo are with the Department of Computer Science, Duke University, NC 27708, USA.
hardware implementations of OTFS modulation and demodulation using software-defined radios (SDRs) [4]–[7], graphics processing units (GPUs) [8], or field-programmable gate arrays (FPGAs) [9]–[13]. These implementations typically focus on demonstrating the functional feasibility of OTFS systems. Despite extensive work on OTFS algorithms and system prototypes, the scalability and real-time processing perspectives of OTFS systems remain largely underexplored. Most prior evaluations consider small-to-moderate DD grid sizes [4], [5], [7], [14], [15], ranging from tens to a few hundred bins in the delay (M ) and Doppler (N ) domains, primarily due to the high complexity of OTFS detection algorithms. As a result, OTFS signal processing latency grows rapidly and becomes prohibitive with large grid sizes (e.g., M or N on the order of several to more than ten thousand). In addition, while several proof-of-concept prototypes have demonstrated functional OTFS systems [4]–[7], [14], performing the full OTFS receiver processing in real time at frame rate–i.e., the DSP of a frame should be completed within the frame duration with a high satisfaction rate–remains challenging. Achieving such real-time processing is challenging for large DD grids, where operations such as channel estimation and equalization involve high-dimensional arithmetic computations. These computational characteristics make GPUs a promising platform for real-time OTFS DSP. Modern GPU architectures are optimized for high-throughput parallel computation and provide specialized hardware and software support for operations such as fast Fourier transforms, matrixvector multiplications (MVM), and general matrix multiplication (GEMM). Importantly, GPUs also provide extremely high on-package memory bandwidth, which is critical when manipulating large DD tensors and channel operators during equalization. As the OTFS grid size increases, the receiver workload increasingly resembles large-scale structured linear algebra, a regime in which GPUs significantly outperform traditional CPUs due to their massive parallelism and memory throughput. At the same time, GPU programmability enables rapid exploration of evolving OTFS receiver algorithms, which is particularly valuable given the active research landscape. However, fully exploiting GPU-based OTFS processing requires careful hardware-algorithm co-design. In this paper, we design and implement a scalable, GPUbased real-time Zak-OTFS signal processing system, guided by three key insights arising from the unique coupling between the Zak-OTFS receiver DSP and GPU processing. First, many OTFS receiver DSP stages involve compact matrix operations with coefficients that depend solely on system parameters such as grid size and index mappings. Examples include phase factors associated with the Zak transform and channel estimation in the DD domain. Precomputing these coefficients eliminates repeated transcendental or indexing-heavy computations, enabling efficient matrix operations during runtime.
2
TABLE I R ECENT OTFS PROCESSING SYSTEMS WITH FULL RECEIVER PIPELINE . Compute Platform OTFS Formulation DD Grid Size, (M, N ) Bandwidth, B (MHz) Frame Duration, T (ms) Equalizer Data Rate Real-Time Processing
[4] (2019) CPU MC-OTFS (32, 32) 25 1.28 MP 3.20 Mbps No
[8] (2023) GPU (A100) MC-OTFS (156, 60) 5 1.000 MP 6.00 Mbps Yes
[9] (2025) FPGA (Alveo U250) MC-OTFS (2048, 8) 40 0.461 MRC 280.32 Mbps Yes
Second, the wireless channel typically exhibits a limited number of paths, resulting in strong sparsity in the DD domain effective channel, heff ∈ CM ×N . While the DD domain channel matrix, Hdd , can be written as an M N × M N matrix, it is directly derived from heff with inherent structured sparsity. Efficient OTFS receiver implementations should therefore avoid explicit construction of dense matrices and instead exploit representations that are aware of such structured sparsity (SS). By leveraging MVM with sparse matrices for equalization, the computational complexity and memory footprint can be significantly reduced while aligning with GPU computation. Third, the computational efficiency of a GPU is significantly undermined by branch conditions, e.g., if-else checks. However, many low-complexity equalizers rely on iteration and conditional termination. Hence, a branchless iteration strategy, rather than dynamic early termination, is necessary. To address this, our implementation determines the optimal number of iterations through comprehensive profiling. By replacing runtime branching with this profile-based approach, the receiver ensures a seamless processing pipeline with peak utilization while maintaining the desired bit-error rate (BER). We implement the real-time Zak-OTFS processing system in Python with PyTorch and customized Triton kernels for GPU, along with a C++ baseline for CPU. Through comprehensive evaluations across four NVIDIA GPU (Jetson Orin, RTX 6000 Ada, A100, and H200) and an Intel CPU (Xeon 6348) platforms, our design meets the 99.9-th percentile processing latency deadline of 2.13 ms for the DD domain grid size of up to (16384, 32), yielding a throughput of 906.52 Mbps with 16QAM modulation and 245.76 MHz bandwidth. In summary, this paper makes the following key contributions: • We leverage compact matrix operation, exploit the sparsity of the DD domain channel, and propose a branchless low-complexity equalizer for efficient Zak-OTFS receiver processing on GPUs; • We implement a practical, scalable computing system for real-time OTFS signal processing that supports an OTFS DD grid size up to (16384, 32); • We conduct a comprehensive evaluation in terms of end-to-end processing latency and BER across various computing platforms and wireless channel conditions. To the best of our knowledge, this is the first implementation that (i) investigates and combines compute hardware- and channel sparsity-awareness for OTFS baseband processing; and (ii) achieves real-time OTFS processing for DD domain grid sizes of up to (M, N ) = (16384, 32). The code, bench-
[10] (2025) CPU MC-OTFS (2048, 512) 300 3.495 MRC 1.80 Gbps No
[14] (2025) CPU Zak-OTFS (31, 37) 0.93 1.233 LMMSE 1.86 Mbps No
[15] (2026) CPU Zak-OTFS (31, 37) 0.93 1.233 CGA 1.86 Mbps No
This Work GPU (H200) Zak-OTFS (16384, 32) 245.76 1.067 SS-CGA 906.52 Mbps Yes
marks, and profiling results are open-sourced at [16]. II. R ELATED W ORK Complexity of OTFS processing. Since OTFS operates on a 2D resource grid in the DD domain and each impulse in the grid is cross-dependent, channel estimation and equalization are of high complexity. Point pilots [3] and spread pilots [17] are two pilot designs in Zak-OTFS, where channel estimation based on the point pilot involves simple elementwise computation in O(M N ), and that on the spread pilot relies on the cross-ambiguity function in O(M 2 N 2 ). Reducing the complexity of equalization has received significant attention [18]–[20]. These methods perform equalization in near-linear complexity (O(M N log M N )) by exploiting the block circulant structure of the modulation. Iterative equalizers such as message-passing (MP) [21], maximum ratio combining (MRC) [22], and conjugate gradient algorithm (CGA) [15] can further reduce complexity to linear per-iteration. Practical OTFS system implementations. Table I compares this work against the related literature across different compute platforms, link parameters, the equalizer algorithm, system throughput, and the real-time capability. The CPU has been used for algorithm prototyping and complexity-reduction studies [4], [14], [15], [22]–[24], e.g., the CGA equalizer [15] is evaluated on a (31, 37) DD domain grid on the CPU. Similarly, the MP equalizer has been implemented on a (32, 32) grid [4], while the MRC equalizer has been evaluated on a larger (2048, 512) grid [10]. However, these CPU-based implementations do not report their real-time processing capability. As the theory of OTFS matures, efforts to implement OTFS processing across a variety of computing platforms aim to enable practical real-time operation under different bandwidth and grid configurations [8]–[10], [14]. Many implementations leverage the FPGA platforms [9], [11]–[13] for their ability to provide deterministic latency and energy-efficient hardware acceleration. For example, a real-time OTFS receiver on an Alveo U250 FPGA with an MRC equalizer achieves a throughput of 280.32 Mbps with a (2048, 8) grid [9]. In contrast, GPU-based real-time OTFS systems remain relatively limited despite their potential for large-scale parallel processing, with only a proof-of-concept demonstration on an NVIDIA A100 using an MP equalizer for a (156, 60) grid [8]. Hence, further investigation leveraging GPU’s massive parallelism for realtime processing of high-bandwidth OTFS signal is in demand. Virtualized RAN (vRAN) systems leveraging GPUs. Recent literature reveals interest in using GPUs to support vRAN
3
Transmitter Processing
Delay-Doppler Domain Pilot Frame
Bits Generation
Channel
Channel Estimation " eff 𝐡 (13)
Data Frame
Modulation
Pilot Generation
IDZT (6)
Pilot
Symbol Mapping
IDZT (6)
Data
Packet (Concat.)
Veh-A + AWGN
Depacket
Pilot
DZT (7)
Data
DZT (7)
$ dd 𝐇 (15)
Equal.
Receiver Processing
Demod. (18)
BER Cal.
Fig. 1. Zak-OTFS signal processing pipeline based on discrete Zak transform (DZT). A point-pilot frame and a data frame are concatenated to form a packet. The simulated channel considers Veh-A [30] for paths with delay and Doppler shifts, and AWGN for SNR adjustment. The receiver pipeline stages include DZT, channel estimation, equalization, and demodulation, which will be discussed in Section III-B.
systems for compute-intensive baseband processing of large matrices. For example, while Agora [25] provides a MIMO baseband processing framework for 5G using CPU and Savannah [26] accelerates it to stringent frequency range 2 (FR2) deadlines, MegaStation [27] demonstrates a reduced tail latency and improved throughput on massive MIMO by migrating the system onto a supercomputer with pooled GPUs. In addition, DecodeX [28] benchmarks LDPC decoding workloads, which are considered the most compute-intensive stage of the receiver pipeline, across CPU/GPU/ASIC, and demonstrates the potential benefit of a GPU-resident pipeline for efficient resource sharing in vRAN. Similarly, X5G [29] integrates NVIDIA Aerial with OpenAirInterface in an end-toend private 5G O-RAN testbed, where GPUs accelerate PHYlayer functions and achieve multi-Gbps downlink throughput. The success of GPU-based OFDM processing in 5G systems reveals its potential to accelerate the processing pipeline for advanced waveforms such as OTFS in NextG networks. III. P RELIMINARIES In this paper, we consider the Zak-OTFS system, shown in Fig. 1, which adopts a twisted convolution representation in the DD domain. Key receiver DSP stages include the Zak transform, effective channel computation (heff ), DD channel construction (Hdd ), equalization, and demodulation, many of which introduce large structured linear algebra workloads. A. Zak-OTFS Modulation and Parameters In Zak-OTFS, M N information symbols, Xdd [k, l], are multiplexed over a frame consisting of an M ×N -dimensional DD grid with M delay bins and N Doppler bins. Let M ≜ {0, 1, . . . , M −1} and N ≜ {0, 1, . . . , N −1}, this DD domain information grid can be denoted by (k · ∆τ, l · ∆ν) , k ∈ M, l ∈ N ,
(1)
with a delay and Doppler period of M · ∆τ and N · ∆ν, respectively, where ∆τ and ∆ν are the delay and Doppler resolutions. ∆τ and ∆ν are given by ∆τ = 1/B and ∆ν = 1/T,
(2)
where B is the signal bandwidth (which determines the sampling rate) and T is the frame duration, B = M · ∆f and T = N/∆f ,
(3)
with ∆f the frequency spacing. Hence, B = M (N · ∆ν) and T = N (M · ∆τ ).
(4)
The DD domain grid size (M × N ) as well as the delay and Doppler resolutions (∆τ and ∆ν) are designed so that (M ∆τ ) · (N ∆ν) = 1 and BT = M N.
(5)
B. Zak-OTFS Channel and Processing Pipeline Fig. 1 illustrates the Zak-OTFS TX/RX processing pipeline and channel simulation considered in this work. This section provides the formulation of each critical DSP stage, with a focus on the receiver pipeline and stage complexity. After data bits are generated and modulated into the delay– Doppler (DD) domain, each symbol of a data packet will be mapped onto a 2D grid and transformed into the time domain using the inverse discrete Zak transform (IDZT) (6). In contrast, the pilot frame is generated directly in the DD domain and converted to the time domain using IDZT. A packet contains one pilot frame and a number of data frames, depending on the channel coherence time [3]. Vehicular-A (Veh-A) [30] model is applied to simulate a channel with multiple path delays and Doppler shifts, and we use AWGN on top of Veh-A as the channel noise (details in Section V). The receiver pipeline consists of both pilot- and dataframe processing, including the discrete Zak transform (DZT) (7), channel estimation including obtaining the channel rebeff (13) and constructing the channel matrix H b dd , sponse grid h followed by equalization (with LMMSE, MRC, and CGA equalizers) and demodulation (18). Finally, the demodulated bits are compared with the ground-truth bits to calculate BER. 1) Zak Transform: OTFS waveform constructs the resource grid in the DD domain, and its conversion to (from) the timedomain waveform is via the Zak (inverse Zak) transform. Inverse discrete Zak transform (IDZT). After the pilot and data frames are generated as DD grids, the transmitter converts the DD domain signal matrix Xdd ∈ CM ×N to a time-domain signal vector x ∈ CM N using IDZT [31, Eq. (13)], N −1 ⌊i/M ⌋ 1 X (6) x[i] = √ Xdd [⟨i⟩M , l] · ej2π N l , ∀i ∈ Q, N l=0 where ⟨·⟩M denotes the modulo M operation and Q ≜ {0, 1, . . . , M N − 1}. x is transmitted through the channel (discussed in Section III-B2), and results in the received timedomain signal vector y ∈ CM N . Discrete Zak transform (DZT). The receiver converts y ∈ CM N back into the DD domain as matrix Ydd ∈ CM ×N using
4
DZT [31, Eq. (1)] for all k ∈ M and l ∈ N : N −1 il 1 X Ydd [k, l] = √ y[k + iM ] · e−j2π N . N i=0
beff [3, Eq. (38)] using from h (7)
Note that both DZT and IDZT reshape the signal vector to/from matrix when transforming the signal domain. 2) Effective Channel in the DD Domain: Zak-OTFS models the time-varying channel as the DD domain response P −1 X h(τ, ν) = hp δ(τ − τp )δ(ν − νp ), (8) p=0
where P denotes the number of resolvable paths in the DD domain, hp is the complex-valued channel gain of the p-th path, δ(·) denotes the Dirac delta function, and τp and νp are the delay and Doppler shift associated with the p-th path, respectively. By associating τp and νp with the DD grid τp = kp /B = kp ∆τ ,
νp = lp /T = lp ∆ν,
(9)
the discrete form of the channel model (8) is denoted as the effective channel, heff = [heff [k, l]] ∈ CM ×N , where heff [k, l] =
P −1 X
hp δ[k − kp ]δ[l − lp ].
(10)
p=0
where δ[·] denotes the Kronecker delta function, and kp and lp are the delay and Doppler shift indices on the DD grid for the p-th path, respectively. The receiver samples the continuous signal undergoing (10) at each time point t = m B , m ∈ M, and results in the discrete time-domain signal y that P MX N −1 X i j2πτp B −τp y[i] = hp x[n]δ[ i−n − τ ]e + w[i], (11) p B p=1
n=0
where δ[·] denotes the Kronecker delta function and w[·] is the additive noise. Converting (11) to DD domain, each element Ydd [k, l] of the received OTFS frame Ydd is expressed as a discrete twisted convolution ∗σd 1 between the quasi-periodic extension of Xdd [k, l]2 and the effective channel heff [k, l], i.e., Ydd [k, l] = heff [k, l] ∗σd Xdd [k, l] + Wdd [k, l], ∀k, l,
(12)
where Wdd is the DD domain additive noise. 3) Channel Estimation: To capture the input/output relationship of the channel, a frame containing a point pilot, i.e., N a single impulse located at (K0 , L0 ) ≜ M , in the DD 2 2 domain, is transmitted. The estimated effective DD domain beff , can be obtained from the received pilot frame: channel, h −jπl b heff [k, l] = Ydd k + M , l + N · e N , (13) 2
2
M N N where − M 2 ≤ k < 2 and − 2 ≤ l < 2 . The DD domain channel relationship (12) can also be described as
ydd = Hdd xdd + wdd ,
(14)
with xdd , ydd , and wdd the vectorized Xdd , Ydd , and Wdd . In b dd ∈ CM N ×M N is obtained this case, the channel matrix H 1 a[k, l] ∗ σd b[k, l] ≜
P
2 X [k + nM, l + mN ] dd
m, n ∈ Z
j2π
k′ ,l′ ∈Z
′
′
a[k − k′ , l − l′ ]b[k′ , l′ ]e M N k (l−l ) 2π
= Xdd [k, l]ej N nl , with k ∈ M, l ∈ N ,
b dd [k ′ N + l′ , kN + l] = H
X
ej2π
(l′ −l−mN )(k+nM ) MN
(15)
n,m∈Z ′
′
×b heff [k − k − nM, l − l − mN ]e ′
j2π N nl
,
′
where k, k ∈ M and l, l ∈ N . Note that only one pair of (k ′ − k − nM, l′ − l − mN ) falls within the valid indices of beff . The channel H b dd depicts the pair-wise input/output relah tionship across the DD domain resource grid points, leading to a large matrix sized M N × M N . 4) Equalization: OTFS systems typically rely on LMMSE equalization to jointly mitigate DD interference and noise b dd , which is computationby inverting the channel matrix, H ally expensive as the DD grid size (M, N ) increases. Lowcomplexity equalizers have been proposed to avoid direct matrix inversion, such as maximum ratio combining (MRC) [22] and the iterative conjugate gradient algorithm (CGA) [15]. Linear minimum mean square error (LMMSE). The LMMSE equalizer solves (14) with direct matrix inversion, −1 H H b b b dd ydd ≜ HLMMSE ydd , (16) bdd = Hdd Hdd + Rn x H where Rn = (SNRlinear )−1 I denotes the noise covariance matrix. It is well-known that LMMSE has a computational complexity of O(M 3 N 3 ). In practice, HLMMSE can also be calculated by solving H bH b b H (17) dd = Hdd Hdd + Rn HLMMSE using factorization techniques such as the Cholesky or lowerupper (LU) decomposition at the same complexity. Maximum ratio combining (MRC). The MRC equalizer [22] operates on symbol vectors, where each vector consists of symbols sharing the same delay index across all Doppler indices. Each transmitted symbol vector contributes to multiple received vectors due to DD coupling. MRC extracts these components across delay branches, cancels interference using current estimates, and coherently combines them to maximize SNR. The updated symbol vector is then mapped to the nearest constellation point and refined iteratively in a decisionfeedback manner. This process requires hard demodulation at each iteration, introducing additional computational overhead. The overall complexity is O(ΞP M N ), where Ξ and P denote beff , respectively. the iteration count and the number of paths in h Conjugate gradient algorithm (CGA). To avoid explicit matrix inversion on Hdd , CGA [15] solves (16) in an iterative bdd = 0, CGA manner. Starting from an initial estimate of x updates the solution using residual-based search directions in each step and terminates when the residual norm falls below a predefined threshold or a maximum number of iterations is reached. Since each iteration in CGA requires only MVMs [15, Algo. 1], it has a computational complexity of O(ΞP M N ) but with lower coefficient compared to MRC. 5) Demodulation: We consider hard demodulation given by x̆dd = arg min ∥xdd − s∥2 , s∈S
∀xdd ∈ xdd ,
(18)
and demap x̆dd into their corresponding binary sequences, from which the BER can be obtained.
5
CPU: Xeon 6348 (16 threads)
GPU: A100
103
2
2
10
101 100 10−1 10−2 10−3 10
−4
Theory: (Nd3) Measured: ≈ Nd2
21 23 25 27 29 211 213 Matrix Dimension, Nd
(a) Matrix Inversion
Latency (ms)
10 Latency (ms)
GPU: RTX 6000
103
101
GPU: H200
(Nd2)
100
matrix inversion on GPU platforms. In particular, MVM is at least one order of magnitude faster than matrix inversion, and up to two orders of magnitude faster with sufficiently large matrices (e.g., Nd ≥ 256).
10−1 10−2
IV. S YSTEM D ESIGN
10−3 10
−4
21 23 25 27 29 211 213 Matrix Dimension, Nd
(b) MVM
Fig. 2. Measured compute latency and scalability on A ∈ CNd ×Nd across CPUand GPU platforms, including (a) matrix inversion, and (b) matrix-vector multiplication (MVM), with varying matrix dimension Nd . GPU-based matrix operations achieve lower latency and better scalability for Nd > 256.
C. Scalability of Matrix Operations in Practical Systems While low-complexity algorithms, such as MRC [22] and CGA [15] equalizers, are a key to enabling real-time and scalable OTFS signal processing, the available computing resources and their varying capability are often not considered in prior work. Specifically, the arithmetic operations may exhibit scalability that differs from theory due to hardware specifications and implementations. To understand how real systems scale with respect to matrix dimensions, we profile the execution latency of matrix operations using a highly optimized BLAS (Basic Linear Algebra Subprograms) library on both CPU and GPU platforms. Fig. 2 illustrates the execution latency of two operations with varying matrix dimensions, Nd : (i) Matrix inversion of A ∈ CNd ×Nd , A−1 , and (ii) Matrix-vector multiplication (MVM) of y = Ax with A ∈ CNd ×Nd and x, y ∈ CNd . We consider two compute platforms (for more details see Section V): (i) server-grade Intel CPU (Xeon 6348) using Intel Math Kernel Library (MKL), and (ii) different NVIDIA GPUs, including RTX 6000 Ada, A100, and H200, using cuBLAS. We configure MKL to use 16 threads to enable parallel processing on the CPU, and each profiling is executed 10K times Fig. 2(a) shows the mean latency of matrix inversion on CPU/GPU platforms with varying Nd . The measured complexity of inverting A ∈ CNd ×Nd scales in O(Nd 2 ) for both CPU and GPU, with the GPU having a larger overhead when Nd < 128, despite the theoretical complexity of O(Nd 3 ). The result suggests that the processing pipeline design should consider practical scalability within the computing system. Fig. 2(b) shows the mean latency of MVM on CPU/GPU platforms with varying Nd . The measurements demonstrates that GPU-based processing maintains nearly O(Nd ) scalability over MVM for Nd ≤ 1024, whereas the CPU latency turns upwards at Nd = 32. On the other hand, GPU has a larger constant overhead than CPU in their flat regions, leading to an intersection around 32 < Nd < 64. Even at moderate bandwidths in modern wireless systems, the resulting problem size already falls into an operating region where GPU-based MVM becomes advantageous; for example, M = 1024 corresponds to a 30.72 MHz bandwidth with ∆f = 30 kHz. The latency comparison between matrix inversion and MVM further supports using low-complexity algorithms (e.g., MRC and CGA) that rely on iterative MVM to avoid explicit
Drawing on these insights, the proposed Zak-OTFS processing pipeline leverages compute-aware optimizations to meet real-time latency requirements while enhancing scalability. Table II summarizes the five major DSP stages in the ZakOTFS receiver which, except for demodulation, require insight into the compute-aware arithmetic design. Our design accelerates these stages by (i) using compact matrix operation, (ii) leveraging channel sparsity, and (iii) designing a computeaware algorithm, as detailed next in this section. A. Compact Matrix Operations In the Zak-OTFS signal processing pipeline (see Section III-B), exponential factors are constantly used in each stage, and are invariant across the frame. Hence, precomputing those terms offline and applying them on the fly reduces the overhead of repeated runtime computation. For example, DZT (7) and estimating heff (13) can be decoupled into offline computed assistant terms and runtime compact matrix operations suitable for highly optimized GPU libraries. GEMM for DZT. Specifically, for DZT, we first construct an auxiliary matrix of FFT kernels EZak ∈ CN ×N , given by l′ l
l′ l
EZak [l′ , l] = e−j2π N · ejπl = (−1)l · e−j2π N .
(19)
Then, we reshape the received time-domain signal vector y ∈ CM N into the matrix Y ∈ CM ×N in column-major order. Under these compact matrix representations, DZT can then be performed in one GEMM with complexity O(M N 2 ) given by Ydd = YEZak .
(20)
Although the complexity is higher than that of the FFT-based expression of O(M N log N ), (20) outperforms in speed in the profiled system, as modern GPU arithmetic units are highly optimized for GEMM. beff . Similarly, we create an auxiliary Hadamard product for h M ×N matrix Etwist ∈ C to store the coefficients required for the twist convolution in (13), Etwist [k, l] = e−j2π
K0 (l−L0 ) MN
,
(21)
and the estimation of heff can be performed via the Hadamard (element-wise) product given by beff = Ydd ⊙ Etwist . h (22) with a complexity of only O(M N ). B. Efficient Matrix Operation via Structured Channel Sparsity In OTFS processing, the DD grid dimension of M × N , results in a channel matrix Hdd ∈ CM N ×M N . This dimension explosion prevents existing OTFS processing pipelines from exploring larger grid sizes, as the massive matrix scales incur prohibitive costs in both computational intensity and memory constraints, including bandwidth and storage requirements.
6
TABLE II A RITHMETIC AND IMPLEMENTATION INSIGHTS FOR THE PROPOSED Z AK -OTFS RECEIVER PROCESSING PIPELINE .
15
0
3 7 15 10 Column index
N Doppler bins
Pilot Location (K0, L0)
N Blocks of M Col's
0 7 15
0 1 Doppler l
Memory O(M N + N 2 ) O(M N ) O(P M N ) O(P M N ) O(M N )
N Blocks of M Rows
10 2
0 1 2 3 4 5 6 7
Complexity O(M N 2 ) O(M N ) O(P M N ) O(P M N ) O(M N )
Row Index q
10 1
7
Arithmetic GEMM Hadamard Product Matrix Construction MVM Element-wise Op.
M Delay bins
100
0
Delay k
1.0 0.8 0.6 0.4 0.2 0.0
Magnitude of Hdd
Magnitude
1.0 0.5 0.0 1 bin 012 3 4 r e l Delay 5 6 7 0 pp bin Do
Primary Insight Compact Matrix Operation (Section IV-A, (20)) Compact Matrix Operation (Section IV-A, (22)) Structured Sparsity (SS) (Section IV-B) SS and Compute-Aware CGA (Section IV-C, Algo. 1) Parallel Decision (Trivial)
Row index
Formulation (7) (13) (15) [15] (18)
Magnitude of heff
OTFS RX Stage DZT b heff Calculation b dd Construction H Equalization Demodulation
0
7 15 Column Index r
b dd Magnitudes b dd (a) b heff Magnitudes (b) H (c) Binarized b heff (d) Binarized H b dd magnitudes in (M, N ) = (8, 2). The binarized versions depict the dominant entries using a threshold of θ = 0.12, where any Fig. 3. The b heff and H b dd entry with magnitude below the threshold is set to zero. (a) The significant bins in b heff correspond to the dominant channel paths. (b) The channel matrix H is divided into N × N smaller sub-matrices of dimension M × M . (c) Thresholding allows isolation of the most significant paths in b heff , with each path at (kp , lp ). In this example, two paths are located at (4, 1) in blue and at (0, 0) in orange. The blue path is associated with zero delay and zero Doppler shift, b dd , where the colors associate while the orange path corresponds to a delay shift of 0.017 ms and a Doppler shift of 15 kHz. (d) The resultant patterns in H b dd . For each row q in H b dd , the dominant path p maps to column rp (q) as a function of (kp , lp ) (33). the dominant path in binarized b heff and H
Note that from (10), b heff [k, l] is given by the sum of P terms, which is naturally sparse in the DD domain. Physically, these terms correspond to the discrete delay and Doppler shifts of the P dominant paths of the wireless propagation channel. beff , which reflects Fig. 3(a) shows an example visualization of h b dd is obtained from two dominant paths in the channel. Since H b heff [k, l] using (15) to expand the dimension with no additional b dd is both sparse and structured. information, the resultant H b dd , we propose an By exploiting the structured sparsity of H efficient matrix representation that simultaneously minimizes storage footprint and accelerates MVM. 1) Structure of Channel Matrix: While the strict blockcirculant property of Hdd theoretically holds only under integer b dd retains a quasiDoppler shifts, our observations reveal that H block-circulant structure even in the more general scenarios. b dd are deterministically comThis is because the entries of H b puted by elements in heff , allowing for a predictable mapping of dominant paths. Despite the presence of fractional Doppler b dd consistently leakage, the primary gain distribution of H follows the indexing pattern of a block-circulant structure. Consequently, we can exploit this structured sparsity to conb dd in low complexity. struct and operate H b dd exhibit a block-wise Fig. 3(b) demonstrates how H beff in Fig. 3(a) diagonal-like pattern even if the corresponding h is irregular. Each physical propagation path in the channel b dd with circular shift in induces an independent pattern in H both row and column directions. For example, each dominant path (distinguished by colors) located in Fig. 3(c) leads to a pattern in Fig. 3(d). This block-circulant structure yields N × N blocks, each of size M × M . Within each block, successive rows are cyclic shifts of one another with linearly increasing offsets; a similar pattern holds column-wise.
2) Exploiting Structured Sparsity: Note that the p-th active path corresponding to a pulse at heff [kp , lp ] induces a circular shift of the transmitted symbol grid in Hdd . We use rp (q) to denote the column index of a significant value in the qb dd induced by the dominant path p. Hence, each the row of H row q exists P significant values where each is mapped to b dd p ∈ P = {1, 2, . . . , P }. With the knowledge of rp (q) for H b dd ∈ (to be discussed in Section IV-B4), the MVM using H CM N ×M N and a vector v ∈ CM N b dd v, u ∈ CM N u=H (23) can be approximated via dominant terms by P MX N −1 X b dd [q, r] · v[r] ≈ Dp,q · v[rp (q)], ∀q ∈ Q, H u[q] = r=0
p=1
(24) where Q = {0, 1, . . . , M N − 1} and coefficients are given by b dd [q, rp (q)]. Dp,q ≜ H (25) Hence, instead of summing over all M N columns, one only needs to sum over the P path-induced dominant terms for each row (e.g., the colored entries shown in Fig. 3(d)). The arithmetic complexity is thus reduced from O(M 2 N 2 ) to O(P M N ), which is highly advantageous for P ≪ M N . In the example shown in Fig. 3, the number of multiplications is reduced by a factor of 8× (from 256 to 32). H
b dd requires an inverse mapping that identifies Similarly, H the row index corresponding to a given output column r and M N ×M N bH path p. Multiplying H by a vector v ∈ CM N dd ∈ C can be written as P X ∗ u[r] ≈ Dp,q · v[qp (r)], ∀r ∈ Q. (26) p (r) p=1
7
Moreover, the reduced MVMs in (24) and (26) have the gather-multiply-reduce computational pattern: (i) gather the required input vector entries v[·] using known indices rp (q) and qp (r); (ii) multiply each gathered vector entry by the ∗ corresponding path coefficient Dp,q or Dp,q ; (iii) reduce p (r) across the path dimension p. This structure avoids irregular sparse traversal and eliminates the need for atomic updates, making it particularly suitable for efficient GPU execution. 3) Efficient Data Structure: Section IV-B2 reveals that to b dd (24) or H bH perform MVM using H dd (26) , storing the entire b dd and M N × M N dense matrix in unnecessary. Instead, H H b Hdd can be respectively stored in a structured-sparse form as oP n P ∗ , q (r) , ∀r. (27) {Dp,q , rp (q)}p=1 , ∀q, and Dp,q p p (r) p=1
This structured-sparsity representation requires memory storage with a size of O(P M N ) elements, which is much smaller than two dense matrices in size O(M N ×M N ) as P ≪ M N . 4) Constructing Hdd in Structured-Sparse Representation: We now describe how to construct the structured-sparse (SS) b dd in (27) directly from the channel response h beff without H M N ×M N explicitly forming the dense matrix in C . Since the b dd is a 2D circular shift with phase rotation, construction of H the position, magnitude, and phase shift of the dominant b dd can be computed analytically. elements in H Dominant path determination. In typical wireless channels, the number of dominant propagation paths is small, hence heff is sparse. We identify the active propagation paths by thresholding the channel magnitude: b = {(k, l) : |heff [k, l]| > θ} , Ω (28) where θ is a small threshold used to suppress noise. We represents the delay-bin and Doppler-bin location of the pth path in the DD grid as (kp , lp ) with p ∈ P, and the corresponding complex channel gains are defined as hp ≜ heff [kp , lp ].
(29)
Analytical indices mapping. Eq. (15) depicts the construcb dd from h beff , and a pulse at heff [kp , lp ] caused by tion of H a dominant path p induces a circular shift of the transmitted b dd , the psymbol grid in Hdd . Consequently, for row q of H th active path contributes only one dominant matrix entry at column rp (q) based on (kp , lp ). b dd by a pair of DD coordinates We index the q-th row in H (k, l) as q = lM + k, k ∈ M, l ∈ N . Equivalently, kq = ⟨q⟩M and lq = ⌊q/M ⌋ for any row index q ∈ Q. We denote b dd induced the column index r of the dominant element in H by the p-th dominant element at row q as rp (q), which can similarly be expressed as rp (q) = lM + k with another set of DD coordinates (k, l) for k ∈ M and l ∈ N . Similarly, the column index rp (q) can be expressed with the input DD grid coordinate (kp (q), lp (q)) as rp (q) = lp (q)M + kp (q),
(30)
the output coordinate (lq , kq ) for q = lq M + kq contributes to the coordinate transform as kp (q) = ⟨K0 + kq − kp ⟩M ,
lp (q) = ⟨L0 + lq − lp ⟩N , (31)
where (K0 , L0 ) denotes the pilot reference location in the DD
grid. Note that the delay bin offset dk (p) and the Doppler bin offset dl (p) only depend on the pulse location (kp , lp ) of p-th beff : path and the pilot reference location (K0 , L0 ) in h dk (p) ≜ K0 − kp ,
dl (p) ≜ L0 − lp ,
(32)
so that rp (q) = ⟨lq + dl (p)⟩N · M + ⟨kq + dk (p)⟩M .
(33)
Fig. 3 provides an example with grid dimension (M, N ) = (8, 2) and pilot location (K0 , L0 ) = (4, 1). The channel has two paths located (kp=0 , lp=0 ) = (4, 1) and (kp=1 , lp=1 ) = (0, 0). In the row q = 7 = 0M + 7, we have lq=7 = 0 and kq=7 = 7. Hence, rp=0 (q = 7) = ⟨lq=7 + L0 − lp=0 ⟩N · M + ⟨kq=7 + K0 − kp=0 ⟩M = ⟨0 + 1 − 1⟩2 · 8 + ⟨7 + 4 − 4⟩8 = 7. and rp=1 (q = 7) = ⟨0 + 1 − 0⟩2 · 8 + ⟨7 + 4 − 0⟩8 = 11. bH For MVM using matrix H dd , we instead need, for each output index r, the row index whose forward mapping equals r. We therefore define the inverse mapping qp (r) such that rp (qp (r)) = r.
(34)
Let kr = ⟨r⟩M and lr = ⌊r/M ⌋ denote the DD coordinates corresponding to index r. Furthermore, define the pathdependent offsets δk (p) ≜ ⟨rp (0)⟩M = kp (0),
r (0)
δl (p) ≜ ⌊ pM ⌋ = lp (0), (35)
which correspond to the delay and Doppler shifts induced by path p. The inverse mapping is given by qp (r) = ⟨lr − dl (p)⟩N · M + ⟨kr − dk (p)⟩M .
(36)
Phase correction and coefficient synthesis. In addition to the index mapping and the channel gain, each matrix entry must include the appropriate OTFS phase rotation, i.e., a complex phase compensation factor γp,q which contributes to Dp,q as γp,q = ejϕp,q ,
Dp,q = hp γp,q
(37)
The quasi-periodic phase term ϕp,q for path p and row q in b dd comes from H ϕp,q = M2πN [(lp − L0 )Ap,q + np,q lp,q M ] ,
(38)
where Ap,q is the raw delay index offset and np,q is the wrap count due to circular shift Ap,q ≜ K0 + kq − kp ,
np,q ≜ ⌊Ap,q /M ⌋ .
(39)
5) Discussions on Complexity and Memory Requirements: The proposed structured-sparse construction avoids explicitly b dd ∈ CM N ×M N , forming the dense channel matrix H 2 which requires storing (M N ) complex-valued entries, and a dense MVM has a computational complexity of O((M N )2 ). Instead, following Section IV-B4 requires only O(P M N ) memory and time to construct the matrix in the sparse form (27). Similarly, (24) and (26) have arithmetic complexity O(P M N ) to perform MVM. For example, with (M, N ) = (48, 32), M N = 1,536. Assuming P = 5 paths, the dense representation requires ≈2.4 M complex values, whereas the diagonal-sparse representation stores only 7,680, resulting in a 307.2× reduction in memory requirement. As we scale to larger Zak-OTFS grids, e.g., (M, N ) = (8192, 32), storing the dense channel matrix alone would require 8 × 81922 × 322 Bytes ≈ 549.8 GB memory per frame, which needs to be updated every 1.07 ms and quickly
8
Algorithm 1 Compute-aware Conjugate Gradient Algorithm (CGA)
SNR=0dB
SNR=15dB 100
b dd and H b Hdd , received vector Require: Structured-sparse channel H ydd , noise covariance matrix Rn , iterations Ξ = 10 bdd Ensure: Estimated transmit vector x b Hdd ydd , x bdd ← 0, c ← b, p ← b, cnorm ← ∥c∥22 1: b ← H 2: for ξ = 1 to Ξ do b Hdd (H b dd p) 3: H′dd ← H ▷ Break into two MVMs 4: ap ← H′dd + λp ▷ Rn = λI 5: α ← cnorm /(pH ap ) bdd ← x bdd + αp, c ← c − αap 6: x 7: cnorm ′ ← ∥c∥22 ▷ No conditional early termination 8: p ← c + (cnorm ′ /cnorm )p, cnorm ← cnorm ′ 9: end for bdd 10: return x
102
10 1
101
10 2
becomes a significant memory bottleneck. The key system-level insight is, therefore, that the OTFS channel matrix is not merely sparse but structurally sparse. This structure allows the original dense operator to be replaced by a compact path-indexed lookup representation, preserving the low arithmetic complexity of sparse processing while enabling GPU-friendly memory access patterns that are much more regular than those of a generic sparse matrix format. C. Compute-Aware Conjugate Gradient-based Equalizer Motivated by the practical scalability of the computation system discussed in Section III-C, equalizers using MVM instead of matrix inversion are preferred for GPU-based realtime OTFS processing, e.g., CGA proposed in [15] (see also Section III-B4). However, compute-aware adoption is necessary to translate operator-level scalability in Fig. 2 to the system level. In addition, unlike the original CGA [15], which performs equalization in the frequency domain, we carry out conjugate gradient-based equalization directly in the DD domain. First, transforming between domains (e.g., DD to frequency domain and back) incurs additional computational overhead and latency, as it requires multiplications involving dense transform matrices applied to the vectorized DD transmit signal xdd . In addition, as shown in Section IV-B, b dd exhibits inherent sparsity in the the effective channel H DD domain due to the limited number of propagation paths. Consequently, operating in the DD domain enables us to exploit the channel sparsity while applying CGA efficiently without incurring transformation overhead. Algo. 1 shows our proposed compute-aware CGA with three key features: (i) substitute of GEMM with equivalent MVM, (ii) dimension reduction of operator, and (iii) fixed termination based on offline profiling. Leveraging scalable operators. In the original CGA [15], the bH b Gram matrix H dd · Hdd is computed once ahead of the iterative processes. This incurs a complexity of O(M N 2 ), which can be reduced to O(P 2 M N ) leveraging the structured sparsity in Hdd . On the other hand, our implementation requires only two MVMs (Line 3, Algo. 1), therefore further reducing the computational complexity to O(P M N ). Dimension reduction of noise covariance. Similarly, for the noise covariance matrix Rn = λ · I, computing Rn explicitly
10 3
100 10 1 0
SNR=30dB
BER
103 cnorm
Equalization.
5 10 15 20 25 30 CGA Iteration
(a) cnorm vs. CGA Iterations
10 4 0
5 10 15 20 25 30 CGA Iteration
(b) BER vs. CGA Iterations
Fig. 4. BER and cnorm vs. iterations in CGA for (M, N ) = (128, 32), 16QAM. Across iterations, cnorm exhibits monotonic decay with a diminishing slope, whereas BER behaves more variably, including flat (SNR = 0 dB) and rebound (SNR = 30 dB). In addition, cnorm is exponentially proportional to SNR, while BER does not have this behavior.
incurs dimension increase as Rn ∈ CM N ×M N . If unified noise variance across the DD bins is assumed, Rn can be reduced to a scalar λ = 1/(SNRlinear ) (Line 4, Algo. 1). Profiling-based fixed termination. While the original CGA benefits from loss-based early termination to reduce iteration counts, this strategy is not well-suited to GPU execution. This is because implementing early termination requires synchronization across the entire GPU and often involves transferring intermediate results to the host for convergence evaluation, i.e., checking the squared residual norm to stop the iteration when cnorm < ϵ2 , with ϵ being the threshold. The introduced kernel launch overhead and CPU-GPU synchronization in each iteration can stall the execution pipeline and offset the benefits of reduced arithmetic operations. Instead, we employ an offline profiling-based fixed iterations for CGA based on the OTFS parameters and link conditions (e.g, DD domain grid size, modulation order, and SNR). BER-driven iteration selection. To determine the iteration budget, we adopt a BER-driven offline profiling approach. We define the convergence of BER by monitoring the signed relative BER change r between consecutive CGA iterations: BERξ−1 − BERξ , (40) rξ = max(BERξ−1 , ζ) where ζ > 0 is a small constant to ensure numerical stability. The BER is considered converged at the first iteration ξ ≥ 2 such that rξ < η, where η is a small tolerance. This criterion captures both cases where the BER begins to rise and where the BER improvement becomes negligible. We set η = 10−2 and ζ = 10−12 in all experiments. The iteration budget is then determined as the maximum convergence iteration across SNRs and grid configurations. To justify this BER-driven criterion, we examine the relationship between cnorm and BER across CGA iterations. Fig. 4 plots the cnorm values and BER across the CGA iterations under different SNR scenarios. As expected, cnorm decreases monotonically with the iteration count for all SNR levels, reflecting the convergence behavior of the conjugate gradient solver toward the least-squares solution. However, the BER curves exhibit significantly different trends: it does not necessarily decrease monotonically with iterations even though cnorm continues to decrease. For example, at high SNR (e.g., 30 dB), the BER first decreases and then slightly increases after reaching a minimum, whereas at lower SNRs (e.g., 0 dB
9
TABLE III H ARDWARE P LATFORMS U SED FOR E VALUATION Type CPU GPU
Platform Intel Xeon Gold 6348 NVIDIA Jetson AGX Orin NVIDIA RTX 6000 Ada NVIDIA A100 NVIDIA H200
Memory DDR4 256 GB LPDDR5 64 GB GDDR6 48 GB HBM2e 40 GB HBM3e 141 GB
and 15 dB), the BER remains largely unchanged across iterations despite the substantial reduction in cnorm . Moreover, the same cnorm value can correspond to substantially different BER levels under different SNR conditions. For example, for cnorm around 10−1 , BER is 10−2 under 30 dB SNR and is larger than 10−1 under 15 dB SNR. These observations indicate that cnorm cannot reliably predict BER across operating regimes, making it unsuitable as a universal convergence indicator. The iteration budget is determined offline using the BERdriven criterion described above, and a fixed-iteration execution model is adopted during runtime. While condition-based early termination is common in CPU implementations, it is inherently inefficient for the GPU’s SIMT architecture due to the costly global synchronization and host-device interactions to decide whether to launch the subsequent demodulation kernel. Predetermined termination avoids unnecessary runtime monitoring while enabling predictable performance with efficient hardware utilization: executing a fixed number of iterations allows the computation to be organized into a fully GPU-resident pipeline with minimal host intervention. For example, for a DD domain grid size of (128, 32) with 10 conjugate gradient iterations, eliminating CPU-GPU synchronization reduces the equalization time from 0.65 ms to 0.54 ms on average, yielding a 17% speedup. V. I MPLEMENTATION Software implementation and hardware platforms. We implement the GPU-based Zak-OTFS processing pipeline in approximately 3,200 lines of Python code. PyTorch is used as the primary backend for general tensor computation, and OpenAI Triton [32] is used to develop custom GPU kernels for largescale matrix operations and performance-critical components, such as MVMs and GEMMs. We evaluate our implementation of the end-to-end Zak-OTFS receiver processing primarily on an NVIDIA H200 GPU. We also consider multiple NVIDIA GPU platforms, including Jetson Orin, RTX 6000 Ada, and A100, to provide a comprehensive performance comparison. To understand the impact of the hardware architecture, we also implement the same Zak-OTFS processing pipeline in C++ using the Armadillo library [33] and evaluate it on a Dell PowerEdge R750 server, equipped with two Intel Xeon Gold 6348 processors, each with 28 physical cores (56 threads) running at 2.6 GHz base frequency. Table III provides the detailed specifications of the considered hardware platforms. Channel model. We verify the proposed system using the ITU Vehicular-A (Veh-A) multipath channel model [30] to model a realistic high-mobility propagation environment, which was also adopted in a number of previous works [14], [15], [24],
TABLE IV P OWER - DELAY P ROFILE OF V EH -A C HANNEL M ODEL Path Index, p Delay τp (µs) Relative Power (dB)
1 0.00 0
2 0.31 −1
3 0.71 −9
4 1.09 −10
5 1.73 −15
6 2.51 −20
[34], [35]. Table IV summarizes the Veh-A channel profile, including six discrete paths with excess delays and their average power levels. The Doppler shift of the p-th path is drawn as νp = νmax cos(2πUp ), where Up ∼ U(0, 1) (uniform distribution), with νmax being the maximum Doppler frequency determined by the carrier frequency and vehicle speed. Note that both the delay and Doppler are fractional. Additive white Gaussian noise (AWGN) is then applied to the p received signal with the noise standard deviation set as σ = ρs /γ, where ρs is the mean signal power and γ is the SNR in linear scale. VI. E VALUATION In this section, we present a comprehensive evaluation of the proposed Zak-OTFS processing system in terms of realtime capability, latency distribution, scalability, and BER with varying compute platforms, equalizers, and data structures. A. Experiment Setup and Evaluation Metrics Experiment setup. We consider a Zak-OTFS system with M delay bins and N = 32 Doppler bins, a frequency spacing of ∆f = N ·∆ν = 30 kHz, corresponding to a Doppler resolution of ∆ν = 937.50 Hz and frame duration of T = 1.07 ms (3). The delay resolution ∆τ depends on the bandwidth B = M ∆f , given by (2)–(3). We set the value of the sampling rate equal to that of B. Quadrature amplitude modulation (QAM) is employed, with QPSK and 16-QAM considered. The channel coherence time Tc is inversely proportional to the Doppler spread, i.e., Tc ≈ 0.423/νmax . As the Doppler spread increases, the channel varies more rapidly within a frame, reducing temporal correlation and making channel estimation more challenging. We set νmax = 100 Hz in the Veh-A channel simulation for all latency-related experiments, unless otherwise specified. This corresponds to a coherence time on the order of Tc ≈ 4.23 ms, which is significantly larger than the OTFS frame duration T = 1.07 ms. This ensures that the channel remains approximately constant over a frame, allowing us to focus on computational performance and system latency without being dominated by channel aging effects. b dd (28) The estimated (sparse) DD domain channel matrix H is constructed using a threshold θ = 0.08. The equalizers, including CGA (Algo. 1) and MRC [22], are executed for Ξ = 10 iterations for a fair comparison, and the effectiveness of this iteration count is discussed in Section VI-D. Real-time processing latency deadline. A high deadline satisfaction rate is necessary for reliable communication services [25], as deadline violations can lead to packet drop and/or backlog accumulation [36]. Hence, we define real-time processing capability as meeting 99.9% latency within the deadline, consistent with prior literature on 5G vRAN processing for OFDM signals [26]. Note that in some evaluations,
10
B. End-to-end Performance A real-time system must process incoming frames within a bounded latency while providing strict timing guarantees. In this section, we evaluate whether our design satisfies these real-time requirements and quantify the achievable complexity reduction while maintaining BER. End-to-end processing latency. The observed practical scalability in Section III-C suggests that MVM demonstrates a constant execution time on GPU within an operation region. This observation reveals the potential of designing a processing with similar scalability. Fig. 5 shows the end-toend median latency of the OTFS processing across various compute platforms (CPU and GPU), equalizers (LMMSE, MRC, and CGA), and matrix data structure (with or without b dd ). The end-to-end latency leveraging structured sparsity of H is measured by transmitting an OTFS packet with a point-pilot frame, followed by a data frame, so the packet deadline will be 2.13 ms. For a fair comparison, we let both CGA and MRC
GPU LMMSE CPU LMMSE
GPU SS-MRC CPU SS-MRC
GPU CGA CPU CGA
GPU SS-CGA CPU SS-CGA
Median Latency (ms)
102
101 Real-time deadline 2.13 ms
100
10−1
2
8
32
8
2
48 192 768 20 8 32 Number of Delay Bins, M 12
51
Fig. 5. Median end-to-end processing latency across hardware platforms (CPU, GPU), equalizers (LMMSE, MRC, and CGA), and structured-sparsity (SS) awareness, with N = 32. The pilot/data frame duration of T = 1.067 ms and processing deadline is 2T = 2.134 ms. CPU 32 cores (M=32, N=32)
GPU (M=32, N=32)
GPU (M=16384, N=32)
100 Complementary CDF
we report median latency because the prohibitive execution latencies of the baseline methods prevent comprehensive data collection. We define the processing deadline of a packet of one pilot and one data frame as two frame durations (2T ): the combined pilot and data frame must be processed before the next packet is fully received. This setup assumes the worstcase scenario that the channel is estimated every other frame, whereas in the real-world scenario, the channel coherence time is longer. Unless otherwise specified, we consider pilot/data frame duration of T = 1.067 ms with N = 32 and ∆f = 30 kHz so that the packet deadline is 2T = 2.13 ms. BER calculation and throughput. We consider hard demodulation and the demodulated bits are compared with the ground truth to calculate frame-level BER. The BER is averaged across multiple frames. The throughput, η, is given by 1 (41) η = · B · bmod · (1 − BER), 2 where B equals the sample rate with a factor of 12 due to the interleaving pilot/data frame arrangement, bmod is the number of bits per DD domain symbol based on the modulation scheme. With QPSK modulation, our system achieves a data rate of 491.44 Mbps with M = 16,384, bmod = 2, and BER = 0.015% at 25 dB SNR. With 16QAM modulation, our system further pushed the data rate over 906.52 Mbps with bmod = 4 and BER = 7.78% under the same configuration. Baseline equalizers. We consider the LMMSE and MRC equalizers introduced in Section III-B4, and refer to the proposed compute-aware CGA in Section IV-C as CGA in the b dd is using the evaluation section. When the channel matrix H structured sparsity (SS) representation (27) in Section IV-B, we refer to the equalizer as SS-aware equalizers, such as SSMRC and SS-CGA. Note that LMMSE inherently requires a b dd for the matrix inversion, but SS is dense channel matrix H designed for efficient MVM (24) only. For both MRC and b dd is constructed in the pilot frame, and the iteraCGA, H tive equalizer is performed for each data frame. In contrast, b dd in the pilot frame and LMMSE constructs and inverts H performs one-shot GEMM equalization for each data frame.
10−1 10−2 p99.9
10−3 10−4 2.13 ms
10−5 0
1
2
3
4
5 6 Latency (ms)
7
8
9
10
Fig. 6. Latency CCDF across 10K packets, where GPU exhibits lower tail latency than the CPU, and increasing the grid size does not significantly affect the distribution. Overall, the GPU platform can meet the 99.9-th percentile (p99.9) processing latency for DD grid dimension of up to (16384, 32).
b dd perform 10 iterations and a threshold θ = 0.08 for sparse H construction in this work, unless otherwise specified. Fig. 5 demonstrates that our system (GPU SS-CGA) can constantly execute the OTFS receiver pipeline within the realtime deadline across the DD domain grid dimension M × N . By leveraging the structured-sparsity data representation, the low-complexity CGA equalizer, and the GPU platform, our system successfully translates the observed unit-test-level scalability of MVM into an end-to-end processing system while meeting the real-time deadline. The LMMSE-based solutions (16) suffer from high complexity incurred by GEMM and matrix inversions. Similarly, the CPU platform is limited by the available hardware parallelism, despite lower latency at smaller grid sizes, leading to suboptimal scalability. Note that even with low-complexity equalizers (e.g., CGA) and a GPU platform, a data structure exploiting the structured sparsity of b dd is the key enabler of constant latency with grid scaling, H as it reduces MVM complexity by leveraging the insight that the number of dominant paths is independent of the grid size. 99.9-th percentile processing latency. In addition to median processing time, the deadline satisfaction rate is a critical metric for evaluating the reliability of a real-time communication link. Fig. 6 presents the tail latency distribution of the OTFS processing pipeline on both CPU and GPU platforms, evaluated using 10K pilot and data frames. Due to significantly
11
10−1 10−2
21 23 25 27 29 211 213 215 Number of Delay Bins, M
1.07 ms
100
1.07 ms
100 10−1 10−2
21 23 25 27 29 211 213 215 Number of Delay Bins, M
(c) CGA (Dense)
0
10−2
21 23 25 27 29 211 213 215 Number of Delay Bins, M
(b) SS-MRC Median Latency (ms)
Median Latency (ms)
101
̂ Est./Inv. HDD Equal. Demod.
2.13 ms
SS-CGA
10−1
(a) LMMSE (Dense) 102
̂ heff Cal.
101
102
1
2
Latency (ms)
3 32
Fig. 8. Breakdown of the 99.9th percentile of the combined one pilot and one data frames with (M, N ) = (128, 32) on GPU. The LMMSE equalizer b dd ; SS-CGA suffers from high overhead due to the large matrix inversion in H trades off equalization time while keeping both the pilot and data frames within the per-frame deadline of T = 1.07 ms if viewed separately.
101 1.07 ms
LMMSE
SS-MRC 100
100
100
10−1
10 1
10 1
10 2
10 2
10−2
21 23 25 27 29 211 213 215 Number of Delay Bins, M
(d) SS-CGA
10 3
10 4 0 5 10 15 20 25 30 SNR (dB)
10 4 0 200 400 600 800 1000 max (Hz)
Processing time breakdown. Fig. 7 shows the latency of processing pilot and data frames, respectively, for each equalizer in Fig. 5 under the same configurations. The pilot and data frames inherently contain different DSP stages and thus different complexity. For example, an LMMSE equalizer performs the entire channel estimation in the pilot frame via matrix inversion and requires only one GEMM per data frame, resulting in unbalanced processing time. In contrast, b dd in the the MRC and CGA equalizers only construct H pilot frame and iteratively equalize the data frame with more complex algorithms. Note that even with structured sparsity, GPU SS-MRC exceeds the deadline for data frames, while GPU SS-CGA meets the deadline for both pilot and data frame constantly for 8 ≤ M ≤ 16384. Fig. 8 further compares LMMSE and SS-CGA for the latency of each OTFS processing stage on the GPU in a packet of one pilot and one data frame. SS-CGA perform Algo. 1 with MVM in each iteration, leading to 0.86 ms equalization time compared to the one-shot GEMM of LMMSE in 0.04 ms. However, the matrix inversion in LMMSE incurs significant overhead in pilotb dd inversion), preventing frame processing (i.e., 31.87 ms for H practical deployment. In contrast, the SS-CGA constructs (27)
(b) BER vs. νmax (QPSK)
(a) BER vs. SNR (QPSK)
100
100
10 1
10 1
10 2
10 2
10 3
10 3
10 4 0 5 10 15 20 25 30 SNR (dB)
10 4 0 200 400 600 800 1000 max (Hz)
(c) BER vs. SNR (16QAM)
(d) BER vs. νmax (16QAM)
Fig. 9. BER vs. SNR at νmax = 100 Hz and BER vs. νmax at 25 dB SNR across equalizers and QAM scheme. (M , N ) = (128, 32), and θ=0.08. SSCGA outperforms MRC and matches LMMSE in BER.
100 10 1 10 2 10 3 10 4 QPSK 10 5 16-QAM 6 10 20 22 24 26 28 210 212 214 Number of Delay Bins, M
BER
higher execution time of CPU processing with large grid sizes, we restrict CPU-based evaluation to (M, N ) = (32, 32). On the other hand, GPU-based processing is evaluated under both small and large grid sizes of (M, N ) = (32, 32) and (16384, 32). The results show that GPU processing not only scales efficiently with increasing grid sizes, but also exhibits a tighter latency distribution compared to CPU processing. Specifically, GPU processing with (M, N ) = (16384, 32) can meet the target deadline of 2.13 ms at 99.984-th percentile, with a 99.9-th percentile latency of 1.69 ms. In contrast, CPU processing suffers from a larger 99.9-th percentile latency of 3.72 ms, approximately 1.75× higher than the target deadline.
SS-CGA
10 3
BER
Fig. 7. The pilot and data frame latency breakdown of Fig. 5 for each equalizers. Among the two low-complexity equalizers (MRC and CGA), only the proposed SS-CGA constantly meets the latency requirement of 1.07 ms for processing both the pilot and the data frame, which is stricter than the packet-level amortized deadline of 2.13 ms, on GPU.
33
BER
100
DZT
LMMSE
BER
1.07 ms
GPU Data
(a) BER vs. M at νmax = 100 Hz
100 10 1 10 2 10 3 10 4 10 5 10 6 0
BER
101
GPU Pilot
102
BER
CPU Data
Median Latency (ms)
Median Latency (ms)
CPU Pilot
102
M=64 M=1024 M=16384
250 500 750 1000 max (Hz)
(b) BER vs. νmax (QPSK)
Fig. 10. BER of SS-CGA with 25 dB SNR. BER decreases when M increases, with QPSK having lower BER. Across multiple M values, the BER increases with νmax but SS-CGA still maintains BER lower than 1%.
in 0.34 ms. Note that, despite the longer equalization, SS-CGA meets the frame deadline for both the pilot and data frames, demonstrating a preferable design trade-off, and maintains a good scalability due to structured sparsity, as shown in Fig. 5. BER. We also evaluate BER versus SNR and νmax across combinations of equalizer design (LMMSE, SS-MRC, and SS-CGA) and QAM modulation (QPSK and 16QAM) using the Veh-A channel (Section V). SNR is swept from 0 to 30 dB, and each operating point is averaged over 1,000 randomly generated bitstreams and independent Veh-A channel realizations with AWGN. For the Veh-A channel model, the maximum Doppler spread νmax ranges from 0 to 1000 Hz in 100 Hz increments. Fig. 9 shows that the proposed SS-CGA achieves similar or better BER than LMMSE, while constantly
100
101
10−2
10−1 −5 10
̂ Construction Hdd
10−4
Equalization
10−4
10−1 10−3 10−2 Sparsity Threshold, θ
100
Iterations
0.1000 0.0100 0.0010 0.0001 0.0000
BER Convergence Lowest BER 35 30 25 20 15 10 5 0 0 5 10 15 20 25 30 5 10 15 20 25 30 SNR (dB) SNR (dB)
(a) M = 128, N = 32
10−6
b dd elements (top) and the execution Fig. 11. The magnitude distribution of H time/BER vs. sparsity thresholds θ (bottom) with (M, N ) = (128, 32), QPSK, SNR = 30 dB. An operating region is observed for 0.03 ≤ θ ≤ 0.08.
outperforms SS-MRC. These results verify that the proposed structured-sparsity representation using dominant paths with an iterative low-complexity equalizer does not sacrifice communication quality. Overall, BER decreases in all cases when SNR increases or νmax decreases, as higher noise power and higher Doppler shift degrade the detection capability of the receiver pipeline. As expected, 16QAM exhibits higher BER across all equalizers due to its greater sensitivity to noise. Fig. 10 reports the BER for SS-CGA with larger OTFS grid sizes by varying the number of delay bins M under 25 dB SNR. Fig. 10(a) shows BER over M under νmax = 100 Hz for QPSK and 16QAM. BER first decreases when M increases since a larger M provides better delay resolution ∆τ . After all the paths are separable around M = 64, increasing M starts to have diminishing gain on BER; BER saturates around 10−4 for QPSK. The 16QAM trend aligns with QPSK but with a higher BER due to its greater sensitivity to noise. Fig. 10(b) shows BER over a range of Doppler spread νmax from 0–1000 Hz for QPSK across three M values. As expected, BER degrades with larger νmax since the number of Doppler bins N does not change in the experiment. However, SS-CGA still maintains BER below 1% across the three cases. C. BER, Sparsity, and Processing Latency Trade-offs In Section IV-B, we discuss how our system leverages the b dd from the P dominant paths, which can be sparsity of H beff (28). We hereby provide obtained via a threshold θ on h an example-based numerical analysis for the three-way tradeb dd sparsity, the processing latency, and the off among the H final BER, as functions of θ. Fig. 11 reports (i) the absolute b dd , (ii) the average (100 magnitude CCDF of elements in H b trials) processing latency of both Hdd construction and equalization of our system on GPU, and (iii) the final resultant BER, by sweeping the threshold θ. The magnitude of elements in b dd is normalized to one, and two clusters around 0.02 and 1 H are presented. Hence, the thresholds between that range (i.e., 0.03 ≤ θ ≤ 0.08) lead to a desired tradeoff. For example, at θ = 0.08, compared to the unthresholded case, (i) 99.98% b dd elements can be pruned, so that (ii) the H b dd of the H construction speeds up from 2.43 ms to 0.34 ms (7.18×) and
(b) M = N = 32
Fig. 12. The iterations in Algo. 1 across SNR at the original cnorm -based termination [15], the lowest BER, and the BER convergence (40). The lowest BER and BER convergence lines are generally close. In contrast, cnorm -based termination constantly exceeds the designed maximum of 30 iterations and requires many more iterations after the lowest BER is achieved.
cnorm convergence 100
BER Convergence 100
10 1
10 1
10 2
10 2
Lowest BER
BER
100
35 30 25 20 15 10 5 0 0
Iterations
10−4 10−1 10−3 10−2 ̂ |Hdd | (Entry Magnitude)
102
100
cnorm convergence
1.0000
BER
100 10−1 10−2 10−3 10−4 10−5 −5 10
BER
Latency (ms)
CCDF
12
10 3
10 3
10 4 0 5 10 15 20 25 30 SNR (dB)
10 4 0 5 10 15 20 25 30 SNR (dB)
(a) M = 128, N = 32
(b) M = N = 32
Fig. 13. The achieved BER at the iteration of the original cnorm -based termination [15], the measured lowest BER, and the converged BER (40) in Algo. 1 across SNR. cnorm -based termination constantly yields higher BER.
equalization from 32.92 ms to 0.78 ms (297.31×), while (iii) the BER is significantly improved from 0.067% to 0.001%. This example suggests that a proper threshold θ benefits both the system processing latency and the BER. This threshold θ can be empirically selected based on the SNR estimated from the pilot frame, inferred from the number of reflectors in the channel, if known, or based on the desired latency requirement. D. BER Convergence in CGA Equalizer The unconditional termination (Section IV-C) requires offline profiled iteration counts to guarantee the desired algorithm convergence. In Fig. 4, the trend of cnorm decreasing over iteration is shown. This section compares the behavior and BER of cnorm -monitoring and the offline-profiling methods using the unconditional CGA (Algo. 1). In each experiment, we run Algo. 1 up to 30 iterations and record the iterations and BER at which (i) the cnorm converges based on a fixed threshold ϵ2 [15], (ii) the lowest BER is achieved, and (iii) the BER convergence iteration defined in Section IV-C. Iterations of BER and CGA convergence. Fig. 12 shows the recorded iterations for cnorm convergence, BER convergence, and the lowest (best) BER on (M, N ) = (128, 32) and (32, 32), varying SNR. The results reveal that cnorm cannot fall below its ϵ2 = 10−3 threshold [15] and reaches the maximum iterations in all cases. This is because cnorm scales linearly with M N in Algo. 1 (line 7) and experiences diminishing slope despite monotonic decay over iterations (Fig. 4). In contrast, both the BER convergence (40) and the lowest BER appear in much earlier iterations, indicating that the additional iterations yield diminishing BER gain if the equalizer uses cnorm as the
Median Latency (ms)
Median Latency (ms)
13
102 Deadline (2N/30 ms)
101 100 10−1
A100 H200
RTX 6000 Jetson Orin
2 8 32 128 512 2048 Number of Doppler Bins, N
(a) M =1024 across GPUs
102 101 Deadline (2N/30 ms)
100 10−1
M=32 M=512 M=1024
M=2048 M=8192
2 8 32 128 512 2048 Number of Doppler Bins, N
(b) Jetson Orin across M
Fig. 14. Median processing latency across GPUs and DD grid sizes (M, N ), with (a) varying N (which affects the processing deadline), and (b) varying both M and N on the Jetson Orin platform. A100
H200
RTX 6000
Jetson Orin
CCDF: P(X > x)
100 10−1 10−2 p99.9
10−3 10−4
2.13 ms 10−5 (N=32) 0 2
4.27 ms (N=64) 4
6
17.07 ms (N=256) 8 10 Latency (ms)
12
14
16
Fig. 15. Complementary CDF of the processing latency CCDF across four GPU platforms. Each GPU is evaluated with M = 1024 and the largest N value that satisfies the 99.9% latency deadline: N = 32 for RTX 6000 and H200, N = 64 for A100, and N = 256 for Jetson Orin.
termination indicator. Instead, our offline profiling empirically suggests that CGA converges fast in early iterations. BER. Fig. 13 shows the BER corresponding to the recorded iterations in Fig. 12 and that achieving the converged BER and the lowest BER requires consistently lower iterations than the maximum iterations (which is also where the cnorm -based termination), indicating that additional iterations after BER convergence are unnecessary without BER improvement. This observation demonstrates the effectiveness of the proposed BER convergence criterion (40) since the converged BER overlaps with the lowest BER. Hence, we set the iteration count Ξ = 10 in Algo. 1 for all SNRs in our evaluation. E. Latency Requirement and GPU Platform Capability The evaluations so far assumed a fixed N as it corresponds to the frame duration T and thus the processing deadline (3). However, GPU platforms exhibit slightly different latencies in unit tests of matrix operations in Section III-C. In this section, we evaluate how different GPU platforms affect the processing latency under varying values of N . Varying N does not change the system throughput when M and ∆f remain constant (3). Fig. 14(a) shows the median end-to-end latency across 10K packets across four GPU platforms (NVIDIA Jetson Orin, RTX 6000 Ada, A100, H200), with M = 1024 and varying N . Given ∆f = 30 kHz, each pilot/data frame has a duration N of 2T , with T = ∆f = 30NkHz (s) = 30NHz (ms). GPU platforms have individual constant overhead and upward turning points, resulting in intersections of latency trends to the real-time processing deadline 2T . These intersections mark their ability to satisfy the tighter deadline constraints with a smaller N . The upward turning points mark where the size of the OTFS grid (M, N ) starts to exceed the parallelization capacity of a GPU and leads to latency scaling. Hence, Jetson Orin,
which is designed for edge computing with lower hardware specifications, exhibits higher constant overhead and an earlier upward turning point. Note that as M increases, the latency lines may not have intersections with the processing deadline line, implying deadline violation for all possible N values. Fig. 14(b) shows the median latency across 10K packets of Jetson Orin for multiple M values with varying N . As M increases, the upward turning points move leftward, as product M × N remains constant. Hence, when M is sufficiently large (M ≥ 2048 for Jetson Orin), a GPU will never be able to meet the real-time deadline, marking its capability limit. Recall that a higher M corresponds to a higher input sample rate when N and ∆f remain constant for a fixed deadline (3). Fig. 15 shows the CCDF of the processing latency across GPUs based on 10K packets, highlighting the differences in their ability to handle smaller frame sizes. As the latency deadline increases with larger N , each GPU is evaluated at the largest N value for which its 99.9% processing latency satisfies the real-time deadline, i.e., the intersection in Fig. 14(a), where N = 32 for RTX 6000 Ada and H200, N = 64 for A100, and N = 256 for Jetson Orin, all with M = 1024. Note that all configurations share the same bandwidth and I/Q sampling rate; hence, the selection of N is based on the channel Doppler spread νmax (and thus the coherence time Tc , Section VI-A) and each GPU’s capabilities. VII. C ONCLUSIONS In this paper, we presented a scalable and practical real-time Zak-OTFS processing system on GPU platforms. Specifically, we exploited compact matrix operations for Zak transform and channel estimation, and a branchless iterative equalizer. In addition, the structured sparsity of the channel matrix was derived and leveraged to significantly reduce the complexity of the channel estimation and equalization. With the insights on hardware-algorithm co-design, our system enabled realtime processing on GPUs that scales up to the grid size of (16384, 32) with a deadline of 2.13 ms, yielding 491.44 Mbps and 906.52 Mbps for QPSK and 16QAM, with BER = 0.015% and 7.78% at 25 dB SNR, respectively. R EFERENCES [1] R. Hadani, S. Rakib, M. Tsatsanis, A. Monk, A. J. Goldsmith, A. F. Molisch, and R. Calderbank, “Orthogonal time frequency space modulation,” in Proc. IEEE WCNC, 2017. [2] S. K. Mohammed, R. Hadani, A. Chockalingam, and R. Calderbank, “OTFS–a mathematical foundation for communication and radar sensing in the delay–Doppler domain,” IEEE BITS Inf. Theory Mag., vol. 2, no. 2, pp. 36–55, 2022. [3] ——, “OTFS–predictability in the delay–Doppler domain and its value to communication and radar sensing,” IEEE BITS Inf. Theory Mag., vol. 3, no. 2, pp. 7–31, 2023. [4] T. Thaj and E. Viterbo, “OTFS modem SDR implementation and experimental study of receiver impairment effects,” in Proc. IEEE ICC Workshops, 2019. [5] X. Wei, W. Yuan, K. Zhang, and F. Liu, “OTFS-assisted ISAC system: Delay Doppler channel estimation and SDR-based implementation,” IEEE Trans. Mob. Comput., 2025. [6] W. Lin, W. Wang, K. Wang, Z. Deng, D. Wan, Y. Deng, Z. Feng, Z. Fan, and S. Deng, “Enabling OTFS for hypersonic aircraft over LEO satellite: A 3D system model and differential delay–Doppler pilot,” IEEE Internet Things J., 2026.
14
[7] M. Nauman, L. Lopacinski, N. Maletic, M. Scheide, J. Gutiérrez, M. Krstic, and E. Grass, “OTFS sensing with SDR: Experimental results and analysis,” in Proc. IEEE VTC–Spring, 2025. [8] Y. X. Yap, N. Bhushan, O. Dizdar, A. Sattarzadeh, D. Redgate, V. Battula, and S. Wang, “A proof of concept for OTFS resilience in doublyselective channels by GPU-enabled real-time SDR,” in Proc. IEEE GLOBECOM, 2023. [9] E. C.-H. Jou, B.-T. Yeh, and T.-D. Chiueh, “Design and implementation of a MIMO-OTFS receiver for LEO satellite communications with overthe-air verification,” IEEE Trans. Circuits Syst. I, Reg. Papers., pp. 1–14, 2025. [10] Y. Wang, M. Ouyang, T. Zhang, and F. Gao, “An RFSoC-based scalable OTFS prototyping platform for integrated sensing and communications,” IEEE Trans. Wireless Commun., vol. 24, no. 9, pp. 7470–7482, 2025. [11] M. Isik, M. Nkomo, A. Das, and K. R. Dandekar, “FPGA implementation of OTFS modulation for 6G communication systems,” in Proc. IEEE FNWF, 2023. [12] A. R. Shadangi, S. S. Das, and I. Chakrabarti, “VLSI architecture for implementing OTFS,” Wireless Pers. Commun., vol. 137, no. 3, pp. 1931–1944, 2024. [13] S. K. Dora, H. B. Mishra, and M. Sahoo, “Low complexity implementation of OTFS transmitter using fully parallel and pipelined hardware architecture,” J. Signal Process. Syst., vol. 95, no. 8, pp. 955–964, 2023. [14] J. Zheng, V. Khammammetti, B. Dabak, S. R. Mattu, T. Chen, and R. Calderbank, “Zak-OTFS with spread pilot in sub-6 GHz: Implementation and over-the-air experimentation,” in Proc. IEEE MILCOM, 2025. [15] S. R. Mattu, N. Mehrotra, S. K. Mohammed, V. Khammammetti, and R. Calderbank, “Low-complexity equalization of Zak-OTFS in the frequency domain,” npj Wireless Technology, vol. 2, no. 9, 2026. [16] “Scalable and real-time Zak-OTFS receivers on GPUs,” https://github. com/functions-lab/zak-otfs-rt, 2026. [17] M. Ubadah, S. K. Mohammed, R. Hadani, S. Kons, A. Chockalingam, and R. Calderbank, “Zak-OTFS to integrate sensing the I/O relation and data communication,” arXiv preprint arXiv:2404.04182, 2024. [18] G. D. Surabhi and A. Chockalingam, “Low-Complexity Linear Equalization for OTFS Modulation,” IEEE Commun. Lett., vol. 24, no. 2, pp. 330–334, 2020. [19] R. Chong, S. Li, Z. Wei, M. Matthaiou, D. W. K. Ng, and G. Caire, “Cross-domain iterative detection for OTFS transmission with frequency domain equalization,” IEEE Trans. Commun., 2025. [20] H. B. Mishra, P. Singh, A. K. Prasad, and R. Budhiraja, “OTFS Channel Estimation and Data Detection Designs with Superimposed Pilots,” IEEE Trans. Wireless Commun., vol. 21, no. 4, pp. 2258–2274, 2021. [21] P. Raviteja, K. T. Phan, Y. Hong, and E. Viterbo, “Interference cancellation and iterative detection for orthogonal time frequency space modulation,” IEEE Trans. Wireless Commun., vol. 17, no. 10, pp. 6501– 6515, 2018. [22] 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, 2020. [23] H. Zhang, X. Huang, and J. A. Zhang, “Zak-OTFS with time-domain offset gradient descent equalization,” IEEE Trans. Veh. Technol., pp. 1– 10, 2025. [24] S. R. Mattu, N. Mehrotra, and R. Calderbank, “Differential communication in channels with mobility and delay spread using Zak-OTFS,” IEEE Wireless Commun. Lett., 2025. [25] J. Ding, R. Doost-Mohammady, A. Kalia, and L. Zhong, “Agora: Realtime massive MIMO baseband processing in software,” in Proc. ACM CoNEXT, 2020. [26] Z. Qi, C.-H. Tung, A. Kalia, and T. Chen, “Savannah: Efficient mmWave baseband processing with minimal and heterogeneous resources,” in Proc. ACM MobiCom, 2024. [27] X. Xie, W. Hou, Z. Guo, and M. Liu, “Building massive MIMO baseband processing on a single-node supercomputer,” in Proc. USENIX NSDI, 2025. [28] Z. Qi, Y. Yao, Y. Li, C.-H. Tung, J. Zheng, D. Zhuo, and T. Chen, “DecodeX: Exploring and benchmarking of LDPC decoding across CPU, GPU, and ASIC platforms,” in Proc. ACM HotMobile, 2026. [29] D. Villa, I. Khan, F. Kaltenberger, N. Hedberg, R. S. da Silva, and S. e. a. Maxenti, “X5G: An open, programmable, multi-vendor, end-to-end, private 5G O-RAN testbed with NVIDIA ARC and OpenAirInterface,” IEEE Trans. Mob. Comput., vol. 24, no. 11, pp. 11 305–11 322, 2025. [30] ITU–R, “Guidelines for evaluation of radio transmission technologies for IMT-2000,” International Telecommunication Union Radiocommunication, Tech. Rep. Recommendation ITU-R M.1225, 1997. [31] F. Lampel, A. Avarado, and F. M. Willems, “On OTFS using the discrete Zak transform,” in Proc. IEEE ICC Workshops, 2022.
[32] P. Tillet, H.-T. Kung, and D. Cox, “Triton: an intermediate language and compiler for tiled neural network computations,” in Proc. ACM MAPL, 2019. [33] C. Sanderson and R. Curtin, “Armadillo: an efficient framework for numerical linear algebra,” in Proc. IEEE ICCAE, 2025. [34] N. Mehrotra, S. R. Mattu, and R. Calderbank, “Zak-OTFS with spread carrier waveforms,” IEEE Wireless Commun. Lett., 2025. [35] B. Dabak, V. Khammammetti, S. K. Mohammed, and R. Calderbank, “Zak-OTFS and LDPC codes,” in Proc. IEEE ICC, 2024. [36] C.-H. Tung, Z. Qi, and T. Chen, “RISE: Real-time image processing for spectral energy detection and localization,” Proc. IEEE DySPAN, 2026.