W przypadku gdy, funkcja (zjawisko), którą się zajmujemy jest okresowa czyli $$g(y+t)=g(y)$$ dla wszystkich $y$, do interpolacji najlepiej użyć funkcji okresowych. Dokonując zamiany zmiennych ($x=\frac{2\pi}{t}$) można rozpatrywać funkcje okresowe o okresie $2\pi$.
Zadanie interpolacji trygonometrycznej polega na znalezieniu dla danej funkcji $f$: $$t_n(x) = \sum_{j=0}^n c_j e^{ijx}$$ ($i=\sqrt{-1}$). Oczekujemy, że wielomian ten przyjmie w $n+1$ punktach $x_k$ z przedziału $[0,2\pi]$ te same wartości co interpolowana funkcja, tzn.: $$t_n(x_k)=f(x_k)$$ gdzie $x_k \ne x_l$, gdy $k \ne l$.
Zadanie te rozwiązuje się podobnie jak poprzednie rozwiązując układ $n+1$ równań liniowych z $n+1$ niewiadomymi $c_o, c_1,\ldots,c_n$ $$\sum_{j=0}^n c_j z_k^j=f(x_k)$$ gdzie $z_k=e^{ix_k}$.
Załóżmy, że węzły interpolacji są równoodległe czyli $x_k=\frac{2k\pi}{n+1}; (k=0,1,\ldots,n)$.
Ze względu na to, że funkcje $e^{ijx}$ są ortogonalne (w sensie iloczynu skalarnego) zagadnienie można potraktować jako rzut dowolnego wektora w przestrzeni na osie. Można pokazać, że współczynniki $c_j$ można wyznaczyć w sposób następujący: $$c_j=\dfrac{(f,e^{ijx})}{n+1}=\dfrac{1}{n+1}\sum_{k=0}^n f(x_k) e^{-ijx_k}$$ (zapis $(\bullet,\bullet)$ oznacza iloczyn skalarny.
Jak wiadomo, wielomian $t_n(x) = \sum_{j=0}^n c_j e^{ijx}$ przedstawiony może być w postaci alternatywnej: $$t_n(x)=\frac{1}{2}a_0 + \sum_{j=1}^m (a_j\cos jx + b_j \sin jx) + \delta \frac{1}{2}a_{m+1} \cos (m+1)x$$ gdzie $\delta=0$, a $m=\frac{1}{2}n$, gdy $n$ parzyste oraz $\delta=1$, a $m=\frac{1}{2}(n-1)$, gdy $n$ nieparzyste. Współczynniki $a_j$ oraz $b_j$ równe są odpowiednio:
$$\begin{split} a_j & =\frac{2}{n+1}\sum_{k=0}^n f(x_k) \cos jx_k\ b_j & =\frac{2}{n+1}\sum_{k=0}^n f(x_k) \sin jx_k \end{split}$$
Powyższy wzór nie jest wykorzystywany do obliczeń. Zazwyczaj stosuje się Szybką Transformatę Fouriera (FFT). Algorytm opracowany przez Cooleya i Tookeya w 1965 pominę. Koszt podejścia klasycznego wymaga $(n+1)^2$ operacji. Natomiast, w przypadku algorytmu FFT, gdy $n+1$ może być przedstawione (rozłożone) jako iloczyn $r_1r_2\ldots r_p$ będzie rzędu $(n+1)(r_1 + r_2 + \cdots + r_p)$. Algorytm FFT najczęściej stosowany jest, gdy $n+1=2^k$, wymaga on wówczas $n\log_2n$ działań.
Pewnym problemem stosowania interpolacji trygonometrycznej jest to, że funkcja jest okresowa. To znaczy zakładamy, że $t_n(0)=t_n(x_{n+1})$. Dokonując pomiarów zazwyczaj tak nie jest. W dalszym ciągu „można" stosować FFT, funkcja będzie okresowa, ale… Aby uwolnić się od tego problemu, nakłada się zmierzone dane funkcję okna. Typowa funkcja okna dla $x=0$ i $x=2\pi$ przyjmuje wartość zero co powoduje zniekształcenie badanego przebiegu.
Mathematica
Interpolacja trygonometryczna
Zajmę się głównie szybką transformatą Fouriera. Funkcja Fourier służy do wyliczenia FFT z zadanego przebiegu danych. Najprostsze jeż użycie będzie takie:
$\pmb{a=\text{Fourier}[{1,1,2,2,1,1,0,0}]}$
${2.82843, +0. i,-0.5+1.20711 i,0., +0. i,0.5, -0.207107 i,0., +0. i,0.5, +0.207107 i,0., +0. i,-0.5-1.20711 i}$
Odwrotna transformata Fouriera:
$\pmb{\text{InverseFourier}[a]}$
${1.,1.,2.,2.,1.,1.,0.,-\text{1.570092458683775$\grave{ }$*${}^{\wedge}$-16}}$ W każdym razie działa. Ósma wartość to „prawie” zero.
Teraz jakiś przebieg okresowy.
$\pmb{\text{err}=.0001;}$
Parametr err to zakres zaburzeń.
$\pmb{\text{fftdata}=\text{Fourier}[\text{data}];}$ $\pmb{\text{ListPlot}[\text{Abs}[\text{fftdata}]]}$
Jak zinterpretować ten wykres? Czemu są cztery ekstrema (piki)?
Można również skorzystać z funkcji Periodogram.
$\pmb{\text{Periodogram}[\text{data}]}$
Jak zinterpretować ten rezultat?
Jak na podstawie powyższych wykresów zidentyfikować częstości przebiegu?
Hint1: Wygenerować jeden okres przebiegu o częstości 1 Hz i pokazać jego transformatę Fouriera oraz periodogram.
Hint2: Funkcja periodogram ma dodatkowy parametr SampleRate. Mówi on ile próbek na jednostkę czasu wykonano. Używany jest do skalowania osi $X$.
Matlab
Interpolacja trygonometryczna
Podstawowym narzędziem jest funkcja fft:
Fs = 1000; % Sampling frequency
T = 1/Fs; % Sample time
L = 1000; % Length of signal
t = (0:L-1)*T; % Time vector
% Sum of a~50 Hz sinusoid and a~120 Hz sinusoid
x = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);
y = x + 2*randn(size(t)); % Sinusoids plus noise
plot(Fs*t(1:50),y(1:50))
title('Signal Corrupted with Zero-Mean Random Noise')
xlabel('time (milliseconds)')
NFFT = 2^nextpow2(L); % Next power of 2 from length of y
Y = fft(y,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);
% Plot single-sided amplitude spectrum.
plot(f,2*abs(Y(1:NFFT/2+1)))
title('Single-Sided Amplitude Spectrum of y(t)')
xlabel('Frequency (Hz)')
ylabel('|Y(f)|')
Zwracam uwagę, że pokazane rozwiązanie sugeruje, żeby długość ciągu danych była potęgą dwójki (polecenie: NFFT = 2^nextpow2(L);) czego nie wymaga Mathematica.
Zadania
Użyć interpolacji trygonometrycznej do znalezienia okresu zmienności tygodniowego (lub miesięcznego) wykresu temperatury. Wyglądać może on jakoś tak:
Przed rozpoczęciem tego zadania sugeruję zapoznanie się z informacjami zawartymi w rozdziale o Jupyterze oraz przeanalizowanie przykładowych notatników na temat wczytywania danych, eliminacji błędnych danych i szybkiej transformaty Fouriera.
Mając transformatę Fouriera warto spróbować przeprowadzić „naiwną, ręczną filtrację danych", która polegać będzie na wyzerowani nieistotnych składowych transformaty i po zostawieniu tylko tych które wnoszą istotny wkład w przebieg oraz porównaniu przebiegów.
Opisać od czego zależy rozdzielczość (zdolność do identyfikowania w przebiegu złożonym składowych o bardzo bliskich częstotliwościach) szybkiej transformaty Fouriera. Można posiłkować się notatnikiem Jupyter.
Zwrócić uwagę na znaczenie twierdzenia Shannona (o próbkowaniu) i wytłumaczyć na co mają wpływ:
- częstotliwość próbkowania,
- czas próbkowania.