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 \(N^2 \) operações, necessárias para computar a TDF, para \(N log_2 N \) 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,

\(N = 2 ^ M\)

Como apresentado em outro post, a transformada discreta de Fourier é definida como
\(F_k = \sum_{j=0}^{N} f_j e^{-i \frac{2 \pi}{N} j k} \)

para \(k = 0, 1, \ldots, N \) .

Supondo que dividimos cada amostra na metade, e decompondo a equação acima em termos de \(\frac{N}{2} \) pontos, temos que:

\(F_k = \sum_{j=0}^{(N/2) – 1} f_j e^{\frac{-i 2 \pi k j}{N}} + \sum_{j=N/2}^{N} f_j e^{\frac{-i 2 \pi k j}{N}}\)

Se criarmos uma nova variável \(m = j – \frac{N}{2} \) , observamos que a segunda somatória possui o mesmo número de elementos da primeira:

\(F_k = \sum_{j=0}^{(N/2) – 1} f_j e^{\frac{-i 2 \pi k j}{N}} + \sum_{m=0}^{(N/2) – 2} f_{m+N} e^{-i (2 \pi / N) k (m + N/2)}\)

com isso, podemos agrupar novamente os termos em um único somatório

\(F_k = \sum_{j=0}^{(N/2) – 1} (f_j + e^{-i \pi k} f_{n+N/2}) e^{\frac{-i 2 \pi k j}{N}}\)

Ao utilizarmos a propriedade

\(e^{-i \pi k} = (-1) ^ k\)

notamos que os pontos são sempre projetados. Isto é, para \(k \) par, esse fator é igual à \(1 \) , para \(k \) ímpar é \(-1 \) . Com isso, decompomos a transformada em duas, uma relativa às pares

\(F_{2k} = \sum_{j=0}^{(N/2)-1} (f_j + f_{j+N/2}) e^{\frac{-i 2 \pi k j}{N/2}}\)

e outra relativa às ímpares

\(F_{2k+1} = \sum_{j=0}^{(N/2)-1} (f_j – f_{j+N/2}) e^{\frac{-i 2 \pi k j}{N/2}} e^{\frac{-i2 \pi j}{N}}\)

para \(k = 0, 1, \ldots, N \).

Se definirmos

\(W = e ^{-i \frac{2 \pi}{N}}\)

podemos reescrever a TDF como

\(F_k = \sum_{j=0}^{N} f_n W^{jk}\)

e as equações decompostas como

\(F_{2k} = \sum_{j=0}^{(N/2)-1} (f_j + f_{j+N/2}) W^{2jk}\)

e
\(F_{2K+1} = \sum_{j=0}^{(N/2)-1} (f_j – f_{j+N/2}) W^j W^{2jk}\)

Dessas expressões, podemos notar que elas se relacionam por

\(g_j = f_j + f_{j+N/2}\)

e

\(h_j = (f_j + f_{j+N/2}) W^j\)

portanto,

\(F_{2k} = G_k \) e \(F_{2k+1} = H_k \)

Ou seja, em vez de um utilizarmos \(N \) pontos, utilizamos dois cálculos para \(N/2\) pontos.

A TDF é calculada formando-se inicialmente as sequências \(g_n\) e \(h_n\) e então calculando-se as TDFs para \(\frac{N}{2} \) pontos para obtermos as transformadas ímpares e pares, independentemente. Os pesos \(W^j\) são frequentemente chamados de twiddle factors na literatura.

O processo de dividir a transformada em duas transformadas independentes, tratando apenas \(\frac{N}{2} \) pontos, é repetida novamente para cada novo intervalo seccionado. Ou seja, calculamos as TDFs para \(\frac{N}{4}\) pontos das quatro sequências de comprimento \(\frac{N}{4} \) compostas dos primeiros e últimos pontos. Assim, repetindo-se esta estratégia até a solução, teremos aproximadamente \(N \log_2 N \) 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 \(N \log_2 N\) 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/.

References

[1] Anthony Ralston and Philip Rabinowitz, A First Course in Numerical Analysis (2nd ed.), McGraw-Hill and Dover, (2001).
[2] N.B Franco, Cálculo Numérico, Pearson Prentice Hall (2006).