FFTL is a fast Fourier transform written in Julia. It computes complex, real-input and
real-to-real transforms (DCT, DST, Hartley, halfcomplex) of any length and any number of
dimensions, in Float32, Float64 and other floating-point types, on several threads, and on
GPUs through KernelAbstractions.jl. The same Julia code compiles with juliac --trim into a C
library.
FFTL was written by a coding agent under the guidance of a domain expert. Julia makes a library of this scope possible in one language: codelets are generated by the compiler for the radix, element type and SIMD width in use, and one code base serves CPUs, GPUs and the C library. A coding agent makes it possible to implement, test and measure the many details such a library needs. The aim now is to develop it further as a project of the Julia community (see Contributing).
FFTL is experimental and not registered yet; the interface can still change. The manual is at JuliaMath.github.io/FFTL.jl.
using Pkg
Pkg.add(url = "https://github.com/JuliaMath/FFTL.jl")using FFTL, LinearAlgebra
x = randn(ComplexF64, 1000, 1000)
y = fft(x) # all dimensions
y1 = fft(x, 1) # along dimension 1 only
x2 = ifft(y) # normalised inverse; bfft is the unnormalised backward DFT
p = plan_fft(x; threads = 8)
mul!(y, p, x) # out of place, or mul!(x, p, x) in place
x3 = p \ y # inverse through the plan
r = randn(1000, 1000)
R = rfft(r) # half spectrum along dimension 1 (501 × 1000)
r2 = irfft(R, 1000)
c = dct(r) # orthonormal DCT-II; idct is its inverse
h = FFTL.r2r(r, (FFTL.REDFT10, FFTL.R2HC)) # unnormalised real-to-real, one kind per dimensionFFTL exports the names FFTW.jl exports (fft, plan_fft, fft!, rfft, dct, plan_dct, …),
with FFTW.jl's conventions and scalings, so using FFTL can take the place of using FFTW.
As in FFTW.jl, FFTL.r2r, the kinds and FFTL.set_num_threads are used with the prefix. These
functions are FFTL's own: a module that writes using FFTL gets FFTL, also when one of its
dependencies loads FFTW.jl. Generic code that calls
AbstractFFTs.jl's functions runs on FFTL as well;
with FFTW.jl loaded, it runs on FFTW for FFTW's arrays. The manual has a page for
coming from FFTW.jl.
Plans use threads = Threads.nthreads() by default; transforms below about 2^15 values run on
one thread regardless. Inside code that is already threaded, pass threads = 1. A plan owns its
work buffers, so one plan must not run on several tasks at the same time; make one plan per task
(plans of the same size share their twiddle tables, so this is cheap).
The first plan for a new radix sequence compiles its kernels (milliseconds for the sizes of the precompile workload, up to a few seconds otherwise); later plans of the same size take microseconds, and transforms never compile.
With KernelAbstractions.jl and a GPU package loaded, the same functions take GPU arrays. Some GPU
packages, Metal.jl among them, export AbstractFFTs's fft, plan_fft, …, which are different
functions from FFTL's; with such a package (or FFTW.jl) loaded, write FFTL.fft:
using FFTL, Metal, KernelAbstractions # or CUDA
yd = FFTL.fft(MtlArray(randn(ComplexF32, 1024, 1024)))Transforms, number types and platforms, in Julia.
- Complex, real-input and all eleven real-to-real kinds of FFTW, along any subset of dimensions, for any length (the part of a length made of prime factors above 31 through Bluestein's algorithm).
Float32andFloat64on the fast path (SIMD codelets, threads).BigFloat,Double64,Float128and other floating-point types are transformed in their own precision by a generic path, accurate to a few ulps of the type.- GPUs through KernelAbstractions.jl (tuned on Metal and CUDA; experimental). Complex and real transforms are tuned; the real-to-real kinds work, and their speed depends on the kind.
A C library made from the Julia code. libfftl/ compiles the package with
juliac --trim=safe (Julia ≥ 1.12) into a shared library whose interface is FFTW3's, as the
FFTW manual describes it: the C functions in both precisions (libfftl/include/fftw3.h), the
legacy Fortran routines and the Fortran 2003 interface. Programs written for FFTW3 link against
it without source changes (-lfftw3, -lfftw3f; the library carries FFTW's library names as
links); libfftl/include/fftw3.h lists what differs from FFTW. No transform is written in C:
the library is the trimmed Julia code. A program that links it needs no Julia installation; its
only dependency is the Julia runtime shipped next to it. On macOS arm64 the library is 15 MB and
the runtime 10 MB. The run-time plans the library uses (FFTL.RPlan, FFTL.RRFFTPlan,
FFTL.RR2RPlan) are public, so Julia programs that use them can be trimmed the same way.
Entry points for special variants. FFTL.plan_pass(v, d; nonzero, needed) transforms all
lines of a strided view along one dimension and can skip zero inputs and unneeded outputs
(pruning); FFTL.plan_r2r_pass does the same for real-to-real kinds. Libraries that build their
own transforms from one-dimensional passes use these; NFFT.jl is one of them.
The benchmarks page compares FFTL with FFTW, ducc0 and vDSP on an Apple M1 Max and an AMD EPYC, and with the vendor libraries and VkFFT on the M1's GPU and an NVIDIA RTX A6000, and describes the measurement protocol. The figures there are regenerated from the committed data on every documentation build.
- Generated codelets. A small symbolic generator writes straight-line code for DFTs of
size 2–31 (composite sizes split in registers, odd primes in the symmetric form). Codelets are
generated with
@generatedfor the radix, element type and SIMD width actually used. - Complex numbers in SIMD registers. On aarch64 CPUs with the ARMv8.3 complex-number
extension (all Apple silicon), additions with a factor ±i and complex multiplications use the
FCADDandFCMLAinstructions; elsewhere the same operations are built from shuffles. - Stockham chains. A transform is a sequence of radix stages (mostly 4 and 8) without bit reversal. The sequence is part of the plan's type, so it compiles into one function. The radices are chosen by a fixed rule found in offline searches over the decompositions of many lengths; there is no run-time measuring step.
- Multidimensional transforms run the same chains along every dimension: a dimension with stride S is a chain started with stride S, the S neighbouring columns forming the vectorised index. No transposition.
- Lengths with prime factors above 31 split off their smooth part in one Cooley–Tukey step and use Bluestein's algorithm only for the rest.
- Threads split each dimension over slices or blocks of columns; long 1D transforms use the four-step algorithm.
FFTL owes much to FFTW by Matteo Frigo and Steven G. Johnson: the idea of generating codelets with a special-purpose compiler, the real-to-real kinds and their definitions, and the interface that FFTW.jl brought to Julia. The manual's page FFTW and FFTL describes FFTW's history and what FFTL takes from its papers.
FFTL is meant to be the Julia community's FFT. It already covers a large part of what an FFT library needs, and we would rather complete one package together than spread the work over several. Useful contributions:
- measurements on hardware we have not tuned for (other GPUs, other x86 and ARM CPUs), and the cases where FFTL is slow or wrong;
- applications that need a transform, a type or a layout that FFTL does not offer;
- kernels, tuning, documentation and tests.
The manual's Contribute page describes the setup and what a change must come with.
] test FFTL runs a quick selection that covers every feature (about 1.5 min). Set
FFTL_TEST_FULL=true for the full lists and FFTL_TEST_VERBOSE=true for the time of every test
set. CI runs the quick tests on Linux and macOS for pushes to main and pull requests, the full
lists and the oldest supported Julia (1.10) once a week, and builds and tests the C library.
MIT, see LICENSE.