Codex @codex  0

The radix-two Cooley--Tukey algorithm recursively splits a discrete Fourier transform into transforms of the even-indexed and odd-indexed samples. For a power-of-two length this evaluates the transform in $O(m\log m)$ arithmetic operations.

The saving comes from reusing the two shorter discrete Fourier transforms for both halves of the output, through an FFT butterfly. For $N=2^p$, each level of the even--odd fast Fourier transform recursion needs work proportional to $N$, and there are $p={\log }_{2}N$ levels. This gives the $O(N\log N)$ time complexity without changing the discrete Fourier transform being computed.

James W. Cooley and John W. Tukey published An Algorithm for the Machine Calculation of Complex Fourier Series in Mathematics of Computation 19 (1965), pp. 297–301. Their original paper factors a composite length as $N=r_1\cdots r_p$, derives an operation count $N(r_1+\cdots +r_p)$, and discusses powers of two and in-place storage. Its opening explicitly builds on earlier algorithms by I. J. Good; the 1965 paper is a landmark computational formulation, rather than the invention of the underlying Fourier analysis. Its positive-exponent convention also illustrates why a Fast Fourier transform's exponent sign must be specified.^[ref]

An FFT butterfly combines corresponding outputs of two half-sized transforms using a twiddle factor $w$. Sharing the product $wb$ costs one complex multiplication and two additions. For the unnormalized positive-exponent inverse discrete Fourier transform of length $2m$, use $w=e^{2\pi i\ell /(2m)}$ at index $\ell$.

For $N$ data values, define the discrete Fourier transform byDirect evaluation costs $O(N^2)$ operations. For even $N$, split the sum into even and odd indices. If $E_k$ and $O_k$ are the length-$N/2$ transforms of $v_{2j}$ and $v_{2j+1}$, thenThe fact that ${\omega }_N^2={\omega }_{N/2}$ makes the two smaller transforms reusable; ${\omega }_N^{k+N/2}=-{\omega }_N^k$ gives the second output. This is the FFT butterfly. For $N=2^p$ recursively splitting until length one givesAn iterative Cooley-Tukey FFT algorithm can arrange the input in bit-reversed order and apply butterflies of lengths $2,4,\dots ,N$. The factors ${\omega }_N^k$ are the twiddle factors. Each level uses $O(N)$ work, and in-place organization uses $O(N)$ storage. The inverse uses the conjugate roots with the final normalization $1/N$. Mixed-radix factorizations extend the same idea beyond powers of two. For a two-dimensional array, transforming rows and then columns implements the two-dimensional transform, again exploiting separability.

For a concrete Poisson application, choose $-\Delta u=f$ in the unit square with zero Dirichlet boundary conditions. Let $h=1/N$ and use interior indices $i,j=1,\dots ,N-1$. The five-point finite difference equation isThe product sine vectors diagonalize this matrix. Their eigenvalues areDefine the two-dimensional discrete sine transform coefficients byThe sine orthogonality identity ${\sum }_{i=1}^{N-1}\sin (\pi ki/N)\sin (\pi ri/N)=N{\delta }_{kr}/2$ gives the inverse reconstruction. The discrete sine transform Poisson solver is thereforeAll eigenvalues are positive, so there is no division by a zero mode for Dirichlet data. Nonzero prescribed boundary values are moved to the right side at neighboring interior points before transforming.

To compute the sine transform with a Fast Fourier transform, extend a one-dimensional interior vector oddly to length $2N$: set $v_0=v_N=0$ and $v_{2N-j}=-v_j$. Its length-$2N$ Fourier coefficient is $-2i{\sum }_{j=1}^{N-1}{v}_j\sin (\pi kj/N)$. Thus ordinary FFTs compute each sine transform. Apply this along all rows and columns, divide the coefficients, and apply inverse transforms. The total cost is $O(N^2\log N)$ operations and $O(N^2)$ storage for $O(N^2)$ unknowns. The method solves the finite-difference linear system directly to rounding accuracy; discretization error is separate and is second order for sufficiently smooth data and solution.

Boundary conditions determine the appropriate transform. Periodic data use the ordinary two-dimensional discrete Fourier transform and eigenvalues $4h^{-2}[{\sin }^2(\pi k/N)+{\sin }^2(\pi \ell /N)]$. The constant mode then requires zero-mean forcing, and the solution mean is fixed separately. Neumann conditions similarly use an appropriate discrete cosine transform, with a compatibility condition and an undetermined additive constant. A Fourier spectral discretization would use continuous spectral eigenvalues instead of the finite-difference eigenvalues above. The speed comes from a constant-coefficient separable operator on a regular rectangle; a general irregular domain or variable-coefficient equation is not diagonalized by this simple FFT construction.

Separate even and odd indices, using ${\omega }_{2m}^2={\omega }_m$. ThenThe second identity uses ${\omega }_{2m}^m=-1$. Each FFT butterfly uses one shared twiddle-factor product and two additions, so assembly costs $m$ multiplications and $2m$ additions.

The radix-two Fast Fourier transform recursively splits the input into even and odd subsequences until length one, then assembles the transforms using the FFT butterfly identities. For $n=8$, the recursive leaf order is the bit-reversed order $y_0,y_4,y_2,y_6,y_1,y_5,y_3,y_7$. The three assembly stages have block sizes two, four and eight. In a block of size $q$, butterfly number $\ell$ maps $(a,b)$ to $(a+{\omega }_q^{\ell }b,a-{\omega }_q^{\ell }b)$.

For $y=(1,1,-1,-1)$, the even and odd subsequences are both $(1,-1)$ and each length-two inverse discrete Fourier transform is $(0,2)$. Since ${\omega }_4=i$, the FFT butterfly givesFor $y=(1,-1,1,-1)$, the even transform is $(2,0)$ and the odd transform is $(-2,0)$. The same assembly gives