FFTs in Python einfach berechnen & visualisieren mit UliEngineering

UliEngineering ist eine gemischte Datenanalyse-Bibliothek in Python – eines der Tools, die sie bietet, ist ein einfach zu verwendendes Paket zur Berechnung von FFTs. Im Gegensatz zu anderen Paketen ist diese Bibliothek auf praktische Anwendungsfälle ausgerichtet und ermöglicht dir, die FFT in nur einer Codezeile durchzuführen! Keine Mathekenntnisse erforderlich. Installiere zunächst UliEngineering.

Erste Schritte

Zuerst generieren wir einige Testdaten. Siehe diesen vorherigen Post für weitere Details zur Generierung von Sinus-Testdaten:

from UliEngineering.SignalProcessing.Simulation import *

# Testdaten generieren: 100 Hz + 400 Hz Ton
data = sine_wave(frequency=100.0, samplerate=1000, amplitude=1.) \
       + sine_wave(frequency=400.0, samplerate=1000, amplitude=0.5)

Die Testdaten bestehen aus einer 100 Hz-Sinuswelle plus einer 400 Hz-Sinuswelle (mit halber Amplitude). Dieses Signal wird mit einer Abtastrate von 1000 Hz abgetastet.

Nun können wir die FFT mit matplotlib berechnen & visualisieren:

# FFT berechnen. Gleiche Abtastrate verwenden
# HINWEIS: Fenster ist standardmäßig "blackman"!
from UliEngineering.SignalProcessing.FFT import compute_fft
fft = compute_fft(data, samplerate=1e3)

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

plt.gcf().set_size_inches(10, 5) # Verwende (20, 10) für einen größeren Plot
plt.plot(fft.frequencies, fft.amplitudes)
plt.xlabel("Frequency")
plt.ylabel("Amplitude")

compute_fft(data, samplerate=1e3) gibt ein FFT-Objekt zurück, das Felder wie frequencies, amplitude und phase enthält. Es führt eine FFT der Größe des Inputs durch (d.h. da data ein Array der Länge 1000 ist, ist die FFT der Größe 1000).

fft.frequencies ist ein Array von Frequenzen (in Hz), entsprechend den Werten in fft.amplitudes. Du kannst auch fft.angles verwenden, um relative Winkel in Grad zu erhalten, was in diesem Blogpost jedoch nicht behandelt wird.

Wie du im oben gezeigten Plot sehen kannst, ist die maximale Frequenz, die erkannt werden kann, immer die Hälfte der Abtastrate, d.h. für unsere Abtastrate von $f_s = 1000,\text{Hz}$ ist es $500,\text{Hz}$. Siehe diese math-zentrierte FFT-Erklärung wenn du mehr Details wissen möchtest.

Intern führt compute_fft() diese Berechnung durch:

• $2$ ist ein Korrekturfaktor, der berücksichtigt, dass wir die hintere Hälfte des rohen FFT-Ergebnisses verwerfen (da wir eine reelle FFT durchführen)
• $\frac{1}{\text{len(data)}}$ normalisiert die FFT-Ergebnisse, sodass sie unabhängig von der Datenlänge sind (d.h. wenn du ein längeres Sample derselben Sinuswelle übergibst, erhältst du dennoch dasselbe Ergebnis)
• $\text{abs}\left(\cdots\right)$ wandelt das komplexe phasenauswertende Ergebnis der FFT in ein Spektrum um, das leichter zu lesen & visualisieren ist.
• Window (standardmäßig blackman) ist das Fenster, das auf die Daten angewendet wird, um einige mathematische Effekte am Anfang und am Ende des Datensatzes zu mildern. Siehe Wikipedia zu Fensterfunktionen für weitere Details. UliEngineering bietet derzeit diese Liste von Fensterfunktionen:
  • blackman
  • bartlett
  • hamming
  • hanning
  • kaiser (Parameter ist fest auf 2.0 gesetzt)
  • none

Frequenzbereiche auswählen

Mit der UliEngineering-API ist die Auswahl eines Frequenzbereichs der FFT trivial einfach: Verwende einfach fft[lowfreq:highfreq]. Du kannst fft[lowfreq:] verwenden, um alles ab lowfreq auszuwählen, oder fft[:highfreq] verwenden, um alles bis highfreq auszuwählen.

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

# Frequenzbereich auswählen: 50 bis 200 Hz
fft = fft[50.0:200.0]

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

plt.gcf().set_size_inches(10, 5) # Verwende (20, 10) für einen größeren Plot
plt.plot(fft.frequencies, fft.amplitudes)
plt.xlabel("Frequency")
plt.ylabel("Amplitude")
plt.savefig("/ram/fft-frequency-range.svg")

Frequenzfenster

Amplitude & Winkel bei einer bestimmten Frequenz extrahieren

Durch die Verwendung von [frequency], d.h. dem Getitem-Operator mit einem einzelnen Wert, erhältst du ein FFTPoint()-Objekt, das die Frequenz, Amplitude und den relativen Winkel für eine gegebene Frequenz enthält. Die Bibliothek wählt automatisch den nächsten FFT-Bucket aus, sodass du auch dann sinnvolle Ergebnisse erhältst, wenn deine FFT keinen Bucket für genau diese Frequenz hat.

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

# Wert bei einer bestimmten Frequenz anzeigen
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)

Kurze FFTs für lange Daten

Mit compute_fft() bedeutet ein extrem langes Datenarray, dass wir eine extrem lange FFT berechnen. In vielen Fällen ist dies nicht wünschenswert und man möchte eine FFT fester Größe berechnen (meistens eine FFT mit Zweierpotenz, z.B. 1024, 2048, 4096 usw.).

Lass uns einige lange Testdaten generieren und annehmen, dass wir eine FFT der Größe 1024 darauf berechnen möchten:

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

Das Plotten der Daten wie oben ergibt

was fast genau so aussieht wie unser compute_fft()-Plot zuvor – genau wie erwartet.

simple_serial_fft_reduce() übernimmt die gesamte Magie und Normalisierung für uns, einschließlich der Partitionierung der Daten in überlappende Chunks, dem Aufaddieren der FFTs und der korrekten Normalisierung des Ergebnisses.

Die Namenskonvention ist hier von Bedeutung:

• simple_...._reduce bedeutet, dass dies eine Variante mit sinnvollen Defaults (simple) zur Verwendung einer reduction-Funktion (Standard: sum) auf mehreren FFTs ist.
• serial bedeutet, dass die einzelnen FFTs

FFTs parallelisieren

Wenn du einen riesigen Datensatz hast, kannst du simple_parallel_fft_reduce() identisch zu simple_serial_fft_reduce() verwenden:

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

In den meisten Fällen möchte man den Executor jedoch manuell initialisieren, um ihn später wiederverwenden zu können:

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

Wir können einen ThreadPoolExecutor() verwenden, da scipy.fftpack (das UliEngineering für die schwere Mathematik verwendet) das Python-GIL freigibt.

Beachte, dass aufgrund der Notwendigkeit vieler Verwaltungsaufgaben simple_parallel_fft_reduce() viel langsamer ist als simple_serial_fft_reduce(), wenn du einen so kleinen Datensatz hast, dass Parallelisierung nicht effektiv ist. Meine anfängliche Empfehlung ist, die parallele Variante in Betracht zu ziehen, wenn die Gesamtausführungszeit der seriellen Variante größer als $0.5s$ ist.