Exercício 26.3. Implemente e teste a transformada rápida de Fourier no Julia (apenas para vetores de dimensão 2^k ).

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 é n=2^k. 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 xjj=0^n-1 é dada por
Xk ;=;sumj=0^n-1 xj ,e^-2pi i,j k/n, k=0,1,dots,n-1
No algoritmo de Cooley–Tuk dividimos o somatório em termos com j par e j ímpar. Definindo
x^(e)m = x2m e x^(o)m = x2m+1 para m=0,dots, n2-1, e
omega = e^-2pi i/n,
chegamos à recursão
Xk = Ek + omega^k,Ok, Xk+tfrac n2 = Ek - omega^k,Ok,
onde Ek e Ok são as FFTs de tamanho n/2 de x^(e) e x^(o), respectivamente.

Passo 3

Vamos traduzir isso em código Julia. A ideia é:
. Testar se e retornar o valor base.
. Calcular recursivamente fft(even) e fft(odd) para metades do vetor.
. Construir o vetor resultado combinando com potências de .
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.

Exercícios de Livros Relacionados