zldsp::fft

zldsp::fft is a header-only C++ Fast Fourier Transform (FFT) library built with Google Highway. zldsp::fft is licensed under Apache-2.0 license. You can obtain the source code at GitHub. zldsp::fft supports SSE2/SSE4/AVX2/NEON SIMD targets and both AoS/SoA layouts for complex numbers.

Usage

1. C++ Standard: C++20 or higher
2. Google Highway: You must include and link Google Highway in your project. The headers in this library expect the following includes to be resolvable:

#include <hwy/aligned_allocator.h>
#include <hwy/highway.h>

Compiler Flags

To activate SIMD, you must set the correct compilation flags. Highway will automatically detect the target architecture from these flags.

API

Include the relevant headers from the /src directory.

Both CFFT and RFFT are templated on the floating-point type (e.g., float or double) and are instantiated using the base-2 order of the FFT (where size = 1 << order).

Complex FFT (CFFT)

#include "src/zldsp_fft_cfft.hpp"
#include <vector>
#include <complex>

// a CFFT of size 1024 (2^10) using float
constexpr size_t order = 10; 
zldsp::fft::CFFT<float> cfft(order);

std::vector<std::complex<float>> in_buffer(1 << order);
std::vector<std::complex<float>> out_buffer(1 << order);

// forward transform (AoS to AoS)
cfft.forward(in_buffer.data(), out_buffer.data());

// backward transform (AoS to AoS)
cfft.backward(out_buffer.data(), in_buffer.data());

Real FFT (RFFT)

#include "src/zldsp_fft_rfft.hpp"
#include <vector>
#include <complex>

// a RFFT of size 1024 (2^10) using float
constexpr size_t order = 10;
zldsp::fft::RFFT<float> rfft(order);

std::vector<float> real_in(1 << order);
std::vector<std::complex<float>> complex_out((1 << order) / 2 + 1);
std::vector<float> sqr_mag_out((1 << order) / 2 + 1);

// forward transform (real to AoS)
rfft.forward(real_in.data(), complex_out.data());

// backward transform (AoS to real)
rfft.backward(complex_out.data(), real_in.data());

// forward transform (real to squared magnitude)
rfft.forward_sqr_mag(real_in.data(), sqr_mag_out.data())

Data Layouts (AoS/SoA)

Both CFFT and RFFT support AoS/SoA for complex numbers:

AoS stores interleaved real/imaginary values in one array:

std::vector<std::complex<float>> out_buffer(1 << order);

SoA stores continuous real/imaginary values in two arrays:

std::vector<float> out_real(1 << order);
std::vector<float> out_imag(1 << order);

Example using SoA:

#include "src/zldsp_fft_cfft.hpp"
#include <vector>
#include <complex>

// a CFFT of size 1024 (2^10) using float
constexpr size_t order = 10; 
zldsp::fft::CFFT<float> cfft(order);

std::vector<std::complex<float>> in_buffer(1 << order);
std::vector<float> out_real(1 << order);
std::vector<float> out_imag(1 << order);

// forward transform (AoS to SoA)
cfft.forward(in_buffer.data(), {out_real.data(), out_imag.data()});

// backward transform (SoA to AoS)
cfft.backward({out_real.data(), out_imag.data()}, in_buffer.data());

Benchmark

Here are some benchmarks results. The benchmark code is available at GitHub. I have tried my best to configure other libraries correctly. All libraries are built with LLVM/Clang or Apple-Clang. Google Benchmark is used for benchmark.

Included libraries are (they are subjected to their own licenses):

• IPP: Intel® Integrated Performance Primitives
• vDSP: Apple vDSP
• ArmPL: Arm Performance Libraries
• FFTW3: FFTW 3.3.10 Mirror with NEON support, take the better one from FFTW_MEASURE and FFTW_ESTIMATE
• KFR: KFR 7.0.0
• PFFFT: PFFFT 1.1.0
• zldsp: zldsp::fft

The benchmarks shown below are Real FFT benchmarks from order 5 to order 25.

Apple M4 Pro 14-core (local)

Intel Core i7-8850H (local)

AMD EPYC 7763 (GitHub Actions runner images)

AMD AMD EPYC 9V74 (GitHub Actions runner images)

Ampere Altra Q80-30 (Oracle Cloud VM)

Design

Low Order

For low orders, 0 <= order <= 5, zldsp_fft uses specialized CFFT kernels. Each transform size therefore has a fixed implementation without general stage iteration or stage dispatch to keep values in SIMD registers where practical.

Medium Order

For medium orders, 6 <= order <= switch_order - 1, zldsp_fft uses out-of-place Stockham DIT design. It is scheduled as follows:

• an even order begins with Radix-4 and following by Radix-4 stages
• an odd order begins with Radix-8 and following by Radix-4 stages

Intermediate data uses an SIMD-friendly AoSoA layout. The first/final stage is fused with conversion from/to the requested AoS or SoA input/output.

High Order

For high orders, switch_order <= order, zldsp_fft uses a hybrid Cooley–Tukey & Stockham design. It uses Radix-4 Cooley–Tukey DIF macro stages to split the transform into micro CFFTs sized to fit approximately within L1 cache. Then, each micro CFFT is completed with the medium-order Stockham DIT design. The resulting matrix is placed in natural order using a reusable terminal buffer and SIMD register-transpose tiles.

Switch Order

The boundary between medium and high orders is derived from the detected L1 and L2 cache sizes. The implementation estimates the largest transform orders that fit its working sets and selects:

switch_order = max(maximum_L1_order + 4, maximum_L2_order)

Reference

• Van Loan, Charles. Computational frameworks for the fast Fourier transform. Society for Industrial and Applied Mathematics, 1992.
• Notes on FFTs: for implementers
• OTFFT documentation

• ←Previous
• Next→