Exercício 08 - Introdução Transformada Discreta de Fourier

Autor: Wagner Machado do Amaral
Data: 11/05/2009

Enunciado: ia636-2009:exercicio8

A função iacos implementa uma equação semelhante à equação da Onda Senoidal.

Equação da Onda Senoidal:

  • (y) Alongamento: posição da onda no tempo t
  • (A) Amplitude: alongamento máximo
  • (w) Freqüência angular:
  • (phi) Fase: Posição inicial da onda
  • (D) offset vertical: Posição central da onda
  • (f) Freqüência: Número de voltas realizadas em uma unidade de tempo
  • (T) Período: Tempo decorrido em uma volta completa

Sendo:

A equação da Onda Senoidal pode ser escrita também da seguinte forma:

A equação acima pode ser implementada da seguite forma:

1 phi=0.
2 f=1/50.
3 A=100.
4 D = 125
5 cols, rows = 256, 128
6 
7 t = arange(cols)
8 g = A*sin( (2*pi*f)*t + phi ) + D
9 mmplot([[g]],['set yrange ['+str(D-A)+':'+str(D+A)+']','set xrange [0:'+str(cols)+']'], title='Gráfico da onda senoidal')

Gráfico da onda senoidal

Considerando que a amplitude da onda não ultrapasse o limite de 0 a 255, a onda senoidal pode ser exibida em uma imagem em tons de cinza. Conforme abaixo:

1 x, y = iameshgrid(range(cols),range(rows))
2 g = A*sin( (2*pi*f)*x + phi ) + D
3 adshow(g,title='Onda senoidal em imagem')

Onda senoidal em imagem

Equação implementada por iacos:

Essa equação utiliza um parâmetro theta que determina um grau de inclinação na onda. Dessa forma a equação não é em função apenas de x, mas também em y.

Nos exemplos de utilização da função iacos os dados são normalizados entre 0 e 255, dessa forma não é necessário informar a amplitude A e o offset vertical D. Porém se os dados não forem normalizados a onda não ficara visível.

Na implementação de iacos a onda gerada é uma onda cosseno, embora a equação seja de uma onda seno. O resultado é o mesmo exceto por uma diferença no tempo, semelhante a alterar phi.

 1 theta=pi/4
 2 phi=0.
 3 T=50.
 4 cols, rows = 256, 128
 5 
 6 x, y = iameshgrid(range(cols),range(rows))
 7 fx = cos(theta)/T
 8 fy = sin(theta)/T
 9 g_sin = sin(2*pi*(fx*x + fy*y) + phi)
10 g_cos = cos(2*pi*(fx*x + fy*y) + phi)
11 adshow(ianormalize(g_sin, [0,255]),title='Onda seno em imagem')
12 adshow(g_sin,title='Onda seno em imagem não normalizada')
13 adshow(ianormalize(g_cos, [0,255]),title='Onda cosseno em imagem')

Onda seno em imagem

Onda seno em imagem não normalizada

Onda cosseno em imagem


Com as funções iadft e iadftview é possível calcular a Transformada de Fourier e exibir o espectro através de uma imagem em função da intensidade. Conforme exemplos abaixo:

 1 f = 255 * iacircle([256,256], 10, [129,129])
 2 iashow(f)
 3 F = iadft(f)
 4 Fv = iadftview(F)
 5 iashow(Fv)
 6 
 7 f = iarectangle([256,256], [10,10], [129,129])
 8 iashow(f)
 9 F = iadft(f)
10 Fv = iadftview(F)
11 iashow(Fv)
12 
13 f = iarectangle([256,256], [1,100], [129,129])
14 iashow(f)
15 F = iadft(f)
16 Fv = iadftview(F)
17 iashow(Fv)
18 
19 f = iarectangle([256,256], [100,1], [129,129])
20 iashow(f)
21 F = iadft(f)
22 Fv = iadftview(F)
23 iashow(Fv)
Warning: downcasting image from int32 to uint16 (may lose precision)


No exemplo abaixo é calculada a Transformada de Fourier de uma função utilizando a função iadft.

Para melhor visualização as componentes de freqüência zero são colocadas no centro do espectro através da função iafftshift.

Finalmente a função é reconstruída pela operação inversa, utilizando a função iaidft.

 1 f = arange(10)*0
 2 f[3:6]=1
 3 
 4 print"\n Funcao: \n", f
 5 mmplot([[f]],['set yrange [-2:2]'],ptitle='Funcao')
 6 
 7 F = iadft(f)
 8 
 9 print"\n Expectro da Transformada de Fourier: \n", F
10 mmplot([[iafftshift(F)]], ptitle='Expectro da Transformada de Fourier')
11 
12 f2 = iaidft(F)
13 print"\n Funcao reconstruida: \n", f2
14 mmplot([[f2]],['set yrange [-2:2]'],ptitle='Funcao reconstruida')
Funcao: 
[0 0 0 1 1 1 0 0 0 0]

 Expectro da Transformada de Fourier: 
[[ 3.00000000 +0.00000000e+00j]
 [-2.11803399 -1.53884177e+00j]
 [ 0.50000000 +1.53884177e+00j]
 [ 0.11803399 -3.63271264e-01j]
 [ 0.50000000 -3.63271264e-01j]
 [-1.00000000 -2.77555756e-16j]
 [ 0.50000000 +3.63271264e-01j]
 [ 0.11803399 +3.63271264e-01j]
 [ 0.50000000 -1.53884177e+00j]
 [-2.11803399 +1.53884177e+00j]]

 Funcao reconstruida: 
[[  2.08616256e-08 -1.75541673e-17j]
 [  2.04326230e-08 +7.02166694e-17j]
 [  2.21428221e-09 +2.10650008e-16j]
 [  9.99999980e-01 +2.45758343e-16j]
 [  9.99999970e-01 -7.92555042e-17j]
 [  9.99999979e-01 -6.31950024e-16j]
 [ -3.09924210e-09 -8.42600032e-16j]
 [  7.33157063e-09 -3.86191682e-16j]
 [  1.00018097e-08 -2.10650008e-16j]
 [  1.26450945e-08 +9.56470983e-17j]]

Veja abaixo a reconstrução passo a passo da função acima, iniciando das senoides de menor frequencia:

 1 F_array = [
 2 [F[0],[0],[0],[0],[0],[0],[0],[0],[0],[0]],
 3 [F[0],F[1],[0],[0],[0],[0],[0],[0],[0],F[9]],
 4 [F[0],F[1],F[2],[0],[0],[0],[0],[0],F[8],F[9]],
 5 [F[0],F[1],F[2],F[3],[0],[0],[0],F[7],F[8],F[9]],
 6 [F[0],F[1],F[2],F[3],F[4],[0],F[6],F[7],F[8],F[9]],
 7 [F[0],F[1],F[2],F[3],F[4],F[5],F[6],F[7],F[8],F[9]]
 8 ]
 9 
10 for i in range(6):
11   f2 = iaidft(F_array[i])
12   mmplot([[f2]],['set yrange [-2:2]'],ptitle='Reconstruindo a funcao (Passo: '+str(i+1)+') ')