W (sporym uproszczeniu) polega na znalezieniu funkcji $f(x)$ takiej, że jeżeli mamy zbiór par punktów $$\{x_i, y_i\};\quad i=1,2,\ldots N$$ $$f(x_i) = y_i\quad \forall\; 1\le i\le N$$
W ogólnym przypadku zadanie jest trudne (sprowadza się do rozwiązania układu $N$ równań nieliniowych), ale jeżeli przyjąć, że funkcja $f$ jest specjalnej postaci — można rozwiązać je łatwiej.
Szczególnym (ale niezbyt dobrym) przypadkiem funkcji $f$ może być wielomian. Jak wiadomo, dla wyrysowania prostej ($y=ax+b$) znaleźć trzeba wartość dwu parametrów, czyli jak wiadomo z geometrii, potrzeba dwu punktów, przez które prosta ma przejść.
W ogólnym przypadku wielomianu stpnia $m$ potrzeba będzie $m+1$ punktów, żeby przeprowadzić przez nie wielomian.
Jeszcze innym szczególny przypadkiem może być wielomian trygonometryczny czyli funkcja o postaci: $$ y = f(x) = \sum_{k=0}^{m} a_{k} \sin(k x) + b_{k} \cos(k x) $$ Układ równań $$ y_i = f(x_i)\quad i=1, 2,\ldots, N $$ w tym przypadku jest bardzo specyficzny. Można go rozwiązywać bezpośrednio, ale zauważono, że w przypadku gdy $m$ jest postaci $2^k$ można go rozwiazać znacznie szybciej (korzystając z różnych regularności funkcji trygonometrycznych).
Tak powstałą 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:
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.
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)
np.asarray()
pozwala skonwertować dane wprowadzone w jakimś formacie do postaci tablicy typu rzeczywistego (float
).x.shape[0]
zwraca liczbę wierszy obiektu x
.np.arange(N)
tworzy tablicę (wiersz) o kolejnych wartościach od zera do N..reshape((N,1))
nadaje obiekt owi inną postać; w tym wypadku konwertuje wiersz w kolumnę o dziesięciu wierszach.M = np.exp(
…) wylicza „na raz” wszystkie elementy tablicy.Funkcja zwraca iloczyn wektorowy M
razy x
.
Wygenerujemy zestaw losowych wartości o 1024 elementach:
x = np.random.random(1024)
DFT_slow(x)
array([ 5.13756711e+02+0.j , 1.70134671e+00+9.70960396j, -4.92046996e+00+7.36887985j, ..., -4.31419041e-01+0.22796975j, -4.92046996e+00-7.36887985j, 1.70134671e+00-9.70960396j])
Funkcja np.fft.fft()
wylicza transformatę Fouriera metodą szybką.
np.fft.fft(x)
array([ 5.13756711e+02+0.j , 1.70134671e+00+9.70960396j, -4.92046996e+00+7.36887985j, ..., -4.31419041e-01+0.22796975j, -4.92046996e+00-7.36887985j, 1.70134671e+00-9.70960396j])
Na pierwszy rzut oka dane są podobne, ale żeby sprawdzić czy wyniki są zbliżone, mozemy użyć funkcji numpy.allclose()
np.allclose(DFT_slow(x), np.fft.fft(x))
True
Czyli nie jest źle.
Żeby sprawdzić które obliczenia są szybsze — użyjemy polecenia timeit
%timeit DFT_slow(x)
%timeit np.fft.fft(x)
82.8 ms ± 4.01 ms per loop (mean ± std. dev. of 7 runs, 10 loops each) 68.7 µs ± 1.21 µs per loop (mean ± std. dev. of 7 runs, 10,000 loops each)
Czyli widać różnicę!
Szybka transformata Fouriera jest najszybsza, gdy N jest potęgą liczby 2. W innym przypadku…
x = np.random.random(1021)
%timeit DFT_slow(x)
%timeit np.fft.fft(x)
68.7 ms ± 3.12 ms per loop (mean ± std. dev. of 7 runs, 10 loops each) 284 µs ± 680 ns per loop (mean ± std. dev. of 7 runs, 1,000 loops each)
x = np.random.random(1023)
%timeit DFT_slow(x)
%timeit np.fft.fft(x)
70.5 ms ± 4.36 ms per loop (mean ± std. dev. of 7 runs, 10 loops each) 132 µs ± 87 ns per loop (mean ± std. dev. of 7 runs, 10,000 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.
x = np.random.random(1026)
%timeit DFT_slow(x)
%timeit np.fft.fft(x)
78.9 ms ± 4.88 ms per loop (mean ± std. dev. of 7 runs, 10 loops each) 96.5 µs ± 193 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)
1026 to $2 \times 3 \times 3 \times 3 \times 19$; im więcej czynników — tym obliczenia szybsze.
Wygenerujemy prostą funkcję okresową i obejrzymy jej transformatę Fouriera.
%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:
plt.plot(t, signal);
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)]
plt.plot(frequencies, magnitude_only, 'r')
[<matplotlib.lines.Line2D at 0x7f6389636100>]
plt.plot(magnitude_only)
#plt.show()
[<matplotlib.lines.Line2D at 0x7f63895a70a0>]
Najpierw przygotujemy wartości „zmiennej niezależnej”.
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.
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]
plt.plot(t, sampled_f)
#plt.show()
[<matplotlib.lines.Line2D at 0x7f6389501f40>]
Wyliczmy FFT dla tego zestawu punktów. Skorzystamy z funkcji fft.fft()
.
fft_output = np.fft.fft(sampled_f)
plt.plot(fft_output);
/opt/anaconda3/lib/python3.9/site-packages/matplotlib/cbook/__init__.py:1298: ComplexWarning: Casting complex values to real discards the imaginary part return np.asarray(x, float)
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()
.
rfft_output = np.fft.rfft(sampled_f)
plt.plot(rfft_output)
/opt/anaconda3/lib/python3.9/site-packages/matplotlib/cbook/__init__.py:1298: ComplexWarning: Casting complex values to real discards the imaginary part return np.asarray(x, float)
[<matplotlib.lines.Line2D at 0x7f638947ddf0>]
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…
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:
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()
.
np.fft.rfftfreq(num_samples, d=1./sample_rate)
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. ])
rfreqs
[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?
plt.plot(rfreqs,rfft_output)
/opt/anaconda3/lib/python3.9/site-packages/matplotlib/cbook/__init__.py:1298: ComplexWarning: Casting complex values to real discards the imaginary part return np.asarray(x, float)
[<matplotlib.lines.Line2D at 0x7f63893f5280>]
Jak widać jest już lepiej.
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}$$
rfft_mag = [np.sqrt(i.real**2 + i.imag**2)/len(rfft_output) for i in rfft_output]
plt.plot(rfreqs, rfft_mag)
[<matplotlib.lines.Line2D at 0x7f638934f340>]
Teraz wygląda to znacznie lepiej.
Osobna kwestia do dyskusji jest taka: czemu amplituda nie jest dokładnie równa 5?
max(rfft_mag)
4.521871837445084
Do jej wyliczania służy funkcja fft.irfft()
(lub fft.ifft
).
irfft_output = np.fft.irfft(rfft_output)
plt.plot(t, irfft_output)
[<matplotlib.lines.Line2D at 0x7f638932b880>]
Sprawdźmy jaka jest różnica między przebiegiem oryginalnym, a odtworzonym:
irfft_output-sampled_f
array([ 2.36847579e-16, -1.77635684e-15, 8.88178420e-16, -8.88178420e-16, 4.44089210e-16, 7.77156117e-16, -4.44089210e-16, 1.77635684e-15, -1.77635684e-15, 0.00000000e+00, -1.11022302e-15, -1.77635684e-15, 8.88178420e-16, -1.77635684e-15, -6.66133815e-16, 6.66133815e-16, 8.88178420e-16, 2.66453526e-15, -1.33226763e-15, 8.88178420e-16, 1.77635684e-15, -1.77635684e-15, 0.00000000e+00, -1.77635684e-15, 3.33066907e-16, 0.00000000e+00, -3.55271368e-15, 1.77635684e-15, 4.44089210e-16, 8.31769454e-16])