This is the implementation guide for the Fast Fourier Transform (FFT). For theoretical background, including the DFT definition, frequency bins, magnitude and phase, see the article on Fourier analysis.
All code examples are written in C for maximum performance and portability. The implementation below follows the Cooley–Tukey decimation‑in‑time (DIT) algorithm.
Direct Implementation — Iterative FFT
The following function computes the FFT of a complex input array using an in‑place iterative algorithm. It uses bit‑reversal permutation followed by iterative butterfly stages.
#include <math.h>
#include <stdlib.h>
#include <string.h>
/**
* Complex number structure.
*/
typedef struct {
double re;
double im;
} complex_t;
/**
* Bit‑reversal permutation — rearranges input array for iterative FFT.
*/
static void bit_reverse(complex_t *x, int N) {
int i, j, k;
for (i = 1, j = 0; i < N; i++) {
int bit = N >> 1;
for (; j & bit; bit >>= 1) {
j ^= bit;
}
j ^= bit;
if (i < j) {
complex_t temp = x[i];
x[i] = x[j];
x[j] = temp;
}
}
}
/**
* Cooley–Tukey iterative FFT (in‑place).
*
* @param x Complex input array (length N) — overwritten with FFT output
* @param N Length of the array (must be a power of two)
* @param invert 0 for forward FFT, 1 for inverse FFT
*/
void fft(complex_t *x, int N, int invert) {
bit_reverse(x, N);
for (int len = 2; len <= N; len <<= 1) {
double angle = 2.0 * M_PI / len;
if (invert) {
angle = -angle;
}
complex_t wlen = { cos(angle), sin(angle) };
for (int i = 0; i < N; i += len) {
complex_t w = { 1.0, 0.0 };
for (int j = 0; j < len / 2; j++) {
complex_t u = x[i + j];
complex_t v = {
w.re * x[i + j + len/2].re - w.im * x[i + j + len/2].im,
w.re * x[i + j + len/2].im + w.im * x[i + j + len/2].re
};
x[i + j].re = u.re + v.re;
x[i + j].im = u.im + v.im;
x[i + j + len/2].re = u.re - v.re;
x[i + j + len/2].im = u.im - v.im;
/* Update twiddle factor */
double next_re = w.re * wlen.re - w.im * wlen.im;
double next_im = w.re * wlen.im + w.im * wlen.re;
w.re = next_re;
w.im = next_im;
}
}
}
/* Normalise for inverse FFT */
if (invert) {
for (int i = 0; i < N; i++) {
x[i].re /= N;
x[i].im /= N;
}
}
}
FFT Helper Functions
For convenience, the following helpers convert real signals to complex arrays and compute the magnitude spectrum.
/**
* Convert real signal to complex array (imaginary part set to zero).
*/
void real_to_complex(const double *x, complex_t *out, int N) {
for (int i = 0; i < N; i++) {
out[i].re = x[i];
out[i].im = 0.0;
}
}
/**
* Compute magnitude spectrum from FFT output (positive frequencies only).
*/
void magnitude_spectrum(const complex_t *x, int N, double *mag) {
int half = N / 2;
for (int i = 0; i < half; i++) {
mag[i] = sqrt(x[i].re * x[i].re + x[i].im * x[i].im);
}
}
Usage Example
The following example demonstrates how to compute the FFT of a real signal, extract its magnitude spectrum, and perform an inverse FFT to recover the original.
#include <stdio.h>
#include <math.h>
#define N 8
int main() {
double signal[N] = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0};
complex_t x[N];
double mag[N/2];
/* Convert to complex */
real_to_complex(signal, x, N);
/* Forward FFT */
fft(x, N, 0);
/* Compute magnitude spectrum */
magnitude_spectrum(x, N, mag);
printf("Magnitude spectrum:\n");
for (int i = 0; i < N/2; i++) {
printf("mag[%d] = %.6f\n", i, mag[i]);
}
/* Inverse FFT to recover original */
fft(x, N, 1);
printf("Recovered signal (real part):\n");
for (int i = 0; i < N; i++) {
printf("x[%d] = %.6f\n", i, x[i].re);
}
return 0;
}
Expected Output
Magnitude spectrum: mag[0] = 36.000000 mag[1] = 9.061284 mag[2] = 4.000000 mag[3] = 2.204964 Recovered signal (real part): x[0] = 1.000000 x[1] = 2.000000 x[2] = 3.000000 x[3] = 4.000000 x[4] = 5.000000 x[5] = 6.000000 x[6] = 7.000000 x[7] = 8.000000
Implementation Notes
- Bit‑reversal: Required for in‑place computation. The input is rearranged before the butterfly stages.
- Twiddle factors: Precomputed per stage using the recurrence w = w · wlen to avoid repeated trigonometric calls.
- In‑place: The algorithm overwrites the input array — memory efficient.
- Power of two: N must be a power of two. For other lengths, zero‑pad or use a general‑purpose DFT.
- Normalisation: The inverse FFT divides by N to recover the original amplitude.
Optimisation: Precomputed Bit‑Reversal
For repeated FFTs on the same length, bit‑reversal indices can be precomputed and reused.
static int *bit_rev_table = NULL;
/**
* Generate bit‑reversal table for a given N.
*/
void bit_reverse_table_init(int N) {
if (bit_rev_table) free(bit_rev_table);
bit_rev_table = malloc(N * sizeof(int));
int bits = log2(N);
for (int i = 0; i < N; i++) {
int rev = 0;
for (int b = 0; b < bits; b++) {
if (i & (1 << b)) {
rev |= (1 << (bits - 1 - b));
}
}
bit_rev_table[i] = rev;
}
}
/**
* Faster bit‑reversal using precomputed table.
*/
static void bit_reverse_table(complex_t *x, int N) {
for (int i = 0; i < N; i++) {
int j = bit_rev_table[i];
if (i < j) {
complex_t temp = x[i];
x[i] = x[j];
x[j] = temp;
}
}
}
Alternative: Recursive FFT
For clarity, the recursive Cooley–Tukey algorithm is often easier to understand, though less efficient due to function call overhead.
/**
* Recursive FFT (not in‑place).
*/
void fft_recursive(const complex_t *x, int N, complex_t *out, int invert) {
if (N == 1) {
out[0] = x[0];
return;
}
complex_t even[N/2], odd[N/2];
for (int i = 0; i < N/2; i++) {
even[i] = x[2 * i];
odd[i] = x[2 * i + 1];
}
complex_t even_out[N/2], odd_out[N/2];
fft_recursive(even, N/2, even_out, invert);
fft_recursive(odd, N/2, odd_out, invert);
double sign = invert ? 1.0 : -1.0;
for (int k = 0; k < N/2; k++) {
double angle = sign * 2.0 * M_PI * k / N;
complex_t twiddle = { cos(angle), sin(angle) };
out[k].re = even_out[k].re + twiddle.re * odd_out[k].re - twiddle.im * odd_out[k].im;
out[k].im = even_out[k].im + twiddle.re * odd_out[k].im + twiddle.im * odd_out[k].re;
out[k + N/2].re = even_out[k].re - twiddle.re * odd_out[k].re + twiddle.im * odd_out[k].im;
out[k + N/2].im = even_out[k].im - twiddle.re * odd_out[k].im - twiddle.im * odd_out[k].re;
}
}
Comparison of Methods
Method Complexity Memory Best for Iterative O(N log N) In‑place General use (fastest) Recursive O(N log N) O(N) Education, clarity Precomputed O(N log N) In‑place Repeated FFTs of same N
See also: Fourier Analysis | Convolution | Convolution Implementation