MN08

Szybka transformacja Fouriera (FFT)

Algorytm FFT dla dyskretnej transformacji Fouriera (DFT) i jego krewniacy dla wyznaczania dyskretnej transformacji cosinusów (DCT i MDCT), choć dotyczą pozornie dość abstrakcyjnych zadań, zrewolucjonizowały wiele dziedzin życia. Między innymi, wykorzystuje się je w

• kompresji obrazów w formacie JPEG (DCT)
• kompresji dźwięku w formacie MP3 i pokrewnych (MDCT)
• rozwiązywaniu ważnych równań różniczkowych cząstkowych (DFT)

a także do

• fitrowania szumów (DFT)
• szybkiego mnożenia wielomianów (DFT)

W niniejszym wykładzie ograniczymy się do przedstawienia szybkiego algorytmu rozwiązywania zadania DFT, algorytmy rozwiązywania zadań pokrewnych (DCT, MDCT, itp.) opierają się na podobnych zasadach.

Zacznijmy od postawienia zadania obliczeniowego DFT. Jest ono (pozornie!) jednocześnie banalne i wydumane:

Dla danego zestawu $N$ liczb $f_{\alpha }\in C$, $\alpha =0,\dots ,N-1$, wyznaczyć $N$ wartości

$c_j=\sum _{\alpha =0}^{N-1}{f}_{\alpha }{\omega }_N^{j\alpha },$

dla $j=0,\dots ,N-1$, przy czym ${\omega }_N=\exp (-i\frac{2\mathrm{\Pi }}{N})$. Jak pamiętamy, jednostka urojona $i$ spełnia $i=\sqrt{-1}$. Taką operację nazywamy dyskretną transformacją Fouriera, DFT.

Ponieważ dane $f$ są wektorem $f=[f_0,\dots ,f_{N-1}{]}^T$, wynik $c$ też jest wektorem, a zadanie jest liniowe, możemy wszystko zapisać macierzowo:

$c=F_{N}f,$

gdzie

$F_N=({\omega }_N^{j\alpha }{)}_{j,\alpha =0,\dots ,N-1}=\left(\begin{array}{ccccc}1& 1& 1& \cdots & 1\\ 1& {\omega }_N& {\omega }_N^2& \cdots & {\omega }_N^{N-1}\\ 1& {\omega }_N^2& {\omega }_N^4& \cdots & {\omega }_N^{2(N-1)}\\ & & & \cdots & \\ 1& {\omega }_N^{N-1}& {\omega }_N^{2(N_1)}& \cdots & {\omega }_N^{(N-1{)}^2}\end{array}\right).$

Tak więc, gdyby naiwnie podejść do zadania, moglibyśmy je rozwiązać brutalnie, tworząc na początek macierz $F_N$, a następnie wyznaczając iloczyn macierzy przez wektor, co łatwo zrobić kosztem $O(N^2)$ operacji. Dlatego istotne jest, że algorytm FFT, który za chwilę omówimy, będzie działać kosztem $O(N\log N)$, czyli praktycznie liniowym (w obu przypadkach stałe proporcjonalności są równe 2) w dodatku będzie miał znacznie mniejsze wymagania pamięciowe. Już nawet gdy $N=8$ (przypadek kompresji JPEG), mamy $N^2=512$, gdy tymczasem $N\log N=24$, więc zysk jest ponaddwudziestokrotny, co ma kapitalne znaczenie, bo zamiast czekać na zapis zdjęcia (powiedzmy) pół sekundy, czekalibyśmy długie 10 (sic!) sekund... albo nasza(?) cyfrówka kosztowałaby majątek...

Fakt

Macierz $F$ jest symetryczna oraz $\frac{1}{N}{F}_N^{*}{F}_N=I$.

Zauważmy, że nasza macierz ma jeszcze więcej specjalnej struktury, powtarza się w niej bardzo wiele takich samych współczynników (sprawdź dla $N=4$, w ogólności ma tylko $N$ różnych wyrazów), gdyż ${\omega }_N^N=1$ (dla $j=0,\dots ,N-1$, ${\omega }_N^j$ to nic innego jak wszystkie zespolone pierwiastki z jedynki), a jak wiadomo, mnożenie przez 1 nic nie kosztuje --- o ile tylko się wie, że gdzie jest ta jedynka. Z grubsza na tym właśnie będzie polegać siła algorytmu FFT!

W wyprowadzeniu algorytmu szybkiej transformacji Fouriera, Fast Fourier Transform, FFT oprzemy się poraz kolejny na regule "dziel i rządź". Dla uproszczenia analizy przyjmiemy, że $N$ jest naturalną potęgą dwójki, w szczególności $N=2m$ dla pewnego naturalnego $m$.

Rzeczywiście, rozbijając naszą sumę na sumę po indeksach parzystych i sumę po indeksach nieparzystych, mamy

Parser nie mógł rozpoznać (nieznana funkcja „\aligned”): {\displaystyle \displaystyle \aligned c_j &= f_0 \omega_N^{0\cdot j} + f_2 \omega_N^{2\cdot j} + \ldots + f_{N-2} \omega_N^{(N-2)\cdot j}\\ && + f_1 \omega_N^{1\cdot j} + f_3 \omega_N^{3\cdot j} + \ldots + f_{N-1} \omega_N^{(N-1)\cdot j}. \endaligned}

Suma po indeksach parzystych da się zapisać w postaci

$\sum _{k=0}^{m-1}{f}_{2k}({\omega }_m{)}^{jk}$

i analogicznie suma po indeksach nieparzystych da się zapisać

${\omega }_N^j\cdot \sum _{k=0}^{m-1}{f}_{2k+1}({\omega }_m{)}^{jk}$

gdyż ${\omega }_N^2={\omega }_m$, a więc wygląda na to, że nasze zadanie wyznaczenia dyskretnej transformacji Fouriera wymiaru $N$ da się sprowadzić do analogicznych zadań mniejszego rozmiaru.

Rzeczywiście, korzystając z tego, że $j=\beta m+\tilde{j}$, gdzie $\tilde{j}=0,\dots ,m-1$, oraz $\beta =0$ lub $1$, mamy ${\omega }_m^j={\omega }_m^{\tilde{j}}$.

Oznaczając

Parser nie mógł rozpoznać (nieznana funkcja „\aligned”): {\displaystyle \displaystyle \aligned \Phi_{\widetilde{j}} &= \sum_{k=0}^{m-1} f_{2k} \omega_m^{\widetilde{j}k}, \qquad \widetilde{j} = 0,\ldots, m-1,\\ \Psi_{\widetilde{j}} &= \sum_{k=0}^{m-1} f_{2k+1} \omega_m^{\widetilde{j}k}, \endaligned}

(jak widać są to DFT dla dwa razy krótszych wektorów, złożonych z tylko parzystych lub tylko nieparzystych współrzędnych $f$) dostajemy ostatecznie

Parser nie mógł rozpoznać (nieznana funkcja „\aligned”): {\displaystyle \displaystyle \aligned c_{\widetilde{j}} = \Phi_{\widetilde{j}} + \omega_N^{\widetilde{j}} \Psi_{\widetilde{j}},\\ c_{\widetilde{j}+m} = \Phi_{\widetilde{j}} - \omega_N^{\widetilde{j}} \Psi_{\widetilde{j}}. \endaligned}

Tym samym, wyznaczenie DFT wymiaru $N$ udało się rzeczywiście sprowadzić do wyznaczenia dwóch DFT wymiaru $m=N/2$ oraz drobnej manipulacji na ich wynikach zgodnie z powyższym wzorem. Oczywiście, te mniejsze transformacje można wyznaczyć takim samym sposobem, co prowadzi do zależności rekurencyjnej, która kończy się na wektorach długości 1, na których DFT to po prostu identyczność.

Proste sprawdzenie pokazuje, że koszt takiego algorytmu jest rzędu $O(N\log N)$, a nie $O(N^2)$, jak dla naiwnego algorytmu mnożenia wektora przez gęstą macierz. Zysk jest więc, nawet dla niewielkich $N$, istotny.

Algorytm Prosta wersja algorytmu FFT

---

function y = fft(x)
	N = length(f);
	if N == 1 
		y = x;
	else
		<math>\displaystyle \omega</math> = <math>\displaystyle \exp(-\frac{2\Pi}{N}i)</math>;
		<math>\displaystyle \omega_k</math> = <math>\displaystyle \omega^{N/2-1}</math>;
		
		u = fft( f[0:2:N-2] );
		v = fft( f[1:2:N-1] );
		
		v = v * <math>\displaystyle \omega_k</math>;
		
		y = [ u+v ; u-v ];
	end
end

Efektywna implementacja FFT

Jak już zdążyliśmy się przyzwyczaić, w algorytmach numerycznych unikamy stosowania jawnej rekurencji, gdy tylko to możliwe. W przypadku FFT można jej również uniknąć, wyznaczając zawczasu --- korzystając z tzw. odwrócenia bitów --- porządek, w którym należy składać 1-wymiarowe DFT w coraz dłuższe wektory zgodnie ze wzorem powyżej tak, by na końcu dostać pożądany wektor $c=F_{N}f$.

Ponadto, istnieją warianty algorytmu FFT, które np. działają na danych rzeczywistych. Na analogicznych zasadach co FFT oparte są również algorytmy wykonujące tzw. dyskretną transformację cosinusów (DCT) i jej wariant MDCT stosowany w kodekach audio takich jak MP3, AAC, czy Ogg-Vorbis..

Interpolacja trygonometryczna

Jednym z wielu zadań, w których daje się zastosować algorytm FFT, jest zadanie interpolacji trygonometrycznej:

Dla danych węzłów $x_k=\frac{2\pi k}{N}$, $k=0,\dots ,N-1$, znaleźć wielomian $P$ (zmiennej rzeczywistej, o wartościach zespolonych) stopnia $\le N-1$ postaci

$P(x)=\sum _{k=0}^{N-1}{\beta }_k{\overline{\omega }}^k,$

gdzie $\overline{\omega }=e^{ix}=\cos (x)+i\sin (x)$ (stąd nazwa: trygonometryczny) taki, że

$P(x_k)=y_k,k=0,\dots ,N-1,$

dla zadanych wartości $y_k$.

W języku macierzy DFT możemy zapisać zadanie interpolacji trygonometrycznej jako zadanie wyznaczania wektora $\beta$ takiego, że

$\overline{F_N}\beta =y$

dla zadanego wektora $y$.

Twierdzenie

Współczynniki $\beta =[{\beta }_0,\dots ,{\beta }_{N-1}{]}^T$ poszukiwanego wielomianu trygonometrycznego wyrażają się wzorem

$\beta =\frac{1}{N}{F}_{N}y.$

Dowód

Jak wiemy, $\frac{1}{N}{F}_N^{*}{F}_N=I$, więc, z symetrii macierzy $F_N$ wynika, że $(\overline{F_N}{)}^{-1}=\frac{1}{N}{F}_N$.

Można pokazać, że gdy dane są rzeczywiste, zadanie interpolacji trygonometrycznej można wyrazić korzystając wyłącznie z liczb rzeczywistych. Jeśli $y_k\in R$, $k=0,\dots ,N-1$, to wtedy (rzeczywisty) wielomian trygonometryczny

$p(x)=\sum _{k=0}^{N-1}(a_k\cos (\frac{2k\pi }{N}x)-b_k\sin (\frac{2k\pi }{N}x)),$

gdzie $a_k=Re({\beta }_k)$, $b_k=Im({\beta }_k)$, interpoluje $y_k$ w węzłach $x_k$. Oczywiście, powyższa formuła w rzeczywistości ma o połowę mniej wyrazów, ze względu na własności funkcji trygonometrycznych.

Splot

W niektórych zastosowaniach potrzebne jest wyznaczenie splotu dwóch wektorów, to znaczy, wyznaczenie wyrażeń postaci

$z_k=\sum _{j=0}^{N-1}{u}_j{v}_{k-j},k=0,\dots ,N-1.$

Zapisując to macierzowo, szukamy iloczynu wektora $u$ z cykliczną macierzą Toeplitza (macierz Toeplitza ma stałe wyrazy wzdłuż diagonali) wyznaczoną przez $v$,

$z=\left(\begin{array}{cccc}{v}_0& v_1& \cdots & v_{N-1}\\ v_{N-1}& v_0& \cdots & v_{N-2}\\ ⋮& \ddots & \ddots & ⋮\\ v_1& \cdots & v_{N-1}& v_0\end{array}\right)\cdot \left(\begin{array}{c}{u}_0\\ u_1\\ ⋮\\ u_{N-1}\end{array}\right),$

tak więc pozornie zadanie wyznaczenia splotu powinno kosztować tyle, co mnożenie macierzy przez wektor, a więc $O(N^2)$ operacji. Tymczasem prosty rachunek pozwala sprawdzić, że odpowiednie transformacje Fouriera, $\hat{z}=F_{N}z$, $\hat{u}=F_{N}u$, $\hat{v}=F_{N}v$ spełniają równanie z macierzą diagonalną!

$\hat{z}=N\cdot \left(\begin{array}{ccc}{\hat{v}}_0& & \\ & \ddots & \\ & & {\hat{v}}_{N-1}\end{array}\right)\cdot \hat{u}=N\cdot \left(\begin{array}{c}{\hat{v}}_0{\hat{u}}_0\\ ⋮\\ {\hat{v}}_{N-1}{\hat{u}}_{N-1}\end{array}\right),$

a to mnożenie daje się wykonać kosztem liniowym, tym samym całe zadanie daje się policzyć kosztem $O(N\log N)$.

Biblioteki

Najpopularniejszą obecnie biblioteką implementującą algorytm FFT dla DFT, DCT i innych pokrewnych (bez ograniczenia, że wymiar jest naturalną potęgą dwójki), jest FFTW której nazwa jest skrótowcem od niezbyt skromnie brzmiącego ang. The Fastest Fourier Transform in the West. Z tej biblioteki korzystają m.in. funkcje MATLABa i Octave'a fft oraz ifft dla transformacji DFT (to znaczy mnożenia przez $F_N$) i, odpowiednio, transformacji odwrotnej, $\frac{1}{N}{\overline{F}}_N$.

FFTW jest napisana w C i w dużym stopniu wykorzystuje możliwości współczesnych procesorów, takie jak potokowanie i instrukcje wektorowe SSE2 i SSE3. Poniżej pokazujemy przykładowy prościutki kod w C, realizujący DFT na pojedynczym wektorze zespolonym.

 
/* Kompilacja:
	gcc -o dft dft.c -lfftw3 -lm */
#include <complex.h> /* rozszerzenie GCC dające operacje na typach zespolonych */
#include <fftw3.h>
#include <math.h>
#define N 8

int main(void)
{
fftw_complex *F, *C;
fftw_plan Plan;
double normfactor = 1.0/N;
int i;

F = fftw_malloc( N * sizeof(fftw_complex) );
C = fftw_malloc( N * sizeof(fftw_complex) );

for( i = 0; i < N; i++ ) /* inicjalizacja wartości tablicy F */
{
	F[i] = i*M_PI*(1-0.5*I); 
}

Plan = fftw_plan_dft_1d( N, F, C, FFTW_FORWARD, FFTW_ESTIMATE );

fftw_execute( Plan );

for( i = 0; i < N; i++ ) /* normalizacja wyznaczonego C */
	C[i] *= normfactor; 

for( i = 0; i < N; i++ )
{
	printf("F[\%d] = \%8.3lf + \%8.3lfi | C[\%d] = \%8.3lf + \%8.3lfi\n", 
	i, creal(F[i]), cimag(F[i]), i, creal(C[i]), cimag(C[i])); 
}

/* .... teraz moglibyśmy zmienić wartości F i ponownie wyznaczyć C,
korzystając z tego samego Planu! */

/* sprzątamy */

fftw_destroy_plan( Plan );
fftw_free( F ); 
fftw_free( C ); 

return(0);
}

Zwrócmy uwagę na linię Plan=fftw_plan_dft_1d(...). To tutaj dokonywane są ustalenia, w jakiej kolejności mają być prowadzone obliczenia. Jest to operacja dość kosztowna, dlatego jeśli mamy wyznaczyć wiele takich samych DFT, tylko na różnych danych, należy przed pierwszą DFT taki plan zachować, a potem już gotowy wykorzystać, gdy chcemy aplikować DFT dla następnych danych.

Interpolacja splajnowa

Interpolacja wielomianami interpolacyjnymi, chociaż korzysta z funkcji gładkich i łatwo reprezentowalnych w komputerze, ma jednak również pewne wady. Zauważmy, że błąd interpolacji może być bardzo duży (zjawisko Rungego), a poza tym, interpolacja jest nielokalna: nawet mała zmiana warości funkcji w pojedynczym węźle może powodować dużą zmianę zachowania całego wielomianu interpolacyjnego. Czasem więc lepiej jest zastosować innego rodzaju interpolację, np. posługując się funkcjami sklejanymi.

Co to są funkcje sklejane?

W ogólności przez funkcję sklejaną rozumie się każdą funkcję przedziałami wielomianową. Nas będą jednak interesować szczególne funkcje tego typu i dlatego termin funkcje sklejane zarezerwujemy dla funkcji przedziałami wielomianowych i posiadających dodatkowe własności, które teraz określimy.

Niech dany będzie przedział skończony $[a,b]$ i węzły

$a=x_0<x_1<\cdots <x_n=b,$

przy czym $n\ge 1$.

Definicja

Funkcję Parser nie mógł rozpoznać (nieznana funkcja „\toR”): {\displaystyle \displaystyle s:R\toR} nazywamy funkcją sklejaną rzędu $r$ ($r\ge 1$) odpowiadającą węzłom $x_j$, $0\le j\le n$, jeśli spełnione są następujące dwa warunki:

(i)

$s$ jest wielomianem stopnia co najwyżej $2r-1$ na każdym

z przedziałów $[x_{j-1},x_j]$, tzn. $s_{|[x_{j-1},x_j]}\in {\mathrm{\Pi }}_{2r-1}$, $1\le j\le n$,

(ii)

$s$ jest $(2r-2)$-krotnie różniczkowalna w sposób

ciągły na całej prostej, tzn. $s\in C^{(2r-2)}(R)$.

Jeśli ponadto

(iii)

$s$ jest wielomianem stopnia co najwyżej $r-1$ poza

$(a,b)$, tzn. $s_{|(-\mathrm{\infty },a]},s_{|[b,+\mathrm{\infty })}\in {\mathrm{\Pi }}_{r-1}$,

to $s$ jest naturalną funkcją sklejaną.

Klasę naturalnych funkcji sklejanych rzędu $r$ opartych na węzłach $x_j$ będziemy oznaczać przez ${{𝒮}}_r(x_0,\dots ,x_n)$, albo po prostu ${{𝒮}}_r$, jeśli węzły są ustalone.

Na przykład, funkcją sklejaną rzędu pierwszego ($r=1$) jest funkcja ciągła i liniowa na poszczególnych przedziałach $[x_{j-1},x_j]$. Jest ona naturalna, gdy poza $(a,b)$ jest funkcją stała. Tego typu funkcje nazywamy liniowymi funkcjami sklejanymi.

Najważniejszymi a praktycznego punktu widzenia są jednak funkcje sklejane rzędu drugiego odpowiadające $r=2$. Są to funkcje, które są na $R$ dwa razy różniczkowalne w sposób ciągły, a na każdym z podprzedziałów są wielomianami stopnia co najwyżej trzeciego. W tym przypadku mówimy o kubicznych funkcjach sklejanych. Funkcja sklejana kubiczna jest naturalna, gdy poza $(a,b)$ jest wielomianem liniowym.

Interpolacja i gładkość

Pokażemy najpierw ważny lemat, który okaże się kluczem do dowodu dalszych twierdzeń.

Niech $W^r(a,b)$ będzie klasą funkcji Parser nie mógł rozpoznać (nieznana funkcja „\toR”): {\displaystyle \displaystyle f:[a,b]\toR} takich, że $f$ jest $(r-1)$ razy różniczkowalna na $[a,b]$ w sposób ciągły oraz $f^{(r)}(x)$ istnieje prawie wszędzie na $[a,b]$ i jest całkowalna z kwadratem, tzn.

Parser nie mógł rozpoznać (nieznana funkcja „\aligned”): {\displaystyle \displaystyle \aligned W^r(a,b) &= \{\,f\in C^{(r-1)}([a,b]):\, f^{(r)}(x) \mbox{ istnieje p.w. na } [a,b] \\ && \qquad\qquad\qquad\qquad \mbox{ oraz } f^{(r)}\in{\cal L}_2(a,b)\,\}. \endaligned}

Oczywiście każda funkcja sklejana rzędu $r$ (niekoniecznie naturalna) należy do $W^r(a,b)$.

Lemat

Niech $f\in W^r(a,b)$ będzie funkcją zerującą się w węzłach, tzn.

$f(x_j)=0,0\le j\le n.$

Wtedy dla dowolnej naturalnej funkcji sklejanej $s\in {{𝒮}}_r$ mamy

$\underset{a}{\overset{b}{\int }}{f}^{(r)}(x)s^{(r)}(x)dx=0.$

Dowód

Dla $r=1$ funkcja $s^{\prime }$ jest przedziałami stała. Oznaczając przez $a_j$ jej wartość na $[x_{j-1},x_j]$ dostajemy

Parser nie mógł rozpoznać (nieznana funkcja „\aligned”): {\displaystyle \displaystyle \aligned \int_a^b f'(x)s'(x)\,dx &= \sum_{j=1}^n\int_{t_{j-1}}^{t_j} f'(x)a_j\,dx \\ &= \sum_{j=1}^n a_j(f(x_j)-f(x_{j-1}))\,=\,0, \endaligned}

ponieważ $f$ zeruje się w $t_j$.

Rozpatrzmy teraz przypadek $r\ge 2$. Ponieważ

$(f^{(r-1)}{s}^{(r)}{)}^{\prime }=f^{(r)}{s}^{(r)}+f^{(r-1)}{s}^{(r+1)},$

to

$\underset{a}{\overset{b}{\int }}{f}^{(r)}(x)s^{(r)}(x)dx=f^{(r-1)}(x)s^{(r)}(x){|}_a^b-\underset{a}{\overset{b}{\int }}{f}^{(r-1)}(x)s^{(r+1)}(x)dx.$

Wobec tego, że $s$ jest poza przedziałem $(a,b)$ wielomianem stopnia co najwyżej $r-1$ oraz $s^{(r)}$ jest ciągła na $R$, mamy $s^{(r)}(a)=0=s^{(r)}(b)$, a stąd

$f^{(r-1)}(x)s^{(r)}(x){|}_a^b=0.$

Postępując podobnie, tzn. całkując przez części $r-1$ razy otrzymujemy w końcu

$\underset{a}{\overset{b}{\int }}{f}^{(r)}(x)s^{(r)}(x)dx=(-1{)}^{r-1}\underset{a}{\overset{b}{\int }}{f}^{\prime }(x)s^{(2r-1)}(x)dx.$

Funkcja $s^{(2r-1)}$ jest przedziałami stała, a więc możemy teraz zastosować ten sam argument jak dla $r=1$, aby pokazać,

że ostatnia całka jest równa zeru.

Funkcje sklejane chcielibyśmy zastosować do interpolacji funkcji. Ważne jest więc, aby odpowiednie zadanie interpolacyjne miało jednoznaczne rozwiązanie.

Twierdzenie

Jeśli $n+1\ge r$, to dla dowolnej funkcji Parser nie mógł rozpoznać (nieznana funkcja „\toR”): {\displaystyle \displaystyle f:[a,b]\toR} istnieje dokładnie jedna naturalna funkcja sklejana $s_f\in {{𝒮}}_r$ interpolująca $f$ w węzłach $x_j$, tzn. taka, że

$s_f(x_j)=f(x_j),0\le j\le n.$

Dowód

Pokażemy najpierw, że jedyną naturalną funkcją sklejaną interpolującą dane zerowe jest funkcja zerowa. Rzeczywiście, jeśli $s(x_j)=0$ dla $0\le j\le n$, to podstawiając w poprzednim Lemacie $f=s$ otrzymujemy

$\underset{a}{\overset{b}{\int }}(s^{(r)}(x){)}^{2}dx=0.$

Stąd $s^{(r)}$ jest funkcją zerową, a więc $s$ jest wielomianem stopnia co najwyżej $r-1$ zerującym się w co najmniej $n+1$ punktach $x_j$. Wobec tego, że $n+1>r-1$, otrzymujemy $s\equiv 0$.

Zauważmy teraz, że problem znalezienia naturalnej funkcji sklejanej $s$ interpolującej $f$ można sprowadzić do rozwiązania kwadratowego układu równań liniowych. Na każdym przedziale $[x_{i-1},x_i]$, $1\le i\le n$, jest ona postaci

$s(x)=w_i(x)=\sum _{j=0}^{2r-1}{a}_{i,j}{x}^j,$

dla pewnych współczynników Parser nie mógł rozpoznać (nieznana funkcja „\inR”): {\displaystyle \displaystyle a_{i,j}\inR} , a na $(-\mathrm{\infty },a]$ i $[b,\mathrm{\infty })$ mamy odpowiednio

$s(x)=w_0(x)=\sum _{j=0}^{r-1}{a}_{0,j}{x}^j$

i

$s(x)=w_{n+1}(x)=\sum _{j=0}^{r-1}{a}_{n+1,j}{x}^j.$

Aby wyznaczyć $s$, musimy więc znaleźć ogółem $2r(n+1)$ współczynników $a_{i,j}$, przy czym są one związane $(2r-1)(n+1)$ warunkami jednorodnymi wynikającymi z gładkości,

$w_i^{(k)}(x_i)=w_{i+1}^{(k)}(x_i)$

dla $0\le i\le n$ i $0\le k\le 2r-2$, oraz $n+1$ niejednorodnymi warunkami interpolacyjnymi,

$w_i(x_i)=f(x_i)$

dla $0\le i\le n$. Otrzymujemy więc układ $2r(n+1)$ równań liniowych ze względu na $2r(n+1)$ niewiadomych $a_{i,j}$.

Naturalna funkcja sklejana interpolująca $f$ jest wyznaczona jednoznacznie wtedy i tylko wtedy, gdy układ ten ma jednoznaczne rozwiązanie. To zaś zachodzi, gdy zero jest jedynym rozwiązaniem układu jednorodnego. Rzeczywiście, układ jednorodny odpowiada zerowym warunkom interpolacyjnym, przy których, jak pokazaliśmy wcześniej, zerowa funkcja sklejana (której odpowiada $a_{i,j}=0$, $\mathrm{\forall }i,j$)

jest jedynym rozwiązaniem zadania interpolacyjnego.

Naturalne funkcje sklejane możemy więc używać do interpolacji funkcji. Pokażemy teraz inną ich własność, która jest powodem dużego praktycznego zainteresowania funkcjami sklejanymi.

Twierdzenie

Niech $f\in W^r(a,b)$ i niech $s_f\in {{𝒮}}_r$ będzie naturalną funkcją sklejaną rzędu $r$ interpolującą $f$ w węzłach $x_j$, $0\le j\le n$. Wtedy

$\underset{a}{\overset{b}{\int }}(f^{(r)}(x){)}^{2}dx\ge \underset{a}{\overset{b}{\int }}(s_f^{(r)}(x){)}^{2}dx.$

Dowód

Jeśli przedstawimy $f$ w postaci $f=s_f+(f-s_f)$, to

Parser nie mógł rozpoznać (nieznana funkcja „\aligned”): {\displaystyle \displaystyle \aligned \int_a^b\Big(f^{(r)}(x)\Big)^2\,dx &= \int_a^b\Big(s_f^{(r)}(x)\Big)^2\,dx\,+\, \int_a^b\Big((f-s_f)^{(r)}(x)\Big)^2\,dx \\ & & \qquad\qquad 2\,\int_a^b s_f^{(r)}(x)(f-s_f)^{(r)}(x)\,dx. \endaligned}

Funkcja $f-s_f$ jest w klasie $W^r(a,b)$ i zeruje się w węzłach $x_j$, $0\le j\le n$. Z Lematu wynika więc, że

$\underset{a}{\overset{b}{\int }}{s}_f^{(r)}(x)(f-s_f{)}^{(r)}(x)dx=0$, a stąd wynika teza.

Wartość całki $\underset{a}{\overset{b}{\int }}(f^{(r)}(x){)}^{2}dx$ może być w ogólności uważana za miarę gładkości funkcji. Dowiedzioną nierówność możemy więc zinterpretować w następujący sposób. Naturalna funkcja sklejana jest w klasie $W^r(a,b)$ najgładszą funkcją spełniającą dane warunki interpolacyjne w wybranych węzłach $x_j$.

Jak już wspomnieliśmy, najczęściej używanymi są kubiczne funkcje sklejane. Dlatego rozpatrzymy je oddzielnie.

Kubiczne funkcje sklejane

Jeśli zdecydowaliśmy się na użycie kubicznych funkcji sklejanych to powstaje problem wyznaczenia $s_f\in {{𝒮}}_2$ interpolującej daną funkcję $f$, tzn. takiej, że $s_f(x_i)=f(x_i)$ dla $0\le i\le n$. W tym celu, na każdym przedziale $[x_i,x_{i+1}]$ przedstawimy $s_f$ w postaci jej rozwinięcia w szereg Taylora w punkcie $x_i$,

$s_f(x)=w_i(x)=a_i+b_i(x-x_i)+c_i\frac{(x-x_i{)}^2}{2}+d_i\frac{(x-x_i{)}^3}{6},$

i podamy algorytm obliczania $a_i,b_i,c_i,d_i$ dla $0\le i\le n-1$.

Warunki brzegowe i warunki ciągłości dla ${s_f}^{″}$ dają nam ${w_0}^{″}(x_0)=0=w_{n^{″}-1}(x_n)$ oraz ${w_i}^{″}(x_{i+1})=w_{i^{″}+1}(x_{i+1})$, czyli

Parser nie mógł rozpoznać (nieznana funkcja „\aligned”): {\displaystyle \displaystyle \aligned c_0 &= 0, \\ c_i+d_ih_i &= c_{i+1}, \qquad 0\le i\le n-2, \\ c_{n-1}+d_{n-1}h_{n-1} &= 0, \endaligned}

gdzie $h_i=x_{i+1}-x_i$. Stąd, przyjmując dodatkowo $c_n=0$, otrzymujemy

$d_i=\frac{c_{i+1}-c_i}{h_i},1\le i\le n-1.$

Z warunków ciągłości dla ${s_f}^{\prime }$ dostajemy z kolei

$b_i+c_i{h}_i+d_i{h}_i^2/2=b_{i+1},0\le i\le n-2,$

oraz

$b_{i+1}=b_i+h_i(c_{i+1}+c_i)/2,0\le i\le n-2.$

Warunki ciągłości $s_f$ dają w końcu

$a_i+b_i{h}_i+c_i{h}_i^2/2+d_i{h}_i^3/6=a_{i+1},0\le i\le n-2.$

Powyższe równania definiują nam na odcinku $[a,b]$ naturalną kubiczną funkcję sklejaną. Ponieważ poszukiwana funkcja sklejana $s_f$ ma interpolować $f$, mamy dodatkowych $n+1$ warunków interpolacyjnych $w_i(x_i)=f(x_i)$, $0\le i\le n-1$, oraz $w_{n-1}(x_n)=f(x_n)$, z których

$a_i=f(x_i),0\le i\le n-1.$

Teraz możemy warunki ciągłości przepisać jako

$f(x_{i+1})=f(x_i)+b_i{h}_i+c_i{h}_i^2+d_i{h}_i^3/6,$

przy czym wzór ten zachodzi również dla $i=n-1$. Po wyrugowaniu $b_i$ i podstawieniu $d_i$ z (Uzupelnic: dei ), mamy

$b_i=f(x_i,x_{i+1})+h_i(c_{i+1}+2c_i)/6,0\le i\le n-1,$

gdzie $f(x_i,x_{i+1})$ jest odpowiednią różnicą dzieloną. Możemy teraz powyższe wyrażenie na $b_i$ podstawić, aby otrzymać

$c_i{h}_i/6+c_{i+1}(h_i+h_{i+1})/3+c_{i+1}{h}_{i+1}/6=f(x_{i+1},x_{i+2})-f(x_i,x_{i+1}).$

Wprowadzając oznaczenie

$c_i^{*}=c_i/6,$

możemy to równanie przepisać jako

$\frac{h_i}{h_i+h_{i+1}}{c}_i^{*}+2c_{i+1}^{*}+\frac{h_{i+1}}{h_i+h_{i+1}}{c}_{i+2}^{*}=f(x_i,x_{i+1},x_{i+2}),$

$0\le i\le n-2$, albo w postaci macierzowej

$\left(\begin{array}{cc}2& w_1\\ u_2& 2& w_2\\ & u_3& 2& w_3\\ & & \ddots & \ddots & \ddots \\ & & & u_{n-2}& 2& w_{n-2}\\ & & & & u_{n-1}& 2\end{array}\right)\left(\begin{array}{c}{c}_1^{*}\\ c_2^{*}\\ c_3^{*}\\ ⋮\\ c_{n-2}^{*}\\ c_{n-1}^{*}\end{array}\right)=\left(\begin{array}{c}{v}_1\\ v_2\\ v_3\\ ⋮\\ v_{n-2}\\ v_{n-1}\end{array}\right),$

gdzie

Parser nie mógł rozpoznać (nieznana funkcja „\aligned”): {\displaystyle \displaystyle \aligned && u_i\,=\,\frac{h_{i-1}}{h_{i-1}+h_i},\qquad w_i\,=\,\frac{h_i}{h_{i-1}+h_i}, \\ && v_i\,=\,f(x_{i-1},x_i,x_{i+1}). \endaligned}

Ostatecznie, aby znaleźć współczynniki $a_i,b_i,c_i,d_i$ należy najpierw rozwiązać układ równań liniowych, a potem zastosować wzory definiujące pozostałe współczynniki.

Zauważmy, że macierz układu równań liniowych jest trójdiagonalna i ma dominującą przekątną. Układ można więc rozwiązać kosztem proporcjonalnym do wymiaru $n$ używając algorytmu przeganiania. Koszt znalezienia wszystkich współczynników kubicznej funkcji sklejanej interpolującej $f$ jest więc też proporcjonalny do $n$.

Na końcu oszacujemy jeszcze błąd interpolacji naturalnymi kubicznymi funkcjami sklejanymi na przedziale $[a,b]$. Będziemy zakładać, że $f$ jest dwa razy różniczkowalna w sposób ciągły.

Twierdzenie

J eśli $f\in F_M^1([a,b])$ to

$\Vert f-s_f{\Vert }_{C([a,b])}\le 5M\underset{1\le i\le n}{\max }(x_i-x_{i-1}{)}^2.$

W szczególności, dla podziału równomiernego Parser nie mógł rozpoznać (błąd składni): {\displaystyle \displaystyle x_i=a+(b-a)i/ń} , $0\le i\le n$, mamy

$\Vert f-s_f{\Vert }_{C([a,b])}\le 5M(\frac{b-a}{n}{)}^2.$

Dowód

Wykorzystamy obliczoną wcześniej postać interpolującej funkcji sklejanej $s_f$. Dla $x\in [x_i,x_{i+1}]$ mamy

Parser nie mógł rozpoznać (nieznana funkcja „\aligned”): {\displaystyle \displaystyle \aligned w_i(x) &= f(x_i)\,+\,\left(\frac{f(x_{i+1})-f(x_i)}{h_i} -h_i(c_{i+1}^*+2c_i^*)\right)(x-x_i) \\ &&\qquad\qquad \,+\, 3c_i^*(x-x_i)^2\,+\,\frac{c_{i+1}^*-c_i^*}{h_i}(x-x_i)^3. \endaligned}

Z rozwinięcia $f$ w szereg Taylora w punkcie $x_i$ dostajemy $f(x)=f(x_i)+f^{\prime }(x_i)(x-x_i)+f^{″}({\xi }_1)(x-x_i{)}^2/2$ oraz $(f(x_{i+1})-f(x_i))/h_i=f^{\prime }(x_i)+h_i{f}^{″}({\xi }_2)/2$. Stąd

Parser nie mógł rozpoznać (nieznana funkcja „\aligned”): {\displaystyle \displaystyle \aligned \lefteqn{f(x)-s_f(x) \,=\, f(x)-w_i(x)} \\ &= \frac{f''(\xi_1)}2(x-x_i)^2-\left(\frac{f''(\xi_2)}2 -(c_{i+1}^*+2c_i^*)\right)h_i(x-x_i) \\ & & \qquad\qquad\qquad -3c_i^*(x-x_i)^2 -\frac{c_{i+1}^*-c_i^*}{h_i}(x-x_i)^3, \endaligned}

oraz

Parser nie mógł rozpoznać (nieznana funkcja „\aligned”): {\displaystyle \displaystyle \aligned |f(x)-s_f(x)| &\le &(M+2|c_{i+1}^*|+6|c_i^*|)h_i^2 \\ &= (M+8\max_{1\le i\le n-1}|c_i^*|)h_i^2. \endaligned}

Niech teraz $\underset{1\le i\le n-1}{\max }|c_i^{*}|=|c_s^{*}|$. Z postaci układu (Uzupelnic: ukltrd ) otrzymujemy

Parser nie mógł rozpoznać (nieznana funkcja „\aligned”): {\displaystyle \displaystyle \aligned |c_s^*| &= 2|c_s^*|-|c_s^*|(u_s+w_s) \,\le\, |u_sc_{s-1}^*+2c_s^*+w_sc_{s+1}| \\ &= |f(x_{s-1},x_s,x_{s+1})|\,\le\, \Big|\frac{f''(\xi_3)}2\Big|\,\le\,\frac 12 M, \endaligned}

a stąd i z (Uzupelnic: psik )

$|f(x)-s_f(x)|\le 5M{h}_i^2,$

co kończy dowód.

Uwaga

Zamiast terminu funkcje sklejane używa się też często terminów splajny (ang. spline-sklejać), albo funkcje gięte.

Uwaga

Niech

$W_M^r(a,b)=\{f\in W^r(a,b):\underset{a}{\overset{b}{\int }}{(f^{(r)}(x))}^{2}dx\le M\}.$

Ustalmy węzły $a=x_0<\cdots <x_n=b$. Dla $f\in W_M^r(a,b)$, niech $s_f$ będzie naturalną funkcją sklejaną interpolującą $f$ w $x_j$, $0\le j\le n$, a $a_f$ dowolną inną aproksymacją korzystającą jedynie z informacji o wartościach $f$ w tych węzłach , tzn.

$a_f=\varphi (f(x_0),\dots ,f(x_n)).$

Załóżmy, że błąd aproksymacji mierzymy nie w normie Czebyszewa, ale w normie średniokwadratowej, zdefiniowanej jako

$\Vert g{\Vert }_{{{ℒ}}_2(a,b)}=\sqrt{\underset{a}{\overset{b}{\int }}(g(x){)}^{2}dx}.$

Wtedy

$\underset{f\in W_M^r(a,b)}{\sup }\Vert f-s_f{\Vert }_{{{ℒ}}_2(a,b)}\le \underset{f\in W_M^r(a,b)}{\sup }\Vert f-a_f{\Vert }_{{{ℒ}}_2(a,b)}.$

Aproksymacja naturalnymi funkcjami sklejanymi jest więc optymalna w klasie $W_M^r(a,b)$.

Można również pokazać, że interpolacja $s_f^{*}$ naturalnymi funkcjami sklejanymi na węzłach równoodległych Parser nie mógł rozpoznać (błąd składni): {\displaystyle \displaystyle x_j=a+(b-a)j/ń} , $0\le j\le n$, jest optymalna co do rzędu w klasie $W_M^r(a,b)$, wśród wszystkich aproksymacji korzystających jedynie z informacji o wartościach funkcji w $n+1$ dowolnych punktach, oraz

$\underset{f\in W_M^r(a,b)}{\max }\Vert f-s_f^{*}{\Vert }_{{{ℒ}}_2(a,b)}\asymp n^{-r}.$

Uwaga

Tak jak wielomiany, naturalne funkcje sklejane interpolujące dane funkcje można reprezentować przez ich współczynniki w różnych bazach. Do tego celu można na przykład użyć bazy kanonicznej $K_j$, $0\le j\le n$, zdefiniowanej równościami

$K_j(x_i)=\{\begin{array}{ll}0& i\ne j,\\ 1& i=j,\end{array}$

przy której $s_f(x)=\sum _{j=0}^{n}f(x_j)K_j(x)$. Baza kanoniczna jest jednak niewygodna w użyciu, bo funkcje $K_j$ w ogólności nie zerują się na żadnym podprzedziale, a tym samym manipulowanie nimi jest utrudnione.

Częściej używa się bazy B-sklejanej $B_j$, $0\le j\le n$. W przypadku funkcji kubicznych, $r=2$, jest ona zdefiniowana przez następujące warunki:

Parser nie mógł rozpoznać (nieznana funkcja „\aligned”): {\displaystyle \displaystyle \aligned B_j(x_j) &= 1, \qquad \mbox{ dla } 0\le j\le n, \\ B_j(x) &= 0,\qquad \mbox{ dla } x\le x_{j-2}, j\ge 2, \mbox{ oraz dla } x\ge x_{j+2}, j\le n-2. \endaligned}

Dla $B_0$ i $B_1$ dodatkowo żądamy, aby

${B_0}^{″}(x_0)=0={B_1}^{″}(x_0),B_1(x_0)=0,$

a dla $B_{n-1}$ i $B_n$ podobnie

$B_{n^{″}-1}(x_n)=0={B_n}^{″}(x_n),B_{n-1}(x_{n-1})=0.$

Wtedy $B_j$ nie zeruje się tylko na przedziale $(x_{j-2},x_{j+2})$, Wyznaczenie współczynników rozwinięcia w bazie $\{B_i{\}}_{i=0}^n$ funkcji sklejanej interpolującej $f$ wymaga rozwiązania układu liniowego z macierzą trójdiagonalną $\{B_j(x_i){\}}_{i,j=0}^n$, a więc koszt obliczenia tych współczynników jest proporcjonalny do $n$.

Uwaga

Oprócz naturalnych funkcji sklejanych często rozpatruje się też okresowe funkcje sklejane. Są to funkcje Parser nie mógł rozpoznać (nieznana funkcja „\toR”): {\displaystyle \displaystyle \tilde s:R\toR} spełniające warunki (i), (ii) z Rozdziału Uzupelnic: warfs , oraz warunek:

(iii)'

${\tilde{s}}^{(i)}$ jest dla $0\le i\le r-1$

funkcją okresową o okresie $(b-a)$, tzn. ${\tilde{s}}^{(i)}(x)={\tilde{s}}^{(i)}(x+(b-a))$, $\mathrm{\forall }x$.

Klasę okresowych funkcji sklejanych rzędu $r$ oznaczymy przez ${\widetilde{{𝒮}}}_r$. Funkcje te mają podobne własności jak naturalne funkcje sklejane. Dokładniej, niech

${\tilde{W}}^r(a,b)=\{f\in W^r(a,b):f^{(i)}(a)=f^{(i)}(b),0\le i\le r-1\},$

tzn. ${\tilde{W}}^r(a,b)$ jest klasą funkcji z $W^r(a,b)$, które można przedłużyć do funkcji, krórych wszystkie pochodne do rzędu $r-1$ włącznie są $(b-a)$-okresowe na $R$. Wtedy dla dowolnej funkcji $f\in {\tilde{W}}^r(a,b)$ zerującej się w węzłach $x_j$, oraz dla dowolnej $\tilde{s}\in {\widetilde{{𝒮}}}_r$ mamy

$\underset{a}{\overset{b}{\int }}{f}^{(r)}(x){\tilde{s}}^{(r)}(x)dx=0.$

Jest to odpowiednik Lematu Uzupelnic: bwazny w przypadku okresowym. W szczególności wynika z niego jednoznaczność rozwiązania zadania interpolacyjnego dla okresowych funkcji $f$ (tzn. takich, że $f(a)=f(b)$), jak również odpowiednia własność minimalizacyjna okresowych funkcji sklejanych. Dokładniej, jeśli $f\in {\tilde{W}}^r(a,b)$ oraz ${\tilde{s}}_f\in {\widetilde{{𝒮}}}_r$ interpoluje $f$ w węzłach $x_j$, $0\le j\le n$, to

$\underset{a}{\overset{b}{\int }}(f^{(r)}(x){)}^{2}dx\ge \underset{a}{\overset{b}{\int }}({\tilde{s}}_f^{(r)}(x){)}^{2}dx.$

Uwaga

Klasyczne zadanie aproksymacyjne w przestrzeniach funkcji definiuje się w następujący sposób.

Niech $F$ będzie pewną przestrzenią liniową funkcji Parser nie mógł rozpoznać (nieznana funkcja „\toR”): {\displaystyle \displaystyle f:[a,b]\toR} , w której określona została norma $\Vert \cdot \Vert$. Niech $V_n\subset F$ będzie podprzestrzenią w $F$ wymiaru $n$. Dla danej $f\in F$, należy znaleźć funkcję $v_f\in F$ taką, że

$\Vert f-v_f\Vert =\underset{v\in V_n}{\min }\Vert f-v\Vert .$

Okazuje się, że tak postawione zadanie ma rozwiązanie $v_f$, choć nie zawsze jest ono wyznaczone jednoznacznie, zob. ćw. Uzupelnic: rkl .

Jako przykład, rozpatrzmy $F=W^r(a,b)$. Utożsamiając funkcje $f_1,f_2\in W^r(a,b)$ takie, że $f_1(x)-f_2(x)\in {\mathrm{\Pi }}_{r-1}$, zdefiniujemy w $W^r(a,b)$ normę

$\Vert f\Vert =\sqrt{\underset{a}{\overset{b}{\int }}{(f^{(r)}(x))}^{2}dx}.$

Dla ustalonych węzłów $a=x_0<\cdots <x_n=b$, niech

$V_{n+1}={{𝒮}}_r$

będzie podprzestrzenią w $W^r(a,b)$ naturalnych funkcji sklejanych rzędu $r$ opartych węzłach $x_j$, $0\le j\le n$. Oczywiście $\text{dim}{{𝒮}}_r=n+1$, co wynika z jednoznaczności rozwiązania w ${{𝒮}}_r$ zadania interpolacji. Okazuje się, że wtedy optymalną dla $f\in W^r(a,b)$ jest naturalna funkcja sklejana $s_f$ interpolująca $f$ w węzłach $x_j$, tzn.

$\Vert f-s_f\Vert =\underset{s\in {{𝒮}}_r}{\min }\Vert f-s\Vert .$

Rzeczywiście, ponieważ norma w przestrzeni $W^r(a,b)$ generowana jest przez iloczyn skalarny

$(f_1,f_2)=\underset{a}{\overset{b}{\int }}{f}_1^{(r)}(x)f_2^{(r)}(x)dx,$

jest to przestrzeń unitarna. Znane twierdzenie mówi, że w przestrzeni unitarnej najbliższą danej $f$ funkcją w dowolnej domkniętej podprzestrzeni $V$ jest rzut prostopadły $f$ na $V$, albo równoważnie, taka funkcja $v_f\in V_{n+1}$, że iloczyn skalarny

$(f-v_f,v)=0,\mathrm{\forall }v\in V.$

W naszym przypadku, ostatnia równość jest równoważna

$\underset{a}{\overset{b}{\int }}(f-v_f{)}^{(r)}(x)s^{(r)}(x)dx=0,\mathrm{\forall }s\in {{𝒮}}_r.$

To zaś jest na mocy Lematu Uzupelnic: bwazny prawdą gdy $v_f$ interpoluje $f$ w punktach $x_j$, czyli $v_f=s_f$.

Dodajmy jeszcze, że nie zawsze interpolacja daje najlepszą aproksymację w sensie klasycznym, zob. ćw. Uzupelnic: intkla .