Szybka transformata Fouriera

FFT (Fast Fourier Transformation) to jedno z najważniejszych narzędzi w przekształcaniu i analizowaniu sygnałów.

W wersji ciągłej zdefiniowana jest w sposób następujący: $$X(f) = \int_{-\infty}^\infty x(t) e^{-i 2 \pi f t} dt$$ Natomiast wersja dyskretna określona jest wzorem: $$X_k = \sum_{n=0}^{N-1} x_n \cdot e^{-i~2\pi~n~\frac{k}{N}}$$ gdzie:

  • $N$ to liczba próbek,
  • $x_i$ to $i$-ta wartość funkcji w dmenie czasu,
  • $X_k$ to (zespolona, na ogół) wartość transformaty.

Jeżeli $X_k$ i $x_k$ potraktować jako składowe wektora, transformatę można zapisać w sposób wektorowy: $$\vec{X} = M \cdot \vec{x}$$ gdzie elementy macierzy $M$ dane są wzorem: $$M_{kn} = e^{-i~2\pi~k~\frac{n}{N}}.$$ Jest to zapis, tak zwanej Wolnej Transformaty Fouriera. Spróbujmy ją policzyć, czy może raczej zaproponować funkcję, która będzie ją liczyła.

In [1]:
import numpy as np
def DFT_slow(x):
    """wylicza Transformatę Fouriera dla jednowymiarowej macierzy x"""
    x = np.asarray(x, dtype=float)
    N = x.shape[0]
    n = np.arange(N)
    k = n.reshape((N, 1))
    M = np.exp(-2j * np.pi * k * n / N)
    return np.dot(M, x)
  1. Funkcja np.asarray() pozwala skonwertować dane wprowadzone w jakimś formacie do postaci tablicy typu rzeczywistego (float).
  2. x.shape[0] zwraca liczbę wierszy obiektu x.
  3. Funkcja np.arange(N) tworzy tablicę (wiersz) o kolejnych wartościach od zera do N.
  4. Metoda obiekt.reshape((N,1)) nadaje obiekt owi inną postać; w tym wypadku konwertuje wiersz w kolumnę o dziesięciu wierszach.
  5. Kolejne polecenie (M = np.exp(…) wylicza „na raz” wszystkie elementy tablicy.

Funkcja zwraca iloczyn wektorowy M razy x.

Wygenerujemy zestaw losowych wartości o 1024 elementach:

In [2]:
x = np.random.random(1024)
In [3]:
DFT_slow(x)
Out[3]:
array([ 514.41243660 +0.j        ,   -4.79884555+11.70083685j,
         -6.96962199-11.06828923j, ...,   12.25441607 -3.87391701j,
         -6.96962199+11.06828923j,   -4.79884555-11.70083685j])

Funkcja np.fft.fft() wylicza transformatę Fouriera metodą szybką.

In [4]:
np.fft.fft(x)
Out[4]:
array([ 514.41243660 +0.j        ,   -4.79884555+11.70083685j,
         -6.96962199-11.06828923j, ...,   12.25441607 -3.87391701j,
         -6.96962199+11.06828923j,   -4.79884555-11.70083685j])

Na pierwszy rzut oka dane są podobne, ale żeby sprawdzić czy wyniki są zbliżone, mozemy użyć funkcji numpy.allclose()

In [5]:
np.allclose(DFT_slow(x), np.fft.fft(x))
Out[5]:
True

Czyli nie jest źle.

Żeby sprawdzić które obliczenia są szybsze — użyjemy polecenia timeit

In [6]:
%timeit DFT_slow(x)
%timeit np.fft.fft(x)
191 ms ± 2.37 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
56.9 µs ± 239 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

Czyli widać różnicę!

Szybka transformata Fouriera jest najszybsza, gdy N jest potęgą liczby 2. W innym przypadku…

In [7]:
x = np.random.random(1021)
In [8]:
%timeit DFT_slow(x)
%timeit np.fft.fft(x)
170 ms ± 3.72 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
2.58 ms ± 315 ns per loop (mean ± std. dev. of 7 runs, 100 loops each)
In [9]:
x = np.random.random(1023)
In [10]:
%timeit DFT_slow(x)
%timeit np.fft.fft(x)
182 ms ± 2.33 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
118 µs ± 113 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

Liczba 1021 jest pierwsza. I o ile wyniki „wolnej” transformaty (zaprogramowanej w naiwny sposób) nie zmieniają się specjalnie, to różnice w przypadku transformaty szybkiej (zakładam, że zaprogramowanej w sposób „optymalny”) już się różnią). 1023 to $3 \times 11 \times 31$ i nawet taki rozkład na „czynniki pierwsze” pozwala przyśpieszyć obliczenia.

In [11]:
x = np.random.random(1026)
In [12]:
%timeit DFT_slow(x)
%timeit np.fft.fft(x)
186 ms ± 1.01 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
90.2 µs ± 119 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

1026 to $2 \times 3 \times 3 \times 3 \times 19$; im więcej czynników — tym obliczenia szybsze.

Proste rzeczy najpierw

Wygenerujemy prostą funkcję okresową i obejrzymy jej transformatę Fouriera.

In [13]:
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import scipy

freq = 1 #Hz - liczba okresów na sekundę
amplitude = 3
time_to_plot = 2 # czas (długość danych wejściowych)
sample_rate = 100 # liczba próbek na sekundę
num_samples = sample_rate * time_to_plot

t = np.linspace(0, time_to_plot, num_samples)
signal = [amplitude * np.sin(freq * i * 2*np.pi) for i in t] 

Nasz sygnał wygląda tak:

In [14]:
plt.plot(t, signal);
In [15]:
fft_output = np.fft.rfft(signal)
magnitude_only = [np.sqrt(i.real**2 + i.imag**2)/len(fft_output) for i in fft_output]
frequencies = [(i*1.0/num_samples)*sample_rate for i in range(num_samples//2+1)]
In [16]:
plt.plot(frequencies, magnitude_only, 'r')
Out[16]:
[<matplotlib.lines.Line2D at 0x7f015987a9b0>]
In [17]:
plt.plot(magnitude_only)
#plt.show()
Out[17]:
[<matplotlib.lines.Line2D at 0x7f01597e5780>]

Prosty przykład

Najpierw przygotujemy wartości „zmiennej niezależnej”.

In [18]:
sample_rate = 10 # próbek na sekundę
total_sampling_time = 3
num_samples = sample_rate * total_sampling_time

t = np.linspace(0, total_sampling_time, num_samples)

Zdefiniujemy funkcję harmoniczną o częstości 1 Hz i amplitudzie 5. $$F(t)=A\sin(2\pi f t)$$

f = lambda x: wave_amplitude * np.sin(frequency_in_hz * 2*np.pi*x)

to konstrukcja pozwalająca stworzyć funkcję „anonimową”, nie związaną z nazwą. Python pozwala zdefiniować funkcję w sposób następujący:

def f (x): return x**2

Tworzac funkcję lambda, jej definicja będzie następująca:

g = lambda x: x**2

Podstawowa różnica polega na tym, że nie ma polecenia return. Zatem w najprostszych zastosowaniach różnicy nie ma. Ale można wymyśleć zastosowania zaawansowane, w któych konstrukcja taka będzie przydatna.

In [19]:
frequency_in_hz = 1
wave_amplitude = 5
f = lambda x: wave_amplitude * np.sin(frequency_in_hz * 2*np.pi*x)
    
sampled_f = [f(i) for i in t]
In [20]:
plt.plot(t, sampled_f)
#plt.show()
Out[20]:
[<matplotlib.lines.Line2D at 0x7f01597ca5c0>]

Wyliczmy FFT dla tego zestawu punktów. Skorzystamy z funkcji fft.fft().

In [21]:
fft_output = np.fft.fft(sampled_f)
In [22]:
plt.plot(fft_output);
/usr/local/lib/python3.5/dist-packages/numpy/core/numeric.py:531: ComplexWarning: Casting complex values to real discards the imaginary part
  return array(a, dtype, copy=False, order=order)

Jak wiadomo, dla rzeczywistego przebiegu wejściowegp transformata Fouriera będzie funkcją zespoloną i symetryczną. Jeżeli nie będziemy transformować wartości zespolonych można ograniczyć się do stosowania funkcji fft.rfft().

In [23]:
rfft_output = np.fft.rfft(sampled_f)
plt.plot(rfft_output)
/usr/local/lib/python3.5/dist-packages/numpy/core/numeric.py:531: ComplexWarning: Casting complex values to real discards the imaginary part
  return array(a, dtype, copy=False, order=order)
Out[23]:
[<matplotlib.lines.Line2D at 0x7f01596bfd68>]

Skale na osiach

Jak widać nie bardzo wiadomo jaki jest związek pomiędzy opisem osi $x$ a częstością przebiegu i jak się ma opis na osi $y$ do amplitudy sygnału. Można poeksperymentować zmieniając odpowiednie wartości w notatniku i przeliczając go ponownie…

Oś $x$

Jeżeli przyjrzymy się zaprezentowanemu na początku równaniu: $$X_k = \sum_{n=0}^{N-1} x_n \cdot e^{-i~2\pi~n~\frac{k}{N}}$$ $X$ jest, tak na prawdę, funkcją $k$ i $N$, a dokładniej $\frac{k}{N}$.

Dodatkowo zawsze interesuje nas tylko połowa wyników (ze względu na symetrię). Numery współczynników transformaty Fouriera można przeliczyć na częstości:

In [24]:
rfreqs = [(i*1.0/num_samples)*sample_rate for i in range(num_samples//2+1)]

Zamiast powyższego wyliczeni amożna również użyć funkcji fft.rfftfreq().

In [25]:
np.fft.rfftfreq(num_samples, d=1./sample_rate)
Out[25]:
array([ 0.        ,  0.33333333,  0.66666667,  1.        ,  1.33333333,
        1.66666667,  2.        ,  2.33333333,  2.66666667,  3.        ,
        3.33333333,  3.66666667,  4.        ,  4.33333333,  4.66666667,  5.        ])
In [26]:
rfreqs
Out[26]:
[0.0,
 0.3333333333333333,
 0.6666666666666666,
 1.0,
 1.3333333333333333,
 1.6666666666666665,
 2.0,
 2.3333333333333335,
 2.6666666666666665,
 3.0,
 3.333333333333333,
 3.6666666666666665,
 4.0,
 4.333333333333334,
 4.666666666666667,
 5.0]

Na pierwszy rzut oka wyniki są identyczne?

In [27]:
plt.plot(rfreqs,rfft_output)
/usr/local/lib/python3.5/dist-packages/numpy/core/numeric.py:531: ComplexWarning: Casting complex values to real discards the imaginary part
  return array(a, dtype, copy=False, order=order)
Out[27]:
[<matplotlib.lines.Line2D at 0x7f015962c7b8>]

Jak widać jest już lepiej.

Oś $y$

Na osi $y$ pojawiają się również wartości ujemne! Chcielibyśmy żeby wszystko było wyskalowane w tych samych jednostkach co ampliktida sygnału.

Moduł każdego współczynnika Fouriera można wyliczyć ze wzoru: $$magnitude(i) = \frac{\sqrt{i.real^2 + i.imag^2}}{N}$$

In [28]:
rfft_mag = [np.sqrt(i.real**2 + i.imag**2)/len(rfft_output) for i in rfft_output]
In [29]:
plt.plot(rfreqs, rfft_mag)
Out[29]:
[<matplotlib.lines.Line2D at 0x7f01596101d0>]

Teraz wygląda to znacznie lepiej.

Osobna kwestia do dyskusji jest taka: czemu amplituda nie jest dokładnie równa 5?

In [30]:
max(rfft_mag)
Out[30]:
4.5218718374450848

Odwrotna transformata Fouriera

Do jej wyliczania służy funkcja fft.irfft() (lub fft.ifft).

In [31]:
irfft_output = np.fft.irfft(rfft_output)
In [32]:
plt.plot(t, irfft_output)
Out[32]:
[<matplotlib.lines.Line2D at 0x7f0159565e80>]

Sprawdźmy jaka jest różnica między przebiegiem oryginalnym, a odtworzonym:

In [33]:
irfft_output-sampled_f
Out[33]:
array([ -2.36847579e-16,   0.00000000e+00,   2.66453526e-15,
         3.55271368e-15,   0.00000000e+00,  -3.33066907e-15,
         4.44089210e-16,   3.55271368e-15,   8.88178420e-16,
         0.00000000e+00,   4.88498131e-15,   4.88498131e-15,
         8.88178420e-16,  -3.55271368e-15,  -1.77635684e-15,
        -3.99680289e-15,  -6.21724894e-15,  -4.44089210e-15,
        -4.88498131e-15,  -4.21884749e-15,  -5.77315973e-15,
         8.88178420e-16,   3.55271368e-15,   6.21724894e-15,
         2.99760217e-15,   3.99680289e-15,   2.66453526e-15,
        -1.77635684e-15,  -1.33226763e-15,  -1.77355391e-15])
In [ ]: