0.10 Implementing ffts in practice  (Page 3/21)

where $0\le k<n$ and ${\omega }_n=exp(-2\pi i/n)$ is a primitive root of unity. Implemented directly, [link] would require $\Theta \left(n^2\right)$ operations; fast Fourier transforms are $O(nlogn)$ algorithms to compute the same result. The most important FFT (and the one primarily used in FFTW) is known as the“Cooley-Tukey” algorithm, after the two authors who rediscovered and popularized it in 1965 [link] , although it had been previously known as early as 1805 by Gauss as well as by laterre-inventors [link] . The basic idea behind this FFT is that a DFT of a composite size $n=n_1{n}_2$ can be re-expressed in terms of smaller DFTs of sizes $n_1$ and $n_2$ —essentially, as a two-dimensional DFT of size $n_1\times n_2$ where the output is transposed . The choices of factorizations of $n$ , combined with the many different ways to implement the data re-orderings of thetranspositions, have led to numerous implementation strategies for the Cooley-Tukey FFT, with many variants distinguished by their ownnames [link] , [link] . FFTW implements a space of many such variants, as described in "Adaptive Composition of FFT Algorithms" , but here we derive the basic algorithm, identify its key features, and outline some important historical variations and their relation to FFTW.

The Cooley-Tukey algorithm can be derived as follows. If $n$ can be factored into $n=n_1{n}_2$ , [link] can be rewritten by letting $\ell ={\ell }_1{n}_2+{\ell }_2$ and $k=k_1+k_2{n}_1$ . We then have:

where $k_{1,2}=0,...,n_{1,2}-1$ . Thus, the algorithm computes $n_2$ DFTs of size $n_1$ (the inner sum), multiplies the result by the so-called [link] twiddle factors ${\omega }_n^{{\ell }_2{k}_1}$ , and finally computes $n_1$ DFTs of size $n_2$ (the outer sum). This decomposition is then continued recursively. The literature uses the term radix to describe an $n_1$ or $n_2$ that is bounded (often constant); the small DFT of the radix is traditionally called a butterfly .

Many well-known variations are distinguished by the radix alone. A decimation in time ( DIT ) algorithm uses $n_2$ as the radix, while a decimation in frequency ( DIF ) algorithm uses $n_1$ as the radix. If multiple radices are used, e.g. for $n$ composite but not a prime power, the algorithm is called mixed radix . A peculiar blending of radix 2 and 4 is called split radix , which was proposed to minimize the count of arithmeticoperations [link] , [link] , [link] , [link] , [link] although it has been superseded in this regard [link] , [link] . FFTW implements both DIT and DIF, is mixed-radix with radices that are adapted to the hardware, and often uses much larger radices (e.g. radix 32) than wereonce common. On the other end of the scale, a “radix” of roughly $\sqrt{n}$ has been called a four-step FFT algorithm (or six-step , depending on how many transposes one performs) [link] ; see "FFTs and the Memory Hierarchy" for some theoretical and practical discussion of this algorithm.

A key difficulty in implementing the Cooley-Tukey FFT is that the $n_1$ dimension corresponds to discontiguous inputs ${\ell }_1$ in $\mathbf{X}$ but contiguous outputs $k_1$ in $\mathbf{Y}$ , and vice-versa for $n_2$ . This is a matrix transpose for a single decomposition stage, and the compositionof all such transpositions is a (mixed-base) digit-reversal permutation (or bit-reversal , for radix 2). The resulting necessity of discontiguous memory access and data re-ordering hindersefficient use of hierarchical memory architectures (e.g., caches), so that the optimal execution order of an FFT for given hardware isnon-obvious, and various approaches have been proposed.