In [2]:
import numpy as np
import matplotlib.pyplot as plt
In [48]:
#funções
def padFFT(x):
x = np.array(x, dtype=np.float)
len_x = len(x)
N = 1
if (len_x % 2 != 0):
while N < len_x:
N *= 2
x_new = np.zeros((N, 1))
for i in range(0, len_x):
x_new[i] = x[i]
x = x_new
return x
#transformação discreta de Fourier no modo convencional
def DFT(x):
N = x.shape[0] #quantidade de elementos
n = np.arange(N)
k = n.reshape((N, 1)) #transposição de n
M = np.exp(-2j * np.pi * k * n / N)
return np.dot(M, x)
#def FFT(x):
#fonte https://jakevdp.github.io/blog/2013/08/28/understanding-the-fft/
def FFT(x):
"""A vectorized, non-recursive version of the Cooley-Tukey FFT"""
x = np.asarray(x, dtype=float)
N = x.shape[0]
if np.log2(N) % 1 > 0:
x = padFFT(x)
N = x.shape[0]
# N_min here is equivalent to the stopping condition above,
# and should be a power of 2
N_min = min(N, 32)
# Perform an O[N^2] DFT on all length-N_min sub-problems at once
n = np.arange(N_min)
k = n[:, None]
M = np.exp(-2j * np.pi * n * k / N_min)
X = np.dot(M, x.reshape((N_min, -1)))
# build-up each level of the recursive calculation all at once
while X.shape[0] < N:
X_even = X[:, :X.shape[1] / 2]
X_odd = X[:, X.shape[1] / 2:]
factor = np.exp(-1j * np.pi * np.arange(X.shape[0])
/ X.shape[0])[:, None]
X = np.vstack([X_even + factor * X_odd,
X_even - factor * X_odd])
return X.ravel()
#calcular as séries de fourier
def serieFourier(A_n, n, x):
len_x = len(x)
L = (x[len_x - 1] - x[0]) / 4
y = np.zeros((len_x, 1))
for i in range(0, len_x):
y[i] = np.sum(A_n * np.sin(n * x[i] * np.pi / L))
return [x, y]
def transfFourier(x, y, fft=1):
dx = x[1] - x[0]
len_x = len(x)
p = np.arange(0, len_x)
len_p = len(p)
k_p = 2 * np.pi * p / (dx * len_x)
Y = np.zeros((len_p, 1), dtype=np.complex_)
if fft == 0:
Y = FFT(y)
k_p = padFFT(k_p)
elif fft == 1:
Y = DFT(y)
elif fft == 2:
Y = np.fft.rfft(y)
return [k_p, Y]
def plotSerieFourier(x, y):
plt.plot(x, y, '-+')
plt.xlabel('x')
plt.ylabel('y')
plt.title('Séries de Fourier')
plt.show()
def plotTransfFourier(k_p, Y):
Ysq = np.square(np.abs(Y))
plt.plot(k_p, Ysq, '-+')
plt.xlabel('$k_p$')
plt.ylabel('$|Y(k_p)^2|$')
plt.title("Transformada de Fourier discreta")
plt.show()
In [49]:
#condições iniciais
x0 = 0
xf = 2.00
dx = 0.01
n = np.array([1])
A_n = np.array([1.0])
#---
L = xf - x0
x = np.arange(x0, xf+dx, dx)
len_x = len(x)
ser_fourier = serieFourier(A_n, n, x)
trans_fourier = transfFourier(*ser_fourier, fft=1)
In [50]:
(x, y) = ser_fourier
In [51]:
(k_p, Y) = trans_fourier
In [52]:
plt.plot(x, y)
plt.show()
In [53]:
plt.plot(k_p, Y)
plt.xlim(0, 200)
plt.show()
In [54]:
plt.plot(k_p, np.fft.fft(y))
plt.xlim(0, 200)
plt.show()
In [57]:
x = np.linspace(0, 2*np.pi, 1e2)
y = np.sin(x)
In [ ]:
plt.plot(np.fft.fft(y))
plt.xlim(0, 200)
plt.show()
In [ ]: