Fast Fourier Transform (FFT): Theory and Python Implementation

Introduction

In signal processing and data analysis, understanding the frequency composition of a signal is critically important. The Fourier Transform decomposes signals into frequency components, and the Discrete Fourier Transform (DFT) applies this concept to discrete data.

However, computing the DFT naively requires \(O(N^2)\) operations, making it impractical for large datasets. The Fast Fourier Transform (FFT) solves this problem by reducing the computational complexity to \(O(N \log N)\).

This article covers the mathematical definition of the DFT, the mechanics of the FFT algorithm, a from-scratch Python implementation, and practical frequency analysis using NumPy and SciPy.

Discrete Fourier Transform (DFT)

Definition of the DFT

For a discrete signal \(x[n]\) of length \(N\) (\(n = 0, 1, \ldots, N-1\)), the DFT is defined as:

where:

• \(X[k]\): the spectrum at frequency index \(k\) (complex-valued)
• \(j\): the imaginary unit (\(j^2 = -1\))
• \(e^{-j2\pi kn/N}\): the twiddle factor

The magnitude \(|X[k]|\) represents the amplitude at each frequency, and the argument \(\angle X[k]\) represents the phase.

Inverse DFT (IDFT)

The inverse transform that recovers the time-domain signal \(x[n]\) from the spectrum \(X[k]\) is:

The DFT and IDFT are invertible, allowing lossless conversion between the time and frequency domains.

Computational Complexity of the DFT

Computing equation \((1)\) directly requires \(N\) complex multiplications and additions for each of the \(N\) values of \(k\), resulting in \(O(N^2)\) complexity. For \(N = 10^6\) (one million samples), this means \(10^{12}\) operations, which is computationally infeasible.

The FFT Algorithm (Cooley-Tukey Method)

Core Idea

The FFT algorithm, published by Cooley and Tukey in 1965, uses a divide-and-conquer strategy to compute the DFT efficiently. The key idea is to split a length-\(N\) DFT into two length-\(N/2\) DFTs based on even and odd indices.

Separating the input signal \(x[n]\) into even- and odd-indexed samples:

The DFT can then be decomposed as:

Here \(E[k]\) and \(O[k]\) are length-\(N/2\) DFTs of the even and odd subsequences, respectively. Using the twiddle factor \(W_N^k = e^{-j2\pi k/N}\):

Equation \((6)\) follows from the symmetry property \(W_N^{k+N/2} = -W_N^k\).

The Butterfly Operation

The pair of equations \((5)\) and \((6)\) is called a butterfly operation. It produces two outputs \((X[k], X[k+N/2])\) from one pair of inputs \((E[k], O[k])\), named for the shape of its signal flow diagram.

Key properties of the butterfly operation:

• Each butterfly requires one complex multiplication and two complex additions/subtractions
• \(N/2\) butterfly operations form one stage
• Recursive subdivision yields \(\log_2 N\) stages in total

Complexity Reduction

Splitting a length-\(N\) DFT into two length-\(N/2\) DFTs recursively gives the recurrence:

Solving this yields \(T(N) = O(N \log N)\). For \(N = 10^6\), this reduces from \(O(N^2) = 10^{12}\) to \(O(N \log N) \approx 2 \times 10^7\), a speedup of roughly 50,000 times.

Recursive FFT Implementation in Python

To build intuition for the algorithm, here is a recursive implementation in Python. This implementation assumes \(N\) is a power of 2.

import numpy as np

def fft_recursive(x):
    """Recursive FFT implementation (Cooley-Tukey radix-2)"""
    N = len(x)
    if N == 1:
        return x

    # Split into even and odd indices, recurse
    even = fft_recursive(x[0::2])
    odd = fft_recursive(x[1::2])

    # Compute twiddle factors
    T = [np.exp(-2j * np.pi * k / N) * odd[k] for k in range(N // 2)]

    # Butterfly operations
    return [even[k] + T[k] for k in range(N // 2)] + \
           [even[k] - T[k] for k in range(N // 2)]

# Verification: compare with NumPy's FFT
x = np.random.randn(256)
X_custom = np.array(fft_recursive(list(x)))
X_numpy = np.fft.fft(x)
print(f"Max error: {np.max(np.abs(X_custom - X_numpy)):.2e}")
# Example output: Max error: 1.42e-12

While this implementation is useful for understanding the algorithm, the Python list operations introduce significant overhead. For practical use, always use np.fft.fft.

Practical Frequency Analysis with NumPy/SciPy

Basic Usage

NumPy provides the np.fft module for efficient FFT computation.

import numpy as np

# Generate a signal
fs = 1000            # Sampling frequency [Hz]
t = np.arange(0, 1, 1/fs)  # 1 second time axis
signal = np.sin(2 * np.pi * 50 * t)  # 50 Hz sine wave

# Compute FFT
N = len(signal)
X = np.fft.fft(signal)           # Complex spectrum
freqs = np.fft.fftfreq(N, 1/fs)  # Frequency axis [Hz]

np.fft.fftfreq(N, d) returns the frequency values corresponding to each frequency bin. The parameter \(d\) is the sampling interval (\(1/f_s\)).

Frequency Analysis of a Composite Signal

In practice, signals often consist of multiple frequency components with additive noise. The following example visualizes the spectrum of a signal containing 50 Hz and 120 Hz sine waves with Gaussian noise.

import numpy as np
import matplotlib.pyplot as plt

# --- Generate signal ---
fs = 1000  # Sampling frequency [Hz]
t = np.arange(0, 1, 1/fs)
signal = np.sin(2 * np.pi * 50 * t) + 0.5 * np.sin(2 * np.pi * 120 * t)
signal += 0.3 * np.random.randn(len(t))

# --- FFT ---
N = len(signal)
X = np.fft.fft(signal)
freqs = np.fft.fftfreq(N, 1/fs)

# --- Plot amplitude spectrum (positive frequencies only) ---
plt.figure(figsize=(10, 4))
plt.plot(freqs[:N//2], 2/N * np.abs(X[:N//2]))
plt.xlabel('Frequency [Hz]')
plt.ylabel('Amplitude')
plt.title('FFT Amplitude Spectrum')
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

The factor 2/N is applied because the FFT output is a two-sided spectrum; when displaying only positive frequencies, the amplitude must be doubled (except for the DC and Nyquist components).

The spectrum clearly shows peaks at 50 Hz and 120 Hz, with the 120 Hz peak at approximately half the amplitude of the 50 Hz peak.

Windowing

Spectral Leakage

The FFT implicitly assumes that the finite-length signal is periodic. When the signal is discontinuous at the edges of the analysis window, frequency components that do not actually exist appear in the spectrum. This phenomenon is called spectral leakage.

For example, a 50.0 Hz signal sampled at 1000 Hz for exactly 1 second (1000 samples) produces no leakage because the signal completes an integer number of cycles. However, a 50.5 Hz signal creates discontinuities at the window boundaries, causing energy to spread to frequencies other than 50.5 Hz.

Mitigating Leakage with Window Functions

Window functions smoothly taper the signal toward zero at both ends, reducing discontinuities. A commonly used window is the Hann window:

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal.windows import hann

# Signal with spectral leakage (50.5 Hz)
fs = 1000
t = np.arange(0, 1, 1/fs)
signal = np.sin(2 * np.pi * 50.5 * t)
N = len(signal)

# Without window
X_no_win = np.fft.fft(signal)

# With Hann window
window = hann(N)
X_hann = np.fft.fft(signal * window)
freqs = np.fft.fftfreq(N, 1/fs)

# Comparison plot
fig, axes = plt.subplots(2, 1, figsize=(10, 6), sharex=True)

axes[0].plot(freqs[:N//2], 20 * np.log10(np.abs(X_no_win[:N//2]) + 1e-12))
axes[0].set_ylabel('Magnitude [dB]')
axes[0].set_title('Without Window (Rectangular)')
axes[0].set_xlim(0, 200)
axes[0].grid(True, alpha=0.3)

axes[1].plot(freqs[:N//2], 20 * np.log10(np.abs(X_hann[:N//2]) + 1e-12))
axes[1].set_xlabel('Frequency [Hz]')
axes[1].set_ylabel('Magnitude [dB]')
axes[1].set_title('With Hann Window')
axes[1].set_xlim(0, 200)
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

Applying a window function slightly widens the main lobe but significantly suppresses the side lobes, effectively reducing spectral leakage.

Time-Frequency Analysis with STFT

Why STFT Is Needed

A standard FFT computes a single spectrum for the entire signal, losing information about “when” each frequency component was present. For signals whose frequency content changes over time, such as speech or vibration data, the Short-Time Fourier Transform (STFT) is needed.

The STFT segments the signal into short windows and applies the FFT to each segment, producing a time-frequency representation (spectrogram).

where \(L\) is the window length, \(H\) is the hop size (window shift), and \(w[n]\) is the window function.

Chirp Signal Analysis Example

A chirp signal, whose frequency changes linearly over time, is a classic test case for time-frequency analysis.

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import chirp, spectrogram

# Generate a chirp signal (0 Hz -> 500 Hz over 2 seconds)
fs = 2000  # Sampling frequency [Hz]
t = np.arange(0, 2, 1/fs)
signal = chirp(t, f0=0, f1=500, t1=2, method='linear')

# Compute spectrogram
f, t_spec, Sxx = spectrogram(signal, fs=fs, nperseg=256, noverlap=128)

# Visualization
plt.figure(figsize=(10, 5))
plt.pcolormesh(t_spec, f, 10 * np.log10(Sxx + 1e-12), shading='gouraud')
plt.colorbar(label='Power [dB]')
plt.xlabel('Time [s]')
plt.ylabel('Frequency [Hz]')
plt.title('Spectrogram of Linear Chirp Signal')
plt.tight_layout()
plt.show()

The spectrogram clearly shows the linearly increasing frequency as a diagonal line. Increasing nperseg (window length) improves frequency resolution but reduces time resolution, and vice versa. This is a fundamental trade-off governed by the uncertainty principle.

Summary

This article covered the definition of the DFT, the mechanics of the Cooley-Tukey FFT algorithm, and practical frequency analysis using NumPy and SciPy.

• DFT decomposes a signal into frequency components with \(O(N^2)\) complexity
• FFT reduces the complexity to \(O(N \log N)\) via divide-and-conquer
• Window functions mitigate spectral leakage
• STFT enables time-frequency analysis for non-stationary signals

The FFT is a cornerstone algorithm in signal processing, with applications spanning filter design, audio analysis, image processing, and telecommunications.

Related Articles

• Frequency Characteristics of the Exponential Moving Average (EMA) Filter - Derives the EMA filter transfer function using the Z-transform and analyzes its frequency response.
• Low-Pass Filter Design and Comparison - Compares moving average, Butterworth, and Chebyshev filter frequency responses using FFT-based methods.
• Types and Comparison of Moving Average Filters - Explains the differences in frequency characteristics of SMA, WMA, and EMA.
• Fundamentals of Filtering Methods in Signal Processing - Provides a systematic overview of Kalman filter, EKF, UKF, and particle filter.
• Matplotlib Practical Tips: Creating Publication-Quality Plots - Covers settings for producing higher-quality spectrum plots.

References

• Cooley, J. W., & Tukey, J. W. (1965). “An algorithm for the machine calculation of complex Fourier series”. Mathematics of Computation, 19(90), 297-301.
• Oppenheim, A. V., & Schafer, R. W. (2009). Discrete-Time Signal Processing (3rd ed.). Prentice Hall.
• NumPy FFT documentation
• SciPy Signal Processing documentation