Fast Fourier Transform (FFT) lab report

1. Title and aim

Fast Fourier Transform (FFT) & Polynomial Multiplication: Cooley-Tukey Radix-2 Decimation, Butterfly Operations, and Frequency Domain Convolution. To investigate the algorithmic principles of the Discrete Fourier Transform (DFT), analyze the Cooley-Tukey Radix-2 Decimation-in-Time (DIT) algorithm, explore in-place bit-reversal permutations, evaluate complex twiddle factor symmetries and recurrence relations, verify the exact butterfly scaling law B(N) = k_fft · (N log₂ N) + B_0 = 0.5000 · (N log₂ N) + 0.0000 butterflies across transform sizes N ∈ [64, 4096], and apply the Convolution Theorem to achieve O(N log N) asymptotic polynomial multiplication.

2. Abstract

Seven benchmark trials were generated across power-of-two transform sizes N = 64 to N = 4096 points (total radix-2 butterfly operations B from 192 to 24576 butterflies). An ordinary least-squares regression of butterfly operations against the composite complexity coordinate x = N log₂(N) yielded an empirical multiplier of 0.5000 butterflies / (N log₂ N) and an intercept of 0.0000 butterflies, demonstrating exact agreement with Cooley and Tukey's decimation-in-time invariant B(N) = 0.5 N log₂(N) (difference = 0.0000%). The empirical telemetry confirms that by recursively bisecting an N-point sequence into even and odd sub-transforms and combining them across log₂(N) stages, each stage executes exactly N/2 dual-output butterfly arithmetic operations.

3. Apparatus and fixed settings

• Cooley-Tukey Radix-2 Decimation-in-Time Engine: Executes recursive and iterative power-of-two discrete Fourier transform decompositions by partitioning input sequences into even and odd indices.
• Bit-Reversal Permutation Permuter: Performs in-place index permutation by reversing the binary bit representation of array indices, enabling non-recursive bottom-up butterfly execution.
• Trigonometric Twiddle Factor Look-Up Table: Precomputes and caches complex unit circle roots of unity W_N^k = exp(-2πi k / N) = cos(2πk/N) - i sin(2πk/N) to eliminate runtime transcendental math overhead.
• Dual-Path Butterfly Arithmetic Kernel Unit: Executes the core butterfly computation: u = A[k], t = W · A[k+m], A[k] = u + t, A[k+m] = u - t using fused multiply-add (FMA) complex arithmetic.
• Convolution & Polynomial Multiplier Subsystem: Implements the Convolution Theorem to perform O(N log N) polynomial coefficient multiplication via forward FFT, pointwise spectral product, and inverse FFT.
• Microarchitectural Profiler & Cache Monitor: Tracks processor memory hierarchy operations, L1/L2 data cache misses across butterfly stride transitions, and verifies exact regression against the theoretical bound 0.5 N log₂ N.
• Power-of-two transform sizes N ∈ [64, 128, 256, 512, 1024, 2048, 4096] (decomposition depths log₂(N) ∈ [6..12]).
• Complex numbers represented in IEEE 754 double precision (8-byte real + 8-byte imaginary, 16 bytes/sample).
• In-place bit-reversal permutation implemented via branchless bit-twiddling lookup.
• Precomputed trigonometric twiddle factor table W_N^k = cos(2πk/N) - i·sin(2πk/N) aligned to 64-byte CPU cache lines.
• Fused multiply-add (FMA) complex butterfly kernel: u = A[k], t = W · A[k+m], A[k] = u + t, A[k+m] = u - t.

4. Theory and graph choice

The Discrete Fourier Transform (DFT) is one of the most fundamental algorithmic transformations in computational science, engineering, and applied mathematics. It converts a discrete finite time-domain or spatial-domain sequence of N complex numbers into its orthogonal frequency-domain representation. Direct evaluation of the DFT via naive matrix-vector multiplication with the dense N×N Fourier Vandermonde matrix requires Θ(N²) complex arithmetic operations. For high-throughput digital signal processing, telecommunications (OFDM/5G), computer vision, audio analysis, and high-precision scientific computing, quadratic scaling becomes completely prohibitive: an N = 65,536 sample buffer demands over 4.29 billion complex multiplications. In 1965, James W. Cooley and John W. Tukey revolutionized computational mathematics by publishing the Fast Fourier Transform (FFT), reducing the computational complexity from Θ(N²) to strictly Θ(N log₂ N). For N = 65,536, Cooley-Tukey evaluates the transform using only ~524,000 operations—a computational acceleration exceeding 8,000×! The radix-2 Decimation-in-Time (DIT) formulation achieves this breakthrough by recursively partitioning an N-point sequence into two (N/2)-point sequences (one containing even-indexed samples and the other odd-indexed samples). Exploiting the rotational symmetries and half-period anti-symmetry of the complex twiddle factors (W_N^(k + N/2) = -W_N^k where W_N = exp(-2πi/N)), the sub-transform outputs are recombined across log₂(N) stages using dual-output butterfly arithmetic kernels. Each butterfly executes 1 complex multiplication and 2 complex additions. Across the entire transform, each of the log₂(N) stages executes exactly N/2 butterflies, establishing the exact theoretical invariant B(N) = 0.5000 · N log₂(N) butterfly operations. In-place evaluation is enabled by the Bit-Reversal Permutation, which reorders the input array into recursive leaf order, allowing non-recursive bottom-up execution with O(1) auxiliary memory. Furthermore, under the Convolution Theorem, cyclic convolution in the time domain is isomorphic to element-wise pointwise multiplication in the frequency domain. This allows multiplying two degree-d polynomials in strictly O(N log N) time, dropping from billions of operations in classical Cauchy multiplication to mere milliseconds. In this laboratory investigation, students benchmark butterfly execution across transform sizes N ∈ [64, 4096], verify the linear scaling law B(N) = 0.5000 · (N log₂ N) via ordinary least-squares regression, explore bit-reversal permutation tables, analyze complex twiddle factor errors, and implement frequency-domain polynomial convolution.

The mathematical regression fits the linear model y = mx + b, where x = N · log₂(N) represents the composite algorithmic complexity coordinate in work units, y = B represents total radix-2 butterfly operations, m represents the butterfly scaling multiplier k_fft, and b represents the baseline arithmetic offset (theoretically B_0 = 0.0000 butterflies). The gradient has physical units butterflies / (N log₂ N).

5. Method

1. Initialize the discrete input sequence or polynomial coefficient vector across benchmark sizes N ∈ [64, 128, 256, 512, 1024, 2048, 4096]. If the original data length is not a power of two, zero-pad trailing elements to the next integral power of two N = 2^M to ensure complete binary tree decomposition.
2. Apply the in-place bit-reversal permutation algorithm. For each index i from 0 to N-1, calculate its M-bit reversed index rev_i. If i < rev_i, swap data elements at indices i and rev_i, placing the array into the exact leaf order required for iterative bottom-up butterfly evaluations.
3. Generate the complex twiddle factor lookup array for the current transform size N. Compute W_N^k = cos(2πk/N) - i·sin(2πk/N) for each sub-problem branch, verifying half-period anti-symmetry W_N^(k + N/2) = -W_N^k and quarter-period quadrature relations.
4. Execute the iterative Cooley-Tukey butterfly stages for s = 1, 2, ..., log₂(N). At each stage s, iterate across blocks of length 2^s with half-length step m = 2^(s-1). For each pair of elements within the block, apply the butterfly kernel and increment the verified butterfly operation counter by exactly 1.
5. For convolution and polynomial multiplication benchmarks, evaluate forward FFTs on both input operands, compute the complex pointwise vector product C_freq = A_freq ⊙ B_freq in linear O(N) operations, and execute the inverse FFT using complex conjugated twiddle factors normalized by 1/N to extract exact product coefficients.
6. Plot the observed butterfly operations B(N) against the composite coordinate x = N·log₂(N). Compute ordinary least-squares regression to verify the theoretical multiplier k_fft = 0.5000 butterflies/(N log₂ N), intercept B_0 = 0.0000 butterflies, and correlation coefficient R² = 1.00000.

6. Observations and calculation table

All benchmark observations use strictly controlled memory buffers. Values are recorded deterministically from the algorithmic comparator register. Residuals are computed relative to the ordinary least-squares regression line.

Download the sample observations (CSV) ↓

7. Worked calculation for one observation

8. Graph, gradient and final calculation

Circular markers represent empirical benchmark observations. The solid line represents the ordinary least-squares asymptotic linear fit.

9. Uncertainty and comparison

Percentage difference = |estimate − theoretical value| / |theoretical value| × 100 = 0.0000%.

For n = 7 benchmark observations, Σresidual² = 0.00000; s² = Σresidual²/(n − 2). Gradient standard error = √[s²/Σ(x − x̄)²] = 2.2934e-11. Intercept standard error = √[s²(1/n + x̄²/Σ(x − x̄)²)] = 4.7909e-07.

Propagating the fit parameters yields a fit standard error of 2.2934e-11 butterflies / (N log₂ N), with an ideal coefficient of determination R² = 1.000000.

In discrete Fourier transform analysis and computational algebra, the butterfly operation count B(N) = 0.5000 · N log₂(N) + 0.0000 butterflies is a strict deterministic mathematical invariant for Radix-2 Cooley-Tukey decompositions. Unlike physical sensor systems, there are no continuous measurement noise, thermal perturbations, or instrument calibration drifts in counting discrete arithmetic butterfly executions. However, in high-performance computing hardware and DSP microprocessors, wall-clock execution time and numerical accuracy are modulated by computer microarchitecture and floating-point dynamics: (1) Floating-point round-off error propagation: in IEEE 754 single and double precision floating-point arithmetic, evaluating complex twiddle factor recurrences across log₂(N) stages accumulates rounding errors with an RMS relative norm bounded by O(log₂(N) · ε_machine), where machine epsilon ε ≈ 1.11 × 10⁻¹⁶ for 64-bit floats; (2) CPU L1/L2/L3 data cache misses and memory stride thrashing: as transform size N increases beyond the L1 data cache capacity (typically 32–64 KB), the power-of-two memory strides (strides of 2^(s-1)) in large butterfly stages cause severe cache line conflicts and TLB stalls unless transpose-based Four-Step or cache-oblivious algorithms are employed; (3) Twiddle factor look-up table latency vs. runtime trigonometric evaluation: dynamically calling transcendental sin() and cos() functions in inner loops introduces 50–100 clock cycles of latency per butterfly, whereas precomputed look-up tables eliminate trigonometric calls at the expense of memory cache occupancy; and (4) SIMD register width and vectorization efficiency: modern AVX2 and AVX-512 vector units process 4 to 8 double-precision complex multiply-add (FMA) instructions simultaneously, accelerating butterflies by up to 6× over scalar implementations. Ordinary least-squares linear regression of observed butterfly executions against the composite complexity coordinate x = N log₂(N) across sample sizes from N = 64 to N = 4096 yields an empirical multiplier k_fft = 0.5000 butterflies/(N log₂ N) and baseline intercept B_0 = 0.0000 butterflies with zero residual variance (R² = 1.00000), verifying the Cooley-Tukey asymptotic complexity bound Θ(N log N) with absolute mathematical certitude.

10. Result and discussion

The empirical benchmark observations provide rigorous quantitative verification of Fourier Transform theory, discrete orthogonal projections, and divide-and-conquer asymptotic complexity. The Discrete Fourier Transform (DFT) converts a discrete complex time-domain sequence of length N into its orthogonal frequency-domain spectrum. Direct matrix-vector evaluation against the dense N×N Vandermonde Fourier matrix scales quadratically as Θ(N²), requiring N² complex multiplications and N(N-1) complex additions. For real-time audio processing, telecommunications (OFDM/5G), and scientific simulation, quadratic scaling is computationally prohibitive: an N = 65,536 sample buffer demands over 4.29 billion multiplications. In 1965, James W. Cooley and John W. Tukey published their Radix-2 Decimation-in-Time (DIT) algorithm, reducing complexity from Θ(N²) to strictly Θ(N log₂ N). For N = 65,536, Cooley-Tukey evaluates the transform using only ~524,000 operations—a computational acceleration exceeding 8,000×! The algorithm recursively partitions an N-point sequence into two (N/2)-point sequences (one containing even-indexed samples and the other odd-indexed samples). Exploiting the rotational periodicity (W_N^(k+N) = W_N^k) and half-period anti-symmetry (W_N^(k+N/2) = -W_N^k) of complex twiddle factors (W_N = exp(-2πi/N)), the sub-transform outputs are recombined across log₂(N) stages using dual-output 'butterfly' arithmetic kernels. Each butterfly executes 1 complex multiplication and 2 complex additions. Across the entire transform, each of the log₂(N) stages executes exactly N/2 butterflies, establishing the exact theoretical invariant B(N) = 0.5000 · N log₂(N) butterfly operations. Across all benchmark trials from N = 64 to N = 4096 points, the empirical regression slope m = 0.5000 butterflies/(N log₂ N) and intercept b = 0.0000 confirm this invariant with zero residual variance (R² = 1.00000). The investigation also verified the profound utility of the Bit-Reversal Permutation, which reorders the input array into recursive leaf order, allowing non-recursive bottom-up execution with O(1) auxiliary memory. Finally, application benchmarks demonstrated the power of the Convolution Theorem: cyclic convolution in the time domain corresponds to element-wise pointwise multiplication in the frequency domain. This reduces polynomial multiplication C(x) = A(x)B(x) from O(N²) in classical Cauchy multiplication down to strictly O(N log N) via forward FFT, pointwise product, and inverse FFT, forming the algorithmic foundation of modern cryptography, computer algebra systems, and large-integer arithmetic.

11. Precautions and limitations

• Always zero-pad input sequences to the nearest power of two N ≥ 2d + 1 when performing polynomial multiplication to prevent circular aliasing from corrupting high-degree terms.
• When precomputing twiddle factors, compute angles using symmetric quarter-period ranges to minimize floating-point round-off error accumulation near π/2 and π.
• In multithreaded FFT implementations, avoid false sharing by ensuring thread partition boundaries align with CPU cache line multiples (64 bytes).
• Verify that the inverse FFT applies the exact 1/N scaling factor and conjugates twiddle factor imaginary signs to maintain energy conservation under Parseval's theorem.

The simulated benchmark environment models idealized double-precision complex arithmetic where radix-2 butterfly operations execute deterministically (B(N) = 0.5 N log₂ N). In physical computing systems and DSP microprocessors, execution throughput is heavily modulated by CPU L1/L2 data cache line capacity and power-of-two memory stride conflicts (cache thrashing at large strides), SIMD vector register availability (AVX-512 vs AVX2 FMA throughput), floating-point round-off accumulation across deep butterfly networks, and operating system scheduling jitter. In the Cooley-Tukey Radix-2 Decimation-in-Time FFT benchmark across transform sizes N in [64, 4096] (where N = 2^M for integer M in [6, 12]), the transform decomposes recursively into M = log₂(N) butterfly stages. At each stage s in [1..log₂(N)], the data vector is divided into contiguous blocks of length 2^s with half-span m = 2^(s-1). Exactly N/2 dual-output butterfly arithmetic operations are executed per stage to combine the lower and upper sub-transform branches with complex twiddle factors W_{2^s}^j = exp(-2πi j / 2^s). Across all log₂(N) stages, the cumulative butterfly operation count is governed by the exact deterministic invariant B(N) = (N/2) · log₂(N) = 0.5000 · (N log₂ N) + 0.0000 butterflies. When plotted against the composite complexity coordinate x = N log₂(N), ordinary least-squares regression yields an empirical slope k_fft = 0.5000 butterflies/(N log₂ N) and baseline intercept B_0 = 0.0000 butterflies with an ideal coefficient of determination R² = 1.000000.

12. Conclusion

This laboratory investigation quantitatively verified the decimation-in-time decomposition, butterfly network mechanics, and asymptotic scaling of the Cooley-Tukey Fast Fourier Transform. Through systematic benchmarking across transform sizes from N = 64 to 4096 points, we verified the empirical butterfly multiplier k_fft = 0.5000 butterflies / (N log₂ N) and baseline intercept B_0 = 0.0000 butterflies with an ideal R² = 1.00000. The investigation established the microarchitectural efficiency and theoretical indispensability of bit-reversal permutations, twiddle factor symmetries, and frequency-domain convolution, proving why the Fast Fourier Transform is universally recognized as one of the most critical algorithmic discoveries of the twentieth century.