Bluestein's FFT algorithm

Remarks

The key identity is $nk = \frac12 \bigl(n^2 + k^2 - (k-n)^2\bigr)$, which gives $e^{-2\pi i nk/N} = e^{-\pi i n^2/N} e^{-\pi i k^2/N} e^{\pi i (k-n)^2/N}$. So the DFT is a pointwise product with a chirp, a convolution with the conjugate chirp, and another pointwise product. The index $k - n$ runs from $-(N-1)$ to $N-1$, and the zero padding to $M \ge 2N-1$ makes the cyclic convolution agree with this linear convolution on $0 \le k < N$. The cyclic convolution of length $M$ equals the inverse DFT of the product of the DFTs of $\alpha$ and $\beta$. Since one can choose $M < 4N$, the constant in $O(N \log N)$ does not depend on the arithmetic of $N$; in particular prime $N$ is allowed. The same identity with $e^{-2\pi i/N}$ replaced by any complex $z$ evaluates the chirp z-transform $\sum_n x_n z^{-nk}$, which is the subject of the source article.

Going deeper

For a student

The algorithm is the standard fallback in FFT libraries for lengths with large prime factors. Compared with Rader's algorithm, it does not need a primitive root and works for every $N$, not only primes. It needs about three FFTs of length $M \approx 2N$ to $4N$, so it is slower than a direct Cooley–Tukey FFT of a highly composite length. The chirp $e^{\pi i m^2/N}$ is the same object as in the fractional Fourier transform and in radar signal processing.