What the Discrete Fourier Transform Computes

Every introduction to the discrete Fourier transform presents the same sum: multiply each of $N$ inputs by a complex exponential of a chosen frequency, add them up, and repeat for $N$ frequencies. The formula is short, which makes it easy to treat as the whole definition. What the sum does not show is why the collection of outputs is invertible, why a periodic input needs no special treatment, and why the direct evaluation costs $N^2$ while practical implementations cost barely $N\log N$. All three follow from one structural fact: the transform is a change of basis in a finite-dimensional complex vector space, and the basis vectors are eigenvectors of a shift.

A fixed basis, not an adaptive analysis

The complex exponentials $e^{j2\pi kn/N}$ for $k = 0 \dots N-1$ form a set of $N$ vectors, each of length $N$. Their inner product with a signal vector $x$ is the DFT coefficient at bin $k$. Because the vectors are mutually orthogonal under the Hermitian inner product, the set is a basis, and each coefficient is the exact projection of $x$ onto one direction.

That framing settles a question that the summation formula leaves open: the transform does not search for structure in the data. It measures the data against $N$ predetermined directions. If the signal is a single exponential at a bin frequency, one coefficient is nonzero and the rest are zero, to numerical precision. If the signal is a sinusoid at a frequency between two bins, every coefficient is nonzero, and the energy spreads across the spectrum in a pattern determined by the signal length. Spectral analysis with a fixed-length transform is therefore always a measurement against a fixed grid, which is why the choice of record length is part of the measurement rather than an implementation detail. That effect is examined in Leakage, Scalloping, and Segment Averaging.

The basis vectors also explain the transform’s symmetry. Exponentials at bin $k$ and at bin $N-k$ are complex conjugates, so for a real input the coefficients at those two bins are conjugates: $X[N-k] = \overline{X[k]}$. Half the output of a real transform is redundant. Libraries that offer a real-input transform exploit this to roughly halve both the arithmetic and the storage, and a full complex transform on real data wastes about half its work.

The transform is a matrix, and the matrix is orthogonal

Writing the transform as a matrix product, $X = W x$ with $W_{kn} = e^{-j2\pi kn/N}$, makes several properties immediate. $W$ is symmetric, $W^T = W$. Its inverse is its conjugate transpose scaled by $1/N$. And its columns are orthogonal with squared norm $N$, so $W^H W = N I$.

Three consequences matter for implementation. First, the inverse transform is the forward transform with the sign of the exponent flipped and the output scaled by $1/N$; a library needs one kernel, not two. Second, the transform preserves energy up to that factor: the sum of squared magnitudes of the outputs equals $N$ times the sum of squared magnitudes of the inputs. Parseval’s relation in this discrete form is what makes power spectral density estimates computable from transform outputs, and it is why normalizing by $N$ is a decision about units rather than a matter of taste. Third, the transform is unitary when scaled by $1/\sqrt{N}$, which means the forward transform does not amplify rounding error. A single evaluation is numerically stable: the error in each output grows only as the square root of the accumulate length, since the terms are added with random signs.

The scale factor placement is a persistent source of confusion when comparing outputs from different tools. Some libraries apply $1/N$ on the inverse, some on the forward transform, some split it as $1/\sqrt{N}$ on both. The magnitude of the mismatch is exactly $N$, and for $N = 1024$ that is a factor of 1024, or about 60 dB. Comparisons between two implementations that disagree by a clean power of two counts of the sample number are almost always a normalization convention rather than a bug.

Why the direct sum is not used

Direct evaluation of $N$ outputs, each a sum of $N$ products, is $N^2$ complex multiply-adds. At $N = 10^6$ that is $10^{12}$ operations, and the transform would be unusable for any interactive purpose. The Cooley-Tukey decomposition observes that the sum can be split by the parity of the input index:

$$X[k] = \sum_{n \text{ even}} x[n] e^{-j2\pi kn/N} + \sum_{n \text{ odd}} x[n] e^{-j2\pi kn/N}$$

The even-indexed terms depend on $n = 2m$ and the exponential becomes $e^{-j2\pi km/(N/2)}$, a length-$N/2$ transform. The odd-indexed terms similarly become a length-$N/2$ transform, multiplied by the twiddle factor $e^{-j2\pi k/N}$. So one length-$N$ transform becomes two length-$N/2$ transforms plus $N/2$ complex multiplications. Recursing gives $\log_2 N$ levels, each of which performs $N/2$ twiddle multiplications, for a total of $(N/2)\log_2 N$ complex multiplications plus $N\log_2 N$ complex additions. In terms of real operations, a complex multiply is four multiplies and two adds, so a radix-2 complex transform on $N = 2^{20}$ costs on the order of $4 \times 10^7$ multiplies instead of $4 \times 10^{12}$ for the direct sum.

import numpy as np

def dft_matrix(x):
    n = x.size
    k = np.arange(n)
    w = np.exp(-2j * np.pi * np.outer(k, k) / n)   # N x N, built explicitly
    return w @ x

x = np.random.default_rng(0).standard_normal(64)
assert np.allclose(dft_matrix(x), np.fft.fft(x), atol=1e-10)

The matrix form is the definition; the recursive split is how the definition is evaluated. The two produce the same numbers, and the loop above is a legitimate correctness oracle for a hand-written transform.

The recursion is a data movement problem

The textbook recursion allocates new arrays for each sub-transform and ends with $O(N\log N)$ total allocation, which dominates the arithmetic on modern hardware. In-place implementations instead reorder the input once, so that the recursion’s base case reads contiguous pairs, and then perform the butterflies level by level over the same buffer. The reordering is a bit-reversal permutation: index $n$ moves to the position whose $\log_2 N$ bit pattern is $n$ reversed. It is a permutation of memory, not an arithmetic step, and combined with the twiddle table it accounts for most of the gap between measured throughput and a pure arithmetic count.

The other reason measured cost differs from the operation count is memory traffic. A length-$2^{20}$ complex transform on double precision touches sixteen megabytes, well beyond any cache, so each butterfly level streams the buffer from main memory. Implementations therefore split the transform into a cache-resident inner transform on blocks and a pass over the blocks, and choose the block size experimentally. Libraries that measure at startup and select among several algorithm variants are choosing a memory-traffic schedule, not a better operation count. The asymptotic exponent is fixed at $N \log N$; the constant is a property of the machine.

For real-valued data, the conjugate symmetry allows a length-$N$ real transform to be computed with a length-$N/2$ complex transform by packing even and odd samples into real and imaginary parts, then separating the two half-spectra afterwards. The unpacking is a fixed set of butterflies, so no extra asymptotic work is added.

What the transform does not determine

The DFT assumes the input is one period of a periodic signal. It extends the record by repetition implicitly, so a record that begins and ends at different values is interpreted as a signal with a discontinuity at every boundary. That discontinuity has broadband content, and the spectrum shows it as energy smeared across all bins. This is a property of the model, not an artifact to be removed by a better algorithm; the remedy is to make the record’s endpoints agree, by windowing or by choosing the record length to contain a whole number of cycles of the components of interest.

The transform also has no notion of time. Each output depends on every input, and the phase of each coefficient encodes the alignment of that frequency component with the record start, which is a reference instant chosen by the analyst. Obtaining a description that has both frequency and time locality requires a different construction, either a short-time transform with overlapping windows or a wavelet basis, and the cost of that locality is a fixed time-frequency trade governed by the uncertainty relation.

Finally, the transform of a finite record gives $N$ coefficients for $N$ samples, so no information is gained in a counting sense. What changes is the representation: structure that is spread over many samples in time becomes concentrated in a few coefficients in frequency. The practical value of the transform is that this concentration is what makes filtering by multiplication, compression by coefficient thresholding, and detection by peak search possible. Filters that operate in the coefficient domain are the subject of FIR Design by Windowing and the Gibbs Ripple.