Algoritmo FFT de base dividida

Algoritmo de transformada rápida de Fourier

La FFT de base dividida es un algoritmo de transformada rápida de Fourier (FFT) para calcular la transformada de Fourier discreta (DFT), y fue descrita por primera vez en un artículo inicialmente poco apreciado por R. Yavne (1968)[1] y posteriormente redescubierta simultáneamente por varios autores en 1984. (El nombre "base dividida" fue acuñado por dos de estos reinventores, P. Duhamel y H. Hollmann.) En particular, la base dividida es una variante del algoritmo Cooley-Tukey FFT que utiliza una combinación de raíces 2 y 4: expresa recursivamente una DFT de longitud N en términos de una DFT más pequeña de longitud N /2 y dos DFT más pequeñas de longitud N /4.

La FFT de base dividida, junto con sus variaciones, tuvo durante mucho tiempo la distinción de lograr el recuento de operaciones aritméticas más bajo publicado (número total exacto de sumas y multiplicaciones reales requeridas ) para calcular una DFT de potencia de dos tamaños N. El recuento aritmético del algoritmo de base dividida original se mejoró en 2004 (con las mejoras iniciales obtenidas en un trabajo inédito de J. Van Buskirk mediante optimización manual para N =64 [2] [3]), pero resulta que uno todavía puede alcanzar el nuevo recuento más bajo mediante una modificación de la base dividida (Johnson y Frigo, 2007). Aunque el número de operaciones aritméticas no es el único factor (ni necesariamente el factor dominante) para determinar el tiempo necesario para calcular una DFT en una computadora , la cuestión del recuento mínimo posible es de interés teórico desde hace mucho tiempo. (Actualmente no se ha demostrado ningún límite inferior estricto en el recuento de operaciones).

El algoritmo de base dividida solo se puede aplicar cuando N es múltiplo de 4, pero dado que divide una DFT en DFT más pequeñas, se puede combinar con cualquier otro algoritmo FFT según se desee.

Descomposición de base dividida

Recuerde que la DFT está definida por la fórmula:

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

donde es un número entero que va desde hasta y denota la raíz primitiva de la unidad : ${\displaystyle k}$ ${\displaystyle 0}$ ${\displaystyle N-1}$ ${\displaystyle \omega _{N}}$

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

y así: . ${\displaystyle \omega _{N}^{N}=1}$

El algoritmo de base dividida funciona expresando esta suma en términos de tres sumas más pequeñas. (Aquí, damos la versión de "diezmado en el tiempo" de la FFT de base dividida; la versión de diezmado dual en frecuencia es esencialmente lo inverso de estos pasos).

Primero, una suma de los índices pares . En segundo lugar, una suma de los índices impares divididos en dos partes: y , según si el índice es 1 o 3 módulo 4. Aquí, denota un índice que va de 0 a . Las sumas resultantes se ven así: ${\displaystyle x_{2n_{2}}}$ ${\displaystyle x_{4n_{4}+1}}$ ${\displaystyle x_{4n_{4}+3}}$ ${\displaystyle n_{m}}$ ${\displaystyle N/m-1}$

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

donde hemos utilizado el hecho de que . Estas tres sumas corresponden a porciones de los pasos Cooley-Tukey de base 2 (tamaño N /2) y base 4 (tamaño N /4), respectivamente. (La idea subyacente es que la subtransformada de índice par de radix-2 no tiene un factor multiplicativo delante, por lo que debe dejarse como está, mientras que la subtransformada de índice impar de radix-2 se beneficia al combinar una segunda subdivisión recursiva .) ${\displaystyle \omega _{N}^{mnk}=\omega _{N/m}^{nk}}$

Estas sumas más pequeñas ahora son exactamente DFT de longitud N /2 y N /4, que pueden realizarse de forma recursiva y luego recombinarse.

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

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

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

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

lo que proporciona todos los resultados si dejamos que el rango de las cuatro expresiones anteriores varíe desde hasta . ${\displaystyle X_{k}}$ ${\displaystyle k}$ ${\displaystyle 0}$ ${\displaystyle N/4-1}$

Observe que estas expresiones están organizadas de modo que necesitamos combinar las distintas salidas DFT mediante pares de sumas y restas, que se conocen como mariposas . Para obtener el recuento mínimo de operaciones para este algoritmo, es necesario tener en cuenta casos especiales para (donde los factores de giro son la unidad) y para (donde los factores de giro son y se pueden multiplicar más rápidamente); véase, por ejemplo, Sorensen et al. (1986). Las multiplicaciones por y normalmente se cuentan como libres (todas las negaciones se pueden absorber convirtiendo sumas en restas o viceversa). ${\displaystyle k=0}$ ${\displaystyle k=N/8}$ ${\displaystyle (1\pm i)/{\sqrt {2}}}$ ${\displaystyle \pm 1}$ ${\displaystyle \pm i}$

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

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

Referencias

• R. Yavne, "Un método económico para calcular la transformada discreta de Fourier", en Proc. Conf. de computadora conjunta de caída de AFIPS. 33 , 115-125 (1968).
• P. Duhamel y H. Hollmann, "Algoritmo FFT de base dividida", Electron. Letón. 20 (1), 14-16 (1984).
• M. Vetterli y HJ Nussbaumer, "Algoritmos simples FFT y DCT con número reducido de operaciones", Signal Processing 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., Habla, procesamiento de señales 32 (4), 750–761 (1984).
• P. Duhamel y M. Vetterli, "Transformadas rápidas de Fourier: una revisión del tutorial y un estado del arte", Signal Processing 19 , 259–299 (1990).
• SG Johnson y M. Frigo, "Una FFT de base dividida modificada con menos operaciones aritméticas", IEEE Trans. Proceso de señal. 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., Habla, procesamiento de señales 34 (1), 152–156 (1986).