[ads2.git] / 9-fft / 9-fft.tex

\input lecnotes.tex
\def\imply{\Rightarrow}
\prednaska{8}{Fourierova transformace}{ \vbox{\hbox{(K.Jakubec, M.Polák
        a~G.Ocsovszky,}\hbox{\ V.Tùma, M.Kozák)}}}

Násobení polynomù mù¾e mnohým pøipadat jako pomìrnì (algoritmicky) snadný
problém. Asi ka¾dého hned napadne \uv{hloupý} algoritmus -- vezmeme
koeficienty prvního polynomu a~vynásobíme ka¾dý se v¹emi koeficienty druhého
polynomu a~pøíslu¹nì u~toho seèteme i~exponenty (stejnì jako to dìláme, kdy¾
násobíme polynomy na~papíøe). Pokud stupeò prvního polynomu je $n$ a~druhého
$m$, strávíme tím èas $\Omega(mn)$. Pro $m=n$ je to kvadraticky pomalé.
Na~první pohled se mù¾e zdát, ¾e rychleji to prostì nejde (pøeci musíme
v¾dy vynásobit \uv{ka¾dý s~ka¾dým}). Ve skuteènosti to ale rychleji fungovat
mù¾e, ale k~tomu je potøeba znát trochu tajemný algoritmus FFT neboli {\I Fast
Fourier Transform}.


\ss{Trochu algebry na~zaèátek:}
\>Libovolný polynom $P$ stupnì $n$ lze reprezentovat dvìma rùznými zpùsoby:

\itemize\ibull
\:svými koeficienty, èili èísly $p_{0}, p_{1}, \ldots ,p_{n}$, nebo
\:svými hodnotami v~$n$ rùzných bodech $x_{0}, x_{1}, \ldots , x_{n}$, èili
èísly $P(x_{0}),$ $P(x_{1}),$ $\ldots , P(x_{n})$.
\endlist

\>Dostateènost $n+1$ hodnot pro urèení druhým zpùsobem lze dokázat
následnovnì: polynom stupnì $n$ má maximálnì $n$ koøenù (indukcí, je-li
$k$ koøenem $P$, pak lze $P$ napsat jako $(x-k)Q$ kde $Q$ je polynom stupnì
o~jedna men¹í, pøitom polynom stupnì 1 má jediný koøen); uvá¾íme-li
dva rùzné polynomy $P$ a~$Q$ stupnì $n$ nabývající v~daných bodech stejných
hodnot, tak $P-Q$ je polynom stupnì maximálnì $n$, ka¾dé
z $x_{0}\ldots x_{n}$ je koøenem tohoto polynomu $\imply$ spor, polynom stupnì
$n$ má $n+1$ koøenù $\imply$ $P-Q$ musí být nulový polynom $\imply$ $P=Q$.
\medskip        

\ss{Konvence:}
Celé polynomy oznaèujeme velkými písmeny, jednotlivé èleny polynomù pøíslu¹nými
malými písmeny (pø.: polynom $W$ stupnì $n$ má koeficienty $w_{0}, w_{1},
w_{2},\ldots, w_{n}$).

Pov¹imnìme si jedné skuteènosti -- máme-li dva polynomy $A$ a~$B$ stupnì $n$
a~body $x_{0}, \ldots, x_{k}$, dále polynom $C=A \cdot B$ (stupnì $2n$), pak
platí $C(x_{k}) = A(x_{k}) \cdot B(x_{k}), k = 0,1,2, \ldots, n.$ Toto èiní
tento druhý zpùsob reprezentace polynomu velice atraktivním pro násobení.
Problémem je, ¾e typicky máme polynom zadaný koeficienty a~ne hodnotami
v~bodech. Tím pádem potøebujeme nìjaký hodnì rychlý algorimtus (tj.
rychlej¹í ne¾ kvadratický, jinak bychom si nepomohli oproti hloupému
algoritmu) na~pøevod polynomu z jedné reprezentace do druhé a~zase zpìt.

Dále bychom si mìli uvìdomit, ¾e stupeò na¹eho výsledného polynomu $C$ bude
$\leq 2n$ (kde $n$ je stupeò výchozích polynomù). Pokud chceme polynom $C$
reprezentovat pomocí jeho hodnot v~bodech, musíme tedy vzít alespoò $2n$
bodù. Tímto konèí malá algebraická vsuvka.

\s{Idea, jak by mìl algoritmus pracovat:}
\algo
\:Vybereme $2n$ bodù $x_{0}, x_{1}, \ldots , x_{2n-1}$.
\:V tìchto bodech vyhodnotíme polynomy $A$ a~$B$.
\:Nyní ji¾ v~lineárním èase získáme hodnoty polynomu $C$ v~tìchto bodech
        (viz vý¹e).
\:Pøevedeme hodnoty polynomu $C$ na~jeho koeficienty.
\endalgo

\>Je vidìt, ¾e klíèové jsou kroky 2 a~4. Vybrání bodù jistì stihneme pohodlnì
v~lineárním èase a~vynásobení samotných hodnot té¾ (máme $2n$ bodù a~$C(x_{k})
= A(x_{k}) \cdot B(x_{k}), k = 0,1,2, \ldots , 2n-1$, tak¾e na~to nepotøebujeme
více ne¾ $2n$ násobení).

Celý trik spoèívá v~chytrém vybrání onìch bodù, ve kterých budeme polynomy
vyhodnocovat. Je na~to potøeba vìdìt pár zajímavostí o~komplexních èíslech,
na~webové stránce pøedná¹ky jsou k dispozici slajdy, zde to bude zapsáno
o~trochu struènìji.

\ss{  Vyhodnocení polynomu metodou Rozdìl a~panuj (algoritmus FFT):}
Mìjme polynom $P$ øádu $n$ a~chtìjme jej vyhodnotit v~$n$ bodech. Vybereme si
body tak, aby byly spárované, èili $\pm x_{0}, \pm x_{1}, \ldots , \pm
x_{n/2-1} $. To nám výpoèet urychlí, proto¾e pak se druhé mocniny $x_{j}$
shodují s~druhými mocninami $-x_{j}$.

Polynom $P$ rozlo¾íme na~dvì èásti, první obsahuje èleny se sudými exponenty,
druhá s~lichými:
$$P(x) = (p_{0}x^{0} + p_{2}x^{2} + \ldots + p_{n-2}x^{n-2}) + (p_{1}x^{1} +
        p_{3}x^{3} + \ldots + p_{n-1}x^{n-1})$$
\>se zavedením znaèení:
$$P_s(x^{2}) = p_{0}x^{0} + p_{2}x^{2} + \ldots + p_{n - 2}x^{n - 2}$$
$$P_l(x^{2}) = p_{1}x^{1} + p_{3}x^{3} + \ldots + p_{n - 1}x^{n - 1}$$

\>Dohromady $P(x) = P_s(x^{2}) + xP_l(x^{2})$ a~$P(-x) = P_s(x^{2}) -
xP_l(x^{2})$. Jinak øeèeno, vyhodnocování $P$ v~$n$ bodech se nám smrskne
na~vyhodnocení $P_s(x)$ a~$P_l(x)$ (oba jsou polynomy stupnì $n/2$ 
a~vyhodnocujeme je nyní v~$x^{2}$) v~$n/2$ bodech (proto¾e $(x_{i})^{2} =
(-x_{i})^{2}$).

\s{Pøíklad:}
$3 + 4x + 6x^{2} + 2x^{3} + x^{4} + 10x^{5} = (3 + 6x^{2} + x^{4}) + x(4 +
2x^{2} + 10x^{4})$.

Teï nám ov¹em vyvstane problém s~oním párováním -- druhá mocina pøece nemù¾e
být záporná a~tím pádem u¾ v~druhé úrovni rekurze body spárované nebudou.
Z~tohoto dùvodu musíme pou¾ít komplexní èísla -- tam druhé mocniny záporné býti
mohou. Jako $x_{0}, \ldots , x_{n-1} $ si zvolíme mocniny
$n$-té primitvní odmocniny z jedné (oznaèíme si ji jako $\omega$). Pøipomeòme:
$\omega$ je $n$-tá odmocnina z $1$ je {\it primitivní} $\equiv$ $(\forall k)
(k<n)\ \omega^k \neq 1$ a~$\omega^n = 1$. Máme $n$ navzájem rùzných odmocnin
z~jednièky, rovnomìrnì rozesetých po jednotkové kru¾nici, BÚNO $n=2^{k}, k \in
N$ (jinak viz slajdy Martina Mare¹e). Jednotlivé $x$ vypadají takto: $1,
\omega, \omega^{2}, \ldots , \omega^{n - 1} $, kde $\omega = e^{2 \pi i/ n}$.

\s{Pìt poznámek:}
\itemize\ibull
\:Pro $k\neq j,\ 0\leq k<j<n$ je $\omega^k \neq \omega^j$, nebo» $\omega^j /
\omega^k = \omega^{j-k}\neq 1$, proto¾e $j-k < n$ a~$\omega$ je $n$-tá
primitivní (a~tedy $\omega^0, \ldots, \omega^{n-1}$ jsou rùzné).
\:$\overline{\omega} = \omega^{-1}$, nebo» $\omega^{-1} = \omega^{-2\pi i/n}$
a~tedy èísla $\omega,\omega^{-1}$ svírají vùèi reálné ose opaèný úhel a~jsou
        komplexnì sdru¾ená.
\:mocniny $\omega$ jsou spárované, èili $\omega^{j} = -\omega^{n/2 + j}$
        (staèí si uvìdomit, ¾e $\omega^{n/2} = -1$ pro $n$ sudé)
\:umocníme-li v¹echny $\omega^{0},\ldots, \omega^{n-1}$ na~druhou, vzniknou nám
        odmnocniny z~$1$ které budou i~nadále spárované: $\omega^n$ je toti¾
        $1$ a~exponenty lze tedy poèítat$\mod n$, $n$ je sudé $\imply$
        dostaneme pùvodní posloupnost ze které zmizí liché mocniny $\omega$, ale $n$
        je mocnina dvojky a~tedy i~$n/2$ je sudé (pro $n=2$ ne, ale tam je to
        triviální), tak¾e sudé mocniny jsou spárované navzájem..
\:$\omega^2$ je $n/2$-tá primitivní odmocnina z 1 -- je toti¾ $(\omega^{-2\pi
        i/n})^2 = \omega^{-2\pi i/(n/2)}$ a~$n$ je mocnina dvojky.
\endlist

\ss{Celý algoritmus bude vypadat takto:}
\>FFT($P$, $ \omega$)

\>{\sl Vstup:} $p_{0}, \ldots , p_{n-1}$, koeficienty polynomu $P$, a~$\omega$,
$n$-tá primitivní odmocina z jedné.

\>{\sl Výstup:} Hodnoty polynomu v~bodech $1, \omega, \omega^{2}, \ldots ,
\omega^{n - 1}$, èili èísla $P(1), P(\omega), P(\omega^{2}),$ $\ldots ,
P(\omega^{n - 1})$.

\algo
\:Pokud $n = 1$, vrátíme $p_{0}$ a~skonèíme.
\:Jinak rozdìlíme $P$ na~sudé a~liché koeficienty a~rekurzivnì zavoláme
        FFT($P_s$, $\omega^{2}$) a~FFT($P_l$, $\omega^{2}$) -- $P_l$ i~$P_s$
        jsou stupnì max. $n/2-1$ a~$\omega^2$ je $n/2$-tá primitivní odmocnina.
\:Pro $j = 0, \ldots , n/2 - 1$ spoèítáme:

\:\qquad $P(\omega^{j}) = P_s(\omega^{2j}) + \omega^{j}\cdot P_l(\omega^{2j})$.

\:\qquad $P(\omega^{j+n/2})=P_s(\omega^{2j})-\omega^{j}\cdot P_l(\omega^{2j})$.

\endalgo


\s{Èasová slo¾itost:}
\>$T(n)=2T(n/2) + \Theta(n) \Rightarrow$ slo¾itost $\Theta(n \log n)$, jako
MergeSort.


Máme tedy algoritmus, který pøevede koeficienty polynomu na~hodnoty tohoto
polynomu v~rùzných bodech. Potøebujeme ale také algoritmus, který doká¾e
reprezentaci polynomu pomocí hodnot pøevést zpìt na~koeficienty polynomu.
K~tomu nám pomù¾e podívat se na~ná¹ algoritmus trochu obecnìji.


\s{Definice:}
\>{\I Diskretní Fourierova transformace} $(DFT)$
je zobrazení $f: { {\bb C} ^n} \rightarrow { {\bb C} ^n}$, kde $$y=f(x) \equiv
\forall j \ y_{j} = \sum \limits ^{n-1}_{k=0} x_{k} \cdot \omega ^{jk}$$
(pøestavme si ji jako funkci vyhodnocující polynom s~koeficienty $x$ v~bodech
$\omega^j$). Takovéto zobrazení je lineární a~tedy popsatelné maticí $\Omega$
s~prvky $\Omega_{jk} = \omega^{jk}$


\s{Jak najít inverzní matici?} Znaème $\overline{\Omega}$ matici, její¾ prvky
jsou komplexnì sdru¾ené odpovídajícím prvkùm $\Omega$, a vyu¾ijeme následující
lemma:

\ss{Lemma:}$\quad \Omega\cdot \overline{\Omega} = n\cdot E$

\proof $$ (\Omega\cdot \overline{\Omega})_{jk} = \sum_{l=0}^{n-1} \omega^{jl}
\cdot \overline{\omega^{lk}} = \sum \omega^{jl} \cdot \overline{\omega}^{lk} =
\sum \omega^{jl} \cdot \omega^{-lk} = \sum \omega^{l(j-k)}$$
\itemize\ibull
\:Pokud $j=k$, pak $ \sum \limits ^{n-1}_{l=0} (\omega ^{0}) ^{l} = n$.

\:Pokud $j\neq k$, pou¾ijeme vzoreèek pro souèet geometrické posloupnosti, kde
$a_{1}=1$ a~$q=\omega ^{(j-k) }$ a~dostaneme ${{\omega^{(j-k)n} -1} \over
{\omega^{(j-k)} -1}} ={1-1 \over r- 1} = {0 \over \neq 0} = 0$, kde $r$
je èíslo rùzné od jednièky (nebo» $\omega$ je $n$-tá primitivní odmoncina).
\endlist
\qed


\>Na¹li jsme inverzi:
$\Omega({1 \over n} \overline{\Omega}) = {1 \over n}\Omega \cdot
\overline{\Omega} = E$, \quad $\Omega^{-1}_{jk} = {1 \over n}\overline{\omega^
{jk}} = {1 \over n}\omega^{-jk} = {1 \over n} {(\omega^{-1})}^{jk}$, \quad
(pøipomínáme, $\omega^{-1}$ je $\overline{\omega}$, $\omega$ je $n$-tá
primitivní odmocnina z jednièky).

Ná¹ algoritmus poèítá tedy i~inverzní transformaci, pouze místo $\omega_n$
pou¾ijeme komplexnì zdru¾ené $\overline{\omega_n}$ a~matici vynásobíme $(1/n)$.
Co¾ je skvìlé -- staèí znát pouze jeden algoritmus u~kterého staèí v~jednom
pøípadì pou¾ít transformovanou matici a~vydìlit $n$.

\s{Výsledek:} Pro $n= 2^k$ lze DFT na~${\bb C}^n$ spoèítat v~èase $\Theta(n
\log n)$ a~DFT$^{-1}$ takté¾.

\s{Dùsledek:}
Polynomy stupnì $n$ lze násobit v~èase $\Theta(n \log n)$:
$\Theta(n \log n)$ pro vyhodnocení, $\Theta(n)$ pro vynásobení a~$\Theta(n
\log n)$ pro pøevedení zpìt.

\s{Pou¾ití:}

\itemize\ibull

\:Zpracování signálu -- rozklad na~siny a~cosiny o~rùzných frekvencích
$\Rightarrow$ spektrální rozklad.
\:komprese dat -- napøíklad formát JPEG.
\:Násobení dlouhých èísel v~èase $\Theta(n \log n)$.
\endlist

\s{Hardwarová implementace FFT}

\figure{img.eps}{Pøíklad prùbìhu algoritmu na~vstupu velikosti 8}{3in}


\>Obrázek ukazuje zapojení kombinaèního obvodu pro DFT pro vstup velikosti~8.
Hladin bude v¾dy $\log_2 n$, tj. v~na¹em pøípadì $\log_2 8 = 3$ hladiny.

\>Podívejme se na~pravou èást obrázku, tedy výstup celého obvodu. Èerná koleèka
pøedstavují podobvody, rovnice vedle nich operaci, kterou provádìjí. Hodnoty
$y_j$ znaèí hodnotu polynomu $P$ v~bodì $\omega^j$ kde $\omega^j$ je $j-tá$
mocnina primitivní $n$-té odmocniny z jednièky. K jejímu spoètení ale
potøebujeme znát hodnoty $s_k$ a~$l_k$ kde $k$ je z intervalu $[0, {n/2} -1]$
a~$s_k$ a~$l_k$ jsou hodnoty polynomu stupnì ${n/2}$ v~bodì $\omega^{2k}$.
V polynomu $s$ jsou sudé koeficienty a~v polynomu $l$ liché koeficienty
polynomu $P$. Vidíme, ¾e se jedná pøesnì o~ná¹ rekurzivní algoritmus pro
poèítání FFT a~tímto zpùsobem postavíme celou sí».
\>Tímto obvodem jsme tedy získali nerekurzivní algoritmus pro poèítání FFT.
V¹imìme si poøadí vstupních hodnot (koeficientù). Èísla jsou v~binárním tvaru
0--7 pøeètená pozpátku. Pro pøedstavu jaké koeficienty polynomu $P$ se objevují
v~rùzných hladinách, na~obrázku jsou naznaèena jejich èísla spolu s~pøíslu¹nými
mocninami primitivní $n$-té odmocniny z jednièky.

\s{Z~toho:}

\itemize\ibull
\:Kombinaèní obvod pro DFT
s~$\Theta(\log n)$ hladinami
a $\Theta(n)$ hradly na~hladinì.
\:Nerekurzivní algoritmus (postupujeme zleva) v~èase $\Theta(n \log n)$.

\endlist
\bigskip

\>Nakonec dodejme, ¾e poèítat lze nejen nad tìlesem komplexních èísel, ale
v podstatì nad ka¾dým komutativním tìlesem s~$n$-tou primitivní odmocninou
(pøíkladem jsou tìlesa øádu $p = 2^k + 1$, $p$ prvoèíslo) èi dokonce
komutativním okruhem s~multiplikativními inversemi pro øád vektoru
a~primitivní odmocninu.

\bye