Algoritmo FFT de Cooley-Tukey

El Cooley-Tukey algoritmo, nombrado por J. W. Cooley y John Tukey, es el algoritmo de transformación más rápido Fourier (FFT). Reexpresa la discreta transformación Fourier (DFT) de un tamaño compuesto arbitrario N=N1N2{displaystyle N=N_{1}N_{2} en términos de N₁ DFTs más pequeños de tamaños N₂, recursivamente, para reducir el tiempo de cálculo a O(N log N) para altamente composite N (números suaves). Debido a la importancia del algoritmo, las variantes específicas y los estilos de implementación se han convertido en conocidos 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 lograr una mayor eficiencia a la hora de 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 Neo- Latín). 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. 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 el algoritmo y describiendo cómo ejecutarlo convenientemente en una computadora.

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. 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. Su artículo citó como inspiración sólo el trabajo de I. J. Good sobre lo que ahora se llama algoritmo FFT de factor primo (PFA); 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 sólo para tamaños que tienen factores relativamente primos y se basa en el teorema del resto chino, a diferencia de el soporte para cualquier tamaño compuesto en Cooley-Tukey).

El caso DIT radix-2

Una FFT radix-2 de diezmado en el tiempo (DIT) es la forma más simple y común del algoritmo Cooley-Tukey, aunque Cooley-Tukey está altamente optimizado. Las implementaciones 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 una etapa recursiva.

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

Xk=.. n=0N− − 1xne− − 2π π iNnk,{displaystyle X_{k}=sum ¿Por qué? {2pi i}{N}nk}

Donde k{displaystyle k} es un entero que va de 0 a 0 N− − 1{displaystyle N-1}.

Radix-2 DIT calcula primero los DFTs de las entradas incluso indexadas ()x2m=x0,x2,...... ,xN− − 2){displaystyle (x_{2m}=x_{0},x_{2},ldotsx_{N-2}}y de las entradas de índice impar ()x2m+1=x1,x3,...... ,xN− − 1){displaystyle (x_{2m+1}=x_{1},x_{3},ldotsx_{N-1}}, y luego combina esos dos resultados para producir el DFT de toda la secuencia. Esta idea se puede realizar recursivamente para reducir el tiempo de ejecución general a O(N log N). Esta forma simplificada supone que N es una potencia de dos; desde el número de puntos de muestra N generalmente se puede elegir libremente por la aplicación (por ejemplo, cambiando la tasa de muestra o la ventana, cero relleno, etcétera), esto a menudo no es una restricción importante.

El algoritmo DIT radix-2 reorganiza el DFT de la función xn{displaystyle x_{n} en dos partes: una suma sobre los índices numerados n=2m{displaystyle No. y una suma sobre los índices impares n=2m+1{displaystyle No.:

Xk=.. m=0N/2− − 1x2me− − 2π π iN()2m)k+.. m=0N/2− − 1x2m+1e− − 2π π iN()2m+1)k{displaystyle {begin{matrix}X_{k} diez= golpesum limits ##{m=0}{N/2-1}x_{2m}e^{-{frac {2pi i}{N}(2m)k}+sum limits ################################################################################################################################################################################################################################################################

Uno puede tener un multiplicador común e− − 2π π iNk{displaystyle e^{-{frac {2ccH00} - Sí. de la segunda suma, como se muestra en la ecuación siguiente. Es entonces claro que las dos sumas son el DFT de la parte de índice uniforme x2m{displaystyle x_{2m} y el DFT de la parte indexada x2m+1{displaystyle x_{2m+1}} de la función xn{displaystyle x_{n}. Denota el DFT del Eentradas de ven-indexadas x2m{displaystyle x_{2m} por Ek{displaystyle E_{k} y el DFT Odd-indexed inputs x2m+1{displaystyle x_{2m+1}} por Ok{displaystyle O... y obtenemos:

Xk=.. m=0N/2− − 1x2me− − 2π π iN/2mk⏟ ⏟ DFTofeven− − indexedpartofxn+e− − 2π π iNk.. m=0N/2− − 1x2m+1e− − 2π π iN/2mk⏟ ⏟ DFTofodd− − indexedpartofxn=Ek+e− − 2π π iNkOkparak=0,...... ,N2− − 1.{displaystyle {begin{matrix}X_{k}=underbrace {sum limits ################################################################################################################################################################################################################################################################ Estoy bien. ##{mathrm {DFT;of;even-indexed;part;of;} x_{n}{}{e^{-{frac {2pi i} {N}k}underbrace {sum limits ################################################################################################################################################################################################################################################################ {fn}mk} _{mhrm {DFT;of;odd-indexed;part;of;} ¿Qué? {fnMicrosoft Sans Serif}

Tenga en cuenta que las igualdades se mantienen k=0,...... ,N− − 1{displaystyle k=0,dotsN-1} pero el crux es eso Ek{displaystyle E_{k} y Ok{displaystyle O... se calculan de esta manera para k=0,...... ,N2− − 1{displaystyle k=0,dots{frac {N}{2}-1} sólo. Gracias a la periodicidad del complejo exponencial, Xk+N2{displaystyle X_{k+{frac {N}{2}}} se obtiene también de Ek{displaystyle E_{k} y Ok{displaystyle O...:

Xk+N2=.. m=0N/2− − 1x2me− − 2π π iN/2m()k+N2)+e− − 2π π iN()k+N2).. m=0N/2− − 1x2m+1e− − 2π π iN/2m()k+N2)=.. m=0N/2− − 1x2me− − 2π π iN/2mke− − 2π π mi+e− − 2π π iNke− − π π i.. m=0N/2− − 1x2m+1e− − 2π π iN/2mke− − 2π π mi=.. m=0N/2− − 1x2me− − 2π π iN/2mk− − e− − 2π π iNk.. m=0N/2− − 1x2m+1e− − 2π π iN/2mk=Ek− − e− − 2π π iNkOk{displaystyle {begin{aligned}X_{k+{frac} {N}{2}} {=sum limits ################################################################################################################################################################################################################################################################ ¿Qué? {2ccH00} {fn} {fn}}sum limits ################################################################################################################################################################################################################################################################ i} {N/2}m(k+{frac {N} {2})}\=sum limits ################################################################################################################################################################################################################################################################ Estoy bien. mi}+e^{-{frac {2pi} ¿Qué? i}sum limits ################################################################################################################################################################################################################################################################ Estoy bien. mi}\=sum limits ################################################################################################################################################################################################################################################################ Estoy bien. {N}k}sum limits ################################################################################################################################################################################################################################################################ Estoy bien. {fn}K}O_{k}end{aligned}}

Podemos reescribir Xk{displaystyle X_{k} y Xk+N2{displaystyle X_{k+{frac {N}{2}}} como:

Xk=Ek+e− − 2π π iNkOkXk+N2=Ek− − e− − 2π π iNkOk{displaystyle {begin{matrix}X_{k} {c}c}+e^{-{frac} {2ccH00} - Sí. {N}{2}} âTMa âTMa âTMa âTMa âTMa âTMa {2pi} {fn} {fn} {fn}} {fnK}}}} {fn}}} {fn}}} {fn}}} {fn}}}}}} {fn}} {fn}}}} {fn}}}} {f} {fn}}} {f}}}}f}}}}}}}}}}}}}}f}}}}}}}}}}f}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}

Este resultado, expresando el DFT de longitud N recursivamente en términos de dos DFT de tamaño N/2, es el núcleo del radio-2 DIT rápido Fourier transformado. El algoritmo gana su velocidad reutilizando los resultados de los cálculos intermedios para calcular múltiples salidas DFT. Observe que los productos finales se obtienen por una combinación +/− Ek{displaystyle E_{k} y Okexp⁡ ⁡ ()− − 2π π ik/N){displaystyle O_{k}exp(-2pi ik/N)}, que es simplemente un tamaño-2 DFT (a veces llamado mariposa en este contexto); cuando esto se generaliza a mayores radios abajo, el tamaño-2 DFT es reemplazado por un DFT más grande (que en sí puede ser evaluado con un FFT).

Este proceso es un ejemplo de la técnica general de 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 un DFT de tamaño N como dos DFT de tamaño N/2 a veces se denomina Danielson-Lanczos. lema, ya que la identidad fue notada por esos dos autores en 1942 (influenciados por el trabajo de Runge de 1903). Aplicaron su lema de forma "al revés" de manera recursiva, duplicando repetidamente el tamaño de DFT hasta que el espectro de transformación convergiera (aunque aparentemente no se dieron cuenta del orden linealítmico [es decir, N log N] complejidad asintótica 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á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 con precisión simple de 36 bits, ~8 dígitos). Reescalando el tiempo por el número de operaciones, esto corresponde aproximadamente a un factor de aceleración de aproximadamente 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:

X0,...N−1 ← ditfft2()x, N, s): DFT de (x0, xs, x2s,... x()N-1)s):
 si N = 1 entonces X0 ← x0 trivial tamaño-1 caso base DFT más X0,...N/2−1 ← ditfft2()x, N/2, 2s) DFT de (x0, x2s, x4s,... x()N-2)s)
 XN/2,...N−1 ← ditfft2()x+s, N/2, 2s) DFT de (xs, xs+2s, xs+4s,... x()N-1)s)
 para k = 0 a N/2−1 do combinar DFTs de dos mitades en el DFT completo:p ← Xkq ← exp(−2πi/N k) Xk+N/2 Xk ← p + q
 Xk+N/2 ← p - q
 final for terminar si

Aquí, ditfft2(x,N,1), calcula X=DFT(x) fuera de lugar por una FFT DIT radix-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 más permutación de inversión de bits; la necesidad mencionada a menudo de una etapa de inversión de bits separada solo surge para ciertos algoritmos locales, como se describe abajo.)

Alto rendimiento Las implementaciones FFT hacen muchas modificaciones a la implementación de tal algoritmo en comparación con este simple pseudocódigo. Por ejemplo, se puede utilizar un caso base más grande que N=1 para amortizar la parte superior de la recursión, los factores de twiddle exp⁡ ⁡ [− − 2π π ik/N]{displaystyle exp[-2pi ik/N]} puede ser precomputado, y los radios más grandes se utilizan a menudo por razones de caché; estas y otras optimizaciones juntas pueden mejorar el rendimiento por un orden de magnitud o más. (En muchas implementaciones del libro de texto se elimina la primera repetición de profundidad a favor de un enfoque no recursivo de amplitud primera, aunque se ha argumentado que la primera recursión de profundidad tiene una mejor localización de memoria.) Varias de estas ideas se describen con más detalle a continuación.

Idea

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

1. Perform N₁ DFTs of size N₂.
2. Multiply por complejas raíces de la unidad (a menudo llamadas los factores twiddle).
3. Perform N₂ DFTs of size 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 denomina algoritmo de diezmado en el tiempo (DIT), mientras que si N ₂ es la base, es diezmado en frecuencia (DIF, también llamado algoritmo 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. (El pequeño DFT de radix a veces se conoce como mariposa, llamado así debido a la forma del diagrama de flujo de datos para el caso de radix-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, normalmente (pero no siempre) empleando el método 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 Rader&# 39;s o algoritmo de 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, aunque las variaciones recientes logran una cuenta uniforme. conteo menor. (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 radices más grandes y/o transformaciones de caso base codificadas de tamaño significativo).

Otra forma de ver el algoritmo de Cooley-Tukey es que reexpresa una DFT unidimensional de tamaño N como un N₁ por N₂ DFT bidimensional (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 seis pasos, dependiendo del número de transposiciones), propuesto inicialmente para mejorar la localidad de la memoria, p.e. para la optimización de la caché o la operación fuera del núcleo, y luego se demostró que era un algoritmo óptimo que no tenía en cuenta la caché.

La factorización general Cooley–Tukey reescribe los índices k y n como k=N2k1+k2{displaystyle K=N_{2}k_{1}+k_{2} y n=N1n2+n1{displaystyle No., respectivamente, donde los índices k_a y n_a huye de 0..N_a-1 (para a de 1 o 2). Es decir, vuelve a indexar la entrada (n) y salida (kcomo N₁ por N₂ arrays bidimensionales en orden columna-major y fila-major, respectivamente; la diferencia entre estos índices es una transposición, como se mencionó anteriormente. Cuando esta re-indexación se sustituye a la fórmula DFT nk, el N1n2N2k1{displaystyle N_{1}n_{2}N_{2}k_{1} el término cruzado desaparece (su exponencial es unidad), y los términos restantes dan

XN2k1+k2=.. n1=0N1− − 1.. n2=0N2− − 1xN1n2+n1e− − 2π π iN1N2⋅ ⋅ ()N1n2+n1)⋅ ⋅ ()N2k1+k2){displaystyle X_{N_{2}k_{1}+k_{2}=sum - ¿Por qué? - ¿Por qué? {2ccH00} i}{N_{2}}cdot (N_{1}n_{2}+n_{1})cdot (N_{2}k_{1}+k_{2}}}}

=.. n1=0N1− − 1[e− − 2π π iN1N2n1k2]().. n2=0N2− − 1xN1n2+n1e− − 2π π iN2n2k2)e− − 2π π iN1n1k1{displaystyle =sum ¿Por qué? {fnfn} {fn}}n_fn}fnfnh}nh}nhnn_ccH3}ccH00}cH00}cH0}cHFF}pcH00}ccH0}cH00}cH0}cH0}cH0}cH0}cH0}cH0}ccH0}ccH00}cH00}cH0}cH00}cH0}ccH0}ccH0}ccccccccccH0}ccH0}ccH00}ccH00}c}c}c}cH00}}ccc}cc}cc}c}ccc}c}c ################################################################################################################################################################################################################################################################ {fnfn} {fnfnK} {fnMicroc {2fnMicroc {2pi} {fn} {fn} {fn}} {fn}}} {fn}} {fn} {fn}} {fn}}}}} {fn}}}}} {fn}}}}}}}} {fn}} {f}}}} {fn}}}}}}}}}} {\\}}}}}}}}}}}}} {\\\\\\n}}}}}}}}}} {\\\\\\\\\\\\\\\\\\}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}}

=.. n1=0N1− − 1().. n2=0N2− − 1xN1n2+n1e− − 2π π iN2n2k2)e− − 2π π iN1N2n1()N2k1+k2){displaystyle =sum - ¿Por qué? ################################################################################################################################################################################################################################################################ {fnfn} {fnfnK} {fnMicroc {2fnMicroc {2pi} ¿Qué?.

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 arbitraria r (así como bases mixtas), como lo demostraron Cooley y Tukey, así como Gauss (quien dio ejemplos de pasos de base 3 y base 6).). Cooley y Tukey asumieron originalmente 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_rN) = 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, la base óptima se encuentra que es 3 (el entero más cercano a e, que minimiza r/log₂r). 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), de ahí la base r en realidad se cancela en la complejidad O(r log(r) N/r log_rN), y el r óptimo está determinado por consideraciones más complicadas. En la práctica, un r bastante grande (32 o 64) es importante para explotar eficazmente, p. la gran cantidad de registros de procesador en los procesadores modernos, 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 grandes N como se mencionó anteriormente.

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 explícitamente bit reversal para algoritmos radix-2 en el lugar. Reversión de bits es la permutación donde los datos en un índice n, escrito en binario con dígitos b₄b₃b₂b₁b₀ (por ejemplo, 5 dígitos para N=32 entradas), se transfiere al índice con dígitos inversos b₀b₁b₂b₃b₄. Considere la última etapa de un algoritmo de radio-2 DIT como el presentado anteriormente, donde la salida está escrita en el lugar sobre la entrada: cuando Ek{displaystyle E_{k} y Ok{displaystyle O... se combinan con un tamaño-2 DFT, esos dos valores son sobrescritos por los productos. Sin embargo, los dos valores de salida deben ir en el primero y segundo mitades de la matriz de salida, correspondiente al más bit significativa b₄ (por N=32); mientras que los dos insumos Ek{displaystyle E_{k} y Ok{displaystyle O... están entrelazados en los elementos uniformes y extraños, correspondientes a los mínimo bit significativa b₀. Así, para obtener la salida en el lugar correcto, b₀ debe ocupar el lugar b₄ y el índice se convierte b₀b₄b₃b₂b₁. Y para la siguiente etapa recurrente, esos 4 bits menos significativos se convertirán en b₁b₄b₃b₂, Si usted incluye todas las etapas recursivas de un algoritmo de radio-2 DIT, Todos los bits deben ser revertidos y por lo tanto uno debe preprocesar la entrada (o post-procesar la salida) con un poco de inversión para obtener la salida en-orden. (Si cada tamaño-N/2 subtransform es operar en datos contiguos, el DIT entrada es pre-procesado por bit-reversal.) Correspondientemente, si usted realiza todos los pasos en orden inverso, usted obtiene un algoritmo de radiox-2 DIF con reversión de bits en post-procesamiento (o pre-procesamiento, 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.

algoritmo iterative-fft es entrada: Array a de n valores complejos donde n es un poder de 2.
 Producto: Array A el DFT de a.

bit-reverse-copy(a, A)
 n ← a.length
 para s = 1 a log(n) do m ← 2s ⋅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← u – t ⋅ ← ⋅ ⋅m retorno A

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

algoritmo bit-reverse-copy(a,A) es entrada: Array a de n valores complejos donde n es un poder de 2.
 Producto: Array A de tamaño n.

 n ← a.length
 para k = 0 a n – 1 do 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.

Sin embargo, muchos usuarios de FFT 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, aunque la inversión de bits se puede realizar en O(N ) tiempo y ha sido objeto de mucha investigación. Además, si bien la permutación es una inversión de bits en el caso de base mixta, en general es 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, 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, equivalentemente, hay disponible una matriz auxiliar del mismo tamaño. La clasificación automática de Stockham El algoritmo 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 en arquitecturas SIMD. Se han propuesto ventajas potenciales aún mayores de SIMD (más accesos consecutivos) para el algoritmo Pease, que también reordena fuera de lugar con cada etapa, pero este método requiere una inversión separada de bits/dígitos y O(N log N) 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 tener beneficios de localidad sin memoria caché en sistemas con memoria jerárquica.

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 de matrices (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 pasa por encima de los datos.