Writing your own FFT implementation

Why do I care about the FFT anyway?

I was first exposed to the DFT (Discrete Fourier Transform) in the Signals class of my electrical engineering undergrad program. I did not realize it at the time, but this single math operation would sit at the core of my professional life for the next seven years.

The DFT, and its optimized algorithmic implementation the FFT, is used everywhere:

Filtering
Visualization of signals via the STFT
Array processing and beamforming via spatial frequency methods
Autocorrelation and cross-correlation
Convolution
Audio pre-processing and the mel spectrogram
Certain communication modulations (frequency division multiplexing)

And that is just the stuff I’ve personally used it for.

Some interesting videos on the history of the FFT, and an explanation of the Fourier transform:

Of course, like any sane engineer, I never implemented one myself. I reached for fft() in MATLAB, numpy.fft in Python, FFTW if I was writing C++, or Nvidia’s cuFFT on the GPU.

I was always using someone else’s implementation behind someone else’s wrapper. But I wanted to know what was actually happening under the hood, and how much performance I have been leaving on the table by picking one implementation over another.

The state of the FFT implementation space on Intel/AMD x86-64

Nearly all production FFT implementations are written in C or C++ and all expose a C API.

Three implementations dominate x86:

Intel MKL/IPP Libraries - the Math Kernel Library dates to the mid-1990s and the Integrated Performance Primitives arrived around 2000. MKL generally focuses on matrix math and batched ahead-of-time computation, while IPP focuses more on signal processing and real time. MKL/IPP has a permissive use license (the Intel Simplified Software License) [1], but is closed source.
FFTW “Fastest Fourier Transform in the West” - The big open source implementation. The current major version FFTW 3 first started shipping in 2003 [3]. It is GPL licensed, though you can purchase a commercial use license [4].
PocketFFT and its successor ducc0.fft by the same author. The core FFT implementation used by numpy and scipy’s fft function. PocketFFT is BSD-3-Clause, and while the wider ducc0 package is GPL, its FFT component (ducc0.fft) is likewise offered under the BSD-3 license [5].

A few others:

RustFFT - A pure Rust SIMD implementation of the FFT algorithm (MIT or Apache-2.0) [6].
KFR - Third party dual licensed GPL (v2/v3) and commercial signal processing library [7].
AOCL-FFTZ - AMD’s vendor library FFT implementation, permissive use license (BSD 3-Clause) [8], parallel to Intel MKL/IPP.

The real question is the tradeoff between ease of use, license cost, performance, and portability.

My personal interest is in low-latency audio and signal processing, so I chose a 256 length FFT, primarily because:

It is a fairly well used power of 2, so libraries that do optimize for it should have good implementations to compare to.
It is small enough to fit entirely in L1 cache, which means that there are fewer memory layout concerns when optimizing it.
It is a nice round number 2**8

To measure this, I built a standardized C++ benchmarking harness using Google Benchmark. The setup:

Calls the FFT operation in sequential order on a 1 MiB data vector.
It runs the FFT on 256 interleaved complex float32 samples, then on the next 256 (a stride of 256).
The FFT calculation is out of place and does not modify the input buffer.
The output buffer is similarly streamed through at a stride of 256.

This models an out-of-place streaming setup, where memory utilization is not a primary concern, but modifying the input isn’t allowed.

These measurements are also specific to my 13th Gen Intel(R) Core(TM) i5-1340P P-Core CPU. On the scales of nanoseconds that these FFT implementations are measuring, the particulars of a microarchitecture matter.

We start with the batch-of-1 case: hand the FFT one buffer of 256 samples, take the result, then hand it the next. This is what you do when you can’t afford to wait — in a low-latency streaming system you process each block the moment it arrives, so you only ever pay the cost of a single transform, never the extra delay of collecting a batch first.

n256typeC2Cdirectionforwardprecisionf32normalizenonelayoutpackedbatch1signal1 MiB streamISAAVX2

One transform at a time (n = 256), across every backend. Hover a row for what it is. Lower is better.

13th Gen Intel(R) Core(TM) i5-1340P · fedora-framework · commit 4a9c0c2

In the graph above, PocketFFT and ducc0 are roughly 7x slower than the leaders. These are the backends behind numpy and scipy. When you call numpy.fft under the hood, you are calling PocketFFT. At the top of the graph, we can see that Intel’s libraries are the fastest (of course this is again on my Intel chip), followed by FFTW, which performs very well, as long as you use a decent runtime kernel optimization setting (‘measure’ vs ‘estimate’).

In the next graph, we hand each backend 16 transforms at once instead of 1. This costs 16x the latency (since we have to wait for 4096 samples to accumulate), but lets an FFT implementation vectorize across transforms rather than just within a single transform. Most of the implementations barely move. Both MKL and FFTW are already heavily optimized within a single transform, and there is not much room left to take advantage of optimizing across multiple transforms. On the other hand PocketFFT and ducc0.fft benefit a lot, since they were not already highly vectorized within a single FFT.

n256typeC2Cdirectionforwardprecisionf32normalizenonelayoutpackedbatch16signal1 MiB streamISAAVX2

A batch of sixteen transforms (n = 256). Same rows, order, and scale as above — so each backend’s bar grows or shrinks in place. Lower is better.

13th Gen Intel(R) Core(TM) i5-1340P · fedora-framework · commit 4a9c0c2

An important aspect when measuring the ‘performance’ of an FFT operation (or any software in general) is that it matters what you are optimizing for. It is rare that you are ever doing an FFT operation in Python where you care about latency on the order of ns. Similarly, most of the time when you are doing sequential FFTs you are doing a lot of them. In this specific domain vectorizing over the FFT operations does not really cost anything.

On the other hand, if you are doing low-latency audio you are probably close to your audio source, which means you are already in a lower-level systems language. There you would reach for a C/C++ library directly, choosing one that has spent effort optimizing for short, in-cache FFTs.

So… Why build an FFT library?

You shouldn’t

The problem with building a general-purpose FFT library is not writing a fast FFT. It is writing a fast FFT for the combinatorial explosion of: {transform size × real/complex × float32/float64 × interleaved/planar × in-place/out-of-place × SSE/AVX2/AVX512/NEON/SVE × microarchitecture} The microarchitecture in this is especially painful, since a different cache layout could cause one algorithmic implementation to outperform another. FFTW solves this explosion by measuring performance of different kernels at runtime, on the end user’s CPU. Intel MKL has a bunch of tuned kernels for its CPUs.

There are two scenarios where writing the FFT yourself can make sense (maybe, if you squint)…

You have a very specific use case, with a very specific processor, that needs a very specific FFT.
Kernel Fusion. This comes from the Matrix Algebra world: you fold in elementwise operations you’d have to do anyway while the data is still in registers, saving a round trip through memory. The operations you fuse to the main operation are called epilogues (if they are after the main operation) or prologues (before the main operation). In an FFT these could be a window function or a scaling factor.

And then there is the real reason I am here:

I just really want to learn how to write a fast FFT.

The Plan

To keep myself sane, my goal is very scoped.

My objective is to write a:

256 Length
Complex float32 Floating Point
Out of place
Complex-to-complex
Batch 1
Single Threaded
Interleaved array-of-structures complex data

FFT which outperforms Intel MKL, on my specific hardware:

13th Gen Intel(R) Core(TM) i5-1340P
Single P-Core
Hyperthreading Disabled
Turbo Disabled

The one gimme I will allow myself is for the FFT to be ahead-of-time compiled. Due to the combinatorial explosion, this is unusual in FFTs, and in specific contrast to Intel.

As part of this I also want to try out Mojo, a new systems-level programming language primarily designed for writing kernels for AI. Since I heard about Mojo a few years ago, in the back of my mind, I have thought Sure, AI is pretty cool, but I think you could use this to write a really good signal processing library. The main nicety, in that respect, is the ability to do extensive compile-time metaprogramming, which allows you to write vectorized code in a much more straightforward fashion.

For a lot of signal processing code, people fall back to writing scalar code, because it is easier, especially when you need to target multiple architectures.

Mojo also has a JIT compiler, and Modular has said that they are open-sourcing their compiler later this summer (2026).

That being said, the Mojo language is new, and does not yet have all the niceties that you would want out of a systems language library (there is no std::thread equivalent as of august 2026). But for math kernels, which get embedded in C code and called externally, it is great.