#!/usr/bin/env python # coding: utf-8 # # Zadanie interpolacji… # # 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. # # 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()`](https://docs.scipy.org/doc/numpy-1.13.0/reference/generated/numpy.asarray.html) 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)`](https://docs.scipy.org/doc/numpy-1.13.0/reference/generated/numpy.arange.html) 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) # Funkcja [`np.fft.fft()`](https://docs.scipy.org/doc/numpy-1.13.0/reference/generated/numpy.fft.fft.html) wylicza transformatę Fouriera metodą szybką. # In[4]: np.fft.fft(x) # Na pierwszy rzut oka dane są podobne, ale żeby sprawdzić czy wyniki są zbliżone, mozemy użyć funkcji [`numpy.allclose()`](https://docs.scipy.org/doc/numpy-1.13.0/reference/generated/numpy.allclose.html) # In[5]: np.allclose(DFT_slow(x), np.fft.fft(x)) # Czyli nie jest źle. # # Żeby sprawdzić które obliczenia są szybsze — użyjemy polecenia `timeit` # In[6]: get_ipython().run_line_magic('timeit', 'DFT_slow(x)') get_ipython().run_line_magic('timeit', 'np.fft.fft(x)') # 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]: get_ipython().run_line_magic('timeit', 'DFT_slow(x)') get_ipython().run_line_magic('timeit', 'np.fft.fft(x)') # In[9]: x = np.random.random(1023) # In[10]: get_ipython().run_line_magic('timeit', 'DFT_slow(x)') get_ipython().run_line_magic('timeit', 'np.fft.fft(x)') # 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]: get_ipython().run_line_magic('timeit', 'DFT_slow(x)') get_ipython().run_line_magic('timeit', 'np.fft.fft(x)') # 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]: get_ipython().run_line_magic('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') # In[17]: plt.plot(magnitude_only) #plt.show() # ## 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() # Wyliczmy FFT dla tego zestawu punktów. Skorzystamy z funkcji [`fft.fft()`](https://docs.scipy.org/doc/numpy-1.13.0/reference/generated/numpy.fft.fft.html). # In[21]: fft_output = np.fft.fft(sampled_f) # In[22]: plt.plot(fft_output); # 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()`](https://docs.scipy.org/doc/numpy-1.13.0/reference/generated/numpy.fft.rfft.html). # In[23]: rfft_output = np.fft.rfft(sampled_f) plt.plot(rfft_output) # ## 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()`](https://docs.scipy.org/doc/numpy-1.13.0/reference/generated/numpy.fft.rfftfreq.html#numpy.fft.rfftfreq). # In[25]: np.fft.rfftfreq(num_samples, d=1./sample_rate) # In[26]: rfreqs # Na pierwszy rzut oka wyniki są identyczne? # In[27]: plt.plot(rfreqs,rfft_output) # 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) # 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) # ## Odwrotna transformata Fouriera # # Do jej wyliczania służy funkcja [`fft.irfft()`](https://docs.scipy.org/doc/numpy-1.13.0/reference/generated/numpy.fft.irfft.html#numpy.fft.irfft) (lub [`fft.ifft`](https://docs.scipy.org/doc/numpy-1.13.0/reference/generated/numpy.fft.ifft.html#numpy.fft.ifft)). # In[31]: irfft_output = np.fft.irfft(rfft_output) # In[32]: plt.plot(t, irfft_output) # Sprawdźmy jaka jest różnica między przebiegiem oryginalnym, a odtworzonym: # In[33]: irfft_output-sampled_f # In[ ]: