FFTArray: Continuous Fourier Transform Automation

• FFTArray is a Python library for discretizing continuous Fourier transforms using FFTs on arbitrarily shifted, multi-dimensional grids, automating phase and scale corrections.
• It abstracts the transition from continuous Fourier integrals to backend FFT calls, enabling accurate pseudo-spectral PDE solvers with optimized lazy factor application.
• Built on the Python Array API Standard, FFTArray interoperates with NumPy, JAX, and PyTorch, supporting GPU acceleration and efficient coordinate-aware operations.

FFTArray is a Python library for implementing discretized continuous Fourier transforms by means of FFTs on arbitrarily shifted, multi-dimensional coordinate grids. Its purpose is not merely to expose a fast discrete transform, but to automate the mathematically correct passage from continuous Fourier integrals on finite sampled domains to backend FFT calls, including the required phase and scale corrections induced by grid spacing and coordinate offsets. Built on the Python Array API Standard, it supports GPU acceleration and interoperates with array backends such as NumPy, JAX, and PyTorch, while aiming to let users write code that follows textbook Fourier formulas rather than backend-specific FFT bookkeeping (Seckmeyer et al., 17 Jul 2025).

1. Scope, design goals, and problem setting

FFTArray is motivated by a common difficulty in pseudo-spectral PDE solvers: one typically needs not an abstract discrete Fourier transform, but a discrete approximation to a specific continuous Fourier integral defined on a finite equidistant grid. That approximation depends on the domain truncation, the coupling between $\Delta x$ and $\Delta f$, the coordinate offsets $x_{\min}$ and $f_{\min}$, and the phase and normalization conventions of the underlying FFT backend. The library is therefore positioned as a separate discretization layer between analytic equations and application code, rather than as a monolithic solver framework (Seckmeyer et al., 17 Jul 2025).

The paper states three goals explicitly: “From formulas to code,” “State-of-the-art performance,” and “Seamless multidimensionality.” In practice, this means that FFTArray is designed so that coordinate construction, space changes, and spectral operators can be written in a form close to their analytic expressions, while the library handles the grid-specific corrections automatically. This design contrasts with software in which Fourier-transform discretization is tightly integrated into a full-stack application, a pattern the paper argues sacrifices generality and encourages code duplication.

A central conceptual distinction is that FFTArray does not redefine the FFT itself. Instead, it formalizes the map

$$\text{continuous Fourier transform on a sampled finite grid} \;\longrightarrow\; \text{correctly phased and scaled backend FFT}.$$

This suggests a broader interpretation of the library: it is a coordinate-aware Fourier-discretization framework whose FFT calls are subordinate to a higher-level numerical representation.

2. Mathematical framework: continuous transforms, sampled grids, and compatibility conditions

FFTArray adopts the continuous Fourier transform pair

$$G(f) = \int_{-\infty}^{\infty} dx \, g(x)\, e^{- 2 \pi i f x}, \qquad g(x) = \int_{-\infty}^{\infty} df\, G(f)\, e^{2 \pi i f x},$$

using linear frequency $f$ rather than angular frequency $\omega$ (Seckmeyer et al., 17 Jul 2025).

To discretize these integrals, the function is sampled on finite equidistant grids

$x_n = x_{\min} + n \Delta x,\quad n=0,\ldots,N-1,$

$f_m = f_{\min} + m \Delta f,\quad m=0,\ldots,N-1,$

with sampled values $\Delta f$0 and $\Delta f$1. The resulting general discretized Fourier transform and inverse are

$\Delta f$2

$\Delta f$3

The position and frequency grids are not independent. FFTArray enforces the sampling-theorem relations

$\Delta f$4

These conditions determine $\Delta f$5 once $\Delta f$6 and $\Delta f$7 are fixed. By contrast, $\Delta f$8 and $\Delta f$9 remain free and enter only through the phase structure of the discretization. The paper emphasizes that aliasing and finite-domain truncation are still the user’s responsibility: FFTArray handles the bookkeeping of the discretization, not the approximation-theoretic choice of domain and resolution.

For real-valued data, the paper notes the usual conjugate symmetry of $x_{\min}$0 and defines the Nyquist frequency as

$x_{\min}$1

This is standard, but in FFTArray it is embedded in a framework where the grids need not begin at the origins assumed by raw FFT calls.

3. Reduction to backend FFTs and the role of phase and scale corrections

Internally, FFTArray uses the NumPy DFT convention

$x_{\min}$2

The library’s main derivation shows how the general discretized transform can be rewritten as a standard FFT surrounded by coordinate-dependent prefactors (Seckmeyer et al., 17 Jul 2025).

For the forward transform, the paper derives

$x_{\min}$3

For the inverse transform, it derives

$x_{\min}$4

These formulas show that a raw FFT is generally insufficient unless the coordinate grid satisfies special offset choices. FFTArray therefore treats three correction classes as first-class numerical objects: a position-space phase factor, a frequency-space phase factor, and a scale factor involving $x_{\min}$5 or $x_{\min}$6. Their application is only $x_{\min}$7, so the overall complexity remains $x_{\min}$8.

The special cases recover familiar FFT idioms. If $x_{\min}$9, the formulas simplify substantially. For a symmetric frequency grid with

$f_{\min}$0

the transform reduces to the standard fftshift representation,

$f_{\min}$1

If both position and frequency grids are chosen symmetrically,

$f_{\min}$2

then both pre- and post-shift corrections collapse into the familiar fftshift/ifftshift pattern. The paper points out that for even $f_{\min}$3 this “symmetric” grid is not perfectly symmetric in the literal sample set, a subtlety that is often glossed over in ad hoc implementations.

The same framework yields textbook spectral identities without hidden assumptions. For example, the derivative identity

$f_{\min}$4

becomes, on sampled grids,

$f_{\min}$5

and in the standard $f_{\min}$6, symmetric-frequency case reduces to the familiar fftfreq formula.

4. Library architecture, constraints, and lazy factor application

FFTArray is organized around two main abstractions: Dimension and Array. A Dimension stores one coordinate axis, including $f_{\min}$7, $f_{\min}$8, $f_{\min}$9, and $$\text{continuous Fourier transform on a sampled finite grid} \;\longrightarrow\; \text{correctly phased and scaled backend FFT}.$$0, while deriving the remaining grid quantities consistently. An Array stores sampled values together with their associated dimensions, the current representation per dimension ("pos" or "freq"), and the state of the correction factors (Seckmeyer et al., 17 Jul 2025).

A notable feature is constraint-based grid construction using the z3 solver. Users may specify any convenient subset of grid constraints, and FFTArray solves for the rest subject to the compatibility condition $$\text{continuous Fourier transform on a sampled finite grid} \;\longrightarrow\; \text{correctly phased and scaled backend FFT}.$$1. When the constraints are ambiguous, it raises NoUniqueSolutionError; when inconsistent, NoSolutionFoundError. The paper also states that FFTArray constrains $$\text{continuous Fourier transform on a sampled finite grid} \;\longrightarrow\; \text{correctly phased and scaled backend FFT}.$$2 to be even or a power of two, and it can relax designated loose_params if exact satisfaction would otherwise require an unsupported size.

The API is representation-centric rather than FFT-centric. Instead of calling fft() or ifft() directly in typical use, one writes $f$0 so that code expresses a change of representation rather than a backend primitive. This is paired with named-dimension broadcasting: arrays align by dimension name, and equally named dimensions must have matching grids and spaces. FFTArray does not automatically transform one operand to match another, because pointwise operations in position and frequency space are mathematically different.

The library’s most distinctive implementation idea is lazy factor application. After a space change, FFTArray may leave the array in an internal state where the FFT or IFFT has been performed but the final phase/scale factors have not yet been applied. This is tracked by the boolean factors_applied, with eager=False by default. The paper introduces internal intermediate states

$$\text{continuous Fourier transform on a sampled finite grid} \;\longrightarrow\; \text{correctly phased and scaled backend FFT}.$$3

to formalize these partially corrected representations, and then derives algebraic rules for addition, multiplication, division, and absolute value when operands have different correction states.

The practical effect is that many workflows avoid explicit phase-factor multiplications entirely. If successive operations would cancel missing corrections anyway, FFTArray skips them. The paper argues that this both improves speed and reduces unnecessary floating-point error. A plausible implication is that the library’s abstraction is not merely semantic; it enables backend-visible optimization by turning phase bookkeeping into symbolic state rather than immediate numerical work.

5. Multi-dimensional formulation and scientific applications

FFTArray’s multi-dimensionality is based on separability: each dimension carries its own coordinate metadata and current space, and multi-dimensional transforms are compositions of the corresponding one-dimensional space changes. Dimensions are labeled rather than addressed by raw axis numbers, so tensor-product constructions can be written in a way closer to the analytic form of the underlying equations (Seckmeyer et al., 17 Jul 2025).

The paper uses this framework in several application classes. A basic example is spectral differentiation. For

$$\text{continuous Fourier transform on a sampled finite grid} \;\longrightarrow\; \text{correctly phased and scaled backend FFT}.$$4

FFTArray implements

$$\text{continuous Fourier transform on a sampled finite grid} \;\longrightarrow\; \text{correctly phased and scaled backend FFT}.$$5

essentially verbatim through a space change, multiplication by the derivative kernel, and a return to position space. The reported result agrees with the analytic derivative to at least 11 decimal digits for the chosen setup.

The main PDE use case is the Schrödinger equation

$$\text{continuous Fourier transform on a sampled finite grid} \;\longrightarrow\; \text{correctly phased and scaled backend FFT}.$$6

Using second-order Trotter splitting,

$$\text{continuous Fourier transform on a sampled finite grid} \;\longrightarrow\; \text{correctly phased and scaled backend FFT}.$$7

and the diagonalization of the Laplacian in frequency space,

$$\text{continuous Fourier transform on a sampled finite grid} \;\longrightarrow\; \text{correctly phased and scaled backend FFT}.$$8

FFTArray allows the split-step routine to be written as alternating multiplications in frequency and position space. The paper emphasizes that in this routine, if the initial state has factors_applied=False, no phase/scale factors need ever be materialized; each into_space call can reduce to a raw backend FFT or IFFT.

The paper also treats Bragg diffraction with Hamiltonian

$$\text{continuous Fourier transform on a sampled finite grid} \;\longrightarrow\; \text{correctly phased and scaled backend FFT}.$$9

imaginary-time propagation for ground-state search, the 2D isotropic harmonic oscillator

$$G(f) = \int_{-\infty}^{\infty} dx \, g(x)\, e^{- 2 \pi i f x}, \qquad g(x) = \int_{-\infty}^{\infty} df\, G(f)\, e^{2 \pi i f x},$$0

and a coupled two-species Gross–Pitaevskii system for $$G(f) = \int_{-\infty}^{\infty} dx \, g(x)\, e^{- 2 \pi i f x}, \qquad g(x) = \int_{-\infty}^{\infty} df\, G(f)\, e^{2 \pi i f x},$$1 and $$G(f) = \int_{-\infty}^{\infty} dx \, g(x)\, e^{- 2 \pi i f x}, \qquad g(x) = \int_{-\infty}^{\infty} df\, G(f)\, e^{2 \pi i f x},$$2. In the latter case, FFTArray is used to express nonlinear effective potentials in position space and kinetic evolution in frequency space within the same split-step framework. The paper specifically notes that the $$G(f) = \int_{-\infty}^{\infty} dx \, g(x)\, e^{- 2 \pi i f x}, \qquad g(x) = \int_{-\infty}^{\infty} df\, G(f)\, e^{2 \pi i f x},$$3 ordering is advantageous because the nonlinearities depend on $$G(f) = \int_{-\infty}^{\infty} dx \, g(x)\, e^{- 2 \pi i f x}, \qquad g(x) = \int_{-\infty}^{\infty} df\, G(f)\, e^{2 \pi i f x},$$4 in position space; using $$G(f) = \int_{-\infty}^{\infty} dx \, g(x)\, e^{- 2 \pi i f x}, \qquad g(x) = \int_{-\infty}^{\infty} df\, G(f)\, e^{2 \pi i f x},$$5 would require extra IFFTs.

These examples collectively show that FFTArray is intended less as an isolated transform library than as reusable infrastructure for pseudo-spectral solvers. This suggests its main scientific significance lies in reducing the cost of changing geometry, coordinate offsets, or solver structure without rewriting Fourier-discretization machinery.

6. Performance, backend integration, and limitations

FFTArray is built on the Python Array API Standard and uses array-api-compat to integrate with NumPy, JAX, and PyTorch, including GPU execution when the backend supports it. The backend namespace is denoted xp, and the paper notes that some detailed behaviors remain backend-defined where the standard is silent. It also implements JAX pytree support via fa.jax_register_pytree_nodes(), with dimension metadata static by default during tracing and an optional dynamically_traced_coords=True mode for traced coordinate parameters except $$G(f) = \int_{-\infty}^{\infty} dx \, g(x)\, e^{- 2 \pi i f x}, \qquad g(x) = \int_{-\infty}^{\infty} df\, G(f)\, e^{2 \pi i f x},$$6, whose fixed shape is required by JAX (Seckmeyer et al., 17 Jul 2025).

The main performance benchmark is a $$G(f) = \int_{-\infty}^{\infty} dx \, g(x)\, e^{- 2 \pi i f x}, \qquad g(x) = \int_{-\infty}^{\infty} df\, G(f)\, e^{2 \pi i f x},$$7 2D isotropic harmonic-oscillator imaginary-time evolution. Three implementations are compared: FFTArray Direct, FFTArray Precomputed, and Raw FFT. Hardware includes AMD Epyc 7543, AMD Ryzen 7950X3D, NVIDIA A100, and NVIDIA RTX 4090; backends include NumPy and JAX. The reported results are that GPUs are about two orders of magnitude faster than CPUs for this workload, that under JAX all three variants perform almost identically on both CPUs and GPUs, and that with NumPy the direct version is slower because propagators are recomputed rather than reused, whereas the precomputed version is comparable to Raw FFT. The paper therefore argues that in the large-array regimes relevant to pseudo-spectral solvers, FFTArray introduces no measurable overhead relative to backend FFT use when the workflow is structured appropriately.

The accuracy study for imaginary-time propagation in the 2D harmonic oscillator reports relative energy errors better than $$G(f) = \int_{-\infty}^{\infty} dx \, g(x)\, e^{- 2 \pi i f x}, \qquad g(x) = \int_{-\infty}^{\infty} df\, G(f)\, e^{2 \pi i f x},$$8 in float64 for both FFTArray and TorchGPE at $$G(f) = \int_{-\infty}^{\infty} dx \, g(x)\, e^{- 2 \pi i f x}, \qquad g(x) = \int_{-\infty}^{\infty} df\, G(f)\, e^{2 \pi i f x},$$9. In float32, energy evaluation becomes a dominant accuracy limiter, and evaluating the energy in float64 after evolving in float32 improves precision by about three orders of magnitude. The paper also states that FFTArray outperforms TorchGPE in the tested scenario, partly because TorchGPE uses doubled frequency-space resolution that cannot be disabled.

The library’s limitations are explicit. It does not remove the need to choose domains and resolutions that control aliasing and truncation error. It does not automatically reconcile operations between position and frequency space. Indexing with a step greater than 1 is unsupported because reducing one-space resolution implies a non-unique modification of the conjugate-space extent. Dynamic traced coordinates in JAX are possible but delicate. The design is based on separable tensor-product dimensions and does not address non-Cartesian transforms. These constraints do not negate the library’s purpose; rather, they define it narrowly as a mathematically explicit layer for discretized Fourier transforms on Cartesian sampled grids.

Taken together, FFTArray represents an attempt to make Fourier-discretization semantics explicit, portable, and optimizable. Its defining claim is that the map from continuous transform formulas to backend FFTs can be automated without obscuring the mathematics and without materially sacrificing performance in the regimes for which pseudo-spectral methods are typically used.