O módulo scipy.fftpack permite calcular transformada rápida de Fourier. Como ilustração, um sinal de entrada (barulhento) pode parecer:
[pastacode lang=”python” message=”” highlight=”” provider=”manual”]
time_step = 0.02
period = 5.
time_vec = np.arange(0, 20, time_step)
sig = np.sin(2 * np.pi / period * time_vec) + \
0.5 * np.random.randn(time_vec.size)[/pastacode]
O observador não sabe a frequência do sinal, apenas o passo do tempo de amostragem do sinal sig. O sinal é suposto vindo de uma função real, de modo que a transformada de Fourier seja simétrica. A função scipy.fftpack.fftfreq() irá gerar as frequências de amostragem e scipy.fftpack.fft() irá calcular a transformada rápida de Fourier:
[pastacode lang=”python” message=”” highlight=”” provider=”manual”]
from scipy import fftpack
sample_freq = fftpack.fftfreq(sig.size, d=time_step)
sig_fft = fftpack.fft(sig)[/pastacode]
Pelo fato de a energia resultante ser simétrica, apenas a parte positiva do espectro deve ser usada para encontrar a frequência:
[pastacode lang=”python” message=”” highlight=”” provider=”manual”]
pidxs = np.where(sample_freq > 0)
freqs = sample_freq[pidxs]
power = np.abs(sig_fft)[pidxs][/pastacode]
Para plotar com o Matplotlib:
[pastacode lang=”python” message=”” highlight=”” provider=”manual”]
pl.figure()
pl.plot(freqs, power)
pl.xlabel('Frequencia [Hz]')
pl.ylabel('Energia')
axes = pl.axes([0.3, 0.3, 0.5, 0.5])
pl.title('Pico de frequencia')
pl.plot(freqs[:8], power[:8])
pl.setp(axes, yticks=[])[/pastacode]
A frequência do sinal pode ser encontrada com:
[pastacode lang=”python” message=”” highlight=”” provider=”manual”]
freq = freqs[power.argmax()]
np.allclose(freq, 1./period) # checa se aquele frequência correta é encontrada[/pastacode]
Agora, o ruído de alta frequência será removido a partir do sinal da transformada de Fourier:
[pastacode lang=”python” message=”” highlight=”” provider=”manual”]
sig_fft[np.abs(sample_freq) > freq] = 0[/pastacode]
O sinal filtrado resultante pode ser calculado pela função scipy.fftpack.ifft():
[pastacode lang=”python” message=”” highlight=”” provider=”manual”]
ain_sig = fftpack.ifft(sig_fft)[/pastacode]
E o resultado pode ser visualizado com:
[pastacode lang=”python” message=”” highlight=”” provider=”manual”]
import pylab as plt
plt.figure()
plt.plot(time_vec, sig)
plt.plot(time_vec, main_sig, linewidth=3)
plt.xlabel('Tempo [s]')
plt.ylabel('Amplitude')[/pastacode]
O Numpy também tem uma implementação de FFT (numpy.fft). No entanto, em geral, deve ser preferido a do scipy, uma vez que utiliza implementações mais eficientes.