El algoritmo Cooley-Tukey , que lleva el nombre de 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 1 DFT más pequeñas de tamaños N 2 , 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, variantes específicas y estilos de implementación se conocen por sus propios nombres, como se describe a continuación.
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 Cooley-Tukey no puede descomponer, o se puede explotar el algoritmo de factores primos para obtener una mayor eficiencia al separar factores relativamente primos .
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.
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 Palas y Juno , pero su trabajo no fue ampliamente reconocido (sólo se publicó póstumamente y en neolatino ). [1] [2] Sin embargo, Gauss no analizó el tiempo computacional asintótico. También se redescubrieron varias formas limitadas 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 convenientemente 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 mediante el empleo de 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 propuesta de prohibición de los ensayos nucleares, 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, a Cooley le dijeron que esto era necesario para determinar las periodicidades de las orientaciones de espín en un cristal tridimensional de helio-3 . Posteriormente, Cooley y Tukey publicaron 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 hubiera 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 citó como inspiración sólo el trabajo de IJ Good sobre lo que ahora se llama algoritmo FFT de factor primo (PFA); [3] Aunque inicialmente se pensó que el algoritmo de Good era equivalente al algoritmo de Cooley-Tukey, rápidamente se dio cuenta de que PFA es un algoritmo bastante diferente (que funciona solo para tamaños que tienen factores relativamente primos y se basa en el teorema del resto chino , a diferencia del soporte para cualquier tamaño compuesto en Cooley-Tukey). [6]
Una FFT de diezmación en el tiempo ( DIT ) de radix-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. Radix-2 DIT divide una DFT de tamaño N en dos DFT entrelazadas (de ahí el nombre "radix-2") de tamaño N /2 con cada etapa recursiva.
La transformada discreta de Fourier (DFT) se define mediante la fórmula:
donde es un número entero que va de 0 a .
Radix-2 DIT primero calcula las DFT de las entradas con índice par y de las entradas con índice impar , y luego combina esos dos resultados para producir la DFT de toda la secuencia. Luego, esta idea se puede realizar de forma recursiva para reducir el tiempo de ejecución general a O ( N log N ). Esta forma simplificada supone que N es una potencia de dos ; Dado que la aplicación normalmente puede elegir libremente el número de puntos de muestreo N (por ejemplo, cambiando la frecuencia o ventana de muestreo, relleno con ceros, etc.), esto no suele ser una restricción importante.
El algoritmo DIT radix-2 reorganiza la DFT de la función en dos partes: una suma de los índices pares y una suma de los índices impares :
Se puede factorizar un multiplicador común a partir de la segunda suma, como se muestra en la siguiente ecuación. Entonces está claro que las dos sumas son la DFT de la parte con índice par y la DFT de la parte con índice impar de la función . Denotamos la DFT de las entradas indexadas pares por y la DFT de las entradas indexadas impares por y obtenemos:
Tenga en cuenta que las igualdades son válidas 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 de y :
Podemos reescribir y como:
Este resultado, que expresa la DFT de longitud N de forma recursiva en términos de dos DFT de tamaño N /2, es el núcleo de la transformada rápida de Fourier radix-2 DIT. El algoritmo gana velocidad al reutilizar los resultados de cálculos intermedios para calcular múltiples salidas DFT. Tenga en cuenta que los resultados finales se obtienen mediante una combinación +/− de y , que es simplemente un DFT de tamaño 2 (a veces llamado mariposa en este contexto); cuando esto se generaliza a raíces más grandes 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).
Este proceso es un ejemplo de la técnica general de los algoritmos de divide y vencerás ; Sin embargo, en muchas implementaciones convencionales se evita la recursividad 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] (influenciado por el trabajo de Runge de 1903 [2 ] ). Aplicaron su lema de forma recursiva "hacia atrás", duplicando repetidamente el tamaño de la DFT hasta que el espectro de transformación convergió (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 es anterior a la disponibilidad generalizada de computadoras mecánicas o electrónicas y requería cálculos manuales (posiblemente con ayudas mecánicas como máquinas sumadoras ); informaron un tiempo de cálculo de 140 minutos para un DFT de tamaño 64 que operaba con entradas reales de 3 a 5 dígitos significativos. El artículo de Cooley y Tukey de 1965 informó un tiempo de ejecución de 0,02 minutos para un DFT complejo 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 en perspectiva el tiempo de cálculo manual, 140 minutos para el tamaño 64 corresponden a una media de como máximo 16 segundos por operación de punto flotante, de los cuales alrededor del 20% son multiplicaciones).
En pseudocódigo , se podría escribir el siguiente procedimiento: [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 tamaño trivial-1 caso base DFT else 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 ) EPS de (x s , x s +2 s , x s +4 s , ..., x ( N -1) s ) para k = 0 a N /2−1 combine las DFT de dos mitades en 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 por 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 x de entrada . 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 permutación de inversión de bits adicional ; la necesidad mencionada a menudo de una etapa de inversión de bits separada solo surge para ciertos algoritmos locales, como se describe a continuación).
Las implementaciones de FFT de alto rendimiento realizan muchas modificaciones en la implementación de dicho algoritmo en comparación con este pseudocódigo simple. Por ejemplo, se puede usar un caso base mayor que N = 1 para amortizar la sobrecarga de la recursividad, los factores de giro se pueden calcular previamente y, a menudo, se usan radices 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 recursividad en profundidad se elimina en favor de un enfoque no recursivo en amplitud , aunque se ha argumentado que la recursividad en profundidad tiene una mejor localidad de memoria . [8] [9] ) Varias de estas ideas se describen con más detalle a continuación.
De manera más general, los algoritmos de Cooley-Tukey reexpresan recursivamente una DFT de un tamaño compuesto N = N 1 N 2 como: [10]
Normalmente, N 1 o N 2 es un factor pequeño ( no necesariamente primo), llamado base (que puede diferir entre las etapas de la recursividad). Si N 1 es la base, se llama algoritmo de diezmado en el tiempo (DIT), mientras que si N 2 es la base, es diezmado en frecuencia (DIF, también llamado algoritmo de Sande-Tukey). La versión presentada anteriormente era un algoritmo DIT radix-2; en la expresión final, la fase que multiplica la transformada impar es el factor de giro, y la combinación +/- ( mariposa ) de las transformadas pares e impares es una DFT de tamaño 2. (La pequeña DFT de la raíz se conoce a veces como mariposa , llamada así debido a la forma del diagrama de flujo de datos para el caso de la raíz-2).
Existen muchas otras variaciones del algoritmo de Cooley-Tukey. Las implementaciones de base mixta manejan tamaños compuestos con una variedad de factores (típicamente pequeños) además de dos, empleando generalmente (pero no siempre) el algoritmo O( N 2 ) para los casos base primos de la recursividad (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 raíces 2 y 4, aprovechando el hecho de que la primera transformación de la base 2 no requiere factor de giro, para lograr lo que durante mucho tiempo fue el recuento de operaciones aritméticas más bajo conocido para potencias de dos tamaños, [10] aunque las variaciones recientes lograr un recuento aún menor. [11] [12] (En las computadoras actuales, el rendimiento está determinado más por consideraciones de canalización de CPU y caché que por recuentos de operaciones estrictos; las implementaciones de FFT bien optimizadas a menudo emplean radices más grandes y/o transformaciones de caso base codificadas de valores significativos). tamaño [13] ).
Otra forma de ver el algoritmo Cooley-Tukey es que reexpresa una DFT unidimensional de tamaño N como una DFT bidimensional N 1 por N 2 (más giros), donde se transpone la matriz de salida . El resultado neto de todas estas transposiciones, para un algoritmo radix-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 matrices de entrada/salida, se denomina algoritmo FFT de cuatro pasos (o de seis pasos , dependiendo del número de transposiciones), propuesto inicialmente. 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 es un algoritmo óptimo que no tiene en cuenta 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 desde 0.. N a -1 (para a de 1 o 2). Es decir, vuelve a indexar la entrada ( n ) y la salida ( k ) como matrices bidimensionales N 1 por N 2 en orden de columna principal y fila principal , 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 por nk , el término cruzado desaparece (su exponencial es la unidad) y los términos restantes dan
donde cada suma interna es una DFT de tamaño N 2 , cada suma externa es una DFT de tamaño N 1 y el término [...] entre corchetes es el factor de giro.
Se puede emplear una base r arbitraria (así como bases mixtas), como lo demostraron Cooley y Tukey [3], así como Gauss (quien dio ejemplos de pasos de base 3 y base 6). [2] Cooley y Tukey originalmente asumieron que la mariposa de base requería trabajo O ( r 2 ) y, por lo tanto, calcularon que la complejidad de una base r era O ( r 2 N / r log r N ) = O ( N log 2 ( N ) r / log2r ) ; a partir del cálculo de los valores de r /log 2 r para valores enteros de r de 2 a 12, se encuentra que la base óptima es 3 (el número entero más cercano a e , que minimiza r /log 2 r ). [3] [17] Sin embargo, este análisis fue erróneo: la raíz-mariposa también es una DFT y se puede realizar mediante 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 el r óptimo está determinado por consideraciones más complicadas. En la práctica, r bastante grandes (32 o 64) son importantes para explotar eficazmente, 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]
Aunque la factorización abstracta de Cooley-Tukey de la DFT anterior 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 sólo almacenamiento auxiliar O(1).
La técnica de reordenamiento más conocida implica la inversión explícita de bits para algoritmos radix-2 in situ. La inversión de bits es la permutación donde los datos en un índice n , escritos en binario con dígitos b 4 b 3 b 2 b 1 b 0 (por ejemplo, 5 dígitos para N = 32 entradas), se transfieren al índice con dígitos invertidos b 0 b 1 segundo 2 segundo 3 segundo 4 . Considere la última etapa de un algoritmo DIT radix-2 como el presentado anteriormente, donde la salida se escribe en el lugar sobre la entrada: cuando y se combinan con un DFT de tamaño 2, las salidas sobrescriben esos dos valores. Sin embargo, los dos valores de salida deben ir en la primera y segunda mitad de la matriz de salida, correspondiente al bit b 4 más significativo (para N =32); mientras que las dos entradas y están entrelazadas en los elementos pares e impares, correspondientes al bit b 0 menos significativo . Por lo tanto, para obtener la salida en el lugar correcto, b 0 debe tomar el lugar de b 4 y el índice se convierte en b 0 b 4 b 3 b 2 b 1 . Y para la siguiente etapa recursiva, esos 4 bits menos significativos se convertirán en b 1 b 4 b 3 b 2. Si incluye todas las etapas recursivas de un algoritmo DIT radix-2, todos los bits deben invertirse y, por lo tanto, se deben pre- procese la entrada ( o posprocese la salida) con una inversión de bits para obtener una salida en orden. (Si cada subtransformada de tamaño N /2 va a operar con datos contiguos, la entrada DIT se preprocesa mediante inversión de bits). En consecuencia, si realiza todos los pasos en orden inverso, obtendrá un algoritmo DIF de base 2. con inversión de bits en el posprocesamiento (o preprocesamiento, respectivamente).
El logaritmo (log) utilizado en este algoritmo es un logaritmo de base 2.
El siguiente es un pseudocódigo para el algoritmo iterativo FFT radix-2 implementado mediante permutación de inversión de bits. [18]
algoritmo iterativo-fft es entrada: Matriz a de n valores complejos donde n es una potencia de 2. salida: Matriz A la DFT de a. copia inversa de bits (a, A) n ← a .length para s = 1 para log( n ) do m ← 2 s ω m ← exp(−2π i / m ) para k = 0 para n -1 por m do ω ← 1 para j = 0 para m / 2 – 1 punto t ← ω A [ k + j + m /2] u ← A [ k + j ] A [ k + j ] ← u + t A [ k + j + m /2] ← u – t ω ← ω ω metro devolver un
El procedimiento de copia inversa de bits se puede implementar de la siguiente manera.
El algoritmo de copia inversa de bits ( a , A ) se ingresa: Matriz a de n valores complejos donde n es una potencia de 2. Salida: Matriz A de tamaño n . n ← a .length para k = 0 an – 1 hacer A [ rev(k)] := a [k]
Alternativamente, algunas aplicaciones (como la convolución) funcionan igualmente bien con datos de bits invertidos, por lo que se pueden realizar transformaciones directas, procesamientos y luego transformaciones inversas, todo sin inversión de bits para producir resultados finales en el orden natural.
Muchos usuarios de FFT, sin embargo, prefieren salidas en 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] a pesar de que 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, si bien la permutación es una inversión de bits en el caso de base mixta, 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 convierten en más complicado de implementar. Además, en muchas arquitecturas de hardware es deseable 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 inversión de bits por separado y/o implican 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, de manera equivalente, hay disponible una matriz auxiliar del mismo tamaño. ElEl algoritmo de ordenación automática de Stockham [22][23]realiza cada etapa de la FFT fuera de lugar, normalmente escribiendo de un lado a otro 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 potenciales aún mayores de SIMD (más accesos consecutivos) para elalgoritmode Pease[25]que también reordena fuera de lugar con cada etapa, pero este método requiere inversión separada de bits/dígitos y O (NlogN) almacenamiento. También se puede aplicar directamente la definición de factorización de Cooley-Tukey con recursividad explícita (primero en profundidad) y raíces pequeñas, lo 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 tenersin memoria cachéen sistemas conmemoria jerárquica.[9][13][26]
Una estrategia típica para algoritmos in situ sin almacenamiento auxiliar y sin pases separados de inversión de dígitos implica pequeñas transposiciones matriciales (que intercambian pares individuales de dígitos) en etapas intermedias, que se pueden combinar con las mariposas de base para reducir el número de pases sobre el datos. [13] [27] [28] [29] [30]
Un algoritmo radix-2 simple y pedagógico en C++
Una implementación simple de Cooley-Tukey de base mixta en C