Split-radix FFT algorithm

The split-radix FFT is a fast Fourier transform (FFT) algorithm for computing the discrete Fourier transform (DFT), and was first described in an initially little-appreciated paper by R. Yavne (1968) and subsequently rediscovered simultaneously by various authors in 1984. (The name "split radix" was coined by two of these reinventors, P. Duhamel and H. Hollmann.) In particular, split radix is a variant of the Cooley–Tukey FFT algorithm that uses a blend of radices 2 and 4: it recursively expresses a DFT of length N in terms of one smaller DFT of length N/2 and two smaller DFTs of length N/4.

The split-radix FFT, along with its variations, long had the distinction of achieving the lowest published arithmetic operation count (total exact number of required real additions and multiplications) to compute a DFT of power-of-two sizes N. The arithmetic count of the original split-radix algorithm was improved upon in 2004 (with the initial gains made in unpublished work by J. Van Buskirk via hand optimization for N=64 ), but it turns out that one can still achieve the new lowest count by a modification of split radix (Johnson and Frigo, 2007). Although the number of arithmetic operations is not the sole factor (or even necessarily the dominant factor) in determining the time required to compute a DFT on a computer, the question of the minimum possible count is of longstanding theoretical interest. (No tight lower bound on the operation count has currently been proven.)

The split-radix algorithm can only be applied when N is a multiple of 4, but since it breaks a DFT into smaller DFTs it can be combined with any other FFT algorithm as desired.

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}}$Estas tres sumas corresponden a porciones de pasos de Cooley-Tukey de base 2 (tamaño N /2) y de base 4 (tamaño N /4), respectivamente. (La idea subyacente es que la subtransformación de índice par de base 2 no tiene ningún factor multiplicativo delante, por lo que debe dejarse tal cual, mientras que la subtransformación de índice impar de base 2 se beneficia al combinarse con una segunda subdivisión recursiva).

Estas sumas más pequeñas son ahora exactamente transformadas discretas de Fourier (DFT) de longitud N /2 y N /4, que se pueden realizar de forma recursiva y luego recombinar.

Más específicamente, dejemos${\displaystyle U_{k}}$denota el resultado de la DFT de longitud N /2 (para${\displaystyle k=0,\ldots ,N/2-1}$), y deja${\displaystyle Z_{k}}$y${\displaystyle Z'_{k}}$denotan los resultados de las DFT de longitud N /4 (para${\displaystyle k=0,\ldots ,N/4-1}$). Luego la salida${\displaystyle X_{k}}$es simplemente:

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

Sin embargo, esto realiza cálculos innecesarios, ya que${\displaystyle k\geq N/4}$resulta que comparte muchos cálculos con${\displaystyle k<N/4}$. En particular, si añadimos N /4 a k , las DFT de tamaño N /4 no cambian (porque son periódicas en N /4), mientras que la DFT de tamaño N /2 no cambia si añadimos N /2 a k . Por lo tanto, lo único que cambia son las${\displaystyle \omega _ {N}^{k}}$y${\displaystyle \omega _ {N}^{3k}}$términos, conocidos como factores de rotación . Aquí, utilizamos las identidades:

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

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

para finalmente llegar a:

${\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),}$

que proporciona todos los resultados${\displaystyle X_{k}}$si dejamos${\displaystyle k}$abarca desde${\displaystyle 0}$a${\displaystyle N/4-1}$en las cuatro expresiones anteriores.

Nótese que estas expresiones están dispuestas de tal manera que necesitamos combinar las distintas salidas de la DFT mediante pares de sumas y restas, que se conocen como mariposas . Para obtener el número mínimo de operaciones para este algoritmo, es necesario tener en cuenta casos especiales para${\displaystyle k=0}$(donde los factores de rotación son la unidad) y para${\displaystyle k=N/8}$(donde los factores de ajuste son${\displaystyle (1\pm i)/{\sqrt {2}}}$y se pueden multiplicar más rápidamente); véase, por ejemplo, Sorensen et al. (1986). Multiplicaciones por${\displaystyle \pm 1}$y${\displaystyle \pm i}$Por lo general, se consideran libres (todas las negaciones pueden absorberse convirtiendo sumas en restas o viceversa).

Esta descomposición se realiza recursivamente cuando N es una potencia de dos. Los casos base de la recursión son N = 1, donde la DFT es simplemente una copia.${\displaystyle X_{0}=x_{0}}$y N =2, donde la DFT es una adición${\displaystyle X_{0}=x_{0}+x_{1}}$y una resta${\displaystyle X_{1}=x_{0}-x_{1}}$.

Estas consideraciones dan como resultado un recuento:${\displaystyle 4N\log _{2}N-6N+8}$sumas y multiplicaciones reales, para N > 1 una potencia de dos. Este recuento supone que, para potencias impares de 2, el factor restante de 2 (después de todos los pasos de base dividida, que dividen N por 4) se maneja directamente mediante la definición de DFT (4 sumas y multiplicaciones reales), o equivalentemente mediante un paso de FFT de Cooley-Tukey de base 2.

• R. Yavne, " Un método económico para calcular la transformada discreta de Fourier ", en Proc. AFIPS Fall Joint Computer Conf. 33 , 115–125 (1968).
• P. Duhamel y H. Hollmann, "Algoritmo FFT de raíz dividida", Electron. Lett. 20 (1), 14–16 (1984).
• M. Vetterli y HJ Nussbaumer, "Algoritmos simples de FFT y DCT con un número reducido de operaciones", Procesamiento de señales 6 (4), 267–278 (1984).
• JB Martens, "Factorización ciclotómica recursiva: un nuevo algoritmo para calcular la transformada discreta de Fourier", IEEE Trans. Acoust., Speech, Signal Processing 32 (4), 750–761 (1984).
• P. Duhamel y M. Vetterli, "Transformadas rápidas de Fourier: una revisión tutorial y el estado del arte", Procesamiento de señales 19 , 259–299 (1990).
• SG Johnson y M. Frigo, " Una FFT de base dividida modificada con menos operaciones aritméticas ", IEEE Trans. Signal Process. 55 (1), 111–119 (2007).
• Douglas L. Jones, " Algoritmos FFT de base dividida ", sitio web de Connexions (2 de noviembre de 2006).
• HV Sorensen, MT Heideman y CS Burrus, "Sobre el cálculo de la FFT de base dividida", IEEE Trans. Acoust., Speech, Signal Processing 34 (1), 152–156 (1986).