Project Nayuki

Free small FFT in multiple languages

Introduction

The fast Fourier transform (FFT) is a versatile tool for digital signal processing (DSP) algorithms and applications. On this page, I provide a free implemen­tation of the FFT in multiple languages, small enough that you can even paste it directly into your application (you don’t need to treat this code as an external library).

Also included is a fast circular convolution function based on the FFT. Note that the FFT, with a bit of pre- and postprocessing, can quickly calculate the discrete cosine transform (DCT), which is used in many multimedia compres­sion algorithms.

My implementation has a reasonable amount of optimization (such as build­ing trigonometric tables), but is not intended to squeeze every last drop of performance. Instead, the implementation is optimized for clarity and conciseness. You’re welcome to add more optimizations on your own if you wish.

Math documentation

Note: In the text below, all input and output vectors are complex, and \(n\) is the length of each vector.

The forward FFT is based on this formula: \( \text{Out}(k) = \displaystyle\sum_{t=0}^{n-1} \, \text{In}(t) \, e^{-2 π i t k / n} \).

The inverse FFT is based on this formula: \( \text{Out}(k) = \displaystyle\sum_{t=0}^{n-1} \, \text{In}(t) \, e^{2 π i t k / n} \).

This FFT does not perform any scaling. So for a vector of length \(n\), after performing a transform and an inverse transform on it, the result will be the original vector multiplied by \(n\) (plus approximation errors).

The circular convolution is based on this formula: \( \text{Out}(j) = \displaystyle\sum_{k=0}^{n-1} \, \text{X}((j - k) \text{ mod } n) \, \text{Y}(k) \).
In other words, \(\text{Out}(j)\) is the sum of all the products \(\text{X}(a) \text{Y}(b)\) where the indices satisfy \(a + b ≡ j \text{ mod } n\).

Source code and API

The test code are runnable main programs that check the FFT and fast convo­lution for correctness against a naive algorithm. They also demonstrate how this FFT code is called. You do not need to include the test code into your own program. Also, feel free to strip down the set of FFT functions to the bare minimum that you need.

The radix-2 FFT code is essential since every function depends on it. The Bluestein FFT code can be omitted if you are not doing FFTs or convolutions of non-power-of-2 sizes. The convolution code can be omitted if you are not computing convolutions and you are not using the Bluestein FFT.

If you are running many FFTs on the same array length, please consider modifying the code to compute the cosine and sine tables once and reuse them. This is especially important when benchmarking the speed of the algorithm.

License summary (MIT License): Please keep the copyright notice and URL intact. You may use, edit, and publish the code for any purpose including commercial use. No warranty is provided by Nayuki. The official license text is included in the source code.

Java (SE 5+) (pairs of reals)

Download: Fft.java, FftTest.java

  • Forward FFT (wrapper) (in-place): void transform(double[] real, double[] imag)
  • Inverse FFT (wrapper) (in-place): void inverseTransform(double[] real, double[] imag)
  • Forward FFT (radix-2) (in-place): void transformRadix2(double[] real, double[] imag)
  • Forward FFT (Bluestein) (in-place): void transformBluestein(double[] real, double[] imag)
  • Circular convolution (real): void convolve(double[] xvec, double[] yvec, double[] outvec)
  • Circular convolution (complex): void convolve(double[] xreal, double[] ximag, double[] yreal, double[] yimag, double[] outreal, double[] outimag)

When calling any of the above functions, all the argument arrays must have the same length.

JavaScript / TypeScript (pairs of reals)

Download: fft.js / fft.ts, fft-test.js, fft-test.xhtml

  • Forward FFT (wrapper) (in-place): transform(real: FloatArray, imag: FloatArray): void
  • Inverse FFT (wrapper) (in-place): inverseTransform(real: FloatArray, imag: FloatArray): void
  • Forward FFT (radix-2) (in-place): transformRadix2(real: FloatArray, imag: FloatArray): void
  • Forward FFT (Bluestein) (in-place): transformBluestein(real: FloatArray, imag: FloatArray): void
  • Circular convolution (real): convolveReal(xvec: FloatArray, yvec: FloatArray, outvec: FloatArray): void
  • Circular convolution (complex): convolveComplex(xreal: FloatArray, ximag: FloatArray, yreal: FloatArray, yimag: FloatArray, outreal: FloatArray, outimag: FloatArray): void

For each function above, type FloatArray = Array<number> | Float64Array;. When calling any of the above functions, all the argument arrays must have the same length. The API is almost the same as the Java port’s.

Python (native complex)

Download: fft.py, fft-test.py

  • Forward FFT (wrapper): transform(list of numeric, False) returning a list of complex
  • Inverse FFT (wrapper): transform(list of numeric, True) returning a list of complex
  • Forward FFT (radix-2): transform_radix2(list of numeric, False) returning a list of complex
  • Inverse FFT (radix-2): transform_radix2(list of numeric, True) returning a list of complex
  • Forward FFT (Bluestein): transform_bluestein(list of numeric, False) returning a list of complex
  • Inverse FFT (Bluestein): transform_bluestein(list of numeric, True) returning a list of complex
  • Circular convolution (real): convolve(list of int/long/float, list of int/long/float, True) returning a list of float
  • Circular convolution (complex): convolve(list of numeric, list of numeric, False) returning a list of complex

The numeric types are int, long, float, and complex. When calling any of the above functions, the returned list has the same length as the argument list(s). When calling convolve(), the two argument lists must have the same length. Note that this Python implementation returns the result instead of storing it back into the argument list(s); this is unlike the implementations for the other languages.

C# (.NET 4.0+) (native complex)

Download: Fft.cs, FftTest.cs

  • Forward FFT (wrapper) (in-place): void Transform(Complex[] vec, false)
  • Inverse FFT (wrapper) (in-place): void Transform(Complex[] vec, true)
  • Forward FFT (radix-2) (in-place): void TransformRadix2(Complex[] vec, false)
  • Inverse FFT (radix-2) (in-place): void TransformRadix2(Complex[] vec, true)
  • Forward FFT (Bluestein) (in-place): void TransformBluestein(Complex[] vec, false)
  • Inverse FFT (Bluestein) (in-place): void TransformBluestein(Complex[] vec, true)
  • Circular convolution (complex): void Convolve(Complex[] xvec, Complex[] yvec, Complex[] outvec)

The convolution function requires all argument arrays to have the same length. Because C# supports complex numbers natively, the codebase is a blend of the Java and Python ports.

C (C99 and above) (native complex)

Download: fft-complex.c, fft-complex.h, fft-complex-test.c

  • Forward FFT (wrapper) (in-place): bool Fft_transform(double complex vec[], size_t n, false)
  • Inverse FFT (wrapper) (in-place): bool Fft_transform(double complex vec[], size_t n, true)
  • Forward FFT (radix-2) (in-place): bool Fft_transformRadix2(double complex vec[], size_t n, false)
  • Inverse FFT (radix-2) (in-place): bool Fft_transformRadix2(double complex vec[], size_t n, true)
  • Forward FFT (Bluestein) (in-place): bool Fft_transformBluestein(double complex vec[], size_t n, false)
  • Forward FFT (Bluestein) (in-place): bool Fft_transformBluestein(double complex vec[], size_t n, true)
  • Circular convolution: bool Fft_convolveComplex(const double complex xvec[], const double complex yvec[], double complex outvec[], size_t n)

When calling any of the above functions, all the argument arrays must have the same length. Each function returns true to indicate success or false to indicate failure.

C (C99 and above) (pairs of reals)

Download: fft-real-pair.c, fft-real-pair.h, fft-real-pair-test.c

  • Forward FFT (wrapper) (in-place): bool Fft_transform(double real[], double imag[], size_t n)
  • Inverse FFT (wrapper) (in-place): bool Fft_inverseTransform(double real[], double imag[], size_t n)
  • Forward FFT (radix-2) (in-place): bool Fft_transformRadix2(double real[], double imag[], size_t n)
  • Forward FFT (Bluestein) (in-place): bool Fft_transformBluestein(double real[], double imag[], size_t n)
  • Circular convolution (real): bool Fft_convolveReal(const double xvec[], const double yvec[], double outvec[], size_t n)
  • Circular convolution (complex): bool Fft_convolveComplex(const double xreal[], const double ximag[], const double yreal[], const double yimag[], double outreal[], double outimag[], size_t n)

When calling any of the above functions, all the argument arrays must have the same length. Each function returns true to indicate success or false to indicate failure.

C++ (C++11 and above) (native complex)

Download: FftComplex.cpp, FftComplex.hpp, FftComplexTest.cpp

  • Forward FFT (wrapper) (in-place): void transform(vector<complex<double> > &vec, false)
  • Inverse FFT (wrapper) (in-place): void transform(vector<complex<double> > &vec, true)
  • Forward FFT (radix-2) (in-place): void transformRadix2(vector<complex<double> > &vec, false)
  • Inverse FFT (radix-2) (in-place): void transformRadix2(vector<complex<double> > &vec, true)
  • Forward FFT (Bluestein) (in-place): void transformBluestein(vector<complex<double> > &vec, false)
  • Inverse FFT (Bluestein) (in-place): void transformBluestein(vector<complex<double> > &vec, true)
  • Circular convolution: vector<complex<double> > convolve(vector<complex<double> > xvec, vector<complex<double> > yvec)

When calling the convolve() function, all the argument vectors must have the same length.

C++ (C++11 and above) (pairs of reals)

Download: FftRealPair.cpp, FftRealPair.hpp, FftRealPairTest.cpp

  • Forward FFT (wrapper) (in-place): void transform(vector<double> &real, vector<double> &imag)
  • Inverse FFT (wrapper) (in-place): void inverseTransform(vector<double> &real, vector<double> &imag)
  • Forward FFT (radix-2) (in-place): void transformRadix2(vector<double> &real, vector<double> &imag)
  • Forward FFT (Bluestein) (in-place): void transformBluestein(vector<double> &real, vector<double> &imag)
  • Circular convolution (real): vector<double> convolve(vector<double> xvec, vector<double> yvec)
  • Circular convolution (complex): pair<vector<double>, vector<double> > convolve(vector<double> xreal, vector<double> ximag, vector<double> yreal, vector<double> yimag)

When calling any of the above functions, all the argument vectors must have the same length.

Rust (pairs of reals)

Download: fft.rs, fft-test.rs

  • Forward FFT (wrapper) (in-place): transform(real: &mut [f64], imag: &mut [f64])
  • Inverse FFT (wrapper) (in-place): inverse_transform(real: &mut [f64], imag: &mut [f64])
  • Forward FFT (radix-2) (in-place): transform_radix2(real: &mut [f64], imag: &mut [f64])
  • Forward FFT (Bluestein) (in-place): transform_bluestein(real: &mut [f64], imag: &mut [f64])
  • Circular convolution (real): convolve_real(xvec: &[f64], yvec: &[f64], outvec: &mut [f64])
  • Circular convolution (complex): convolve_complex(xreal: &[f64], ximag: &[f64], yreal: &[f64], yimag: &[f64], outreal: &mut [f64], outimag: &mut [f64])

When calling any of the above functions, all the argument slices must have the same length. Each function returns nothing (unit/void).

More info