Algoritmo FFT de Cooley–Tukey

Esse algoritmo, incluindo sua aplicação recursiva, foi inventado por volta de 1805 por Carl Friedrich Gauss, que o usou para interpolar as trajetórias dos asteroides Pallas e Juno, mas seu trabalho não foi amplamente reconhecido (sendo publicado apenas postumamente e em neolatim). Gauss não analisou o tempo computacional assintótico, no entanto. Várias formas limitadas também foram redescobertas diversas vezes ao longo do século XIX e início do século XX. FFTs se tornaram populares depois que James Cooley da IBM e John Tukey de Princeton publicaram um artigo em 1965 reinventando o algoritmo e descrevendo como executá-lo convenientemente em um computador. Tukey supostamente teve a ideia durante uma reunião do Comitê Consultivo Científico do Presidente Kennedy discutindo maneiras de detectar testes de armas nucleares na União Soviética empregando sismômetros localizados fora do país. Esses sensores gerariam séries temporais sismológicas. No entanto, a análise desses dados exigiria algoritmos rápidos para calcular DFTs devido ao número de sensores e à duração das séries temporais. Essa tarefa foi crítica para a ratificação da proposta de proibição de testes nucleares para que quaisquer violações pudessem ser detectadas sem a necessidade de visitar instalações soviéticas. Outro participante daquela reunião, Richard Garwin da IBM, reconheceu o potencial do método e colocou Tukey em contato com Cooley. No entanto, Garwin garantiu que Cooley não soubesse o propósito original. Em vez disso, Cooley foi informado de que isso era necessário para determinar as periodicidades das orientações de spin em um cristal 3-D de hélio-3. Cooley e Tukey publicaram posteriormente seu artigo conjunto, e a ampla adoção ocorreu rapidamente devido ao desenvolvimento simultâneo de conversores analógico-digitais capazes de amostrar em taxas de até 300 kHz.

Uma FFT de dizimação no tempo (DIT) de base 2 é a forma mais simples e comum do algoritmo de Cooley-Tukey, embora implementações altamente otimizadas de Cooley-Tukey normalmente utilizem outras formas do algoritmo, conforme descrito abaixo. A DIT de base 2 divide uma DFT de tamanho N em duas DFTs intercaladas (daí o nome "base 2") de tamanho N/2 com cada estágio recursivo. A transformada discreta de Fourier (DFT) é definida pela fórmula: onde n {\displaystyle n} é um número inteiro que varia de 0 a N − 1 {\displaystyle N-1} . A DIT de base 2 primeiro calcula os DFTs das entradas de índices pares ( x 2 m = x 0 , x 2 , … , x N − 2 ) {\displaystyle (x_{2m}=x_{0},x_{2},\ldots ,x_{N-2})} e das entradas de índices ímpares ( x 2 m + 1 = x 1 , x 3 , … , x N − 1 ) {\displaystyle (x_{2m+1}=x_{1},x_{3},\ldots ,x_{N-1})} e, em seguida, combina esses dois resultados para produzir a DFT de toda a sequência. Essa ideia pode ser aplicada recursivamente para reduzir o tempo de execução total para O(N log N). Essa forma simplificada assume que N é uma potência de dois; como o número de pontos de amostra N geralmente pode ser escolhido livremente pela aplicação (por exemplo, ajustando a taxa de amostragem ou a janela, aplicando preenchimento com zeros, etc.), isso geralmente não é uma restrição importante.

De forma mais geral, os algoritmos de Cooley–Tukey reexpressam recursivamente uma DFT de tamanho composto N = N1N 2 como: Normalmente, N1 ou N2 é um fator pequeno (não necessariamente primo), chamado de radix (que pode diferir entre os estágios da recursão). Se N1 for o radix, é chamado de algoritmo de dizimação no tempo (DIT), enquanto se N2 for o radix, é de dizimação na frequência (DIF, também chamado de algoritmo de Sande–Tukey). A versão apresentada acima era um algoritmo radix-2 DIT (DIT de base 2); na expressão final, a fase que multiplica a transformada ímpar é o fator de rotação (twiddle factors), e a combinação +/- (borboleta) das transformadas pares e ímpares é uma DFT de tamanho 2. (A pequena DFT do radix às vezes é conhecida como "borboleta", assim chamada por causa do formato do diagrama de fluxo de dados para o caso radix-2.)

Existem muitas outras variações do algoritmo de Cooley-Tukey. Implementações de base mista (mixed-radix) lidam com tamanhos compostos com uma variedade de fatores (tipicamente pequenos) além de dois, geralmente (mas nem sempre) empregando o algoritmo O(N²) para os casos base primos da recursão (também é possível empregar um algoritmo N log N para os casos base primos, como o algoritmo de Rader ou o algoritmo de Bluestein). A abordagem split radix combina as bases 2 e 4, explorando o fato de que a primeira transformação da base 2 não requer fator de rotação (twiddle factor), a fim de alcançar o que foi por muito tempo a menor contagem de operações aritméticas conhecida para tamanhos de potência de dois, embora variações recentes alcancem uma contagem ainda menor. (Em computadores atuais, o desempenho é determinado mais por considerações de cache e pipeline da CPU do que por contagens de operações estritas; implementações de FFT bem otimizadas frequentemente empregam bases maiores e/ou transformadas de caso base codificadas de tamanho considerável).

Embora a fatoração do tipo Cooley-Tukey da DFT, acima, se aplique de alguma forma a todas as implementações do algoritmo, existe uma diversidade muito maior nas técnicas para ordenar e acessar os dados em cada estágio da FFT. De especial interesse é o problema de construir um algoritmo in-place que sobrescreva sua entrada com seus dados de saída usando apenas armazenamento auxiliar de O(1). A técnica de reordenação mais conhecida envolve a inversão explícita de bits para algoritmos radix-2 in-place. A inversão de bits é a permutação na qual os dados em um índice n, escritos em binário com os digitos b₄b₃b₂b₁b₀ (por exemplo, 5 dígitos para N= 32 entradas), são alocados para o índice com os dígitos invertidos b₀b₁b₂b₃b₄. Considere o último estágio de um algoritmo radix-2 DIT como o apresentado acima, onde a saída é escrita in-place sobre a entrada: quando E k {\displaystyle E_{k}} e O k {\displaystyle O_{k}} são combinados com uma DFT de tamanho 2, esses dois valores são sobrescritos pelas saídas. No entanto, os dois valores de saída devem ir na primeira e na segunda metade do arranjo de saída, correspondendo ao bit mais significativo b4 (para N=32); enquanto as duas entradas E k {\displaystyle E_{k}} e O k {\displaystyle O_{k}} estão intercaladas nos elementos pares e ímpares, correspondendo ao bit menos significativo b0. Assim, para obter a saída no local correto, b0 deve substituir b4 e o índice passa a ser b0b4b3b2b1. E para o próximo estágio recursivo, esses 4 bits menos significativos passarão a ser b1b4b3b2, se você incluir todos os estágios recursivos de um algoritmo radix-2 DIT, todos os bits devem ser invertidos e, portanto, é necessário pré-processar a entrada (ou pós-processar a saída) com uma inversão de bits para obter a saída na ordem correta. (Se cada subtransformação de tamanho N/2 operar em dados contíguos, a entrada DIT é pré-processada por inversão de bits). Correspondentemente, se você executar todas as etapas na ordem inversa, obterá um algoritmo DIF de base 2 com inversão de bits no pós-processamento (ou pré-processamento, respectivamente).