O algoritmo de Cooley–Tukey designa uma família de algoritmos de transformada rápida de Fourier (FFT), nomeada em homenagem a J. W. Cooley e John Tukey, usada para calcular a transformada discreta de Fourier (DFT). Para um comprimento composto , a decomposição geral expressa uma DFT de tamanho por meio de DFTs de tamanho , multiplicações por fatores de rotação e DFTs de tamanho . Quando aplicada recursivamente a comprimentos que podem ser decompostos em fatores pequenos, ela reduz o custo da avaliação direta, de , para . Isso ocorre, em particular, para potências de dois e, de forma mais geral, para números suaves. Essa redução de custo explica a importância do algoritmo em processamento digital de sinais e em outras aplicações da DFT.

A fatoração pode empregar diferentes bases em seus estágios e pode ser combinada com outros métodos. Quando uma subtransformação tem comprimento primo, podem ser usados o algoritmo de Rader ou o algoritmo de Bluestein. Para fatores coprimos, o algoritmo FFT de fator primo oferece uma decomposição alternativa.

Uma fatoração equivalente, incluindo sua aplicação recursiva, foi descrita por Carl Friedrich Gauss por volta de 1805, embora seu trabalho só tenha sido publicado postumamente. Cooley e Tukey redescobriram e popularizaram o método em 1965.

História Uma fatoração equivalente à empregada no algoritmo, incluindo sua aplicação recursiva, foi descrita por volta de 1805 por Carl Friedrich Gauss, que a empregou para interpolar as trajetórias dos asteroides Palas 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. As FFTs se popularizaram 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, que contribuiu para a ampla difusão da FFT.

O fato de Gauss ter descrito o mesmo algoritmo (embora sem analisar seu custo assintótico) não foi percebido até vários anos após o artigo de Cooley e Tukey de 1965. Seu artigo citou como inspiração apenas o trabalho de I. J. Good sobre o que agora é chamado de algoritmo FFT de fator primo (PFA); embora o algoritmo de Good tenha sido inicialmente considerado equivalente ao algoritmo de Cooley–Tukey, rapidamente percebeu-se que o PFA é um algoritmo bastante diferente (funcionando apenas para tamanhos que têm fatores relativamente primos e contando com o teorema do resto chinês, ao contrário do suporte a qualquer tamanho composto presente no algoritmo de Cooley–Tukey).

Decomposição geral Para uma sequência de comprimento , adota-se a definição

Seja e defina-se . Reescrevendo os índices como

obtém-se a fatoração

Nessa expressão, e . A fatoração pode ser interpretada em três etapas:

calcular DFTs de tamanho ; multiplicar os resultados pelos fatores de rotação ; calcular DFTs de tamanho .

direita|miniaturadaimagem|300px|Interpretação matricial de uma etapa do algoritmo de Cooley–Tukey. A sequência de comprimento é organizada como uma matriz; calculam-se DFTs em uma direção, aplicam-se os fatores de rotação e calculam-se DFTs na outra direção. As transposições e permutações dependem da organização adotada para os dados.

O fator pequeno empregado em cada estágio é denominado base (radix). Na dizimação no tempo (DIT), a decomposição começa pela separação dos índices da sequência de entrada; na dizimação na frequência (DIF), começa pela separação dos índices da saída. As duas formas implementam a mesma fatoração, mas diferem na ordem das operações e nas permutações dos dados.

Se e a mesma base fixa for empregada recursivamente, o custo satisfaz uma recorrência da forma , que resulta em .

O caso DIT de base 2 A FFT de dizimação no tempo de base 2 é obtida tomando na decomposição geral. Em cada estágio, a sequência de entrada é separada em duas subsequências intercaladas: uma formada pelas amostras de índices pares e outra pelas amostras de índices ímpares. Calculam-se então duas DFTs de tamanho , cujos resultados são combinados por operações denominadas borboletas. A expressão “base 2” decorre dessa decomposição por um fator 2.

Aplicada recursivamente, essa forma pressupõe que seja uma potência de dois. Para outros comprimentos, podem ser usadas bases mistas e, para fatores primos, os algoritmos de Rader ou de Bluestein. O preenchimento com zeros pode adequar o comprimento a uma implementação e produzir uma amostragem mais densa do espectro, mas não acrescenta informação ao sinal nem aumenta sua resolução espectral intrínseca.

O algoritmo DIT de base 2 reorganiza a DFT da sequência em duas partes: uma soma sobre os índices pares e uma soma sobre os índices ímpares :

Pode-se fatorar um multiplicador comum da segunda soma, conforme mostrado na equação abaixo. Fica então claro que as duas somas correspondem à DFT da parte de índices pares e à DFT da parte de índices ímpares da sequência . Denotando a DFT das entradas de índices pares por e a DFT das entradas de índices ímpares por , obtém-se:

Como e são DFTs de comprimento , basta calculá-las para . Os demais valores de são obtidos pelas relações de periodicidade usadas a seguir. Graças à periodicidade do exponencial complexo, também é obtido de e :

Assim, podem-se reescrever e como:

As saídas são obtidas pela soma e pela diferença de e . Essas duas operações constituem uma DFT de tamanho 2 e são usualmente denominadas operações borboleta. Em decomposições de bases maiores, essa DFT de tamanho 2 é substituída por uma DFT do tamanho correspondente à base.

Este processo é um exemplo de divisão e conquista. Em muitas implementações, evita-se a recursão explícita e processam-se os estágios iterativamente.

A reescrita acima de uma DFT de tamanho N como duas DFTs de tamanho N/2 é às vezes chamada de lema de Danielson–Lanczos, uma vez que a identidade foi observada por esses dois autores em 1942Danielson, G. C., and C. Lanczos, "Some improvements in practical Fourier analysis and their application to X-ray scattering from liquids," J. Franklin Inst. 233, 365–380 and 435–452 (1942). (influenciados pelo trabalho de Runge de 1903). Eles aplicaram seu lema de forma recursiva "reversa", dobrando repetidamente o tamanho da DFT até que o espectro de transformação convergisse (embora aparentemente não tenham reconhecido que haviam obtido uma complexidade assintótica log-linear, de ordem ). O trabalho de Danielson–Lanczos antecedeu a ampla disponibilidade de computadores mecânicos ou eletrônicos e exigiu cálculo manual (possivelmente com auxílios mecânicos como máquinas de somar); eles relataram um tempo de computação de 140 minutos para uma DFT de tamanho 64 operando em entradas reais com 3 a 5 dígitos significativos. O artigo de Cooley e Tukey de 1965 relatou um tempo de execução de 0,02 minuto para uma DFT complexa de tamanho 2048 em um IBM 7094.

Implementação recursiva Em pseudocódigo, o procedimento abaixo poderia ser escrito: X0,...,N−1 ← ditfft2(x, N, s): # DFT de (x0, xs, x2s, ..., x(N-1)s): if N = 1 then X0 ← x0 # caso base trivial da DFT de tamanho 1 else 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) for k = 0 to (N/2)-1 do # combinar as DFTs das duas metades p ← Xk q ← exp(−2πi k/N) Xk+N/2 Xk ← p + q Xk+N/2 ← p − q end for end if Aqui, a chamada ditfft2(x, N, 1) calcula X=DFT(x) em um vetor separado. Nesse procedimento, N é uma potência de 2 e s representa o passo entre amostras do vetor de entrada. A notação x+s representa a subsequência iniciada em xs.

(Os resultados estão na ordem correta em X e nenhuma outra permutação de reversão de bits é necessária; a frequentemente mencionada necessidade de um estágio separado de reversão de bits só surge em certos algoritmos in-place, como descrito abaixo.)

Implementações de FFT de alto desempenho fazem muitas modificações na implementação de tal algoritmo em comparação com este pseudocódigo simples. Por exemplo, pode-se usar um caso base maior que N=1 para amortizar a sobrecarga da recursão, os fatores de rotação (twiddle factors) podem ser pré-calculados, e bases maiores são frequentemente usadas por motivos de cache; essas e outras otimizações, em conjunto, podem melhorar o desempenho em uma ordem de magnitude ou mais. (Em muitas implementações de livros didáticos, a recursão em profundidade é eliminada em favor de uma abordagem não recursiva em largura, embora se argumente que a recursão em profundidade oferece melhor localidade de memória.Johnson, Steven G.; Frigo, Matteo. "Implementing FFTs in practice". In: Burrus, C. S. (ed.). Fast Fourier Transforms. cap. 11. Rice University: Connexions, 2008.) Várias dessas ideias são descritas em mais detalhes abaixo.

Contagem de operações Para , a FFT de base 2 possui estágios, cada qual formado por operações borboleta. Em uma implementação direta, cada borboleta envolve, no máximo, uma multiplicação complexa por um fator de rotação e duas adições complexas. Obtêm-se, assim, no máximo,

multiplicações complexas e

adições complexas. Em comparação, a avaliação direta da DFT requer multiplicações complexas e adições complexas. Alguns fatores de rotação, como , , e , podem ser aplicados sem uma multiplicação complexa geral, de modo que implementações otimizadas podem reduzir essas contagens.

Variações Implementações de base mista (mixed-radix) empregam diferentes bases ao longo da recursão e permitem tratar comprimentos compostos com fatores distintos. Quando um caso-base tem comprimento primo , ele pode ser avaliado diretamente, com custo , ou por algoritmos especializados, como os de Rader e Bluestein.

A variante split-radix combina decomposições de bases 2 e 4 e apresentou, por longo tempo, a menor contagem conhecida de operações aritméticas para comprimentos que são potências de dois. Modificações posteriores reduziram ainda mais essa contagem. Em computadores atuais, o desempenho também depende da hierarquia de memória, da vetorização e da organização dos acessos aos dados; por isso, implementações otimizadas frequentemente empregam bases maiores ou mistas.

Outra maneira de ver o algoritmo de Cooley–Tukey é que ele reexpressa uma DFT unidimensional de tamanho N como uma DFT bidimensional N1 por N2 além dos fatores de rotação, onde a matriz de saída é transposta.