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

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 é exibida abaixo. Nota-se que o autor utilizou o método Leapfrog para integrar as equações de movimento fazendo a seguinte sequência de operações: dar um meio passo no espaço real, tomar a transformada, dar um passo no espaço de momento, tomar a transformada inversa e finalmente dar um passo completo no espaço real.

def time_step():
    global t, psi_mod_x, psi_mod_k, t_max
    psi_mod_x *= x_evolve_half
    for i in range(N_steps):
        psi_mod_k = fft(psi_mod_x)
        psi_mod_k *= k_evolve
        psi_mod_x = ifft(psi_mod_k)
        psi_mod_x *= x_evolve_half*x_evolve_half
    psi_mod_k = fft(psi_mod_x)
    t += dt * N_steps
    t_max = t_max - 1

O restante do código ocupa-se somente de inicialização de valores, definição do potencial e escrita dos resultados.

Potencial Quadrado

O potencial quadrado foi o exemplo escolhido pelo autor do código. Nesse exemplo fica visível que uma porção da onda de probabilidade atravessa o potencial embora classicamente não tenha energia suficiente para tal feito. Essa é uma demonstração numérica do efeito de tunelamento.

Potencial Triangular

O potencial triangular, por ter uma base mais larga, coíbe o efeito de tunelamento, no entanto, uma pequena porção da onda incidente ainda consegue transpor a barreira.

Potencial Batman

Em cursos básicos de quântica normalmente resolve-se apenas o potencial quadrado, às vezes, excepcionalmente, apresenta-se o potencial triangular. Isso se deve a dificuldade em tratar casos de potenciais irregulares analiticamente. Mas em casos reais os potenciais são extremamente irregulares. A solução numérica trata tais casos com naturalidade e sem qualquer alteração. Em particular, um potencial tipo batman é encarado sem problemas:

A parte superior da curva é apenas ilustrativa (pois não faz sentido matemático nem físico)