Tekst źródłowy strony MN10


<!-- 
Konwertowane  z pliku LaTeX przez latex2mediawiki, zob. http://www.ii.uj.edu.pl/&nbsp;pawlik1/latex2mediawiki.php.
Niezb�dne rozszerzenia i modyfikacje oryginalnego latex2mediawiki 
wprowadzi� przykry@mimuw.edu.pl
-->
 
=Szybka transformacja Fouriera (FFT)=

{{powrot |Metody numeryczne | do strony głównej 
przedmiotu <strong>Metody numeryczne</strong>}}

Algorytm FFT (''Fast Fourier Transform'') dla dyskretnej transformacji Fouriera (DFT, ''Discrete Fourier Transform'') i jego krewniacy dla
wyznaczania dyskretnej transformacji cosinusów (DCT, ''Discrete Cosine Transform'' i MDCT, ''Modified Discrete Cosine Transform''), 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 
* filtrowania szumów (DFT)
* [[Zaawansowane_algorytmy_i_struktury_danych/Wykład_4#Mnożenie wielomianów w punktach|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!)
banalne i jednocześnie wydumane:

Dla danego zestawu <math>N</math> liczb <math>f_\alpha\in C</math>, <math>\alpha = 0,\ldots, N-1</math>,
wyznaczyć  <math>N</math> wartości

<center><math>c_j = \sum_{\alpha = 0}^{N-1} f_\alpha \omega_N^{j\alpha}</math>,</center>

dla <math>j=0,\ldots, N-1</math>, przy czym <math>\omega_N = \exp(-i\frac{2\Pi}{N})</math>. Jak
pamiętamy, jednostka urojona <math>i</math> spełnia <math>i = \sqrt{-1}</math>. Taką operację nazywamy
<strong>dyskretną transformacją Fouriera</strong>, DFT.

[[grafika:Fourier.jpg|thumb|right||Jean Baptiste Joseph Fourier<br>  [[Biografia Fourier|Zobacz biografię]]]]

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

<center><math>c = F_N f</math>,</center>

gdzie

<center><math>
F_N = (\omega_N^{j\alpha})_{j,\alpha = 0,\ldots, N-1} = 
\begin{pmatrix} 
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{pmatrix} </math></center>

Tak więc, gdyby naiwnie podejść do zadania, moglibyśmy je rozwiązać brutalnie,
tworząc na początek macierz <math>F_N</math>, a następnie wyznaczając iloczyn macierzy
przez wektor, co łatwo zrobić kosztem <math>O(N^2)</math> operacji. Dlatego istotne jest,
że algorytm FFT, który za chwilę omówimy, będzie działać kosztem <math>O(N\log\,N)</math>,
czyli praktycznie liniowym (w obu przypadkach stałe proporcjonalności są niewielkie) w dodatku będzie miał znacznie mniejsze wymagania pamięciowe. 

==Algorytm FFT==

Aby móc taniej mnożyć przez <math>F_N</math>, musimy odkryć kilka bardzo szczególnych własności tej macierzy.

{{fakt|||
Macierz <math>F</math> jest symetryczna oraz <math>\frac{1}{N} F_N^* F_N = I</math>.
}}

{{dowod|||
Dowód pozostawiamy jako [[MN10LAB|ćwiczenie]].
}}

Zauważmy, że nasza macierz ma jeszcze więcej specjalnej struktury, powtarza się
w niej bardzo wiele takich samych współczynników (sprawdź dla <math>N=4</math>, w ogólności
<math>F_N</math> ma tylko <math>N</math> różnych wyrazów), gdyż <math>\omega_N^N = 1</math> (dla
<math>j=0,\ldots,N-1</math>, <math>\omega_N^j</math>  to nic innego jak kolejne zespolone pierwiastki z jedynki).

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

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

<center><math>\begin{align} 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}.
\end{align}</math></center>

Suma po indeksach parzystych da się zapisać w postaci

<center><math>
\sum_{k=0}^{m-1} f_{2k} (\omega_m) ^ {jk}
</math></center>

i analogicznie suma po indeksach nieparzystych da się zapisać 

<center><math>
\omega_N^j\cdot \sum_{k=0}^{m-1} f_{2k+1} (\omega_m) ^ {jk}
</math></center>

gdyż <math>\omega_N^2 = \omega_m</math>, a więc wygląda na to, że nasze zadanie wyznaczenia
dyskretnej transformacji Fouriera wymiaru <math>N</math> da się sprowadzić do analogicznych zadań mniejszego
rozmiaru.

Rzeczywiście, korzystając z tego, że <math>j = \beta m + \widetilde{j}</math>, gdzie <math>\widetilde{j}
= 0,\ldots, m-1</math>, oraz <math>\beta = 0</math> lub <math>1</math>, mamy <math>\omega_m^j =
\omega_m^{\widetilde{j}}</math>.

Oznaczając 

<center><math>\begin{align} \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}, 
\end{align}</math></center>

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

<center><math>\begin{align} 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}}.
\end{align}</math></center>

Tym samym, wyznaczenie DFT wymiaru <math>N</math> udało się rzeczywiście sprowadzić
do wyznaczenia dwóch DFT wymiaru <math>m = N/2</math> 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
<math>O(N\log\,N)</math>, a nie <math>O(N^2)</math>, jak dla naiwnego algorytmu mnożenia wektora przez
gęstą macierz. Zysk jest więc, nawet dla niewielkich <math>N</math>, istotny.

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

Jak już zdążyliśmy się przyzwyczaić, gdy tylko to możliwe, w algorytmach numerycznych unikamy stosowania jawnej rekurencji. 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 <math>c = F_Nf</math>.

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 [http://www.xiph.org/ogg  OggVorbis].

==Interpolacja trygonometryczna==

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

Dla danych węzłów <math>x_k = \frac{2\pi k }{N}</math>, <math>k = 0,\ldots,N-1</math>, znaleźć
wielomian <math>P</math> (zmiennej rzeczywistej, o wartościach zespolonych) stopnia <math>\leq N-1</math> postaci
<center><math>
P(x) = \sum_{k=0}^{N-1} \beta_k \bar{\omega}^k</math>,</center>

gdzie <math>\bar{\omega} = e^{ix} = \cos(x) + i\, \sin(x)</math> (stąd nazwa: trygonometryczny)
taki, że 
<center><math>
P(x_k) = y_k, \qquad k = 0,\ldots,N-1</math>,</center>

dla zadanych wartości <math>y_k</math>. 

W języku macierzy DFT możemy zapisać zadanie interpolacji trygonometrycznej jako
zadanie wyznaczania wektora <math>\beta</math> takiego, że 
<center><math>
\overline{F_N}\beta  = y
</math></center>

dla zadanego wektora <math>y</math>.

{{twierdzenie|O współczynnikach wielomianu interpolacji trygonometrycznej|O współczynnikach wielomianu interpolacji trygonometrycznej|

Współczynniki <math>\beta = [\beta_0,\ldots,\beta_{N-1}]^T</math> poszukiwanego wielomianu
trygonometrycznego wyrażają się wzorem
<center><math>
\beta = \frac{1}{N} F_N y</math></center>

}}

{{dowod|||
Jak wiemy, <math>\frac{1}{N} F_N^* F_N = I</math>, z symetrii macierzy <math>F_N</math> wynika więc, 
że <math>(\overline{F_N})^{-1} = \frac{1}{N} F_N</math>.
}}

Można pokazać, że  gdy dane są rzeczywiste, zadanie interpolacji
trygonometrycznej możemy wyrazić
korzystając wyłącznie z liczb rzeczywistych.
Jeśli <math>y_k\in R</math>, <math>k = 0,\ldots,N-1</math>, to wtedy (rzeczywisty) wielomian trygonometryczny
<center><math>
p(x) = \sum_{k=0}^{N-1} \left( a_k \cos(\frac{2k\pi}{N} x) - b_k \sin(\frac{2k\pi}{N} x)\right)</math>,</center>

gdzie <math>a_k = Re(\beta_k)</math>, <math>b_k = Im(\beta_k)</math>, interpoluje <math>y_k</math> w węzłach
<math>x_k</math>. 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 <strong>splotu</strong> dwóch
wektorów, to znaczy wyznaczenie wyrażeń postaci
<center><math>
z_k = \sum_{j=0}^{N-1} u_j v_{k-j}, \qquad k = 0, \ldots, N-1</math></center>
(przyjmujemy tu konwencję, że wektory rozszerzamy <math>N</math>-periodyczne tak, że z definicji <math>v_{j+kN} = v_j</math> dla dowolnych całkowitych <math>j, k</math>).

Zapisując to macierzowo, szukamy iloczynu wektora <math>u</math> z cykliczną macierzą Toeplitza
(macierz Toeplitza ma stałe wyrazy wzdłuż diagonali) wyznaczoną przez <math>v</math>,
<center><math>
z = 
\begin{pmatrix} 
v_0 & v_1 & \cdots      & v_{N-1}  \\
v_{N-1} & v_0 & \cdots & v_{N-2}\\
 \vdots   & \ddots & \ddots & \vdots \\
 v_1    &  \cdots    &  v_{N-1} & v_0
\end{pmatrix}  \cdot \begin{pmatrix}  u_0\\u_1\\\vdots\\u_{N-1}\end{pmatrix} . 
</math></center>

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

<center><math>
\hat{z} = N \cdot \begin{pmatrix} \hat{v}_0 & & \\
 & \ddots & \\
 & & \hat{v}_{N-1}
\end{pmatrix}  \cdot \hat{u} = N\cdot \begin{pmatrix} \hat{v}_0 \hat{u}_0\\
 \vdots  \\
 \hat{v}_{N-1}\hat{u}_{N-1}
\end{pmatrix} </math>,</center>

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

==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 biblioteka o niezbyt skromnie brzmiącej nazwie [http://www.fftw.org  FFTW] (''The Fastest Fourier Transform in the West''). Z tej biblioteki korzystają m.in. funkcje MATLABa i Octave'a <code style="color: #006">fft</code> oraz <code style="color: #006">ifft</code> dla transformacji DFT (to znaczy mnożenia przez <math>F_N</math>) i, odpowiednio, transformacji odwrotnej, <math>\frac{1}{N}\overline{F}_N</math>.

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.

 <Source>/* 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);
}
</Source>
 
Zwrócmy uwagę na linię <code>Plan=fftw_plan_dft_1d(...)</code>. 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, ale na
różnych danych <math>f</math>, należy przed pierwszą DFT taki plan zachować, a potem, aplikując DFT dla następnych danych, wykorzystać gotowy.

==Literatura==

W celu dogłębnego zapoznania się z omawianym na wykładzie materiałem, przeczytaj <b>rozdział 6.12--6.13</b> w
* D. Kincaid, W. Cheney <cite>Analiza numeryczna</cite>, Wydawnictwa Naukowo-Techniczne, Warszawa 2006, ISBN 83-204-3078-X.