Grupo4 - FFT

A Transformada rápida de Fourier (em inglês Fast Fourier Transform, ou FFT) é um algoritmo que torna viável o cálculo da Transformada Discreta de Fourier (DFT), que é a forma discretizada da [Transformada de Fourier]. A FFT permite transformar de forma rápida uma série de sinais discretos em uma amostra contendo as frequências desses sinais, desde que satisfaça algumas propriedades.

Transformada Discreta de Fourier

Em muitas aplicações se tem informação sobre um conjunto de dados, ao invés de uma função contínua. A Transformada Discreta de Fourier transforma esse conjunto de dados em um conjunto de tamanho igual com informação sobre as frequências da função que satisfaz o conjunto de dados.

Para um conjunto de dados igualmente espaçados, pode-se, ao considerar os dados como um período de uma função periódica, cujo período normalmente é considerado entre $[-\pi ,\pi ]$ para facilitar o cálculo (e que pode sempre ser transformada em uma função nesse interválo), mostrar que a transformada discreta de Fourier pode ser dada pela equação:

$F_{k}=\sum _{n=0}^{N-1}f_{n}e^{-i2\pi nk/N}$

A sua inversa é, em paralelo ao caso da transformada contínua,

$f_{n}={\frac {1}{N}}\sum _{n=0}^{N-1}F_{k}e^{i2\pi nk/N}$

A transformada também pode ser expressa em forma vetorial, como

${\vec {A}}=\mathbf {W^{nk}} {\vec {a}}~$ onde $~\mathbf {W^{nk}} ~$ é definido como

$\mathbf {W^{nk}} ={\begin{bmatrix}e^{-2\pi i0\cdot 0/N}&e^{-2\pi i0\cdot 1/N}&\cdots &e^{-2\pi i0\cdot k/N}\\e^{-2\pi i1\cdot 0/N}&e^{-2\pi i1\cdot 1/N}&\cdots &e^{-2\pi i1\cdot k/N}\\\vdots &\vdots &\ddots &\vdots \\e^{-2\pi in\cdot 0/N}&e^{-2\pi in\cdot 1/N}&\cdots &e^{-2\pi in\cdot k/N}\end{bmatrix}}$

O cálculo dessa expressão leva em torno de $N^{2}$ passos para o resultado. Uma amostra com 3,000 pontos precisa de 9,000,000 operações para a transformada ser obtida, tornando a DFT inviável para aplicações rápidas.

Transformada Rápida de Fourier

É possível calcular a transformada com $N\log _{2}N$ passos. Para isso se dispõe de um algoritmo chamado Transformada Rápida de Fourier. Considera-se um conjunto de pontos $N=2^{p}$ (com $p$ inteiro, então, da definição da DFT

$F_{k}=\sum _{n=0}^{N-1}f_{n}\cdot e^{-i{\frac {2\pi }{N}}\cdot k\cdot n}$

podemos dividir o somatório em 2:

$F_{k}=\color {red}\sum _{n=0}^{N/2-1}f_{2n}\cdot e^{-i{\frac {2\pi }{N}}\cdot k\cdot 2n}\color {black}+\color {blue}\sum _{n=0}^{N/2-1}f_{2n+1}\cdot e^{-i{\frac {2\pi }{N}}\cdot k\cdot (2n+1)}$

$\color {black}F_{k}=\color {red}\sum _{n=0}^{N/2-1}f_{2n}\cdot e^{-i{\frac {2\pi }{N/2}}\cdot k\cdot n}\color {black}+\color {blue}e^{-i{\frac {2\pi }{N}}k}\sum _{n=0}^{N/2-1}f_{2n+1}\cdot e^{-i{\frac {2\pi }{N/2}}\cdot k\cdot n}$

$\color {black}F_{k}=\color {red}\sum _{n=0}^{N/2-1}f_{2n}\cdot e^{-i{\frac {2\pi }{N/2}}\cdot k\cdot n}\color {black}+\color {blue}~C(k)~\sum _{n=0}^{N/2-1}f_{2n+1}\cdot e^{-i{\frac {2\pi }{N/2}}\cdot k\cdot n}$

onde a soma em vermelho é a parte par e a soma em azul é a parte ímpar da transformada. As duas somas tem o mesmo expoente, que agora é dividido por $N/2$. Desse expoente, é evidente a relação entre o ponto $k$ e o ponto $k+N/2$

$e^{-i{\frac {2\pi }{N/2}}\cdot k\cdot n}=e^{-i{\frac {2\pi }{N/2}}\cdot (k+N/2)\cdot n}=e^{-i{\frac {2\pi }{N/2}}\cdot k\cdot n}\cdot e^{-i{\frac {2\pi }{N/2}}\cdot N/2\cdot n}=e^{-i{\frac {2\pi }{N/2}}\cdot k\cdot n}\cdot 1$

Com essa relação, podemos ver que $F_{k}$ e $F_{k+N/2}$ tem o mesmo expoente e podem ser calculadas ao mesmo tempo. Mais que isso, a nova forma da transformada pode ser sucessivamente dividida, cada vez produzindo somas com limites menores.

Exemplo

Suponha que temos a função sinusoidal $a(t)=sin(2\pi \cdot 1Hz\cdot t)$ e fazemos quatro medidas no intervalo de 1 segundo, resultando em

$a_{0}=0.00$ $a_{1}=1.00$ $a_{2}=0.00$ $a_{3}=-1.00$

Com essas 4 medidas, podemos dividir a soma 2 vezes:

$A_{k}=\sum _{t=0}^{3}a_{t}\cdot e^{-i{\frac {2\pi }{4}}\cdot k\cdot t}$

$A_{k}=\sum _{t_{1}=0}^{1}a_{2t_{1}}\cdot e^{-i{\frac {2\pi }{2}}\cdot k\cdot t_{1}}+C_{k}^{1}\sum _{t_{1}=0}^{1}a_{2t_{1}+1}\cdot e^{-i{\frac {2\pi }{2}}\cdot k\cdot t_{1}}$

$A_{k}=\sum _{t_{2}=0}^{0}a_{4t_{2}}\cdot e^{-i{\frac {2\pi }{1}}\cdot k\cdot t_{2}}+C_{k}^{2}\sum _{t_{2}=0}^{0}a_{4t_{2}+2}\cdot e^{-i{\frac {2\pi }{1}}\cdot k\cdot t_{2}}+C_{k}^{1}\sum _{t_{2}=0}^{0}a_{4t_{2}+1}\cdot e^{-i{\frac {2\pi }{1}}\cdot k\cdot t_{2}}+C_{k}^{3}\sum _{t_{2}=0}^{0}a_{4t_{2}+3}\cdot e^{-i{\frac {2\pi }{1}}\cdot k\cdot t_{2}}$

e como temos $C_{k}^{j}=(e^{-i{\frac {2\pi }{N}}k})^{j}$ podemos calcular

$A_{0}=1.00\cdot C_{0}^{1}-1.00\cdot C_{0}^{3}=0.00+i0.00$

$A_{1}=1.00\cdot C_{1}^{1}-1.00\cdot C_{1}^{3}=0.00-i2.00$

$A_{2}=1.00\cdot C_{2}^{1}-1.00\cdot C_{2}^{3}=0.00+i0.00$

$A_{3}=1.00\cdot C_{3}^{1}-1.00\cdot C_{3}^{3}=0.00+i2.00$

FFT para N diferente de uma potência de 2

Mesmo com a FFT sendo um algoritmo extremamente eficiente para $N=2^{p}$, esse dificilmente é o caso que encotramos. Ainda assim, para $N$ altamente composto ($N=r_{1}\cdot r_{2}\cdot ...\cdot r_{m}$) o algoritmo ainda resulta em uma boa queda no tempo de cálculo.

Para o caso mais simples $N=r_{1}\cdot r_{2}$ a transformada pode ser escrita como

$F(k_{1},k_{0})=\sum _{n_{0}=0}^{r_{2}-1}\left[\sum _{n_{1}=0}^{r_{1}-1}x(n_{1},n_{0})e^{-i2\pi \cdot k\cdot n_{1}\cdot r_{2}}\right]e^{-i2\pi \cdot k\cdot n_{0}}$

onde

$k=k_{1}\cdot r_{1}+k_{0}~~~~~~~~k_{0}=0,1,...,r_{1}-1~~~~~~~~k_{1}=0,1,...,r_{2}-1$

$n=n_{1}\cdot r_{2}+n_{0}~~~~~~~~n_{0}=0,1,...,r_{2}-1~~~~~~~~n_{1}=0,1,...,r_{1}-1$

assim a transformada que antes necessitava de N calculos, agora pode ser vista como $r_{1}+r_{2}$ calculos.

Algorítmo (usando recursão)

$>>\mathbf {function} ~FFT(N,{\vec {a}})$

$>>>~~~~\mathbf {if} ~N=1~\mathbf {then}$

$>>>~~~~~~~~return~{\vec {a}};$

$>>>~~~~\mathbf {else}$

$>>>~~~~~~~~E=FFT(N/2,(a_{0},a_{2},...,a_{N-2}));$

$>>>~~~~~~~~O=FFT(N/2,(a_{1},a_{3},...,a_{N-1}));$

$>>>~~~~~~~~\mathbf {for} ~k=0~to~k<=N/2-1~\mathbf {do}$

$>>>~~~~~~~~~~~~A_{k}=E_{k}+exp(-i2\pi k/N)*O_{k};$

$>>>~~~~~~~~~~~~A_{k+N/2}=E_{k}-exp(-i2\pi k/N)*O_{k};$

$>>>~~~~return~{\vec {A}};$

Implementação

A natureza recursiva do algoritmo da transformada rápida de fourier o torna ideal para implementações em linguagens funcionais como Haskell. Abaixo exibimos as partes mais relevantes do código, onde omitimos inclusões de bibliotecas e a definição das funções auxiliares evens e odds.

--esse termo multiplica o somatório de termos ímpares
f xs n = exp $ -(0:+2*pi)*n / genericLength xs

--caso a função seja chamada com apenas dois termos, temos o caso base
ffti [x,y] n = x + y * (exp $ -(0:+pi)*n)

--caso contrário aplique a função recursivamente separando o somatório em 
-- termos com índice par e com índice ímpar
ffti xs n = ffti (evens xs) n + f xs n * ffti (odds xs) n

--função que calcula os n coeficientes (ffti calcula um coeficiente)
fft xs [] = []
fft xs (y:ys) = (ffti xs y):(fft xs ys)

--exemplo
fft [0, 1, 4, 9] [0,1,2,3]

--saída ( :+ indica número complexo em Haskell)
[14.0 :+ 0.0, -4.0 :+ 8.0,  -6.0 :+ 0.0,  -4.0 :+ -8.0)]

Embora não seja uma linguagem puramente funcional Wolfram Mathematica também se presta a uma implementação direta e clara.

// caso base de uma lista contendo apenas dois termos
fft[{x_, y_}, n_] := x + y Exp[-I Pi n ]
// caso contrário, divida a lista entre termos com índice ímpar e par 
// e aplique a função recursivamente
fft[f_List, n_] :=                fft[Downsample[f, 2, 1], n] + 
    Exp[-((I 2 Pi n )/Length[f])] fft[Downsample[f, 2, 2], n]
//acima calcula-se com um coeficiente abaixo calcula-se todos
fft[f_List] := Table[fft[f, n] // N, {n, 0, Length[f] - 1}]

//exemplo
fft[{0, 1, 4, 9}]
{14., -4. + 8. I, -6., -4. - 8. I}

Exemplos

Alguns exemplos básicos da transformada de fourier utilizando o código acima:

Gaussiana