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

$\overrightarrow{A}={\mathbf{𝐖}}^{\mathbf{𝐧}\mathbf{𝐤}}\overrightarrow{a}$ onde ${\mathbf{𝐖}}^{\mathbf{𝐧}\mathbf{𝐤}}$ é definido como

${\mathbf{𝐖}}^{\mathbf{𝐧}\mathbf{𝐤}}=\left[\begin{array}{cccc}{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}\\ ⋮& ⋮& \ddots & ⋮\\ e^{-2\pi in\cdot 0/N}& e^{-2\pi in\cdot 1/N}& \cdots & e^{-2\pi in\cdot k/N}\end{array}\right]$

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=\sum _{n=0}^{N/2-1}{f}_{2n}\cdot e^{-i\frac{2\pi }{N}\cdot k\cdot 2n}+\sum _{n=0}^{N/2-1}{f}_{2n+1}\cdot e^{-i\frac{2\pi }{N}\cdot k\cdot (2n+1)}$

$F_k=\sum _{n=0}^{N/2-1}{f}_{2n}\cdot e^{-i\frac{2\pi }{N/2}\cdot k\cdot n}+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}$

$F_k=\sum _{n=0}^{N/2-1}{f}_{2n}\cdot e^{-i\frac{2\pi }{N/2}\cdot k\cdot n}+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}[\sum _{n_1=0}^{r_1-1}x(n_1,n_0)e^{-i2\pi \cdot k\cdot n_1\cdot r_2}]e^{-i2\pi \cdot k\cdot n_0}$

onde

$k=k_1\cdot r_1+k_0{k}_0=0,1,...,r_1-1k_1=0,1,...,r_2-1$

$n=n_1\cdot r_2+n_0{n}_0=0,1,...,r_2-1n_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{𝐟}\mathbf{𝐮}\mathbf{𝐧}\mathbf{𝐜}\mathbf{𝐭}\mathbf{𝐢}\mathbf{𝐨}\mathbf{𝐧}FFT(N,\overrightarrow{a})$

$>>>\mathbf{𝐢}\mathbf{𝐟}N=1\mathbf{𝐭}\mathbf{𝐡}\mathbf{𝐞}\mathbf{𝐧}$

$>>>return\overrightarrow{a};$

$>>>\mathbf{𝐞}\mathbf{𝐥}\mathbf{𝐬}\mathbf{𝐞}$

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

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

$>>>\mathbf{𝐟}\mathbf{𝐨}\mathbf{𝐫}k=0tok<=N/2-1\mathbf{𝐝}\mathbf{𝐨}$

$>>>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\overrightarrow{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

Abaixo encontra-se alguns exemplos básicos da transformada de fourier com base no código acima porém com modificações para lidar com o caso em que o intervalo de amostragem é diferente de $[0,2\pi ]$ e cuja taxa de amostragem é qualquer (diferente de $\frac{2\pi k}{N}$ como assumido acima)

Gaussiana

A menos de uma normalização a função gaussiana é definida por: $f(x)=e^{-\frac{(x-<x>{)}^2}{2{\sigma }^2}}$ Aplicando a transformada de fourier na curva gaussiana obtem-se outra gaussiana.

Observa-se que uma curva larga no espaço real corresponde a uma curva estreita no espaço de fourier e vice-versa. Esse é um fato matemático amplamente conhecido e se exprime também na física pelo Princípio da Incerteza de Heisenberg. Abaixo uma curva gaussiana mais larga que a anterior e sua transformada correspondentemente mais estreita.

Uma função par em relação à origem possui uma expansão em fourier com termos ímpares nulos, no entanto, deslocando a gaussiana um pouco para a direita e, portanto, quebrando sua simetria, vemos que a sua transformada possui uma parte complexa:

Oscilador Harmônico Amortecido

A curva que representa um oscilador harmônico amortecido é descrita por $f(x)=sin(ax)e^{-bx}$ onde o termo senoidal, periódico, é amortecido pelo termo exponencial que leva a oscilação rapidamente a zero.

A transformada de fourier dessa curva possui parte real e imaginária, da parte imaginária dessa curva é possível obter a frequência de ressonância do fenômeno representado por ela. No caso abaixo vemos que a frequência de ressonância vale $\pm 10Hz$

Abaixo a parte real da transformada:

Aplicações

Equação de Schrödinger: propagação de onda de probabilidade

O cálculo da evolução temporal da equação de schrödinger é computacionalmente exigente. A cada passo de tempo é necessário recalcular a posição da onda de probabilidade que representa a partícula em cada ponto do domínio. Portanto, torna-se vantajoso o uso da transformada de fourier rápida. Além disso a própria transformada faz parte da solução, pois também deseja-se obter a sua evolução temporal no espaço de momento. O código foi desenvolvido por [https://jakevdp.github.io/blog/2012/09/05/quantum-python/ Jake Vanderplas, no entanto, sob extensas modificações foi possível entender o código e modifica-lo para representar situações diferentes das propostas pelo autor. Em seu site o autor explica a teoria por trás do problema e também como funciona seu algoritmo. A seguir exploramos a interação de uma onda de probabilidade com potenciais de diferentes formatos.

A parte principal do código modificado é exibidade abaixo:

Potencial Quadrado

Square.gif

Potencial Batman

Batman.gif