Split-radix FFT algorithm

Recall that the DFT is defined by the formula:

${\displaystyle X_{k}=\sum _{n=0}^{N-1}x_{n}\omega _{N}^{nk}}$

where ${\displaystyle k}$ is an integer ranging from ${\displaystyle 0}$ to ${\displaystyle N-1}$ and ${\displaystyle \omega _{N}}$ denotes the primitive root of unity:

${\displaystyle \omega _{N}=e^{-{\frac {2\pi i}{N}}},}$

and thus: ${\displaystyle \omega _{N}^{N}=1}$.

The split-radix algorithm works by expressing this summation in terms of three smaller summations. (Here, we give the "decimation in time" version of the split-radix FFT; the dual decimation in frequency version is essentially just the reverse of these steps.)

First, a summation over the even indices ${\displaystyle x_{2n_{2}}}$. Second, a summation over the odd indices broken into two pieces: ${\displaystyle x_{4n_{4}+1}}$ and ${\displaystyle x_{4n_{4}+3}}$, according to whether the index is 1 or 3 modulo 4. Here, ${\displaystyle n_{m}}$ denotes an index that runs from 0 to ${\displaystyle N/m-1}$. The resulting summations look like:

${\displaystyle X_{k}=\sum _{n_{2}=0}^{N/2-1}x_{2n_{2}}\omega _{N/2}^{n_{2}k}+\omega _{N}^{k}\sum _{n_{4}=0}^{N/4-1}x_{4n_{4}+1}\omega _{N/4}^{n_{4}k}+\omega _{N}^{3k}\sum _{n_{4}=0}^{N/4-1}x_{4n_{4}+3}\omega _{N/4}^{n_{4}k}}$

where we have used the fact that ${\displaystyle \omega _{N}^{mnk}=\omega _{N/m}^{nk}}$. These three sums correspond to portions of radix-2 (size N/2) and radix-4 (size N/4) Cooley–Tukey steps, respectively. (The underlying idea is that the even-index subtransform of radix-2 has no multiplicative factor in front of it, so it should be left as-is, while the odd-index subtransform of radix-2 benefits by combining a second recursive subdivision.)

These smaller summations are now exactly DFTs of length N/2 and N/4, which can be performed recursively and then recombined.

More specifically, let ${\displaystyle U_{k}}$ denote the result of the DFT of length N/2 (for ${\displaystyle k=0,\ldots ,N/2-1}$), and let ${\displaystyle Z_{k}}$ and ${\displaystyle Z'_{k}}$ denote the results of the DFTs of length N/4 (for ${\displaystyle k=0,\ldots ,N/4-1}$). Then the output ${\displaystyle X_{k}}$ is simply:

${\displaystyle X_{k}=U_{k}+\omega _{N}^{k}Z_{k}+\omega _{N}^{3k}Z'_{k}.}$

This, however, performs unnecessary calculations, since ${\displaystyle k\geq N/4}$ turn out to share many calculations with ${\displaystyle k<N/4}$. In particular, if we add N/4 to k, the size-N/4 DFTs are not changed (because they are periodic in N/4), while the size-N/2 DFT is unchanged if we add N/2 to k. So, the only things that change are the ${\displaystyle \omega _{N}^{k}}$ and ${\displaystyle \omega _{N}^{3k}}$ terms, known as twiddle factors. Here, we use the identities:

${\displaystyle \omega _{N}^{k+N/4}=-i\omega _{N}^{k}}$

${\displaystyle \omega _{N}^{3(k+N/4)}=i\omega _{N}^{3k}}$

to finally arrive at:

${\displaystyle X_{k}=U_{k}+\left(\omega _{N}^{k}Z_{k}+\omega _{N}^{3k}Z'_{k}\right),}$

${\displaystyle X_{k+N/2}=U_{k}-\left(\omega _{N}^{k}Z_{k}+\omega _{N}^{3k}Z'_{k}\right),}$

${\displaystyle X_{k+N/4}=U_{k+N/4}-i\left(\omega _{N}^{k}Z_{k}-\omega _{N}^{3k}Z'_{k}\right),}$

${\displaystyle X_{k+3N/4}=U_{k+N/4}+i\left(\omega _{N}^{k}Z_{k}-\omega _{N}^{3k}Z'_{k}\right),}$

which gives all of the outputs ${\displaystyle X_{k}}$ if we let ${\displaystyle k}$ range from ${\displaystyle 0}$ to ${\displaystyle N/4-1}$ in the above four expressions.

Notice that these expressions are arranged so that we need to combine the various DFT outputs by pairs of additions and subtractions, which are known as butterflies. In order to obtain the minimal operation count for this algorithm, one needs to take into account special cases for ${\displaystyle k=0}$ (where the twiddle factors are unity) and for ${\displaystyle k=N/8}$ (where the twiddle factors are ${\displaystyle (1\pm i)/{\sqrt {2}}}$ and can be multiplied more quickly); see, e.g. Sorensen et al. (1986). Multiplications by ${\displaystyle \pm 1}$ and ${\displaystyle \pm i}$ are ordinarily counted as free (all negations can be absorbed by converting additions into subtractions or vice versa).

This decomposition is performed recursively when N is a power of two. The base cases of the recursion are N=1, where the DFT is just a copy ${\displaystyle X_{0}=x_{0}}$, and N=2, where the DFT is an addition ${\displaystyle X_{0}=x_{0}+x_{1}}$ and a subtraction ${\displaystyle X_{1}=x_{0}-x_{1}}$.

These considerations result in a count: ${\displaystyle 4N\log _{2}N-6N+8}$ real additions and multiplications, for N>1 a power of two. This count assumes that, for odd powers of 2, the leftover factor of 2 (after all the split-radix steps, which divide N by 4) is handled directly by the DFT definition (4 real additions and multiplications), or equivalently by a radix-2 Cooley–Tukey FFT step.

• R. Yavne, "An economical method for calculating the discrete Fourier transform," in Proc. AFIPS Fall Joint Computer Conf. 33, 115–125 (1968).
• P. Duhamel and H. Hollmann, "Split-radix FFT algorithm," Electron. Lett. 20 (1), 14–16 (1984).
• M. Vetterli and H. J. Nussbaumer, "Simple FFT and DCT algorithms with reduced number of operations," Signal Processing 6 (4), 267–278 (1984).
• J. B. Martens, "Recursive cyclotomic factorization—a new algorithm for calculating the discrete Fourier transform," IEEE Trans. Acoust., Speech, Signal Processing 32 (4), 750–761 (1984).
• P. Duhamel and M. Vetterli, "Fast Fourier transforms: a tutorial review and a state of the art," Signal Processing 19, 259–299 (1990).
• S. G. Johnson and M. Frigo, "A modified split-radix FFT with fewer arithmetic operations," IEEE Trans. Signal Process. 55 (1), 111–119 (2007).
• Douglas L. Jones, "Split-radix FFT algorithms," Connexions web site (Nov. 2, 2006).
• H. V. Sorensen, M. T. Heideman, and C. S. Burrus, "On computing the split-radix FFT", IEEE Trans. Acoust., Speech, Signal Processing 34 (1), 152–156 (1986).