Skip to content

FastFourierTransform.md

Latest commit

 

History

History
132 lines (86 loc) · 6.92 KB

FastFourierTransform.md

File metadata and controls

132 lines (86 loc) · 6.92 KB
Raw

Fast Fourier Transform (FFT)

Overview & Motivation

Many problems in signal processing, communications, and control require knowing the frequency content of a signal. The Discrete Fourier Transform (DFT) answers this question exactly, but its direct computation costs $O(N^2)$ multiplications — prohibitive for real-time embedded systems processing thousands of samples.

The Fast Fourier Transform is a family of algorithms that compute the DFT in $O(N \log N)$ operations by exploiting symmetry and periodicity of the complex exponential. This library implements the Radix-2 Cooley-Tukey variant, which recursively splits an $N$-point DFT into two $N/2$-point DFTs, requiring $N$ to be a power of two.

The core insight is that the DFT matrix has a recursive block structure: half of its entries are shared between the even-indexed and odd-indexed sub-problems, so most of the work can be reused rather than recomputed.

Mathematical Theory

Prerequisites

The Discrete Fourier Transform of a length-$N$ sequence $x[n]$ is defined as:

$$X[k] = \sum_{n=0}^{N-1} x[n] , e^{-j2\pi kn/N}, \quad k = 0, 1, \ldots, N-1$$

The inverse DFT recovers the original sequence:

$$x[n] = \frac{1}{N} \sum_{k=0}^{N-1} X[k] , e^{j2\pi kn/N}$$

We define the twiddle factor $W_N = e^{-j2\pi/N}$, so $X[k] = \sum_{n} x[n] , W_N^{kn}$.

Cooley-Tukey Decomposition

Split the sum into even-indexed and odd-indexed terms:

$$X[k] = \underbrace{\sum_{m=0}^{N/2-1} x[2m] , W_{N/2}^{km}}_{E[k]} ;+; W_N^k \underbrace{\sum_{m=0}^{N/2-1} x[2m+1] , W_{N/2}^{km}}_{O[k]}$$

Because $E[k]$ and $O[k]$ are periodic with period $N/2$:

$$X[k] = E[k] + W_N^k , O[k]$$ $$X[k + N/2] = E[k] - W_N^k , O[k]$$

This pair of equations is the butterfly operation: two outputs from two inputs and one twiddle-factor multiplication.

Bit-Reversal Permutation

The recursive decomposition reorders input by bit-reversing the sample indices. For in-place computation the input array is pre-permuted so that butterfly stages proceed sequentially from the smallest sub-problems to the full $N$-point result.

Complexity Analysis

Case Time Space Notes
All $O(N \log_2 N)$ $O(N)$ $N$ must be a power of 2

Why $O(N \log N)$: There are $\log_2 N$ butterfly stages, each performing $N/2$ butterfly operations (one complex multiply + two complex adds). Total: $\frac{N}{2} \log_2 N$ complex multiplications.

Twiddle factors are pre-computed and stored in a lookup table of size $N/2$, so no trigonometric evaluations occur at runtime.

Step-by-Step Walkthrough

Input: $x = [1, 2, 3, 4]$, $N = 4$

Step 1 — Bit-reversal permutation

Index (binary) Reversed Value
0 (00) 0 (00) 1
1 (01) 2 (10) 3
2 (10) 1 (01) 2
3 (11) 3 (11) 4

Reordered: $[1, 3, 2, 4]$

Step 2 — Stage 1 (size-2 butterflies), twiddle $W_2^0 = 1$

  • Butterfly$(1, 3)$: outputs $(1+3, ; 1-3) = (4, ; -2)$
  • Butterfly$(2, 4)$: outputs $(2+4, ; 2-4) = (6, ; -2)$

Result: $[4, -2, 6, -2]$

Step 3 — Stage 2 (size-4 butterfly), twiddles $W_4^0 = 1$, $W_4^1 = -j$

  • $X[0] = 4 + 1 \cdot 6 = 10$
  • $X[1] = -2 + (-j)(-2) = -2 + 2j$
  • $X[2] = 4 - 1 \cdot 6 = -2$
  • $X[3] = -2 - (-j)(-2) = -2 - 2j$

Output: $X = [10,; -2+2j,; -2,; -2-2j]$

Verification: $X[0] = 1+2+3+4 = 10$ ✓

Pitfalls & Edge Cases

  • Power-of-2 constraint. Non-power-of-2 inputs require zero-padding, which introduces interpolation artifacts in the frequency domain.
  • Fixed-point overflow. Each butterfly stage can double the magnitude. For Q15/Q31 types, scale the input by $1/N$ before the transform or apply per-stage scaling.
  • Spectral leakage. A finite-length DFT implicitly applies a rectangular window. Use a window function before the transform to reduce sidelobes.
  • Aliasing. The Nyquist-Shannon theorem requires $f_s \geq 2 f_{\max}$. Violating this folds high frequencies onto low frequencies irreversibly.
  • Frequency resolution is $\Delta f = f_s / N$. Increasing $N$ improves resolution at the cost of more memory and computation.

Variants & Generalizations

Variant Key Difference
Mixed-radix FFT Supports arbitrary composite $N$ using radix-2, -3, -5, etc. stages
Split-radix FFT Reduces multiplications by ~33 % compared to standard radix-2
Real-valued FFT Exploits conjugate symmetry of real inputs to halve work and memory
Bluestein's algorithm Handles arbitrary $N$ by converting DFT to a circular convolution
Inverse FFT Same structure; negate twiddle exponents and scale by $1/N$

Applications

  • Spectral analysis — Identifying frequency components, harmonics, and noise floors in sensor data.
  • Convolution and filtering — Multiplication in the frequency domain is equivalent to time-domain convolution; enables efficient FIR filtering.
  • Power Spectral Density — Used internally by Welch's method for spectral estimation.
  • Communication systems — OFDM modulation/demodulation is built directly on the FFT/IFFT pair.
  • Vibration and fault diagnosis — Examining spectral peaks to detect mechanical anomalies.
  • Data compression — Foundation for transforms (DCT, MDCT) used in JPEG, MP3, etc.

Connections to Other Algorithms

graph LR
    FFT --> PSD["Power Spectral Density"]
    FFT --> DCT["Discrete Cosine Transform"]
    WIN["Window Functions"] --> FFT
    FFT -.-> DK["Durand-Kerner (polynomial eval)"]
Loading
Algorithm Relationship
Power Spectral Density Calls FFT internally on each windowed segment
Discrete Cosine Transform Computed via FFT with input rearrangement and twiddle-factor post-processing
Window Functions Applied before FFT to reduce spectral leakage

References & Further Reading

  • Cooley, J.W. and Tukey, J.W., "An algorithm for the machine calculation of complex Fourier series", Mathematics of Computation, 19(90), 1965.
  • Oppenheim, A.V. and Schafer, R.W., Discrete-Time Signal Processing, 3rd ed., Pearson, 2009 — Chapters 8–9.
  • Smith, S.W., The Scientist and Engineer's Guide to Digital Signal Processing, California Technical Publishing, 1997 — Chapters 8, 12.