3.2.5 Aproximação de Funções — Análise de Fourier — A Transformada Discreta Rápida de Fourier
1. A Transformada Rápida de Fourier
A transformada rápida de Fourier é um algoritmo desenvolvido para calcular a transformada discreta de Fourier utilizando menos recursos computacionais. Conhecida por FFT (Fast Fourier Transform), Sua eficiência se deve ao fato dela reaproveitar resultados computados anteriormente, reduzindo o número total de operações aritméticas.
Algoritmos que implementam FFTs exploram a simetria e a periodicidade das funções trigonométricas para calcular a transformada. Este tipo de abordagem reduz a quantidade de operações, necessárias para computar a TDF, para operações. Em problemas reais, isto significa um ganho de eficiência considerável, conforme podemos comparar no gráfico abaixo:
Em 1965, J.W. Cooley e J.W. Tukey publicaram o primeiro algoritmo para o cálculo da transformada rápida de Fourier, utilizando-se decimação no domínio do tempo. Posteriormente, muitos algoritmos derivaram-se da técnica de Cooley-Tukey, utilizando-se do princípio de que uma transformada discreta pode ser dividida — decimada — em transformadas menores, uma vez que muitas
operações são redundantes.
Outras abordagens permitem o cálculo da FFT, mas sem melhoria na complexidade do algoritmo. Além das classes de algoritmos derivados de Cooley-Tukey, temos aqueles que implementam decimação no domínio da frequência. Esses algoritmos são também conhecidos como métodos de Sande-Tukey.
1.1. Algoritmo de Sande-Tukey
Para que seja possível calcular a FFT por este algoritmo, devemos ter como entrada uma quantidade par de valores. Ou seja,
Como apresentado em outro post, a transformada discreta de Fourier é definida como
para .
Supondo que dividimos cada amostra na metade, e decompondo a equação acima em termos de pontos, temos que:
Se criarmos uma nova variável , observamos que a segunda somatória possui o mesmo número de elementos da primeira:
com isso, podemos agrupar novamente os termos em um único somatório
Ao utilizarmos a propriedade
notamos que os pontos são sempre projetados. Isto é, para par, esse fator é igual à , para ímpar é . Com isso, decompomos a transformada em duas, uma relativa às pares
e outra relativa às ímpares
para .
Se definirmos
podemos reescrever a TDF como
e as equações decompostas como
e
Dessas expressões, podemos notar que elas se relacionam por
e
portanto,
Ou seja, em vez de um utilizarmos pontos, utilizamos dois cálculos para pontos.
A TDF é calculada formando-se inicialmente as sequências e e então calculando-se as TDFs para pontos para obtermos as transformadas ímpares e pares, independentemente. Os pesos são frequentemente chamados de twiddle factors na literatura.
O processo de dividir a transformada em duas transformadas independentes, tratando apenas pontos, é repetida novamente para cada novo intervalo seccionado. Ou seja, calculamos as TDFs para pontos das quatro sequências de comprimento compostas dos primeiros e últimos pontos. Assim, repetindo-se esta estratégia até a solução, teremos aproximadamente operações.
1.2. Algoritmo de Cooley-Tukey
A abordagem que utiliza decimação no tempo, utiliza-se do inverso do que faz o algoritmo de Sande-Tukey. Embora os dois métodos sejam diferentes na ordem de processamento dos dados, ambos possuem complexidade equivalente, de operações. Neste algoritmo, a amostra original é dividida inicialmente em pontos rotulados com numeração par e ímpar, e as diversas TDFs são aplicadas sobre essas divisões.
2. Implementação
2.1. Algoritmo de Sande-Tukey
def fft_sandetukey(inputf):
"""
Return a list with the coeficients of Fast Fourier Transform
(Sande-Tukey -- Decimation in Frequency algorithm).
F = fft_sandetukey(inputf)
INPUT:
* f: list with discrete points of time-domain
return: list discrete points of frequency-domain
Author: Pedro Garcia [sawp@sawp.com.br]
see: http://www.sawp.com.br
License: Creative Commons
http://creativecommons.org/licenses/by-nc-nd/2.5/br/
Mar 2011
"""
def fft_st(x, n):
if n != 1:
d = [complex(0.0) for i in xrange(n)]
e = [complex(0.0) for i in xrange(n)]
nh = n / 2
for k in xrange(nh, -1, -1):
s = x[k]
t = x[k + nh]
e[k] = s + t
d[k] = s - t
for k in xrange(nh):
theta = k * omega * n
real = cos(theta)
imag = sin(theta)
d[k] = d[k] * complex(real, imag)
# divide and conquer
bige = fft_st(e, nh)
bigd = fft_st(d, nh)
j = 0
for k in xrange(nh):
x[j] = bige[k]
x[j + 1] = bigd[k]
j = j + 2
return x
f = map(complex, inputf)
m = len(f)
omega = -2.0 * pi / m
return fft_st(f, m - 1)
Esta função está disponível em http://www.sawp.com.br/code/interpolation/fft_sandetukey.py assim como um exemplo de como utilizá-la.
2.2. Algoritmo de Cooley-Tukey
def fft_cooleytukey(inputf):
"""
Return a list with the coeficients of Fast Fourier Transform
(Sande-Tukey -- Decimation in Time algorithm).
F = fft_cooleytukey(inputf)
INPUT:
* f: list with discrete points of time-domain
return: list discrete points of frequency-domain
Author: Pedro Garcia [sawp@sawp.com.br]
see: http://www.sawp.com.br
License: Creative Commons
http://creativecommons.org/licenses/by-nc-nd/2.5/br/
Mar 2011
"""
def is_power_of_two(number):
l = log(number, 2)
return (l == floor(l))
def fft_ct(x):
length = len(x)
if length == 1:
return x
even = x[0::2]
odd = x[1::2]
evenFFT = fft_cooleytukey(even)
oddFFT = fft_cooleytukey(odd)
pivot = int(length / 2)
omega = 2 * pi / length
left = []
right = []
for k in xrange(pivot):
angle = k * omega
real = cos(angle)
imag = sin(angle)
base = complex(real, imag)
left = left + [evenFFT[k] + base * oddFFT[k]]
right = right + [evenFFT[k] - base * oddFFT[k]]
F = left + right
return F
f = map(complex, inputf)
n = len(f)
if is_power_of_two(n):
F = fft_ct(f)
else:
F = dft(f)
return F
Esta função está disponível em http://www.sawp.com.br/code/interpolation/fft_cooleytukey.py assim como um exemplo de como utilizá-la.
3. Copyright
Este documento é disponível sob a licença Creative Commons. As regras dos direitos de cópia deste conteúdo estão acessíveis em http://creativecommons.org/licenses/by-nc-nd/2.5/br/.