4.2.3 Transformada de Fourier Discreta (TFD)

Antes de iniciarmos a TFD revisaremos alguns itens de matemática e também de Fortran. Revisaremos um pouco de paridade de funções, operações de números complexos e como trabalhar com números complexos com o Fortran 90.

Paridade de uma Função

Uma função é dita par ou ímpar se:

par: $\displaystyle f(x) = f(-x)$

(4.2.9)

ímpar: $\displaystyle f(x) = -f(-x)$

(4.2.10)

A paridade de uma função define a simetria que é fica melhor representada na figura 4.3a e 4.3b.

**Figura 4.3:** Exemplo de uma (a) função par () e uma (b) função ímpar ().

Essa nomenclatura advem de:

$\displaystyle f(x)= x^k \begin{cases} \text{será par, } & \text{se k é um nú... ... par},\\ \text{será ímpar, } & \text{se k é um número ímpar}. \end{cases}$

(4.2.11)

Algumas propriedades das funções pares e ímpares:

Uma função é para e ímpar ao mesmo tempo de ;
Podemos ter uma função sem paridade (?) (gostaria de colocar um exemplo aqui);
Uma função ímpar definida na origem é nula na origem;
O resultado da soma de duas funções pares é par;
O resultado da some de duas funções ímpares é ímpar;
O resultado da multiplicação de par $\times$ par ou ímpar $\times$ ímpar é par;
O resultado da multiplicação de par $\times$ ímpar é ímpar;
O resultado da derivada de uma função par é ímpar;
O resultado da derivada de uma função ímpar é par;

Números complexos - básico

Os números complexos são elementos do conjunto complexo $\mathbb{C}$ e o conjunto $\mathbb{C}$ está contido no conjunto dos números reais $\mathbb{R}$ ^4.1. Quando se busca encontrar raízes de uma equação e a raiz apresenta um número $\sqrt{-1}$ , esse número é denominado como uma raiz complexa, assim $i=sqsrt{-1}$ e $i^2 = -1$ .

Um número complexo é denominado como:

$\displaystyle K = a + bi$

(4.2.12)

em que $K$ é um números complexo composto pela parte real $a$ e a parte imaginária $bi$ , sendo $b$ um número real.

Algumas operações com números complexos:

soma: ;
multiplicação: $K \cdot L = (a+bi) \cdot (c+di) = (ac-bd)+(bc+ad)i$ ;
multiplicação: $K \cdot \bar{K} = (a+bi) \cdot (a-bi) = a^2 + b^2$ , $\bar{K}$ é o complexo conjugado de ;
módulo: $\vert K\vert = r = \sqrt{a^2 + b^2}$

Na forma polar podemos escrever o número complexo $K$ como:

$\displaystyle K = r (cos\theta + i sen\theta)$

(4.2.13)

Pela igualdade de Euler temos que:

$\displaystyle e^{\pm i \theta} = cos \theta \pm i sen \theta$

(4.2.14)

Assim o número complexo $K$ pode ser escrito na forma:

$\displaystyle K = r e^{\pm i \theta}$

(4.2.15)

Números complexos - no Fortran 90

No Fortran 90 podemos trabalhar com números complexos de uma forma muito parecida como a que trabalhamos no lápis e papel. No programa a seguir temos dois exemplos de números complexos (um vetor complexo e uma variável complexa), como é escrita a parte real e a parte complexa (na verdade são escritas como dois números reais), como atribuir à um número complexo dois números reais, imprimir um número complexo, como utilizar a parte real, imaginária e o módulo do número complexo.

!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
! EXEMPLO DA UTILIZACAO DE NUMEROS COMPLEXOS
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> Outubro/2009
program num_complex2
implicit none
! DECLARACAO DE VARIAVEIS >>> INTEIROS
integer :: i, n
! DECLARACAO DE VARIAVEIS >>> REAIS
real(kind=8), parameter :: pi = 3.1415926535897967d0
real(kind=8) :: a, b, theta, soma_a, soma_b, passo
! DECLARACAO DE VARIAVEIS >>> COMPLEXOS
complex(kind=8), allocatable, dimension(:) :: K1
complex(kind=4) :: K2
n = 18
soma_a = 0.d0
soma_b = 0.d0
passo = 0.d0
allocate(K1(n))
write(*,*)'Vetor Complexo K1(i) >>> ( PARTE REAL , PARTE IMAGINARIA )'
do i = 1, n
        if ( i == 1 ) then
                passo = 0.0
        else
                passo = passo + 5.0
        end if
        theta = passo * pi / 180.0
        a = cos(theta)
        b = sin(theta)
        soma_a = soma_a + a
        soma_b = soma_b + b
        ! *** UM VETOR COMPLEXO E ESCRITO COMO DOIS NUMEROS REAIS ***
        !   *** VETOR_COMPLEXO = (PARTE REAL, PARTE IMAGINARIA) ***  
        ! escrevendo o numero complexo com os dois numeros reais     
        K1(i) = CMPLX(a,b)
        write(*,*)'Vetor Complexo K1(',i,') >>>',K1(i)
end do
soma_a = soma_a / real(n)
soma_b = soma_b / real(n)
! *** UM NUMERO COMPLEXO E ESCRITO COMO DOIS NUMEROS REAIS ***
!   *** NUMERO_COMPLEXO = (PARTE REAL, PARTE IMAGINARIA) ***  
K2 = CMPLX(soma_a,soma_b)
write(*,*)'>>>>>>>>>>>>><<<<<<<<<<<<'
write(*,*)'Numero complexo       K2 =',K2
write(*,*)'Parte Real de         K2 =',REAL(K2)
write(*,*)'Parte Imaginaria de   K2 =',AIMAG(K2)
write(*,*)'Complexo Conjugado de K2 =',CONJG(K2)
write(*,*)'Valor absoluto de     K2 =',CABS(K2) 
write(*,*)'>>>>>>>>>>>>><<<<<<<<<<<<'
write(*,*)'>>> TERMINOUT <<<'
deallocate(K1)
!
stop
end program num_complex2
!

Transformada de Fourier Discreta - TFD

A partir deste ponto iremos implementar como exercício da linguagem Fortran 90 uma TFD. Os princípios e propriedades de uma Transformada de Fourier não estão descritos aqui nesta apostila, para maiores detalhes sobre a parte matemática é importante que seja consultado o livro da referência [1] e as referências citadas pelo autor. Pode ser consultado também os livros de métodos matemáticos aplicados à Física entre outros bons livros disponíveis. Nesta parte simplesmente nos preocuparemos somente com a implementação e tradução do que é visto teoricamente para a uma linguagem computacional, que é o Fortran 90.

Assumiremos que temos uma função $f(t)$ , em que $t$ é uma variável independente e não necessariamente indica o tempo e que a integral 4.2.16 exista.

$\displaystyle g(\omega) = \int \limits_{-\infty}^{+\infty} e^{-i\omega t}f(t) dt$

(4.2.16)

O resultado da integral dada pela $g(\omega)$ é chamada de Transformada de Fourier de $f(t)$ . $\omega$ é a variável independente da Transformada de Fourier chamada de frequência angular. Se $t$ é uma variável espacial a variável $\omega$ será o número de onda $k=2\pi/\lambda$ , sendo $\lambda$ o comprimento de onda.

Entre as transformadas existe uma correspondência bi-unívoca entre as funções $f(t)$ e $g(\omega)$ tal que a equação 4.2.17 exista.

$\displaystyle f(t) = \frac{1}{2\pi} \int \limits_{-\infty}^{+\infty} e^{i\omega t}g(\omega) d\omega$

(4.2.17)

esta equação é chamada de Transformada de Fourier inversa.

Como vimos na seção de integração numérica, podemos resolver a integral que contém $f(t)$ e obter a Transformada de Fourier $g(\omega)$ . Neste caso como estamos discretizando o problema da integral para obter a Transformada de Fourier, dizemos que esta Transformada de Fourier numérica é a TFD dentro de um intervalo qualquer $[a,b]$ em que a $f(t)$ é integrável. Assim a equação 4.2.16 toma a forma da equação 4.2.18.

$\displaystyle g(\omega) = \Delta t \sum \limits_{j=1}^{n} e^{-i\omega_k t_j}f(t_j)$

(4.2.18)

Ao observarem a TFD (ou qualquer outra transformada de fourier) verão que teremos uma dependência de duas variáveis, uma em $t_j$ e outra em $\omega_k$ que estão uma no espaço real e outra no espaço dos complexos, respectivamente. A idéia da transformada de fourier é manter a bi-unicidade entre a $f(t_j)$ e $g(\omega_k)$ e para isso a quantidade $t_j$ deve ser igual a quantidade de $\omega_k$ , ou vice-verso. Essa igualdade de pontos nos diz que: se entre $t_i$ e $t_f$ temos 100 pontos, então entre $\omega_i$ e $\omega_f$ também teremos 100 pontos ou intervalos.

Vamos pegar como exemplo $t_i=-30$ e $t_f=30$ , se $\Delta t = 0.6$ então teremos 100 intervalos entre $t_i$ e $t_f$ . Como $\omega$ varia sempre de $-m\pi$ até $+m\pi$ , então $\Delta \omega = 2m\pi/100$ , onde $m$ é um número inteiro. Deve-se sempre explorar os valores de $\Delta t$ e $\Delta \omega$ para que tenhamos sempre um conjunto de pontos que deixe as curvas $f(t_j)$ e $g(\omega_k)$ suave.

Uma vez definido os intervalos, vemos que é fácil perceber que a TFD é uma multiplicação de matrizes em que $g(\omega_k)$ é uma matriz coluna, o somatório da exponencial é uma matriz bi-dimensional e a $f(t_j)$ é uma matriz coluna, assim definimos a matriz $W_{kj}$ pela expressão 4.2.19.

$\displaystyle W_{kj} = e^{-i\omega_k t_j}$

(4.2.19)

assim a equação 4.2.18 assume a forma matricial:

$\displaystyle g = W * f$

(4.2.20)

$\displaystyle \begin{bmatrix} g(\omega_1) \\ g(\omega_2) \\ \vdots \\ g(\omega_... ...atrix}\times \begin{bmatrix} f(t_1) \\ f(t_2) \\ \ddots \\ f(t_j) \end{bmatrix}$

(4.2.21)

A matriz representada pela equação 4.2.21 possui elementos que são números complexos, assim pela igualdade de Euler podemos escrever a exponencial conforme a equação 4.2.22.

$\displaystyle e^{-i\omega_k t_j} = cos(\omega_k t_j) - i sen(\omega_k t_j)$

(4.2.22)

a equação 4.2.18 e 4.2.21 pode ser escrita na forma da equação 4.2.23

$\displaystyle g(\omega) = \Delta t \sum \limits_{j=1}^{n} \left[ cos(\omega_k t_j) - i sen(\omega_k t_j) \right] f(t_j)$

(4.2.23)

sendo a parte real (eq. 4.2.24) e imaginária (eq. 4.2.25) dada por:

$\displaystyle Real\left( g(\omega) \right) = \Delta t \sum \limits_{j=1}^{n} cos(\omega_k t_j) f(t_j)$

(4.2.24)

$\displaystyle Imaginario\left( g(\omega) \right) = \Delta t \sum \limits_{j=1}^{n} sen(\omega_k t_j) f(t_j)$

(4.2.25)

A nossa atividade será calcular a TFD das equações 4.2.26 e 4.2.27.

$\displaystyle f(t) = e^{ -\frac{(t-t_0)^2}{2\sigma^2} }$

(4.2.26)

A equação 4.2.26 é a de uma gaussiana em que $t_0$ é o centro da gaussiana e $\sigma$ é a largura da gaussiana. A parte muito interessante aqui é alterar os centros e principalmente as larguras das gaussianas para ver o comportamento da sua TFD. Utilize os intervalos $-30 \le t \le 30$ e $-\pi \le \omega \le \pi$ .

$\displaystyle f(t) = sen(10t) \cdot e^{-t}$

(4.2.27)

A equação 4.2.27 é a de um oscilador harmônico amortecido, o intervalo no espaço real é $0 \le t \le 10$ e espaço complexo $-25\pi \le \omega \le 25\pi$ .

Atenção: para as duas equações encontre a TFD e faça 4 gráficos com xmgrace, sendo os gráficos da $f(t)$ , valor absoluto da $g(\omega)$ , da parte real de $g(\omega)$ e da parte imaginária de $g(\omega)$ .

Além da TFD apresentada acima, na literatura temos um outro métodos de calcular a TFD de uma maneira mais rápida e eficiente que é conhecida como a Fast Fourier Transform - FFT. Um dos algoritmos mais conhecidos é o do Cooley-Tukey^4.2 [9]. Para aprender mais sobre a FFT é bom verificar a referência [1] e também a mais completa/gratuita/on-line a referência [10].

A idéia desta FFT é calcular a TFD mais rápida manipulando os índices do somatório saindo de $k,j=1,\ldots,n$ e indo para $k,(2j, 2j+1)=1,\ldots, n/2$ . A equação da TFD (eq. 4.2.18) assume a forma da equação 4.2.28.

$\displaystyle g(\omega_k) = \sum \limits_{j=0}^{(n/2)-1} \left[ e^{-\frac{2\pi ... ...mits_{j=0}^{(n/2)-1} \left[ e^{-\frac{2\pi i}{n}} \right]^{k(2j+1)} f(t_{2j+1})$

(4.2.28)

Além de podermos implementar este código de FFT, podemos também utilizar uma sub-rotina já existente de FFT (que no fundo realiza a TFD) já otimizada e bem testada pelo usuários. Dentro do desenvolvimento computacional chamamos estas sub-rotinas de bibliotecas matemáticas (ou no inglês math library). As bibliotecas de FFT podem ser pegas para livre uso no site http://www.fftw.org^4.3. Para ilustrar a utilização da biblioteca FFTW3 vejamos o programa a seguir:

!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
! COMO COMPILAR: gfortran -o fftw3f90.x fftw3f90.f90 -lfftw3
!
! NA LINHA: include "/usr/include/fftw3.f" DEVE SER ALTERADO O CAMINHO
! ONDE SE ENCONTRA A BIBLIOTECA/SUBROTINA fftw3.f
!
! ESTE PROGRAMA FOI ADAPTADO DA VERSAO FORTRAN 77 DO SITE
! http://www.psc.edu/general/software/packages/fftw/examples/index.php
! CUJA LICENSA E GERIDA PELA GNU General Public License
!
! ATENCAO: VARIAVEIS REAIS EM KIND=8 E ALGUMAS INTEIRAS COM KIND=8 (8 BITS)
!
! PROGRAMADOR: Fernando Sato
! E_MAIL: fsato01@gmail.com
! DATA: Outubro/2009
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
program fftw3f90
implicit none

! INCLUINDO O ARQUIVO QUE CONTEM AS DECLARACOES DA FFTW3
include "/usr/include/fftw3.f"

! DECLARACAO DE VARIAVEIS
! UTILIZAR >>> n2 = n/2 + 1
integer, parameter :: n = 60, n2 = 31
real(kind=8), parameter :: pi = 3.1415926535897967d0
real(kind=8), dimension(n) :: a, b 
real(kind=8) :: sigma, a0, abslt, pr, pim
integer(kind=8) :: plan_forward, plan_backward
real(kind=8), dimension(n)     :: in, out2
complex(kind=8), dimension(n2) :: out
integer :: i

! ARQUIVOS DE SAIDA
open(UNIT=30,FILE="fftw_fx_out.dat",STATUS='REPLACE',ACTION='WRITE')
open(UNIT=31,FILE="fftw_tfd_out.dat",STATUS='REPLACE',ACTION='WRITE')
open(UNIT=32,FILE="fftw_ttfd_out.dat",STATUS='REPLACE',ACTION='WRITE')

! CONTANTES DA GAUSSIANA
sigma = 10.d0
a0 = 0.d0

!>>>>>>>>>>>>>>>>>> ATRIBUINDO VALORES
do i = 1, n

	! ENTRADA X
	a(i) = dfloat(i) - dfloat(n2)
	
	! VALOR DA F(X)
	!b(i) = exp( - (a(i)-a0)**2.d0  / (2.d0*sigma*sigma) )
	b(i) = exp(-a(i)**2.d0)

	! ESCREVENDO X E F(X)
	write(30,500) a(i),b(i)

	! ENTRADA P/ FFT
	in(i) = b(i)
end do

! FAZENDO TRANSFORMADA CALCULANDO g(w)
!>>>>>>>>>>>>>>>>>> PASSANDO VARIAVEIS P/ FFT IN:REAL OUT:COMPLEX
call dfftw_plan_dft_r2c_1d ( plan_forward, n, in, out, FFTW_ESTIMATE )

!>>>>>>>>>>>>>>>>>> EXECUTANDO FFT
call dfftw_execute( plan_forward )

!>>>>>>>>>>>>>>>>>> SAIDA DA FFT OUT:COMPLEX
do i=n2,1,-1
	pr = REAL(out(i))
	pim = AIMAG(out(i))
	abslt = DSQRT(pr*pr + pim*pim)
	write(31,501) -i, pr/dfloat(n2), pim/dfloat(n2), abslt/dfloat(n2)
end do
do i = 1,n2
	pr = REAL(out(i)) 
	pim = AIMAG(out(i))
	abslt = DSQRT(pr*pr + pim*pim)
	write(31,501) i, pr/dfloat(n2), pim/dfloat(n2), abslt/dfloat(n2)
end do
!>>>>>>>>>>>>>>>>>> ENCERRA A TRANSFORMADA DE FOURIER
call dfftw_destroy_plan ( plan_forward )

! FAZENDO TRANSFORMADA DA TRANSFORMADA CALCULANDO f(t)
!>>>>>>>>>>>>>>>>>> PASSANDO VARIAVEIS P/ FFT IN:COMPLEX OUT:REAL
call dfftw_plan_dft_c2r_1d ( plan_backward, n, out, out2, FFTW_ESTIMATE )

!>>>>>>>>>>>>>>>>>> EXECUTANDO A TRANSFORMADA DA TRANSFORMADA DE FOURIER
call dfftw_execute( plan_backward )

!>>>>>>>>>>>>>>>>>> SAIDA: TRANFORMADA DA TRANSFORMADA DE FOURIER
do i = 1,n
	write(32,502) i, out2(i)/dfloat(n)
end do
!>>>>>>>>>>>>>>>>>> ENCERRA A TRANSFORMADA DE FOURIER
call dfftw_destroy_plan ( plan_backward )

close(30)
close(31)
close(32)

!>>>>>>>>>>>>>>>>>> FORMATOS DE ESCRITA
500 format(f18.8,e18.8)
501 format(i6,1x,e18.8,e18.8,e18.8)
502 format(i6,1x,e18.8)
!
stop
end program fftw3f90

Somente lembrando que este código é um exemplo de aplicação das bibliotecas FFTW [11] e não há nenhuma garantia de que este código dará resultados em ambientes de pesquisa ou em qualquer outro ambiente, de qualquer forma o autor não se responsabiliza pelos efeitos colaterais advindos do código acima.