Fast Fourier Transform

A fast Fourier transform (FFT) is an algorithm for efficiently computing the discrete Fourier transform (DFT) or its inverse. Rather than defining a different mathematical transform, an FFT reorganizes the DFT calculation to reuse intermediate results. Standard FFT methods reduce the arithmetic cost of transforming $N$ samples from $O(N^2)$ to $O(N\log N)$, making Fourier analysis practical for large datasets. FFTs are fundamental tools in signal processing and scientific computing. (fftw.org)

Mathematical definition

For a sequence $x_0,x_1,\ldots,x_{N-1}$ of complex numbers, one common convention defines the DFT as

where $i^2=-1$. The inverse is

An FFT evaluates these finite sums more efficiently; in exact arithmetic, it produces the same result as direct evaluation. Sign and normalization conventions differ between implementations: some place the factor $1/N$ on the forward transform, while a unitary convention uses $1/\sqrt N$ in both directions. (numpy.org)

The DFT is a finite-dimensional counterpart of the Fourier transform. For uniformly sampled time-domain data, its outputs describe frequency components on a discrete grid. With sampling frequency $f_s$, adjacent frequency bins are separated by $f_s/N$. In the usual output ordering, the zero-frequency component comes first, followed by positive-frequency components and then components interpreted as negative frequencies. (numpy.org)

Direct evaluation computes $N$ sums containing $N$ terms each. Its quadratic computational complexity contrasts with the near-linear cost of an FFT. The notation $O(N\log N)$, expressed using big-O notation, describes asymptotic growth rather than a fixed execution time or speedup. (sigproc.mit.edu)

The Cooley–Tukey principle

The best-known FFT family is the Cooley–Tukey algorithm. It applies divide and conquer: a transform of composite length is decomposed into smaller transforms, whose outputs are combined using complex phase factors. (fftw.org)

Radix-2 decomposition

For even $N$, separate the input into its even- and odd-indexed samples:

Writing $W_N=e^{-2\pi i/N}$, the original transform becomes

Thus two length-$N/2$ transforms produce all $N$ outputs after only a linear amount of additional work. The paired addition and subtraction is called a butterfly, and $W_N^k$ is a twiddle factor. (sigproc.mit.edu)

When $N=2^m$, repeated recursive subdivision yields $\log_2 N$ stages. The corresponding recurrence relation is

which gives $T(N)=O(N\log N)$. The speedup comes from shared subcomputations and the structure of the complex exponentials, not from discarding terms or approximating the transform. (sigproc.mit.edu)

Other factorizations

Cooley–Tukey is not restricted to powers of two. A factorization $N=ab$ permits a decomposition into transforms of lengths $a$ and $b$, together with twiddle-factor multiplications and rearrangement of indices. Mixed-radix algorithms combine different factors; split-radix algorithms combine decompositions of different sizes, commonly one transform of length $N/2$ with two of length $N/4$. (fftw.org)

Algorithm families and transform lengths

Different FFT methods exploit different arithmetic structures:

• Prime-factor, or Good–Thomas, algorithms use relatively prime factors of $N$ to avoid the interstage twiddle factors required by a corresponding Cooley–Tukey decomposition.
• Rader’s algorithm converts a transform of prime length $p$ into a cyclic convolution of length $p-1$, apart from the zero-frequency term.
• Bluestein’s algorithm rewrites a DFT as a convolution using quadratic phase factors. It handles arbitrary lengths, including primes.
• Split-radix methods seek to reduce arithmetic counts by combining different recursive decompositions. (fftw.org)

Consequently, an FFT does not inherently require a power-of-two input length. Sizes with small prime factors are often convenient, but efficient implementations can also compute prime-length transforms in $O(N\log N)$ operations. Padding to a different length may improve execution time, but it changes the transform being evaluated. (fftw.org)

Real-valued and multidimensional transforms

For real-valued input, the output has conjugate symmetry:

with indices interpreted modulo $N$. Approximately half the frequency-domain values are therefore redundant. Real-input FFT algorithms exploit this symmetry to reduce computation and storage. The zero-frequency coefficient is real, as is the Nyquist-frequency coefficient when $N$ is even. (numpy.org)

A multidimensional DFT is separable: it can be computed by applying one-dimensional transforms along each array dimension. A two-dimensional transform, for example, can be evaluated by transforming every row and then every column. For a fixed-dimensional array containing $M$ elements, this approach gives an arithmetic cost of $O(M\log M)$. Multidimensional FFTs support applications in image processing and computational science. (fftw.org)

Applications

Spectral analysis

FFTs convert sampled signals into frequency-domain representations, allowing analysis of oscillations, amplitude, and phase. They are used in audio analysis and other forms of signal processing. Applying transforms to successive data segments provides a time-dependent view of spectral content rather than a single spectrum for an entire record. (numpy.org)

Fast convolution and filtering

The convolution theorem connects convolution in one domain with multiplication in the other. With the normalization above,

where multiplication is elementwise and $\circledast$ denotes circular convolution. This gives an efficient route to filtering long sequences. (mathworks.com)

Circular convolution is not automatically the same as linear convolution. For sequences of lengths $L$ and $M$, both must be padded to a common transform length of at least $L+M-1$ to recover the full linear convolution without wraparound. The same coefficient-convolution relationship underlies FFT-based multiplication of polynomials. (mathworks.com)

Numerical computation

Fourier-transform methods can also assist in solving partial differential equations. Transforming along spatial dimensions can simplify suitable differential operators, while the choice of transform must match the problem’s boundary conditions. Related fast cosine and sine transforms are useful for different symmetry and boundary requirements. (fftw.org)

Implementation and numerical accuracy

Arithmetic count alone does not determine performance. Memory access, cache behavior, data layout, and processor capabilities can outweigh modest differences in operation counts. Optimized implementations therefore combine small specialized transform kernels with different decompositions and execution strategies. (fftw.org)

Some libraries construct a plan before execution, selecting an implementation for a particular transform size, layout, and machine. The planning cost can be reused across repeated transforms. FFTW is a prominent example of this adaptive approach. (fftw.org)

In floating-point arithmetic, intermediate operations introduce rounding errors. Different FFT algorithms may consequently produce slightly different numerical outputs even though they represent the same exact transform. Well-implemented Cooley–Tukey methods have favorable numerical stability, but accurate twiddle factors are important: poorly chosen recurrences for generating them can accumulate substantial error. (fftw.org)

Interpretation and limitations

An FFT accelerates a calculation; it does not eliminate the limitations of the sampled data. In particular, zero-padding creates a denser sampling of the finite record’s spectrum but does not improve its intrinsic ability to resolve closely spaced frequency components. That ability depends on the original observation length and sampling rate. (mathworks.com)

Normalization and output ordering also matter. Raw coefficients, normalized amplitudes, and power spectra are different representations; a one-sided spectrum for real input requires appropriate treatment of the redundant negative-frequency components. Comparisons between implementations must account for their sign, scaling, and ordering conventions. (numpy.org)

Computing a full FFT is not always necessary when only a few frequency components are required. Specialized or pruned methods can calculate selected outputs, although their savings depend on how many outputs are needed and on implementation overhead. (fftw.org)

History

FFT-like decompositions predate electronic computers. Carl Friedrich Gauss described a method in 1805 that was later recognized as an early form of the Cooley–Tukey approach. This connection was established through subsequent historical research. (ftp.fftw.org)

The modern widespread adoption of the FFT followed James W. Cooley and John W. Tukey’s 1965 paper, An Algorithm for the Machine Calculation of Complex Fourier Series. Their presentation showed how factorizations of the transform length could greatly reduce computational work. Later research developed additional algorithms, real-data and multidimensional methods, and implementations optimized for particular computing architectures. (research.ibm.com)

References

1. FFTW Home Pagefftw.org
2. Discrete Fourier Transform (numpy.fft) — NumPy v2.2 Manualnumpy.org
3. FFTW FAQ - Section 3fftw.org
4. Fast Fourier Transform — 6.300: Signal Processingsigproc.mit.edu
5. The Design and Implementation of FFTW3fftw.org
6. Introduction (FFTW 3.3.11)fftw.org
7. Amplitude Estimation and Zero Paddingmathworks.com
8. Multi-dimensional Transforms (FFTW 3.3.11)fftw.org
9. Equalization, Convolution, and Cyclic Prefix Additionmathworks.com
10. More DFTs of Real Data (FFTW 3.3.11)fftw.org