Calculer et visualiser facilement des FFT en Python avec UliEngineering

UliEngineering est une bibliothèque d’analyse de données mixte en Python — l’un des utilitaires qu’elle fournit est un paquet facile à utiliser pour calculer des FFT. Contrairement à d’autres paquets, cette bibliothèque est orientée vers des cas d’usage pratiques et vous permet de faire la FFT en une seule ligne de code ! Aucune connaissance en mathématiques requise. Tout d’abord, installez UliEngineering.

Pour commencer

D’abord, nous allons générer quelques données de test. Voir ce billet précédent pour plus de détails sur la façon de générer des données de test sinusoïdales :

from UliEngineering.SignalProcessing.Simulation import *

# Générer des données de test : ton 100 Hz + 400 Hz
data = sine_wave(frequency=100.0, samplerate=1000, amplitude=1.) \
       + sine_wave(frequency=400.0, samplerate=1000, amplitude=0.5)

Les données de test consistent en une onde sinusoïdale de 100 Hz plus une onde sinusoïdale de 400 Hz (avec la moitié de l’amplitude). Ce signal est échantillonné avec une fréquence d’échantillonnage de 1000 Hz.

Nous pouvons maintenant calculer et visualiser la FFT avec matplotlib :

# Calculer la FFT. Utiliser la même fréquence d'échantillonnage
# NOTE : la fenêtre par défaut est "blackman" !
from UliEngineering.SignalProcessing.FFT import compute_fft
fft = compute_fft(data, samplerate=1e3)

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

plt.gcf().set_size_inches(10, 5) # Utiliser (20, 10) pour un graphique plus grand
plt.plot(fft.frequencies, fft.amplitudes)
plt.xlabel("Frequency")
plt.ylabel("Amplitude")

compute_fft(data, samplerate=1e3) renvoie un objet FFT contenant des champs tels que frequencies, amplitude et phase. Elle effectue une FFT de la taille de l’entrée (c.-à-d. comme data est un tableau de longueur 1000, la FFT sera de taille 1000).

fft.frequencies est un tableau de fréquences (en Hz), correspondant aux valeurs dans fft.amplitudes. Vous pouvez aussi utiliser fft.angles pour obtenir des angles relatifs en degrés, mais cela n’est pas couvert dans ce billet.

Comme vous pouvez le voir dans le graphique ci-dessus, la fréquence maximale qui peut être détectée est toujours la moitié de la fréquence d’échantillonnage, c.-à-d. pour notre fréquence d’échantillonnage de $f_s = 1000,\text{Hz}$, il s’agit de $500,\text{Hz}$. Voir cette explication des FFT plus axée sur les mathématiques si vous souhaitez connaître plus de détails.

En interne, compute_fft() effectue ce calcul :

• $2$ est un facteur de correction qui prend en compte le fait que nous jetons la seconde moitié du résultat brut de la FFT (puisque nous faisons une FFT réelle)
• $\frac{1}{\text{len(data)}}$ normalise les résultats de la FFT afin qu’ils soient indépendants de la longueur des données (c.-à-d. si vous passez un échantillon plus long de la même onde sinusoïdale, vous obtiendrez toujours le même résultat)
• $\text{abs}\left(\cdots\right)$ convertit le résultat complexe tenant compte de la phase de la FFT en un spectre qui est plus facile à lire et à visualiser.
• Window (qui par défaut est blackman) est la fenêtre appliquée aux données pour atténuer certains effets mathématiques au début et à la fin du jeu de données. Voir Wikipedia sur les fonctions de fenêtrage pour plus de détails. UliEngineering propose actuellement cette liste de fonctions de fenêtrage :
  • blackman
  • bartlett
  • hamming
  • hanning
  • kaiser (le paramètre est fixé à 2.0)
  • none

Sélectionner des plages de fréquences

Avec l’API UliEngineering, sélectionner une plage de fréquences de la FFT est extrêmement facile : utilisez simplement fft[lowfreq:highfreq]. Vous pouvez utiliser fft[lowfreq:] pour sélectionner tout à partir de lowfreq ou utiliser fft[:highfreq] pour sélectionner tout jusqu’à highfreq.

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

# Sélectionner la plage de fréquences : 50 à 200 Hz
fft = fft[50.0:200.0]

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

plt.gcf().set_size_inches(10, 5) # Utiliser (20, 10) pour un graphique plus grand
plt.plot(fft.frequencies, fft.amplitudes)
plt.xlabel("Frequency")
plt.ylabel("Amplitude")
plt.savefig("/ram/fft-frequency-range.svg")

Fenêtre de fréquence

Extraire l’amplitude et l’angle à une certaine fréquence

En utilisant [frequency], c.-à-d. l’opérateur getitem avec une seule valeur, vous obtenez un objet FFTPoint() contenant la fréquence, l’amplitude et l’angle relatif pour une fréquence donnée. La bibliothèque sélectionne automatiquement le bucket FFT le plus proche, donc même si votre FFT n’a pas de bucket pour cette fréquence spécifique, vous obtiendrez des résultats cohérents.

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

# Afficher la valeur à une certaine fréquence
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)

Des FFT courtes pour des données longues

Avec compute_fft(), si nous avons un tableau de données extrêmement long, cela signifie que nous calculerons une FFT extrêmement longue. Dans de nombreux cas, cela n’est pas souhaitable et vous voulez calculer une FFT de taille fixe (le plus souvent une FFT de puissance de deux, p. ex. 1024, 2048, 4096 etc).

Générons quelques données de test longues et supposons que nous voulons calculer une FFT de taille 1024 dessus

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

Tracer les données comme ci-dessus donne

simple_serial_fft_reduce() s’occupe de toute la magie et de la normalisation pour nous, y compris le partitionnement des données en morceaux chevauchants, l’addition des FFT et la normalisation correcte du résultat.

La convention de nommage est significative ici :

• simple_...._reduce signifie qu’il s’agit d’une variante avec des valeurs par défaut sensées (simple) pour utiliser une fonction de reduction (par défaut : sum) sur plusieurs FFT.
• serial signifie les FFT individuelles

Paralléliser les FFT

Si vous avez un jeu de données énorme, vous pouvez utiliser simple_parallel_fft_reduce() de manière identique à simple_serial_fft_reduce() :

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

Cependant, dans la plupart des cas, vous voudrez initialiser l’exécuteur manuellement afin de pouvoir le réutiliser plus tard :

from UliEngineering.SignalProcessing.FFT import *
from concurrent.futures import ThreadPoolExecutor
executor = ThreadPoolExecutor() # Pas d'argument => utiliser num_cpus threads
fft = simple_parallel_fft_reduce(data, samplerate=1e3, fftsize=1024, executor=executor)

Nous pouvons utiliser un ThreadPoolExecutor() car scipy.fftpack (que UliEngineering utilise pour faire les calculs complexes) libère le GIL de Python.

Notez qu’en raison du besoin d’effectuer de nombreuses tâches de gestion, simple_parallel_fft_reduce() est beaucoup plus lent que simple_serial_fft_reduce() si vous avez un jeu de données si petit que la parallélisation n’est pas efficace. Ma recommandation initiale est d’envisager d’utiliser la variante parallèle si le temps d’exécution total de la variante série est supérieur à $0.5s$