Fast Fourier, Cosine and Sine Transforms

Computes discrete Fourier transforms (complex, real-to-complex and complex-to-real) and discrete cosine and sine transforms of types one to four over any axes of an array, with the fast Fourier transform (FFT) code of 'ducc0' (< https://github.com/mreineck/ducc>) bundled in the BSD-licensed subset that 'SciPy' uses. Transforms of any length run in O(n log n) time, real input uses true real transforms, and several threads can share one transform.

zufft

Overview

zufft wraps ducc0, the C++ library behind scipy.fft, to provide fast Fourier transforms (FFTs) and discrete cosine and sine transforms (DCTs and DSTs):

One difference matters when you switch: zufft scales inverse transforms by 1/n, as SciPy does, while stats and fftwtools don’t. Use norm = "none" to get their numbers. vignette("migrating-from-fftwtools") gives the exact arguments for each fftwtools function.

The conventions for axes, normalisation and the DCT and DST definitions follow scipy.fft. Helpers such as fft_next_len(), fft_freq() and fft_shift() mirror SciPy’s next_fast_len(), fftfreq() and fftshift().

Installation

# Install zufft from CRAN:
install.packages("zufft")

# Or the development version from GitHub:
# install.packages("pak")
pak::pak("pedrobtz/zufft")

Usage

library(zufft)

fft_c2c() computes the full spectrum, just like stats::fft():

x <- c(1, 4, 2, 8, 5, 7)
fft_c2c(x)
#> [1]  27+0.000000e+00i  -5+5.196152e+00i   0+4.440892e-16i -11+0.000000e+00i
#> [5]   0-4.440892e-16i  -5-5.196152e+00i

For real input, the spectrum is Hermitian-symmetric: the negative-frequency terms are the complex conjugates of the positive ones. fft_r2c() returns only the n %/% 2 + 1 non-negative frequency terms, and fft_c2r() takes you back. Give n to recover an odd length: by default fft_c2r() assumes an even one, as SciPy’s irfft() does:

X <- fft_r2c(x)
X
#> [1]  27+0.000000e+00i  -5+5.196152e+00i   0+4.440892e-16i -11+0.000000e+00i
fft_c2r(X, n = length(x))
#> [1] 1 4 2 8 5 7

Use axes to choose the dimensions to transform, and workers to spread the work over several threads (fft_workers() sets a default for the session). Here we transform each column of a matrix on two threads:

m <- matrix(rnorm(4096 * 64), 4096)
spec <- fft_r2c(m, axes = 1, workers = 2)
dim(spec)
#> [1] 2049   64

R stores arrays in column-major order, so the first dimension is contiguous in memory. That’s why, with axes = NULL, the real transforms halve the first dimension (to n %/% 2 + 1): it’s fast, and it matches fftwtools. If you give axes explicitly, the last axis listed is halved instead, as in SciPy. With CRAN’s builds, the results are bit-identical whatever the number of threads.

fft_dct() and fft_dst() compute DCTs and DSTs with SciPy’s definitions. The default is type 2, and the inverse takes you back to x:

fft_dct(x, norm = "ortho")
#> [1] 11.0227038 -4.6508890 -1.0000000  0.4082483  0.0000000 -3.8343924
fft_dct(fft_dct(x), inverse = TRUE)
#> [1] 1 4 2 8 5 7

Lengths with only small prime factors transform fastest. fft_next_len() returns the smallest fast length at or above n, so you can zero-pad the input to it:

fft_next_len(1009)
#> [1] 1024

Speed

You don’t need a native build for zufft to keep up with fftwtools. Here are a few cases on a Linux x86-64 runner, using one thread, against Ubuntu’s FFTW (built with SSE2 and AVX SIMD code):

A ratio below 1 means zufft is faster. On one thread, CRAN’s portable build is faster on 30 of the 33 benchmark cases. The three exceptions are single 1-D complex transforms of about 10⁵ to 10⁶ points, at worst 1.26 times slower. A native build is faster on all 33.

Most of the speed comes from planning: zufft caches its plans, while fftwtools creates and destroys an FFTW plan on every call. Threads (workers) add more, which fftwtools doesn’t offer. And every length runs in O(n log n), even a prime one, unlike in stats::fft().

bench/RESULTS.md has all 33 cases, how they were measured, and the numbers for macOS on Apple silicon.

Native build

A native build is only worth it on Linux x86-64. CRAN’s binaries are compiled to run on any CPU of their architecture, and on x86-64 that means SSE2, whose SIMD registers hold two doubles. If your CPU has AVX2 or AVX-512, a build for your own machine runs faster. It takes two steps in R:

1. Set the environment variable ZUFFT_CXXFLAGS with Sys.setenv(). zufft adds it to its compiler flags, and no other package sees it.

2. Install from source with install.packages("zufft", type = "source"). A prebuilt binary is already compiled and ignores the variable.

Sys.setenv(ZUFFT_CXXFLAGS = "-march=native")
install.packages("zufft", type = "source")

You can check the result with zufft_info(). The SIMD width, in doubles, is 4 with AVX2, 8 with AVX-512, and 2 for a portable build:

zufft::zufft_info()$simd_width

How much does it help? On the Linux x86-64 runners of bench/RESULTS.md, where this gives AVX2, the 16384 × 64 complex column case goes from 0.96 to 0.60 of fftwtools’ time. A complex transform of prime length 10 007 goes from 0.81 to 0.49.

Other platforms gain little or nothing:

• On Apple silicon the SIMD width stays at two doubles (NEON), so a native build gains little.

• Windows builds ignore ZUFFT_CXXFLAGS and stay portable, because code built for AVX by the GCC in Rtools crashed in our tests.

Keep in mind that a native build with GCC, on a CPU with fused multiply-add (FMA), can give results that differ in the last bit between thread counts. zufft::zufft_info()$fp_contract reports "fast" when that applies. CRAN’s builds report "off".

Licence

zufft is MIT-licensed. Its bundled ducc0 files are the “BSD-3-Clause OR GPL-2.0-or-later” subset that SciPy vendors, used under BSD-3-Clause. inst/COPYRIGHTS lists every copyright holder. No GPL-only ducc0 code is included.

Getting help

If you encounter a clear bug, please file an issue with a minimal reproducible example on GitHub.