Calcule e visualize facilmente FFTs em Python usando UliEngineering

UliEngineering é uma biblioteca de análise de dados mista em Python - uma das utilidades que ela fornece é um pacote fácil de usar para calcular FFTs. Em contraste com outros pacotes, esta biblioteca é orientada para casos de uso práticos e permite que você faça a FFT em apenas uma linha de código! Nenhum conhecimento de Matemática necessário. Primeiro, instale UliEngineering.

Começando

Primeiro, geraremos alguns dados de teste. Veja este post anterior para mais detalhes sobre como gerar dados de teste sinusoidais:

from UliEngineering.SignalProcessing.Simulation import *

# Gerar dados de teste: tom de 100 Hz + 400 Hz
data = sine_wave(frequency=100.0, samplerate=1000, amplitude=1.) \
       + sine_wave(frequency=400.0, samplerate=1000, amplitude=0.5)

Os dados de teste consistem de uma onda senoidal de 100 Hz mais uma onda senoidal de 400 Hz (com metade da amplitude). Aquele sinal é amostrado com uma taxa de amostragem de 1000 Hz.

Agora podemos calcular & visualizar a FFT usando matplotlib:

# Calcular FFT. Use a mesma taxa de amostragem
# NOTA: Janela padrão é "blackman"!
from UliEngineering.SignalProcessing.FFT import compute_fft
fft = compute_fft(data, samplerate=1e3)

# Plotar
from matplotlib import pyplot as plt
plt.style.use("ggplot")

plt.gcf().set_size_inches(10, 5) # Use (20, 10) para obter um plot maior
plt.plot(fft.frequencies, fft.amplitudes)
plt.xlabel("Frequência")
plt.ylabel("Amplitude")

compute_fft(data, samplerate=1e3) retorna um objeto FFT contendo campos como frequencies, amplitude e phase. Ele realiza uma FFT do tamanho da entrada (ou seja, já que data é um array de comprimento 1000, a FFT será de tamanho 1000).

fft.frequencies é um array de frequências (em Hz), correspondendo aos valores em fft.amplitudes. Você também pode usar fft.angles para obter ângulos relativos em graus, mas isso não está sendo coberto neste post de blog.

Como você pode ver no plot mostrado acima, a frequência máxima que pode ser detectada é sempre metade da taxa de amostragem, ou seja, para nossa taxa de amostragem de $f_s = 1000,\text{Hz}$ é $500,\text{Hz}$. Veja esta explicação de FFT mais centrada em matemática se você quiser saber mais detalhes.

Internamente, compute_fft() realiza este cálculo:

• $2$ é um fator de correção que leva em conta que descartamos a última metade do resultado bruto da FFT (já que estamos fazendo uma FFT real)
• $\frac{1}{\text{len(data)}}$ normaliza os resultados da FFT para que sejam independentes do comprimento dos dados (ou seja, se você passar uma amostra mais longa da mesma onda senoidal você ainda obterá o mesmo resultado
• $\text{abs}\left(\cdots\right)$ Converte o resultado complexo consciente de fase da FFT para um espectro que é mais fácil de ler & visualizar.
• Window (que padrão é blackman) é a janela que é aplicada aos dados para aliviar alguns efeitos matemáticos no início e no final do conjunto de dados. Veja wikipedia sobre funções de janela para mais detalhes. UliEngineering atualmente oferece esta lista de funções de janela:
  • blackman
  • bartlett
  • hamming
  • hanning
  • kaiser (Parâmetro é fixado em 2.0)
  • none

Selecionando faixas de frequência

Usando a API UliEngineering, selecionar uma faixa de frequência da FFT é trivialmente fácil: Apenas use fft[lowfreq:highfreq]. Você pode usar fft[lowfreq:] para selecionar tudo começando de lowfreq ou usar fft[:highfreq] para selecionar tudo até highfreq.

from UliEngineering.SignalProcessing.FFT import compute_fft
fft = compute_fft(data, samplerate=1e3)

# Selecionar faixa de frequência: 50 a 200 Hz
fft = fft[50.0:200.0]

# Plotar
from matplotlib import pyplot as plt
plt.style.use("ggplot")

plt.gcf().set_size_inches(10, 5) # Use (20, 10) para obter um plot maior
plt.plot(fft.frequencies, fft.amplitudes)
plt.xlabel("Frequência")
plt.ylabel("Amplitude")
plt.savefig("/ram/fft-frequency-range.svg")

Janela de frequência

Extrair amplitude & ângulo em uma certa frequência

Usando [frequency], ou seja, operador getitem com um único valor, você obtém um objeto FFTPoint() contendo a frequência, amplitude e ângulo relativo para uma dada frequência. A biblioteca automaticamente seleciona o bucket FFT mais próximo, então mesmo se sua FFT não tem um bucket para aquela frequência específica, você obterá resultados sensatos.

from UliEngineering.SignalProcessing.FFT import compute_fft
fft = compute_fft(data, samplerate=1e3)

# Mostrar valor em uma certa frequência
print(fft[30]) # FFTPoint(frequency=30.0, value=2.91e-08, angle=0.0)
print(fft[100]) # FFTPoint(frequency=100.0, value=0.419, angle=0.0)

FFTs curtas para dados longos

Usando compute_fft(), se temos um array de dados extremamente longo, isso significa que calcularemos uma FFT extremamente longa. Em muitos casos, isso não é desejável e você quer calcular uma FFT de tamanho fixo (mais frequentemente uma FFT potência-de-dois, ex.: 1024, 2048, 4096 etc).

Vamos gerar alguns dados de teste longos e assumir que queremos calcular uma FFT de tamanho 1024 neles

from UliEngineering.SignalProcessing.FFT import *
fft = simple_serial_fft_reduce(data, samplerate=1e3, fftsize=1024)
# fft.frequencies, fft.amplitudes etc

Plotar os dados como acima produz

que parece quase exatamente como nosso plot compute_fft() antes - exatamente como esperávamos.

simple_serial_fft_reduce() cuida de toda a mágica e normalização para nós, incluindo particionar os dados em chunks sobrepostos, adicionar as FFTs e normalizar propriamente o resultado.

A convenção de nomenclatura é significativa aqui:

• simple_...._reduce significa que esta é uma variante com padrões sensatos (simple) para usar uma função de redução (padrão: sum) em múltiplas FFTs.
• serial significa as FFTs individuais

Paralelizando as FFTs

Se você tem um conjunto de dados enorme, você pode usar simple_parallel_fft_reduce() identicamente a simple_serial_fft_reduce():

from UliEngineering.SignalProcessing.FFT import *
fft = simple_parallel_fft_reduce(data, samplerate=1e3, fftsize=1024)
# Use fft.frequencies, fft.amplitudes etc

No entanto na maioria dos casos você quer inicializar o executor manualmente para poder reutilizá-lo depois:

from UliEngineering.SignalProcessing.FFT import *
from concurrent.futures import ThreadPoolExecutor
executor = ThreadPoolExecutor() # Sem argumento => use num_cpus threads
fft = simple_parallel_fft_reduce(data, samplerate=1e3, fftsize=1024, executor=executor)

Podemos usar um ThreadPoolExecutor() já que scipy.fftpack (que UliEngineering usa para fazer a matemática pesada) desbloqueia o GIL Python.

Note que devido à necessidade de fazer muitas tarefas de housekeeping, simple_parallel_fft_reduce() é muito mais lento que simple_serial_fft_reduce() se você tem um conjunto de dados tão pequeno que a paralelização não é efetiva. Minha recomendação inicial é considerar usar a variante paralela se o tempo total de execução da variante serial é maior que $0.5s$