MN10

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)

Plik:XXX.png

Odtwarzacz MP3

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$, wyznacz $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.

Jean Baptiste Joseph Fourier
Zobacz biografię

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$.

Algorytm FFT

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, 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ę po raz 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

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 $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 OggVorbis.

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$, z symetrii macierzy $F_N$ wynika więc, że $(\overline{F_N}{)}^{-1}=\frac{1}{N}{F}_N$.

Można pokazać, że gdy dane są rzeczywiste, zadanie interpolacji trygonometrycznej możemy 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 \PIPEREAD  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.