Discrete Fourier transform (source code)

= Discrete Fourier transform
{title2=$X_k=\sum_{n=0}^{N-1}x_ne^{-2\pi ink/N}$}
{wiki}

= DFT
{c}
{synonym}

The <discrete Fourier transform> turns a finite list of samples into the weights of its <discrete Fourier modes>. The goal is to reveal a simple oscillatory recipe hidden in complicated-looking data, while retaining enough information to reconstruct every sample exactly. In <signal processing>, those weights identify sampled <frequencies> and <wave phases>; in <numerical analysis>, they simplify <circular convolution> and the <discrete sine transform Poisson solver>.

Start with eight real samples, indexed by $n=0,\ldots,7$:
$$
x_n=1+\cos\left(\frac{2\pi n}{8}\right)+\frac12\sin\left(\frac{4\pi n}{8}\right).
$$
Writing $a=\sqrt2/2$, the actual list is
$$
x=\left(2,\frac32+a,1,\frac12-a,0,\frac32-a,1,\frac12+a\right).
$$
The <discrete Fourier transform> recovers the recipe as
$$
X=(8,4,-2i,0,0,0,2i,4).
$$
The constant contributes $X_0=8$; the <cosine> contributes $X_1=X_7=4$; the half-sized <sine> contributes $X_2=-2i$ and $X_6=2i$. All other <coefficients> vanish exactly. The imaginary <coefficients> carry the <sine>'s <wave phase>, which a plot of <absolute values> alone would hide.

\Image[/dft-frequency-decomposition.png]
{title=Eight samples become a constant, a cosine pair, and a sine pair in the discrete Fourier transform}

Precisely, for a positive <integer> $N$ and a <vector> $x\in\mathbb C^N$, the forward <discrete Fourier transform> and <inverse discrete Fourier transform> used here are
$$
X_k=\sum_{n=0}^{N-1}x_n\exp\left(-\frac{2\pi ink}{N}\right),\qquad k=0,\ldots,N-1,
$$
$$
x_n=\frac1N\sum_{k=0}^{N-1}X_k\exp\left(\frac{2\pi ink}{N}\right),\qquad n=0,\ldots,N-1.
$$
Thus the forward exponent has a minus sign and no scale factor; the inverse exponent has a plus sign and the factor $1/N$. The <coefficients> in the reconstruction are $X_k/N$, rather than $X_k$. In the example, $X_0/8=1$, the pair $4/8$ reconstructs the <cosine>, and the pair $-2i/8,2i/8$ reconstructs the half-sized <sine>. These conventions agree with the scaled <Fourier coefficients> in MIT's authored notes.https://ocw.mit.edu/courses/6-02-introduction-to-eecs-ii-digital-communication-systems-fall-2012/06d7972059805d20e6a7db3a52ec23f2_MIT6_02F12_chap13.pdf[Hari Balakrishnan, Christopher Terman, and George Verghese, MIT 6.02 notes, Chapter 13, equations (13.15)–(13.16), 2012]{ref}

The proof idea is <orthogonality of roots of unity>: for <integers> $k,\ell\in\{0,\ldots,N-1\}$,
$$
\sum_{n=0}^{N-1}e^{2\pi in(k-\ell)/N}=\begin{cases}N,&k=\ell,\\0,&k\ne\ell.\end{cases}
$$
Distinct <discrete Fourier modes> cancel in this finite <geometric series>. Substituting the forward formula into the <inverse discrete Fourier transform> leaves only the desired sample. The same <orthogonality> makes the normalized <discrete Fourier modes> an <orthonormal basis> and gives the <discrete Parseval identity>, in the present convention,
$$
\sum_{n=0}^{N-1}|x_n|^2=\frac1N\sum_{k=0}^{N-1}|X_k|^2.
$$
For the eight-sample example, direct squaring gives $\sum_nx_n^2=13$, while $\sum_k|X_k|^2=64+16+4+4+16=104$, so $104/8=13$. This identity preserves the squared <L2 norm> after accounting for normalization. The existing <discrete Parseval identity> instead scales the forward <discrete Fourier transform> by $N^{-1/2}$, producing a <unitary matrix> and equal squared <L2 norms> without the factor $1/N$.

For real input, <discrete Fourier conjugate symmetry> pairs positive and negative <frequencies>. The label $k$ is a residue modulo $N$: $k=7$ when $N=8$ represents the same <discrete Fourier mode> as $k=-1$, since $e^{2\pi i7n/8}=e^{-2\pi in/8}$. It therefore represents negative <frequency> $-1/8$ cycles per sample. If samples are $\Delta t$ seconds apart, the corresponding signed <frequency> is $-1/(8\Delta t)$ cycles per second. Similarly, index $6$ represents $-2$, while index $4$ is its own negative modulo $8$.

The <discrete Fourier transform> is an invertible map between two finite <vectors>. The continuous <Fourier transform> instead describes a <function> on the real line by an <integral> over a continuous <frequency> variable. A finite list of samples does not determine the full continuous <frequency spectrum>: <aliasing of discrete Fourier modes> already supplies distinct waveforms with identical observations. Applying the inverse formula at arbitrary <integer> indices gives the <periodic extension of a finite sequence>, repeating the recorded block every $N$ samples. This is a property of that representation, rather than evidence that the physical <signal> repeats beyond the recording. Laurent Demanet's MIT notes place the <discrete Fourier transform> alongside the continuous <Fourier transform> and <Fourier series> to explain the role of sampling and repetition.https://math.mit.edu/icg/resources/teaching/18.330/18.330-Mar21.pdf[Laurent Demanet, Introduction to Numerical Analysis, MIT 18.330, §§6.2–6.3, pp. 57–60, March 22, 2012]{ref}

The beauty is that a change of <basis> exposes cancellation and structure without losing data. Direct evaluation needs $N$ sums of $N$ terms; a <Fast Fourier transform> computes the same <discrete Fourier transform> with far less work. Computational speed and the mathematical meaning of the <coefficients> are separate questions.