Calcular y visualizar FFTs fácilmente en Python usando UliEngineering

UliEngineering es una biblioteca de análisis de datos mixtos en Python - una de las utilidades que proporciona es un paquete fácil de usar para calcular FFTs. A diferencia de otros paquetes, esta biblioteca está orientada hacia casos de uso prácticos y te permite hacer la FFT en ¡solo una línea de código! No se requiere conocimiento de matemáticas. Primero, instala UliEngineering.

Introducción

Primero, generaremos algunos datos de prueba. Ver este post anterior para más detalles sobre cómo generar datos de prueba sinusoidales:

from UliEngineering.SignalProcessing.Simulation import *

# Generar datos de prueba: tono 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)

Los datos de prueba consisten en una onda seno de 100 Hz más una onda seno de 400 Hz (con la mitad de la amplitud). Esa señal se muestrea con una tasa de muestreo de 1000 Hz.

Ahora podemos calcular y visualizar la FFT usando matplotlib:

# Calcular FFT. Usar la misma tasa de muestreo
# NOTA: ¡La ventana por defecto es "blackman"!
from UliEngineering.SignalProcessing.FFT import compute_fft
fft = compute_fft(data, samplerate=1e3)

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

plt.gcf().set_size_inches(10, 5) # Usar (20, 10) para obtener una gráfica más grande
plt.plot(fft.frequencies, fft.amplitudes)
plt.xlabel("Frecuencia")
plt.ylabel("Amplitud")

compute_fft(data, samplerate=1e3) devuelve un objeto FFT que contiene campos como frequencies, amplitude y phase. Realiza una FFT del tamaño de la entrada (es decir, como data es un array de longitud 1000, la FFT será de tamaño 1000).

fft.frequencies es un array de frecuencias (en Hz), correspondiente a los valores en fft.amplitudes. También puedes usar fft.angles para obtener ángulos relativos en grados, pero eso no se cubre en este blog post.

Como puedes ver en la gráfica mostrada arriba, la frecuencia máxima que se puede detectar es siempre la mitad de la tasa de muestreo, es decir, para nuestra tasa de muestreo de $f_s = 1000,\text{Hz}$ es $500,\text{Hz}$. Ver esta explicación de FFT más centrada en matemáticas si quieres conocer más detalles.

Internamente, compute_fft() realiza este cálculo:

• $2$ es un factor de corrección que tiene en cuenta que descartamos la segunda mitad del resultado bruto de la FFT (ya que estamos haciendo una FFT real)
• $\frac{1}{\text{len(data)}}$ normaliza los resultados de la FFT para que sean independientes de la longitud de los datos (es decir, si pasas una muestra más larga de la misma onda seno seguirás obteniendo el mismo resultado
• $\text{abs}\left(\cdots\right)$ Convierte el resultado complejo consciente de la fase de la FFT a un espectro que es más fácil de leer y visualizar.
• Window (que por defecto es blackman) es la ventana que se aplica a los datos para aliviar algunos efectos matemáticos al principio y al final del conjunto de datos. Ver wikipedia sobre funciones de ventana para más detalles. UliEngineering actualmente ofrece esta lista de funciones de ventana:
  • blackman
  • bartlett
  • hamming
  • hanning
  • kaiser (El parámetro está fijado a 2.0)
  • none

Seleccionar rangos de frecuencia

Usando la API de UliEngineering, seleccionar un rango de frecuencia de la FFT es trivialmente fácil: Solo usa fft[lowfreq:highfreq]. Puedes usar fft[lowfreq:] para seleccionar todo empezando desde lowfreq o usar fft[:highfreq] para seleccionar todo hasta highfreq.

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

# Seleccionar rango de frecuencia: 50 a 200 Hz
fft = fft[50.0:200.0]

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

plt.gcf().set_size_inches(10, 5) # Usar (20, 10) para obtener una gráfica más grande
plt.plot(fft.frequencies, fft.amplitudes)
plt.xlabel("Frecuencia")
plt.ylabel("Amplitud")
plt.savefig("/ram/fft-frequency-range.svg")

Ventana de frecuencia

Extraer amplitud y ángulo a una cierta frecuencia

Usando [frequency], es decir, el operador getitem con un solo valor, obtienes un objeto FFTPoint() que contiene la frecuencia, amplitud y ángulo relativo para una frecuencia dada. La biblioteca selecciona automáticamente el bucket FFT más cercano, así que incluso si tu FFT no tiene un bucket para esa frecuencia específica, obtendrás resultados sensatos.

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

# Mostrar valor a una cierta frecuencia
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 cortas para datos largos

Usando compute_fft(), si tenemos un array de datos extremadamente largo, esto significa que calcularemos una FFT extremadamente larga. En muchos casos, esto no es deseable y quieres calcular una FFT de tamaño fijo (más a menudo una FFT de potencia de dos, p. ej. 1024, 2048, 4096 etc).

Generemos algunos datos de prueba largos y asumamos que queremos calcular una FFT de tamaño 1024 en ellos

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

Trazar los datos como arriba produce

simple_serial_fft_reduce() se encarga de toda la magia y normalización por nosotros, incluyendo la partición de los datos en fragmentos superpuestos, la suma de las FFTs y la normalización adecuada del resultado.

La convención de nomenclatura es significativa aquí:

• simple_...._reduce significa que esta es una variante con valores por defecto sensatos (simple) para usar una función de reduction (por defecto: sum) en múltiples FFTs.
• serial significa las FFTs individuales

Paralelizar las FFTs

Si tienes un conjunto de datos enorme, puedes usar simple_parallel_fft_reduce() idénticamente a simple_serial_fft_reduce():

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

Sin embargo en la mayoría de los casos querrás inicializar el executor manualmente para poder reutilizarlo más tarde:

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

Podemos usar un ThreadPoolExecutor() ya que scipy.fftpack (que UliEngineering usa para hacer los cálculos difíciles) libera el GIL de Python.

Ten en cuenta que debido a la necesidad de hacer muchas tareas de mantenimiento, simple_parallel_fft_reduce() es mucho más lenta que simple_serial_fft_reduce() si tienes un conjunto de datos tan pequeño que la paralelización no es efectiva. Mi recomendación inicial es considerar usar la variante paralela si el tiempo total de ejecución de la variante serial es mayor que $0.5s$