Algoritmo FFT de Cooley-Tukey

Algoritmo de transformada rápida de Fourier

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 ₁ 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, variantes específicas y estilos de implementación se conocen 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 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.

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 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]

El caso del DIT radix-2

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:

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

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 : ${\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 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). ${\displaystyle E_{k}}$ ${\displaystyle O_{k}\exp(-2\pi ik/N)}$

Diagrama de flujo de datos para N =8: una FFT de base 2 de diezmado en el tiempo divide una DFT de longitud N en dos DFT de N /2 de longitud seguida de una etapa de combinación que consta de DFT de N /2 tamaño 2 llamada "mariposa" operaciones (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 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).

Pseudocódigo

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. ${\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 1d como algo así como una DFT 2d. La matriz de entrada 1d de longitud N = N ₁ N ₂ se reinterpreta como una matriz 2d N ₁ × N ₂ almacenada en orden de columna principal . Se realizan DFT 1d 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 giro) y finalmente se realizan DFT 1d a lo largo de la dirección N ₁ . 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. Realice N ₁ DFT de tamaño N ₂ .
2. Multiplicar por raíces complejas de la unidad (a menudo llamadas factores de giro ).
3. Realice N ₂ DFT de tamaño N ₁ .

Normalmente, N ₁ o N ₂ es un factor pequeño ( no necesariamente primo), llamado base (que puede diferir entre las etapas de la recursividad). Si N ₁ es la base, se llama algoritmo de diezmado en el tiempo (DIT), mientras que si N ₂ 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).

Variaciones

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 ² ) 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 casos 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 ₁ 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 ₁ 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 para nk , el término cruzado desaparece (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 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 base de la mariposa requería trabajo O ( r ² ) y, por lo tanto, calcularon que la complejidad de una base r era O ( r ²  N / r  log _r N ) = O ( N  log ₂ ( N )  r / _log2r ) ; a partir del cálculo de los valores de r /log ₂ 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 ₂ 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]}

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

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 ₄ b ₃ b ₂ b ₁ b ₀ (por ejemplo, 5 dígitos para N = 32 entradas), se transfieren al índice con dígitos invertidos b ₀ b ₁ segundo ₂ segundo ₃ segundo ₄ . 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 ₄ 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 ₀ menos significativo . 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 _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). ${\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 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 SIMD potenciales aún mayores (más accesos consecutivos) para elde Pease,^[25]que también reordena fuera de lugar con cada etapa, pero este método requiere inversión de bits/dígitos por separado 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]}

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", revista IEEE ASSP, 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. Computadora. 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) . Transacciones IEEE sobre audio y electroacústica . 15 (2): 76–79. CiteSeerX 10.1.1.467.7209 . doi :10.1109/tau.1967.1161903. 
5. Rockmore, Daniel N., Computadora. Ciencia. Ing. 2 (1), 60 (2000). La FFT: un algoritmo que toda la familia puede utilizar Número especial sobre los "diez mejores algoritmos del siglo" Barry A. Cipra. "Lo mejor del siglo XX: los editores nombran los 10 algoritmos principales" (PDF) . Noticias SIAM . 33 (4). Archivado desde el original (PDF) el 7 de abril de 2009 . Consultado el 31 de marzo de 2009 .
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úm. 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 Transformadas rápidas de Fourier (CS Burrus, ed.), cap. 11, Universidad Rice, Houston TX: Connexions, septiembre de 2008.
9. Singleton, Richard C. (1967). "Sobre el cálculo de la transformada rápida de Fourier". Comunitario. 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 del 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. Proceso de señal. 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. Código Bib : 2005IEEEP..93..216F. CiteSeerX 10.1.1.66.3097 . doi :10.1109/JPROC.2004.840301. S2CID  6644892. 
14. Gentleman WM y G. Sande, "Transformaciones rápidas de Fourier, por diversión y ganancias", 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 ajenos al caché. En Actas del 40º Simposio IEEE sobre fundamentos de la informática (FOCS 99), p.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 monoprocesadores". Revisión SIAM . 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 (n.º de catálogo 98CB36280) . págs. 544–553. CiteSeerX 10.1.1.46.9319 . doi :10.1109/SFCS.1998.743505. ISBN  978-0-8186-9172-0. S2CID  14307262.
21. Rubio, M.; Gómez, P.; Drouiche, K. (2002). "Un nuevo algoritmo de inversión de bits ultrarrápido". Revista internacional de control adaptativo y procesamiento de señales . 16 (10): 703–707. doi :10.1002/acs.718. S2CID  62201722.
22. Atribuido originalmente 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. Por favor, MC (1968). "Una adaptación de la transformada rápida de Fourier para 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 radix-2 in situ y en orden". Proc. ICASSP : 28A.2.1–28A.2.4.
28. Temperton, C. (1991). "Transformación rápida de Fourier autoclasificada in situ". Revista SIAM de 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 de autoclasificación in situ con espacio de trabajo mínimo". Traducción IEEE. ASSP . 52 (10): 2835–2836. Código Bib : 1994ITSP...42.2835Q. doi : 10.1109/78.324749.
30. Hegland, M. (1994). "Un algoritmo de transformada rápida de Fourier autoclasificable in situ adecuado para procesamiento vectorial y paralelo". Matemática numérica . 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 Cooley-Tukey . Artículo. 10. Un algoritmo radix-2 simple y pedagógico en C++
• "BESO". GitHub . 11 de febrero de 2022. Una implementación simple de Cooley-Tukey de base mixta en C
• Dsplib en GitHub
• "Algoritmo FFT de diezmado en el tiempo Radix-2". Archivado desde el original el 31 de octubre de 2017. "Алгоритм БПФ по основанию два с прореживанием по времени" (en ruso).
• "Algoritmo FFT de diezmado Radix-2 en frecuencia". Archivado desde el original el 14 de noviembre de 2017. "Алгоритм БПФ по основанию два с прореживанием по частоте" (en ruso).