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.