Algoritmo de FFT de Cooley-Tukey

Algoritmo de transformada rápida de Fourier

El algoritmo Cooley–Tukey , llamado así por JW Cooley y John Tukey , es el algoritmo de transformada rápida de Fourier (FFT) más común. Reexpresa la transformada discreta de Fourier (DFT) de un tamaño compuesto arbitrario en términos de N ₁ DFT más pequeñas de tamaños N ₂ , de forma recursiva , para reducir el tiempo de cálculo a O( N log N ) para N altamente compuestos ( números suaves ). Debido a la importancia del algoritmo, las variantes específicas y los estilos de implementación se han hecho conocidos por sus propios nombres, como se describe a continuación. ${\displaystyle N=N_{1}N_{2}}$

Debido a que el algoritmo Cooley-Tukey divide la DFT en DFT más pequeñas, se puede combinar arbitrariamente con cualquier otro algoritmo para la DFT. Por ejemplo, el algoritmo de Rader o Bluestein se puede utilizar para manejar factores primos grandes que no se pueden descomponer mediante Cooley-Tukey, o el algoritmo de factores primos se puede aprovechar para lograr una mayor eficiencia en la separación de factores primos relativos .

El algoritmo, junto con su aplicación recursiva, fue inventado por Carl Friedrich Gauss . Cooley y Tukey lo redescubrieron y popularizaron de forma independiente 160 años después.

Historia

Este algoritmo, incluida su aplicación recursiva, fue inventado alrededor de 1805 por Carl Friedrich Gauss , quien lo utilizó para interpolar las trayectorias de los asteroides Pallas y Juno , pero su trabajo no fue ampliamente reconocido (siendo publicado solo póstumamente y en neolatín ). ^{[1] [2]} Sin embargo, Gauss no analizó el tiempo computacional asintótico. Varias formas limitadas también fueron redescubiertas varias veces a lo largo del siglo XIX y principios del XX. ^[2] Las FFT se hicieron populares después de que James Cooley de IBM y John Tukey de Princeton publicaran un artículo en 1965 reinventando ^[2] el algoritmo y describiendo cómo ejecutarlo cómodamente en una computadora. ^[3]

Según se informa, a Tukey se le ocurrió la idea durante una reunión del Comité Asesor Científico del Presidente Kennedy en la que se discutían formas de detectar pruebas de armas nucleares en la Unión Soviética empleando sismómetros ubicados fuera del país. Estos sensores generarían series temporales sismológicas. Sin embargo, el análisis de estos datos requeriría algoritmos rápidos para calcular las DFT debido a la cantidad de sensores y la duración del tiempo. Esta tarea fue fundamental para la ratificación de la prohibición de pruebas nucleares propuesta, de modo que se pudiera detectar cualquier violación sin necesidad de visitar las instalaciones soviéticas. ^{[4] [5]} Otro participante en esa reunión, Richard Garwin de IBM, reconoció el potencial del método y puso a Tukey en contacto con Cooley. Sin embargo, Garwin se aseguró de que Cooley no conociera el propósito original. En cambio, le dijeron a Cooley que esto era necesario para determinar las periodicidades de las orientaciones de espín en un cristal 3-D de helio-3 . Cooley y Tukey publicaron posteriormente su artículo conjunto, y rápidamente se adoptó ampliamente debido al desarrollo simultáneo de convertidores analógicos a digitales capaces de muestrear a velocidades de hasta 300 kHz.

El hecho de que Gauss había descrito el mismo algoritmo (aunque sin analizar su coste asintótico) no se comprendió hasta varios años después del artículo de Cooley y Tukey de 1965. ^[2] Su artículo citaba como inspiración solo el trabajo de IJ Good sobre lo que ahora se llama algoritmo FFT de factores primos (PFA); ^[3] aunque inicialmente se pensó que el algoritmo de Good era equivalente al algoritmo de Cooley-Tukey, rápidamente se comprendió que PFA es un algoritmo bastante diferente (que funciona solo para tamaños que tienen factores primos relativos y se basa en el teorema del resto chino , a diferencia del soporte para cualquier tamaño compuesto en Cooley-Tukey). ^[6]

El caso DIT de base 2

Una FFT de diezmado en el tiempo ( DIT ) de base 2 es la forma más simple y común del algoritmo Cooley-Tukey, aunque las implementaciones de Cooley-Tukey altamente optimizadas suelen utilizar otras formas del algoritmo, como se describe a continuación. La DIT de base 2 divide una DFT de tamaño N en dos DFT intercaladas (de ahí el nombre "de base 2") de tamaño N /2 con cada etapa recursiva.

La transformada de Fourier discreta (DFT) se define mediante la fórmula:

${\displaystyle X_{k}=\sum _{n=0}^{N-1}x_{n}e^{-{\frac {2\pi i}{N}}nk},}$

donde es un número entero que va de 0 a . ${\displaystyle k}$ ${\displaystyle N-1}$

Radix-2 DIT primero calcula las DFT de las entradas de índice par y de las entradas de índice impar , y luego combina esos dos resultados para producir la DFT de toda la secuencia. Esta idea se puede llevar a cabo de forma recursiva para reducir el tiempo de ejecución total a O( N log N ). Esta forma simplificada supone que N es una potencia de dos ; dado que la cantidad de puntos de muestra N generalmente se puede elegir libremente por la aplicación (por ejemplo, cambiando la frecuencia de muestreo o la ventana, rellenando con ceros, etcétera), esto no suele ser una restricción importante. ${\displaystyle (x_{2m}=x_{0},x_{2},\ldots ,x_{N-2})}$ ${\displaystyle (x_{2m+1}=x_{1},x_{3},\ldots ,x_{N-1})}$

El algoritmo DIT de base 2 reorganiza la DFT de la función en dos partes: una suma sobre los índices pares y una suma sobre los índices impares : ${\displaystyle x_{n}}$ ${\displaystyle n={2m}}$ ${\displaystyle n={2m+1}}$

${\displaystyle {\begin{matrix}X_{k}&=&\sum \limits _{m=0}^{N/2-1}x_{2m}e^{-{\frac {2\pi i}{N}}(2m)k}+\sum \limits _{m=0}^{N/2-1}x_{2m+1}e^{-{\frac {2\pi i}{N}}(2m+1)k}\end{matrix}}}$

Se puede factorizar un multiplicador común a partir de la segunda suma, como se muestra en la ecuación siguiente. Entonces queda claro que las dos sumas son la DFT de la parte de índice par y la DFT de la parte de índice impar de la función . Denotemos la DFT de las entradas de índice par por y la DFT de las entradas de índice O por y obtenemos: ${\displaystyle e^{-{\frac {2\pi i}{N}}k}}$ ${\displaystyle x_{2m}}$ ${\displaystyle x_{2m+1}}$ ${\displaystyle x_{n}}$ ${\displaystyle x_{2m}}$ ${\displaystyle E_{k}}$ ${\displaystyle x_{2m+1}}$ ${\displaystyle O_{k}}$

${\displaystyle {\begin{matrix}X_{k}=\underbrace {\sum \limits _{m=0}^{N/2-1}x_{2m}e^{-{\frac {2\pi i}{N/2}}mk}} _{\mathrm {DFT\;of\;even-indexed\;part\;of\;} x_{n}}{}+e^{-{\frac {2\pi i}{N}}k}\underbrace {\sum \limits _{m=0}^{N/2-1}x_{2m+1}e^{-{\frac {2\pi i}{N/2}}mk}} _{\mathrm {DFT\;of\;odd-indexed\;part\;of\;} x_{n}}=E_{k}+e^{-{\frac {2\pi i}{N}}k}O_{k}\qquad {\text{ for }}k=0,\dots ,{\frac {N}{2}}-1.\end{matrix}}}$

Nótese que las igualdades se cumplen para , pero el quid de la cuestión es que y se calculan de esta manera solo para . Gracias a la periodicidad de la exponencial compleja , también se obtiene a partir de y : ${\displaystyle k=0,\dots ,N-1}$ ${\displaystyle E_{k}}$ ${\displaystyle O_{k}}$ ${\displaystyle k=0,\dots ,{\frac {N}{2}}-1}$ ${\displaystyle X_{k+{\frac {N}{2}}}}$ ${\displaystyle E_{k}}$ ${\displaystyle O_{k}}$

${\displaystyle {\begin{aligned}X_{k+{\frac {N}{2}}}&=\sum \limits _{m=0}^{N/2-1}x_{2m}e^{-{\frac {2\pi i}{N/2}}m(k+{\frac {N}{2}})}+e^{-{\frac {2\pi i}{N}}(k+{\frac {N}{2}})}\sum \limits _{m=0}^{N/2-1}x_{2m+1}e^{-{\frac {2\pi i}{N/2}}m(k+{\frac {N}{2}})}\\&=\sum \limits _{m=0}^{N/2-1}x_{2m}e^{-{\frac {2\pi i}{N/2}}mk}e^{-2\pi mi}+e^{-{\frac {2\pi i}{N}}k}e^{-\pi i}\sum \limits _{m=0}^{N/2-1}x_{2m+1}e^{-{\frac {2\pi i}{N/2}}mk}e^{-2\pi mi}\\&=\sum \limits _{m=0}^{N/2-1}x_{2m}e^{-{\frac {2\pi i}{N/2}}mk}-e^{-{\frac {2\pi i}{N}}k}\sum \limits _{m=0}^{N/2-1}x_{2m+1}e^{-{\frac {2\pi i}{N/2}}mk}\\&=E_{k}-e^{-{\frac {2\pi i}{N}}k}O_{k}\end{aligned}}}$

Podemos reescribir y como: ${\displaystyle X_{k}}$ ${\displaystyle X_{k+{\frac {N}{2}}}}$

${\displaystyle {\begin{matrix}X_{k}&=&E_{k}+e^{-{\frac {2\pi i}{N}}{k}}O_{k}\\X_{k+{\frac {N}{2}}}&=&E_{k}-e^{-{\frac {2\pi i}{N}}{k}}O_{k}\end{matrix}}}$

Este resultado, que expresa la DFT de longitud N recursivamente en términos de dos DFT de tamaño N /2, es el núcleo de la transformada rápida de Fourier DIT de base 2. El algoritmo gana velocidad al reutilizar los resultados de los cálculos intermedios para calcular múltiples salidas DFT. Nótese que las salidas finales se obtienen mediante una combinación +/− de y , que es simplemente una DFT de tamaño 2 (a veces llamada mariposa en este contexto); cuando esto se generaliza a bases mayores a continuación, la DFT de tamaño 2 se reemplaza por una DFT más grande (que a su vez se puede evaluar con una FFT). ${\displaystyle E_{k}}$ ${\displaystyle O_{k}\exp(-2\pi ik/N)}$

Diagrama de flujo de datos para N = 8: una FFT de radix 2 con diezmado en el tiempo divide una DFT de longitud N en dos DFT de longitud N /2 seguida de una etapa de combinación que consta de N /2 DFT de tamaño 2 llamadas operaciones "mariposa" (llamadas así debido a la forma de los diagramas de flujo de datos).

Este proceso es un ejemplo de la técnica general de los algoritmos de dividir y vencer ; sin embargo, en muchas implementaciones convencionales se evita la recursión explícita y en su lugar se recorre el árbol computacional en amplitud .

La reexpresión anterior de una DFT de tamaño N como dos DFT de tamaño N /2 a veces se denomina lema de Danielson - Lanczos , ya que la identidad fue notada por esos dos autores en 1942 ^[7] (influenciados por el trabajo de Runge de 1903 ^[2] ). Aplicaron su lema de manera recursiva "hacia atrás", duplicando repetidamente el tamaño de la DFT hasta que el espectro de la transformada convergiera (aunque aparentemente no se dieron cuenta de la complejidad asintótica linealítmica [es decir, orden N  log  N ] que habían logrado). El trabajo de Danielson-Lanczos fue anterior a la disponibilidad generalizada de computadoras mecánicas o electrónicas y requirió cálculo manual (posiblemente con ayudas mecánicas como máquinas sumadoras ); informaron un tiempo de cálculo de 140 minutos para una DFT de tamaño 64 que operaba con entradas reales de 3 a 5 dígitos significativos. El artículo de 1965 de Cooley y Tukey informó un tiempo de ejecución de 0,02 minutos para una DFT compleja de tamaño 2048 en un IBM 7094 (probablemente en precisión simple de 36 bits , ~8 dígitos). ^[3] Reescalando el tiempo por el número de operaciones, esto corresponde aproximadamente a un factor de aceleración de alrededor de 800.000. (Para poner el tiempo para el cálculo manual en perspectiva, 140 minutos para tamaño 64 corresponden a un promedio de como máximo 16 segundos por operación de punto flotante, alrededor del 20% de los cuales son multiplicaciones).

Pseudocódigo

En pseudocódigo , el siguiente procedimiento podría escribirse: ^[8]

X 0,..., N −1 ← ditfft2 ( x , N , s ): DFT de (x 0 , x s , x 2 s , ..., x ( N -1) s ): si  N = 1 entonces  X 0 ← x 0  caso base de DFT trivial de tamaño 1  de lo contrario  X 0,..., N /2−1 ← ditfft2 ( x , N /2, 2 s ) DFT de (x 0 , x 2 s , x 4 s , ..., x ( N -2) s ) X N /2,..., N −1 ← ditfft2 ( x +s, N /2, 2 s ) DFT de (x s , x s +2 s , x s +4 s , ..., x ( N -1) s ) para  k = 0 a N /2−1 combina  las DFT de dos mitades en una DFT completa: p ← X k q ← exp(−2π i / N  k ) X k + N /2  X k ← p + q X k + N /2 ← p − q fin para  fin si

Aquí, ditfft2( x , N ,1), calcula X = DFT( x ) fuera de lugar mediante una FFT DIT de base 2, donde N es una potencia entera de 2 y s = 1 es el paso de la matriz de entrada x . x + s denota la matriz que comienza con x _s .

(Los resultados están en el orden correcto en X y no se requiere ninguna otra permutación de inversión de bits ; la necesidad, a menudo mencionada, de una etapa de inversión de bits separada solo surge para ciertos algoritmos en el lugar, como se describe a continuación).

Las implementaciones de FFT de alto rendimiento hacen muchas modificaciones a la implementación de dicho algoritmo en comparación con este pseudocódigo simple. Por ejemplo, se puede utilizar un caso base más grande que N = 1 para amortizar la sobrecarga de la recursión, los factores de twiddle se pueden calcular previamente y a menudo se utilizan bases más grandes por razones de caché ; estas y otras optimizaciones juntas pueden mejorar el rendimiento en un orden de magnitud o más. ^[8] (En muchas implementaciones de libros de texto, la recursión en profundidad se elimina a favor de un enfoque no recursivo en amplitud , aunque se ha argumentado que la recursión en profundidad tiene una mejor localidad de memoria . ^{[8] [9]} ) Varias de estas ideas se describen con más detalle a continuación. ${\displaystyle \exp[-2\pi ik/N]}$

Idea

El paso básico de la FFT de Cooley–Tukey para factorizaciones generales puede verse como la reinterpretación de una DFT unidimensional como algo así como una DFT bidimensional. La matriz de entrada unidimensional de longitud N = N ₁ N ₂ se reinterpreta como una matriz bidimensional N ₁ × N _{2 almacenada en} orden de columna principal . Se realizan DFT unidimensionales más pequeñas a lo largo de la dirección N ₂ (la dirección no contigua), luego se multiplica por factores de fase (factores de twiddle) y, finalmente, se realizan DFT unidimensionales a lo largo de la dirección N _1. El paso de transposición se puede realizar en el medio, como se muestra aquí, o al principio o al final. Esto se hace de forma recursiva para las transformaciones más pequeñas.

De manera más general, los algoritmos de Cooley-Tukey reexpresan recursivamente una DFT de un tamaño compuesto N = N ₁ N ₂ como: ^[10]

1. Realizar N ₁ DFT de tamaño N ₂ .
2. Multiplica por raíces complejas de la unidad (a menudo llamadas factores de torsión ).
3. Realizar N ₂ DFT de tamaño N ₁ .

Por lo general, N ₁ o N ₂ es un factor pequeño ( no necesariamente primo), llamado base (que puede diferir entre las etapas de la recursión). Si N ₁ es la base, se denomina algoritmo de diezmado en el tiempo (DIT), mientras que si N ₂ es la base, se denomina algoritmo de diezmado en frecuencia (DIF, también llamado algoritmo de Sande-Tukey). La versión presentada anteriormente era un algoritmo DIT de base 2; en la expresión final, la fase que multiplica la transformada impar es el factor twiddle, y la combinación +/- ( mariposa ) de las transformadas par e impar es una DFT de tamaño 2. (La DFT pequeña de la base a veces se conoce como mariposa , llamada así por la forma del diagrama de flujo de datos para el caso de base 2).

Variaciones

Existen muchas otras variaciones del algoritmo Cooley-Tukey. Las implementaciones de base mixta manejan tamaños compuestos con una variedad de factores (normalmente pequeños) además de dos, empleando habitualmente (pero no siempre) el algoritmo O( N ² ) para los casos base primos de la recursión (también es posible emplear un algoritmo N  log  N para los casos base primos, como el algoritmo de Rader o Bluestein ). La base dividida fusiona las bases 2 y 4, explotando el hecho de que la primera transformación de la base 2 no requiere factor de torsión, para lograr lo que durante mucho tiempo fue el recuento de operaciones aritméticas más bajo conocido para tamaños de potencia de dos, ^[10] aunque las variaciones recientes logran un recuento incluso menor. ^{[11] [12]} (En las computadoras actuales, el rendimiento está determinado más por consideraciones de caché y canalización de CPU que por recuentos de operaciones estrictos; las implementaciones de FFT bien optimizadas a menudo emplean bases más grandes y/o transformaciones de caso base codificadas de forma rígida de tamaño significativo. ^[13] ).

Otra forma de ver el algoritmo Cooley–Tukey es que reexpresa una DFT unidimensional de tamaño N como una DFT bidimensional N ₁ por N _{2 (más twiddles), donde la matriz de salida está} transpuesta . El resultado neto de todas estas transposiciones, para un algoritmo de base 2, corresponde a una inversión de bits de los índices de entrada (DIF) o salida (DIT). Si, en lugar de utilizar una base pequeña, se emplea una base de aproximadamente √ N y transposiciones explícitas de la matriz de entrada/salida, se denomina algoritmo FFT de cuatro pasos (o de seis pasos , dependiendo del número de transposiciones), inicialmente propuesto para mejorar la localidad de la memoria, ^{[14] [15]} por ejemplo para la optimización de la caché o la operación fuera del núcleo , y más tarde se demostró que era un algoritmo óptimo que ignora la caché . ^[16]

La factorización general de Cooley-Tukey reescribe los índices k y n como y , respectivamente, donde los índices k _a y n _a van de 0 a N _a -1 (para a de 1 o 2). Es decir, reindexa la entrada ( n ) y la salida ( k ) como N ₁ por N ₂ matrices bidimensionales en orden de columna y fila , respectivamente; la diferencia entre estas indexaciones es una transposición, como se mencionó anteriormente. Cuando esta reindexación se sustituye en la fórmula DFT para nk , el término cruzado se anula (su exponencial es la unidad) y los términos restantes dan ${\displaystyle k=N_{2}k_{1}+k_{2}}$ ${\displaystyle n=N_{1}n_{2}+n_{1}}$ ${\displaystyle N_{1}n_{2}N_{2}k_{1}}$

${\displaystyle X_{N_{2}k_{1}+k_{2}}=\sum _{n_{1}=0}^{N_{1}-1}\sum _{n_{2}=0}^{N_{2}-1}x_{N_{1}n_{2}+n_{1}}e^{-{\frac {2\pi i}{N_{1}N_{2}}}\cdot (N_{1}n_{2}+n_{1})\cdot (N_{2}k_{1}+k_{2})}}$

${\displaystyle =\sum _{n_{1}=0}^{N_{1}-1}\left[e^{-{\frac {2\pi i}{N_{1}N_{2}}}n_{1}k_{2}}\right]\left(\sum _{n_{2}=0}^{N_{2}-1}x_{N_{1}n_{2}+n_{1}}e^{-{\frac {2\pi i}{N_{2}}}n_{2}k_{2}}\right)e^{-{\frac {2\pi i}{N_{1}}}n_{1}k_{1}}}$

${\displaystyle =\sum _{n_{1}=0}^{N_{1}-1}\left(\sum _{n_{2}=0}^{N_{2}-1}x_{N_{1}n_{2}+n_{1}}e^{-{\frac {2\pi i}{N_{2}}}n_{2}k_{2}}\right)e^{-{\frac {2\pi i}{N_{1}N_{2}}}n_{1}(N_{2}k_{1}+k_{2})}}$.

donde cada suma interna es una DFT de tamaño N ₂ , cada suma externa es una DFT de tamaño N ₁ , y el término entre corchetes [...] es el factor de torsión.

Se puede emplear un radio arbitrario r (así como radios mixtos), como lo demostraron tanto Cooley y Tukey ^[3] como Gauss (quien dio ejemplos de pasos de radio 3 y radio 6). ^[2] Cooley y Tukey asumieron originalmente que la mariposa de radio requería trabajo O( r ² ) y por lo tanto calcularon la complejidad para un radio r como O( r ²  N / r  log _r N ) = O( N  log ₂ ( N )  r /log ₂ r ); a partir del cálculo de valores de r /log ₂ r para valores enteros de r de 2 a 12, se descubre que el radio óptimo es 3 (el entero más cercano a e , que minimiza r /log ₂ r ). ^{[3] [17]} Sin embargo, este análisis era erróneo: la mariposa de base también es una DFT y se puede realizar a través de un algoritmo FFT en operaciones O( r log r ), por lo tanto, la base r en realidad se cancela en la complejidad O( r  log( r )  N / r  log _r N ), y la r óptima se determina mediante consideraciones más complicadas. En la práctica, r bastante grandes (32 o 64) son importantes para explotar de manera efectiva, por ejemplo, la gran cantidad de registros de procesador en los procesadores modernos, ^[13] e incluso una base ilimitada r = √ N también logra una complejidad O( N  log  N ) y tiene ventajas teóricas y prácticas para N grande como se mencionó anteriormente. ^{[14] [15] [16]}

Reordenamiento de datos, inversión de bits y algoritmos in situ

Aunque la factorización abstracta de Cooley-Tukey de la DFT, antes mencionada, se aplica de alguna forma a todas las implementaciones del algoritmo, existe una diversidad mucho mayor en las técnicas para ordenar y acceder a los datos en cada etapa de la FFT. De especial interés es el problema de diseñar un algoritmo in situ que sobrescriba su entrada con sus datos de salida utilizando solo almacenamiento auxiliar O(1).

La técnica de reordenamiento más conocida implica la inversión explícita de bits para algoritmos de base 2 en el lugar. La inversión de bits es la permutación donde los datos en un índice n , escritos en binario con dígitos b ₄ b ₃ b ₂ b ₁ b ₀ (por ejemplo, 5 dígitos para N = 32 entradas), se transfieren al índice con dígitos invertidos b ₀ b ₁ b ₂ b ₃ b ₄ . Considere la última etapa de un algoritmo DIT de base 2 como el presentado anteriormente, donde la salida se escribe en el lugar sobre la entrada: cuando y se combinan con una DFT de tamaño 2, esos dos valores se sobrescriben con las salidas. Sin embargo, los dos valores de salida deben ir en la primera y segunda mitad de la matriz de salida, correspondientes al bit más significativo b ₄ (para N = 32); mientras que las dos entradas y se intercalan en los elementos pares e impares, correspondientes al bit menos significativo b ₀ . Por lo tanto, para obtener la salida en el lugar correcto, b ₀ debe tomar el lugar de b ₄ y el índice se convierte en b ₀ b ₄ b ₃ b ₂ b ₁ . Y para la siguiente etapa recursiva, esos 4 bits menos significativos se convertirán en b ₁ b ₄ b ₃ b ₂ . Si incluye todas las etapas recursivas de un algoritmo DIT de base 2, todos los bits deben invertirse y, por lo tanto, uno debe preprocesar la entrada ( o posprocesar la salida) con una inversión de bits para obtener una salida en orden. (Si cada subtransformación de tamaño N /2 debe operar en datos contiguos, la entrada DIT se preprocesa mediante inversión de bits). Correspondientemente, si realiza todos los pasos en orden inverso, obtiene un algoritmo DIF de base 2 con inversión de bits en el posprocesamiento (o preprocesamiento, respectivamente). ${\displaystyle E_{k}}$ ${\displaystyle O_{k}}$ ${\displaystyle E_{k}}$ ${\displaystyle O_{k}}$

El logaritmo (log) utilizado en este algoritmo es un logaritmo de base 2.

El siguiente es un pseudocódigo para un algoritmo FFT iterativo de base 2 implementado mediante permutación de inversión de bits. ^[18]

El algoritmo iterativo-fft tiene  como entrada una matriz a de n valores complejos donde n es una potencia de 2. Salida: la matriz A, la DFT de a.  copia inversa de bit(a, A) n ← a .length para  s = 1 a log( n ) do  m ← 2 s  ω m ← exp(−2π i / m ) para  k = 0 a  n -1 por  m  do  ω ← 1 para  j = 0 a  m / 2 – 1 do  t ← ω  A [ k + j + m /2] u ← A [ k + j ] A [ k + j ] ← u + t  A [ k + j + m /2] ← u – t  ω ← ω  ω m  devolver  A

El procedimiento de copia inversa de bits se puede implementar de la siguiente manera.

El algoritmo bit-reverse-copy( a , A ) tiene  como entrada: una matriz a de n valores complejos donde n es una potencia de 2. Salida: una matriz A de tamaño n . n ← a .longitud para  k = 0 a  n – 1 hacer  A [rev(k)] := a [k]

Como alternativa, algunas aplicaciones (como la convolución) funcionan igualmente bien con datos con inversión de bits, por lo que se pueden realizar transformaciones hacia adelante, procesamiento y luego transformaciones inversas, todo sin inversión de bits para producir resultados finales en el orden natural.

Sin embargo, muchos usuarios de FFT prefieren salidas de orden natural, y una etapa de inversión de bits explícita y separada puede tener un impacto no despreciable en el tiempo de cálculo, ^[13] aunque la inversión de bits se puede realizar en tiempo O( N ) y ha sido objeto de mucha investigación. ^{[19] [20] [21]} Además, mientras que la permutación es una inversión de bits en el caso de base 2, es más generalmente una inversión de dígitos arbitraria (de base mixta) para el caso de base mixta, y los algoritmos de permutación se vuelven más complicados de implementar. Además, es deseable en muchas arquitecturas de hardware reordenar las etapas intermedias del algoritmo FFT para que operen en elementos de datos consecutivos (o al menos más localizados). Con estos fines, se han ideado varios esquemas de implementación alternativos para el algoritmo Cooley-Tukey que no requieren una inversión de bits separada y/o involucran permutaciones adicionales en etapas intermedias.

El problema se simplifica enormemente si está fuera de lugar : la matriz de salida es distinta de la matriz de entrada o, equivalentemente, hay disponible una matriz auxiliar de igual tamaño.El algoritmo de autoordenamiento de Stockham ^[22][23]realiza cada etapa de la FFT fuera de lugar, típicamente escribiendo de ida y vuelta entre dos matrices, transponiendo un "dígito" de los índices con cada etapa, y ha sido especialmente popular enarquitecturasSIMD^[23][24] Se han propuesto ventajas SIMD potenciales aún mayores (más accesos consecutivos) para elPease,^[25]que también reordena fuera de lugar con cada etapa, pero este método requiere inversión de bit/dígito separada y almacenamiento O(NlogN). También se puede aplicar directamente la definición de factorización de Cooley-Tukey con recursión explícita (en profundidad) y raíces pequeñas, que produce una salida fuera de lugar de orden natural sin un paso de permutación separado (como en el pseudocódigo anterior) y se puede argumentar que tienebeneficios de localidadsin tener en cuenta la cachémemoria jerárquica.^[9][13][26]

Una estrategia típica para algoritmos locales sin almacenamiento auxiliar y sin pases de inversión de dígitos separados implica pequeñas transposiciones de matriz (que intercambian pares individuales de dígitos) en etapas intermedias, que se pueden combinar con las mariposas de radix para reducir la cantidad de pases sobre los datos. ^{[13] [27] [28] [29] [30]}

Referencias

1. Gauss, Carl Friedrich (1866). "Theoria interpolationis Methodo nova tractata" [Teoría sobre un nuevo método de interpolación]. Nachlass (manuscrito inédito). Werke (en latín y alemán). vol. 3. Göttingen, Alemania: Königlichen Gesellschaft der Wissenschaften zu Göttingen. págs. 265–303.
2. Heideman, MT, DH Johnson y CS Burrus , "Gauss y la historia de la transformada rápida de Fourier", IEEE ASSP Magazine, 1, (4), 14–21 (1984)
3. Cooley, James W.; Tukey, John W. (1965). "Un algoritmo para el cálculo automático de series complejas de Fourier". Matemáticas. Computación. 19 (90): 297–301. doi : 10.2307/2003354 . JSTOR  2003354.
4. Cooley, James W.; Lewis, Peter AW; Welch, Peter D. (1967). "Notas históricas sobre la transformada rápida de Fourier" (PDF) . IEEE Transactions on Audio and Electroacoustics . 15 (2): 76–79. CiteSeerX 10.1.1.467.7209 . doi :10.1109/tau.1967.1161903. 
5. Rockmore, Daniel N., Comput. Sci. Eng. 2 (1), 60 (2000). La FFT: un algoritmo que puede utilizar toda la familia Número especial sobre los "diez mejores algoritmos del siglo" Barry A. Cipra. "Lo mejor del siglo XX: los editores nombran los 10 mejores algoritmos" (PDF) . SIAM News . 33 (4). Archivado desde el original (PDF) el 2009-04-07 . Consultado el 2009-03-31 .
6. James W. Cooley, Peter AW Lewis y Peter W. Welch, "Notas históricas sobre la transformada rápida de Fourier", Proc. IEEE , vol. 55 (n.º 10), pág. 1675–1677 (1967).
7. Danielson, GC y C. Lanczos, "Algunas mejoras en el análisis práctico de Fourier y su aplicación a la dispersión de rayos X de líquidos", J. Franklin Inst. 233 , 365–380 y 435–452 (1942).
8. SG Johnson y M. Frigo, "Implementación de FFT en la práctica", en Fast Fourier Transforms (CS Burrus, ed.), cap. 11, Rice University, Houston TX: Connexions, septiembre de 2008.
9. Singleton, Richard C. (1967). "Sobre el cálculo de la transformada rápida de Fourier". Commun. ACM . 10 (10): 647–654. doi : 10.1145/363717.363771 . S2CID  6287781.
10. Duhamel, P., y M. Vetterli, "Transformadas rápidas de Fourier: una revisión tutorial y un estado del arte", Signal Processing 19 , 259–299 (1990)
11. Lundy, T. y J. Van Buskirk, "Un nuevo enfoque matricial para FFT reales y convoluciones de longitud 2 ^k ", Computing 80 , 23–45 (2007).
12. Johnson, SG y M. Frigo, "Una FFT de base dividida modificada con menos operaciones aritméticas", IEEE Trans. Signal Process. 55 (1), 111–119 (2007).
13. Frigo, M.; Johnson, SG (2005). "El diseño y la implementación de FFTW3" (PDF) . Actas del IEEE . 93 (2): 216–231. Bibcode :2005IEEEP..93..216F. CiteSeerX 10.1.1.66.3097 . doi :10.1109/JPROC.2004.840301. S2CID  6644892. 
14. Gentleman WM y G. Sande, "Transformadas rápidas de Fourier: por diversión y beneficio", Proc. AFIPS 29 , 563–578 (1966).
15. Bailey, David H., "FFT en memoria externa o jerárquica", J. Supercomputing 4 (1), 23–35 (1990)
16. M. Frigo, CE Leiserson, H. Prokop y S. Ramachandran. Algoritmos que ignoran la memoria caché. En Actas del 40.º Simposio IEEE sobre Fundamentos de la Ciencia de la Computación (FOCS 99), págs. 285-297. 1999. Resumen ampliado en IEEE, en Citeseer.
17. Cooley, JW, P. Lewis y P. Welch, "La transformada rápida de Fourier y sus aplicaciones", IEEE Trans on Education 12 , 1, 28–34 (1969)
18. Cormen, Thomas H.; Leiserson, Charles; Rivest, Ronald; Stein, Clifford (2009). Introducción a los algoritmos (3ª ed.). Cambridge, Massachusetts: MIT Press. págs. 915–918. ISBN 978-0-262-03384-8.
19. Karp, Alan H. (1996). "Inversión de bits en procesadores monoprocesadores". SIAM Review . 38 (1): 1–26. CiteSeerX 10.1.1.24.2913 . doi :10.1137/1038001. JSTOR  2132972. 
20. Carter, Larry; Gatlin, Kang Su (1998). "Hacia un programa óptimo de permutación de inversión de bits". Actas del 39.° Simposio anual sobre fundamentos de la informática (Cat. N.° 98CB36280) . págs. 544–553. CiteSeerX 10.1.1.46.9319 . doi :10.1109/SFCS.1998.743505. ISBN .  978-0-8186-9172-0. Número de identificación del sujeto  14307262.
21. Rubio, M.; Gómez, P.; Drouiche, K. (2002). "Un nuevo algoritmo de inversión de bits superrápido". Revista Internacional de Control Adaptativo y Procesamiento de Señales . 16 (10): 703–707. doi :10.1002/acs.718. S2CID  62201722.
22. Originalmente atribuido a Stockham en WT Cochran et al. , ¿Qué es la transformada rápida de Fourier?, Proc. IEEE vol. 55, 1664–1674 (1967).
23. PN Swarztrauber, Algoritmos FFT para computadoras vectoriales, Parallel Computing vol. 1, 45–63 (1984).
24. Swarztrauber, PN (1982). "Vectorización de las FFT". En Rodrigue, G. (ed.). Cálculos paralelos . Nueva York: Academic Press. págs. 51–83. ISBN 978-0-12-592101-5.
25. Pease, MC (1968). "Una adaptación de la transformada rápida de Fourier para el procesamiento paralelo". J. ACM . 15 (2): 252–264. doi : 10.1145/321450.321457 . S2CID  14610645.
26. Frigo, Mateo; Johnson, Steven G. "FFTW".Una biblioteca C gratuita ( GPL ) para calcular transformadas de Fourier discretas en una o más dimensiones, de tamaño arbitrario, utilizando el algoritmo Cooley-Tukey
27. Johnson, HW; Burrus, CS (1984). "Una FFT de raíz 2 en orden y en el lugar". Proc. ICASSP : 28A.2.1–28A.2.4.
28. Temperton, C. (1991). "Transformada rápida de Fourier autoordenada en el lugar". Revista SIAM sobre computación científica y estadística . 12 (4): 808–823. doi :10.1137/0912043.
29. Qian, Z.; Lu, C.; An, M.; Tolimieri, R. (1994). "Algoritmo FFT in situ de autoordenación con espacio de trabajo mínimo". IEEE Trans. ASSP . 52 (10): 2835–2836. Bibcode :1994ITSP...42.2835Q. doi :10.1109/78.324749.
30. Hegland, M. (1994). "Un algoritmo de transformada rápida de Fourier autoordenable en el lugar adecuado para el procesamiento vectorial y paralelo". Numerische Mathematik . 68 (4): 507–547. CiteSeerX 10.1.1.54.5659 . doi :10.1007/s002110050074. S2CID  121258187. 

Enlaces externos

• "Transformada rápida de Fourier - FFT". Técnica de Cooley-Tukey . Artículo. 10. Un algoritmo de base 2 simple y pedagógico en C++
• "KISSFFT". GitHub . 11 de febrero de 2022. Una implementación sencilla de Cooley–Tukey de base mixta en C
• Dsplib en GitHub
• "Algoritmo de FFT de diezmado en tiempo de base 2". Archivado desde el original el 31 de octubre de 2017. "Алгоритм БПФ по основанию два с прореживанием по времени" (en ruso).
• "Algoritmo de FFT de frecuencia de decimación de base 2". Archivado desde el original el 14 de noviembre de 2017. "Алгоритм БПФ по основанию два с прореживанием по частоте" (en ruso).