MN10

Szybka transformacja Fouriera (FFT)

<<< Powrót do strony głównej przedmiotu Metody numeryczne

Algorytm FFT (\ang.{Fast Fourier Transform}) dla dyskretnej transformacji Fouriera (DFT, \ang.{Discrete Fourier Transform}) i jego krewniacy dla wyznaczania dyskretnej transformacji cosinusów (DCT, \ang.{Discrete Cosine Transform} i MDCT, \ang.{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)

Plik:XXX.png

Odtwarzacz MP3

a także do

• filtrowania 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!) banalne i jednocześnie 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.

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

Algorytm FFT

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

Fakt

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

Dowód

Dowód pozostawiamy jako ć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 $N=4$, w ogólności $F_N$ 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 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 $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 O współczynnikach wielomianu interpolacji trygonometrycznej

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

Mogłoby się zdawać, że 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 biblioteka o niezbyt skromnie brzmiącej nazwie FFTW (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}{\bar{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, ale na różnych danych $f$, 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 rozdział 6.12--6.13 w

• D. Kincaid, W. Cheney Analiza numeryczna, Wydawnictwa Naukowo-Techniczne, Warszawa 2006, ISBN 83-204-3078-X.