6.14. FFT

Transformata Fouriera rozkłada sygnał w dziedzinie czasu na jego składowe częstotliwości. Szybka transformata Fouriera (FFT) to standardowy algorytm jej wydajnego obliczania – powód istnienia tego modułu oraz powód, dla którego transformaty wystarczająco długie, by były użyteczne (1024 punkty, 4096 punktów), są w ogóle wykonalne na mikrokontrolerze.

Na kamerze wejściem jest zazwyczaj bufor równomiernie rozmieszczonych próbek z mikrofonu, osi akcelerometru, czujnika prądu lub sondy drgań; wyjściem jest widmo, które aplikacja następnie bada pod kątem szczytów, energii lub cech dla klasyfikatora.

numpy.fft dostarcza jednowymiarową dyskretną transformatę Fouriera. Dwie funkcje to fft() (prosta) oraz ifft() (odwrotna, znormalizowana przez N).

Długość transformaty musi być potęgą dwójki. Inne długości zgłaszają ValueError. N = 256 i N = 1024 to częste wybory: wystarczająco długie, by zapewnić użyteczną rozdzielczość częstotliwościową przy większości wbudowanych częstotliwości próbkowania, a zarazem na tyle krótkie, by zmieścić się w pamięci RAM i zakończyć obliczenia z dużym zapasem w obrębie jednego okresu ramki.

fft() przyjmuje część rzeczywistą wejścia jako pierwszy argument pozycyjny oraz opcjonalną część urojoną jako drugi, a zwraca 2-krotkę (real, imag) tablic rzeczywistych.

6.14.1. Wejście rzeczywiste, wyjście jako moduł

Gdy aplikacja podaje wyłącznie wejście rzeczywiste i potrzebuje jedynie widma modułu:

real, imag = np.fft.fft(x)
magnitude  = np.sqrt(real * real + imag * imag)

Wynikiem jest rzeczywista tablica ndarray o tej samej długości co wejście. Funkcja pomocnicza spectrogram() robi dokładnie to samo bez pośrednich alokacji – właściwe rozwiązanie wewnątrz pętli strumieniowej.

6.14.2. Transformata odwrotna

ifft() jest znormalizowana tak, że obieg w obie strony zwraca pierwotne wejście z dokładnością do błędu zmiennoprzecinkowego. Oba wywołania przyjmują i zwracają parę tablic rzeczywistych:

real, imag       = np.fft.fft(y)
re_back, im_back = np.fft.ifft(real, imag)

6.14.3. Wyznaczanie dominującej częstotliwości

Częstym wzorcem jest zlokalizowanie największego szczytu w buforze próbek – dźwięku, drgań, modulowanej podczerwieni. Przepis:

  1. Pobierz N próbek, gdzie N jest potęgą dwójki.

  2. Wykonaj FFT.

  3. Znajdź prążek o największym module.

  4. Przelicz indeks prążka na Hz przy użyciu częstotliwości próbkowania.

from ulab import numpy as np

N  = 1024
fs = 8000.0                       # sample rate, Hz

# ... fill ``samples`` with N data points ...

real, imag = np.fft.fft(samples)
spectrum   = np.sqrt(real * real + imag * imag)

half      = spectrum[:N // 2]     # only first half is meaningful
peak_bin  = np.argmax(half)
peak_hz   = peak_bin * fs / N

print("peak at", peak_hz, "Hz")