Exercício 26.3. Implemente e teste a transformada rápida de Fourier no Julia (apenas para vetores de dimensão
).
Passo 1
O enunciado pede para escrever do zero uma rotina em Julia que calcule a transformada rápida de Fourier (FFT) para vetores cujo tamanho é
. O assunto aqui é a FFT pelo método de Cooley–Tuk, usando divisão em índices pares e ímpares e fator de torção. Para isso vamos precisar entender a definição da DFT e como ela se reduz recursivamente em metade do tamanho.
Passo 2
A DFT de um vetor
é dada por
No algoritmo de Cooley–Tuk dividimos o somatório em termos com
par e
ímpar. Definindo
e
para
, e
,
chegamos à recursão
chegamos à recursão
onde
e
são as FFTs de tamanho
de
e
, respectivamente.
Passo 3
Vamos traduzir isso em código Julia. A ideia é:
. Testar se
e retornar o valor base.
. Calcular recursivamente
. Construir o vetor resultado combinando com potências de
.
fft(even)
e fft(odd)
para metades do vetor.function \texttt{minha\_fft}(x::Vector{ComplexF\(64\)})
n = length(x)
if n == 1
return copy(x)
elseif n % 2 != 0
error("mínimo: tamanho deve ser potência de \(2\)")
end
# separar pares e ímpares
even = minha_fft(x[1:2:end])
odd = minha_fft(x[2:2:end])
# fator de torção
X = Vector{ComplexF\(64\)}(undef, n)
for k in 0:(n \div 2 - 1)
\omega_k = exp(-2im * \pi * k / n)
X[k+1] = even[k+1] + \omega_k * odd[k+1]
X[k+1 + n\div 2] = even[k+1] - \omega_k * odd[k+1]
end
return X
end
Passo 4
Agora vamos testar comparando com a FFT built-in de Julia:
using Random, FFTW
# vetor de teste de dimensão \(8\) (\(= 2^3\))
Random.seed!(\(1234\))
x \(=\) randn(\(8\)) .\(+\) im*randn(\(8\))
Y\(1= minha_fft\)(x)
Y\(2=\) fft(x) # FFTW.fft
println("Erro máximo: ", maximum(abs.(Y\(1\) .- Y\(2\))))
Em meu teste o erro máximo fica na ordem de
, o que confirma que nossa implementação está correta para entradas de tamanho
.
Resposta
Coletei acima a função
\(minha_fft\)
em Julia e mostrei um teste de comparação com fft
do pacote FFTW, obtendo erro numérico desprezível.