Snadný výpočet a vizualizace FFT v Pythonu pomocí UliEngineering

UliEngineering je knihovna pro smíšenou analýzu dat v Pythonu — jedním z nástrojů, které poskytuje, je snadno použitelný balíček pro výpočet FFT. Na rozdíl od jiných balíčků je tato knihovna orientována na praktické případy použití a umožňuje vám provést FFT jen jedním řádkem kódu! Nevyžaduje se znalost matematiky. Nejprve nainstalujte UliEngineering.

Začínáme

Nejprve vygenerujeme nějaká testovací data. Viz tento předchozí příspěvek pro více detailů o tom, jak generovat sinusoidní testovací data:

from UliEngineering.SignalProcessing.Simulation import *

# Generovat testovací data: tón 100 Hz + 400 Hz
data = sine_wave(frequency=100.0, samplerate=1000, amplitude=1.) \
       + sine_wave(frequency=400.0, samplerate=1000, amplitude=0.5)

Testovací data se skládají z 100 Hz sinusovky plus 400 Hz sinusovky (s poloviční amplitudou). Tento signál je vzorkován vzorkovací frekvencí 1000 Hz.

Nyní můžeme vypočítat a vizualizovat FFT pomocí matplotlib:

# Vypočítat FFT. Použít stejnou vzorkovací frekvenci
# POZNÁMKA: Okno má výchozí hodnotu "blackman"!
from UliEngineering.SignalProcessing.FFT import compute_fft
fft = compute_fft(data, samplerate=1e3)

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

plt.gcf().set_size_inches(10, 5) # Pro větší graf použijte (20, 10)
plt.plot(fft.frequencies, fft.amplitudes)
plt.xlabel("Frekvence")
plt.ylabel("Amplituda")

compute_fft(data, samplerate=1e3) vrací FFT objekt obsahující pole jako frequencies, amplitude a phase. Provádí FFT o velikosti vstupu (tj. jelikož data jsou pole o délce 1000, FFT bude velikosti 1000).

fft.frequencies je pole frekvencí (v Hz), odpovídající hodnotám v fft.amplitudes. Můžete také použít fft.angles pro získání relativních úhlů ve stupních, ale to není pokryto v tomto blogpostu.

Jak můžete vidět ve výše zobrazeném grafu, maximální frekvence, kterou lze detekovat, je vždy polovina vzorkovací frekvence, tj. pro naši vzorkovací frekvenci $f_s = 1000,\text{Hz}$ je to $500,\text{Hz}$. Viz toto více matematicky zaměřené vysvětlení FFT, pokud chcete vědět více detailů.

Interně compute_fft() provádí tento výpočet:

• $2$ je korekční faktor, který bere v úvahu, že zahazujeme druhou polovinu surového výsledku FFT (jelikož děláme reálnou FFT)
• $\frac{1}{\text{len(data)}}$ normalizuje výsledky FFT, takže jsou nezávislé na délce dat (tj. pokud předáte delší vzorek stejné sinusovky, stále dostanete stejný výsledek)
• $\text{abs}\left(\cdots\right)$ Převede komplexní výsledek FFT citlivý na fázi na spektrum, které je snazší číst a vizualizovat.
• Okno (které má výchozí hodnotu blackman) je okno, které se aplikuje na data pro zmírnění některých matematických efektů na začátku a na konci datasetu. Viz Wikipedia o okenních funkcích pro více detailů. UliEngineering aktuálně nabízí tento seznam okenních funkcí:
  • blackman
  • bartlett
  • hamming
  • hanning
  • kaiser (Parametr je pevně nastaven na 2.0)
  • none

Výběr frekvenčních rozsahů

Pomocí API UliEngineering je výběr frekvenčního rozsahu FFT triviálně snadný: Stačí použít fft[lowfreq:highfreq]. Můžete použít fft[lowfreq:] pro výběr všeho začínajícího od lowfreq nebo použít fft[:highfreq] pro výběr všeho až do highfreq.

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

# Vybrat frekvenční rozsah: 50 až 200 Hz
fft = fft[50.0:200.0]

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

plt.gcf().set_size_inches(10, 5) # Pro větší graf použijte (20, 10)
plt.plot(fft.frequencies, fft.amplitudes)
plt.xlabel("Frekvence")
plt.ylabel("Amplituda")
plt.savefig("/ram/fft-frequency-range.svg")

Frekvenční okno

Extrakce amplitudy a úhlu při určité frekvenci

Pomocí [frequency], tj. operátoru getitem s jednou hodnotou, dostanete objekt FFTPoint() obsahující frekvenci, amplitudu a relativní úhel pro danou frekvenci. Knihovna automaticky vybírá nejbližší FFT bucket, takže i když váš FFT nemá bucket pro tuto specifickou frekvenci, dostanete rozumné výsledky.

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

# Zobrazit hodnotu při určité frekvenci
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)

Krátké FFT pro dlouhá data

Pomocí compute_fft(), pokud máme extrémně dlouhé pole dat, znamená to, že počítáme extrémně dlouhou FFT. V mnoha případech to není žádoucí a chcete počítat FFT s pevnou velikostí (nejčastěji FFT o velikosti mocniny dvojky, např. 1024, 2048, 4096 atd).

Vygenerujme nějaká dlouhá testovací data a předpokládejme, že chceme na nich spočítat FFT o velikosti 1024

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

Vykreslení dat jako výše dává

simple_serial_fft_reduce() se stará o veškerou magii a normalizaci za nás, včetně rozdělení dat do překrývajících se bloků, sčítání FFT a správné normalizace výsledku.

Konvence pojmenování je zde významná:

• simple_...._reduce znamená, že toto je varianta s rozumnými výchozími hodnotami (simple) pro použití reduction funkce (výchozí: sum) na více FFT.
• serial znamená jednotlivé FFT

Paralelizace FFT

Pokud máte obrovský dataset, můžete použít simple_parallel_fft_reduce() identicky k simple_serial_fft_reduce():

from UliEngineering.SignalProcessing.FFT import *
fft = simple_parallel_fft_reduce(data, samplerate=1e3, fftsize=1024)
# Použít fft.frequencies, fft.amplitudes atd.

Nicméně ve většině případů chcete inicializovat executor ručně, abyste ho mohli později znovu použít:

from UliEngineering.SignalProcessing.FFT import *
from concurrent.futures import ThreadPoolExecutor
executor = ThreadPoolExecutor() # Žádný argument => použít num_cpus vláken
fft = simple_parallel_fft_reduce(data, samplerate=1e3, fftsize=1024, executor=executor)

Můžeme použít ThreadPoolExecutor() jelikož scipy.fftpack (který UliEngineering používá pro náročné matematické výpočty) uvolňuje Python GIL.

Pamatujte, že kvůli nutnosti dělat hodně úklidových úkonů je simple_parallel_fft_reduce() mnohem pomalejší než simple_serial_fft_reduce(), pokud máte dataset tak malý, že paralelizace není efektivní. Moje počáteční doporučení je zvážit použití paralelní varianty, pokud celková doba běhu sériové varianty je větší než $0.5s$