From: Martin Mares Date: Wed, 30 Nov 2011 18:40:28 +0000 (+0100) Subject: FFT: Vyhrabana z archivu X-Git-Url: http://mj.ucw.cz/gitweb/?a=commitdiff_plain;h=436ad78e04e86cb79b689538098dded24c46ec5f;p=ads2.git FFT: Vyhrabana z archivu --- diff --git a/7-fft/7-fft.tex b/7-fft/7-fft.tex new file mode 100644 index 0000000..9f7f51c --- /dev/null +++ b/7-fft/7-fft.tex @@ -0,0 +1,391 @@ +\def\scharfs{\char"19} +\input lecnotes.tex +\def\imply{\Rightarrow} +\prednaska{9}{Fourierova transformace}{} + +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} +Celé polynomy oznaèujeme velkými písmeny, jednotlivé èleny polynomù pøíslu¹nými +malými písmeny (pø.: polynom $W$ stupnì $d$ má koeficienty $w_{0}, w_{1}, +w_{2},\ldots, w_{d}$). + +Libovolný polynom $P$ stupnì (nejvý¹e) $d$ lze reprezentovat +jednak jeho koeficienty, tedy èísly $p_{0}, p_{1}, \ldots ,p_{d}$, druhak +i~pomocí hodnot: + +\>{\bf Lemma:} Polynom stupnì nejvý¹e $d$ je jednoznaènì urèeò svými +hodnotami v~$d+1$ rùzných bodech. + +\>{\it Dùkaz:} +Polynom stupnì $d$ má maximálnì $d$ 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ì $d$ nabývající v~daných bodech stejných +hodnot, tak $P-Q$ je polynom stupnì maximálnì $d$, ka¾dé +z $x_{0}\ldots x_{d}$ je koøenem tohoto polynomu $\imply$ spor, polynom stupnì +$d$ má $d+1$ koøenù $\imply$ $P-Q$ musí být nulový polynom $\imply$ $P=Q$. +\qed +\medskip + +Pov¹imnìme si jedné skuteènosti -- máme-li dva polynomy $A$ a~$B$ stupnì $d$ +a~body $x_{0}, \ldots, x_{n}$, dále polynom $C=A \cdot B$ (stupnì $2d$), pak +platí $C(x_{j}) = A(x_{j}) \cdot B(x_{j})$ pro $j = 0,1,2, \ldots, n$. Toto +èiní tento druhý zpùsob reprezentace polynomu velice atraktivním pro násobení +-- máme-li $A$ i $B$ reprezentované hodnotami v $n \geq 2d+1$ bodech, pak +snadno (v $\Theta(n)$) spoèteme takovou reprezentaci $C$. +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ý algoritmus (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. + +\s{Idea, jak by mìl algoritmus pracovat:} +\algo +\:Vybereme $n\geq 2d+1$ bodù $x_{0}, x_{1}, \ldots , x_{n-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: + $C(x_i) = A(x_i)\cdot B(x_i)$ +\:Pøevedeme hodnoty polynomu $C$ na~jeho koeficienty. +\endalgo + +\>Je vidìt, ¾e klíèové jsou kroky 2 a~4. +Celý trik spoèívá v~chytrém vybrání onìch bodù, ve kterých budeme polynomy +vyhodnocovat -- zvolí-li se obecná $x_j$, tak se to rychle neumí, pro speciální +$x_j$ ale uká¾eme, ¾e to rychle jde. + +\ss{Vyhodnocení polynomu metodou Rozdìl a~panuj (algoritmus FFT):} +Mìjme polynom $P$ stupnì $\leq d$ 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_{d-2}x^{d-2}) + (p_{1}x^{1} + + p_{3}x^{3} + \ldots + p_{d-1}x^{d-1})$$ +\>se zavedením znaèení: +$$P_s(t) = p_0t^0 + p_{2}t^{1} + \ldots + p_{d - 2}t^{d - 2\over 2}$$ +$$P_l(t) = p_1t^0 + p_3t^1 + \ldots + p_{d - 1}t^{d - 2\over 2}$$ + +\>bude $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í polynomu $P$ v~$n$ bodech se nám +smrskne na~vyhodnocení $P_s$ a~$P_l$ v~$n/2$ bodech -- oba jsou polynomy +stupnì nejvý¹e $d/2$ a~vyhodnocujeme je v~$x^{2}$ (vyu¾íváme +rovnosti $(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á mocnina 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. + +% komplex +\h{Komplexní intermezzo} +\def\i{{\rm i}} +\def\\{\hfil\break} + +\s{Základní operace} + +\itemize\ibull +\:Definice: ${\bb C} = \{a + b\i \mid a,b \in {\bb R}\}$ + +\:Sèítání: $(a+b\i)\pm(p+q\i) = (a\pm p) + (b\pm q)\i$. \\ +Pro $\alpha\in{\bb R}$ je $\alpha(a+b\i) = \alpha a + \alpha b\i$. + +\:Komplexní sdru¾ení: $\overline{a+b\i} = a-b\i$. \\ +$\overline{\overline x} = x$, $\overline{x\pm y} = \overline{x} \pm +\overline{y}$, $\overline{x\cdot y} = \overline x \cdot \overline y$, $x +\cdot \overline x \in {\bb R}$. + +\:Absolutní hodnota: $\vert x \vert = \sqrt{x\cdot\overline{x}}$, tak¾e $\vert +a+b\i \vert = \sqrt{a^2+b^2}$. \\ +Také $\vert \alpha x \vert = \vert \alpha\vert \cdot \vert x \vert$. + +\:Dìlení: $x/y = (x\cdot \overline{y}) / (y \cdot \overline{y})$. +\endlist + +\s{Gau{\scharfs}ova rovina a goniometrický tvar} + +\itemize\ibull +\:Komplexním èíslùm pøiøadíme body v~${\bb R}^2$: $a+b\i \leftrightarrow (a,b)$. + +\:$\vert x\vert$ je vzdálenost od~bodu $(0,0)$. + +\:$\vert x\vert = 1$ pro èísla le¾ící na~jednotkové kru¾nici ({\I komplexní jednotky}). \\ +Pak platí $x=\cos\varphi + \i\sin\varphi$ pro nìjaké $\varphi\in\left[ +0,2\pi \right)$. + +\:Pro libovolné $x\in{\bb C}$: $x=\vert x \vert \cdot (\cos\varphi(x) + +\i\sin\varphi(x))$. \\ +Èíslu $\varphi(x)\in\left[ 0,2\pi \right)$ øíkáme {\I argument} +èísla~$x$, nìkdy znaèíme $\mathop{\rm arg} x$. + +\:Navíc $\varphi({\overline{x}}) = -\varphi(x)$. +\endlist + +\s{Exponenciální tvar} + +\itemize\ibull +\:Eulerova formule: $e^{i\varphi} = \cos\varphi + \i\sin\varphi$. + +\:Ka¾dé $x\in{\bb C}$ lze tedy zapsat jako $\vert x\vert \cdot e^{\i\cdot +\varphi(x)}$. + +\:Násobení: $xy = \left(\vert x\vert\cdot e^{\i\cdot\varphi(x)}\right) \cdot + \left(\vert y\vert\cdot e^{\i\cdot\varphi(y)}\right) = \vert x\vert + \cdot \vert y\vert \cdot e^{\i\cdot(\varphi(x) + \varphi(y))}$. \\ +(absolutní hodnoty se násobí, argumenty sèítají) + +\:Umocòování: $x^\alpha = \left(\vert x\vert\cdot e^{\i\cdot\varphi(x)}\right)^ + \alpha = {\vert x\vert}^\alpha\cdot e^{\i \alpha \varphi(x)}$. + +\:Odmocòování: $\root n\of x = {\vert x\vert}^{1/n} \cdot e^{\i\cdot +\varphi(x)/n}$. \\ +Pozor -- odmocnina není jednoznaèná: $1^4=(-1)^4=\i^4=(-\i)^4=1$. +\endlist + +\s{Odmocniny z~jednièky} + +\itemize\ibull +\:Je-li nìjaké $x\in{\bb C}$ $n$-tou odmocninou z~jednièky, musí platit: +$\vert x \vert = 1$, tak¾e $x=e^{\i\varphi}$ pro nìjaké~$\varphi$. +Proto $x^n = e^{\i\varphi n} = \cos{\varphi n} + \i\sin\varphi n = 1$. +Platí tedy $\varphi n = 2k\pi$ pro nìjaké $k\in{\bb Z}$. + +\:Z~toho plyne: $\varphi = 2k\pi/n$ \\ +(pro $k=0,\ldots,n-1$ dostáváme rùzné $n$-té odmocniny). + +\:Obecné odmocòování: $\root n \of x = {\vert x\vert}^{1/n} \cdot e^{\i\varphi + (x)/n} \cdot u$, kde $u=\root n\of 1$. + +\:Je-li $x$ odmocninou z 1, pak $\overline{x} = x^{-1}$ -- je toti¾ $1 = \vert +x\cdot \overline{x}\vert = x\cdot \overline{x}$. +\endlist + +\s{Primitivní odmocniny} + +\s{Definice:} $x$ je {\I primitivní} $k$-tá odmocnina z 1 $\equiv x^k=1 \land \forall j: 0Tuto definici splòují napøíklad èísla $\omega = e^{2\pi \i / k}$ a $\overline\omega = e^{-2\pi\i/k}$. +Platí toti¾, ¾e $\omega^j = e^{2\pi\i j/k}$, co¾ je rovno~1 právì tehdy, +je-li $j$ násobkem~$k$ (jednotlivé mocniny èísla~$\omega$ postupnì obíhají +jednotkovou kru¾nici). Analogicky pro~$\overline\omega$. + +\>Uka¾me si nìkolik pozorování fungujících pro libovolné èíslo~$\omega$, +které je primitivní $k$-tou odmocninou z~jednièky (nìkdy budeme potøebovat, +aby navíc $k$ bylo sudé): + +\itemize\ibull +\:Pro $0\leq jFFT($P$, $ \omega$) + +\>{\sl Vstup:} $p_{0}, \ldots , p_{n-1}$, koeficienty polynomu $P$ stupnì +nejvý¹e $n-1$, 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 èleny se sudými a lichými exponenty (jako v pùvodní +my¹lence) 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$, $\omega^2$ je +$n/2$-tá primitivní odmocnina, a mocniny $\omega^2$ jsou stále po dvou +spárované ($n$ je mocnina dvojky, a tedy i $n/2$ je sudé; popø. $n=2$ a je to +zøejmé). +\: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 Diskrétní 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}$$ +(DFT si lze mimo jiné pøedstavit jako funkci vyhodnocující polynom +s~koeficienty $x_k$ v~bodech $\omega^j$). Takovéto zobrazení je lineární +a~tedy popsatelné maticí $\Omega$ s~prvky $\Omega_{jk} = \omega^{jk}$. +Chceme-li umìt pøevádìt z~hodnot polynomu na koeficienty, zajímá nás inverze +této matice. + + +\ss{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¾ijme 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)}\hbox{.}$$ +\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é øady +s kvocientem $\omega ^{(j-k) }$ a~dostaneme ${{\omega^{(j-k)n} -1} \over +{\omega^{(j-k)} -1}} ={1-1 \over \neq 0} = 0$ ( +$\omega^{j-k} - 1$ je jistì $\neq 0$, nebo» $\omega$ je $n$-tá primitivní +odmoncina a $j-kNa¹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}$). + +Vyhodnocení polynomu lze provést vynásobením $\Omega$, pøevod do pùvodní +reprezentace vynásobením $\Omega^{-1}$. My jsme si ale v¹imli chytrého +spárování, a vyhodnocujeme polynom rychleji ne¾ kvadraticky (proto FFT, +jako¾e {\it fast}, ne jako {\it fuj}). Uvìdomíme-li si, ¾e $\overline \omega = +\omega^{-1}$ je také $n$-tá primitivní odmocnina z 1 (má akorát +úhel s opaèným znaménkem), tak mù¾eme stejným trikem vyhodnotit i~zpìtný +pøevod -- nejprve vyhodnotíme $A$ a $B$ v $\omega^j$, poté pronásobíme +hodnoty a~dostaneme tak hodnoty polynomu $C=A\cdot B$, a pustíme na nì +stejný algoritmus s~$\omega^{-1}$ (hodnoty $C$ +vlastnì budou v~algoritmu \uv{koeficienty polynomu}). Nakonec jen získané +hodnoty vydìlíme $n$ a~máme chtìné koeficienty. + +\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í FFT:} + +\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{Paralelní implementace FFT} + +\figure{img.eps}{Pøíklad prùbìhu algoritmu na~vstupu velikosti 8}{3in} + +Zkusme si prùbìh algoritmu FFT znázornit graficky (podobnì, jako jsme kreslili +hradlové sítì). Na~levé stranì obrázku se nachází vstupní vektor $x_0,\ldots,x_{n-1}$ +(v~nìjakém poøadí), na~pravé stranì pak výstupní vektor $y_0,\ldots,y_{n-1}$. +Sledujme chod algoritmu pozpátku: Výstup spoèítáme z~výsledkù \uv{polovièních} +transformací vektorù $x_0,x_2,\ldots,x_{n-2}$ a $x_1,x_3,\ldots,x_{n-1}$. +Èerné krou¾ky pøitom odpovídají výpoètu lineární kombinace $a+\omega^kb$, +kde $a,b$ jsou vstupy krou¾ku a $k$~nìjaké pøirozené èíslo závislé na poloze +krou¾ku. Ka¾dá z~polovièních transformací se poèítá analogicky z~výsledkù +transformace velikosti $n/4$ atd. Celkovì výpoèet probíhá v~$\log_2 n$ vrstvách +po~$\Theta(n)$ operacích. + +Jeliko¾ operace na~ka¾dé vrstvì probíhají na sobì nezávisle, mù¾eme je poèítat +paralelnì. Ná¹ diagram tedy popisuje hradlovou sí» pro paralelní výpoèet FFT +v~èase $\Theta(\log n\cdot T)$ a prostoru $\O(n\cdot S)$, kde $T$ a~$S$ znaèí +èasovou a prostorovou slo¾itost výpoètu lineární kombinace dvou komplexních èísel. + +\s{Cvièení:} Doka¾te, ¾e permutace vektoru $x_0,\ldots,x_{n-1}$ odpovídá bitovému +zrcadlení, tedy ¾e na pozici~$b$ shora se vyskytuje prvek $x_d$, kde~$d$ je +èíslo~$b$ zapsané ve~dvojkové soustavì pozpátku. + +\s{Nerekurzivní FFT} + +Obvod z~pøedchozího obrázku také mù¾eme vyhodnocovat po~hladinách zleva doprava, +èím¾ získáme elegantní nerekurzivní algoritmus pro výpoèet FFT v~èase $\Theta(n\log n)$ +a prostoru $\Theta(n)$: + +\algo +\algin $x_0,\ldots,x_{n-1}$ +\:Pro $j=0,\ldots,n-1$ polo¾íme $y_k\= x_{r(k)}$, kde $r$ je funkce bitového zrcadlení. +\:Pøedpoèítáme tabulku $\omega^0,\omega^1,\ldots,\omega^{n-1}$. +\:$b\= 1$ \cmt{velikost bloku} +\:Dokud $b{\bf Lemma:} Polynom stupnì nejvý¹e $d$ je jednoznaènì urèeò svými -hodnotami v~$d+1$ rùzných bodech. - -\>{\it Dùkaz:} -Polynom stupnì $d$ má maximálnì $d$ 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ì $d$ nabývající v~daných bodech stejných -hodnot, tak $P-Q$ je polynom stupnì maximálnì $d$, ka¾dé -z $x_{0}\ldots x_{d}$ je koøenem tohoto polynomu $\imply$ spor, polynom stupnì -$d$ má $d+1$ koøenù $\imply$ $P-Q$ musí být nulový polynom $\imply$ $P=Q$. -\qed -\medskip - -Pov¹imnìme si jedné skuteènosti -- máme-li dva polynomy $A$ a~$B$ stupnì $d$ -a~body $x_{0}, \ldots, x_{n}$, dále polynom $C=A \cdot B$ (stupnì $2d$), pak -platí $C(x_{j}) = A(x_{j}) \cdot B(x_{j})$ pro $j = 0,1,2, \ldots, n$. Toto -èiní tento druhý zpùsob reprezentace polynomu velice atraktivním pro násobení --- máme-li $A$ i $B$ reprezentované hodnotami v $n \geq 2d+1$ bodech, pak -snadno (v $\Theta(n)$) spoèteme takovou reprezentaci $C$. -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ý algoritmus (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. - -\s{Idea, jak by mìl algoritmus pracovat:} -\algo -\:Vybereme $n\geq 2d+1$ bodù $x_{0}, x_{1}, \ldots , x_{n-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: - $C(x_i) = A(x_i)\cdot B(x_i)$ -\:Pøevedeme hodnoty polynomu $C$ na~jeho koeficienty. -\endalgo - -\>Je vidìt, ¾e klíèové jsou kroky 2 a~4. -Celý trik spoèívá v~chytrém vybrání onìch bodù, ve kterých budeme polynomy -vyhodnocovat -- zvolí-li se obecná $x_j$, tak se to rychle neumí, pro speciální -$x_j$ ale uká¾eme, ¾e to rychle jde. - -\ss{Vyhodnocení polynomu metodou Rozdìl a~panuj (algoritmus FFT):} -Mìjme polynom $P$ stupnì $\leq d$ 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_{d-2}x^{d-2}) + (p_{1}x^{1} + - p_{3}x^{3} + \ldots + p_{d-1}x^{d-1})$$ -\>se zavedením znaèení: -$$P_s(t) = p_0t^0 + p_{2}t^{1} + \ldots + p_{d - 2}t^{d - 2\over 2}$$ -$$P_l(t) = p_1t^0 + p_3t^1 + \ldots + p_{d - 1}t^{d - 2\over 2}$$ - -\>bude $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í polynomu $P$ v~$n$ bodech se nám -smrskne na~vyhodnocení $P_s$ a~$P_l$ v~$n/2$ bodech -- oba jsou polynomy -stupnì nejvý¹e $d/2$ a~vyhodnocujeme je v~$x^{2}$ (vyu¾íváme -rovnosti $(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á mocnina 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. - -% komplex -\h{Komplexní intermezzo} -\def\i{{\rm i}} -\def\\{\hfil\break} - -\s{Základní operace} - -\itemize\ibull -\:Definice: ${\bb C} = \{a + b\i \mid a,b \in {\bb R}\}$ - -\:Sèítání: $(a+b\i)\pm(p+q\i) = (a\pm p) + (b\pm q)\i$. \\ -Pro $\alpha\in{\bb R}$ je $\alpha(a+b\i) = \alpha a + \alpha b\i$. - -\:Komplexní sdru¾ení: $\overline{a+b\i} = a-b\i$. \\ -$\overline{\overline x} = x$, $\overline{x\pm y} = \overline{x} \pm -\overline{y}$, $\overline{x\cdot y} = \overline x \cdot \overline y$, $x -\cdot \overline x \in {\bb R}$. - -\:Absolutní hodnota: $\vert x \vert = \sqrt{x\cdot\overline{x}}$, tak¾e $\vert -a+b\i \vert = \sqrt{a^2+b^2}$. \\ -Také $\vert \alpha x \vert = \vert \alpha\vert \cdot \vert x \vert$. - -\:Dìlení: $x/y = (x\cdot \overline{y}) / (y \cdot \overline{y})$. -\endlist - -\s{Gau{\scharfs}ova rovina a goniometrický tvar} - -\itemize\ibull -\:Komplexním èíslùm pøiøadíme body v~${\bb R}^2$: $a+b\i \leftrightarrow (a,b)$. - -\:$\vert x\vert$ je vzdálenost od~bodu $(0,0)$. - -\:$\vert x\vert = 1$ pro èísla le¾ící na~jednotkové kru¾nici ({\I komplexní jednotky}). \\ -Pak platí $x=\cos\varphi + \i\sin\varphi$ pro nìjaké $\varphi\in\left[ -0,2\pi \right)$. - -\:Pro libovolné $x\in{\bb C}$: $x=\vert x \vert \cdot (\cos\varphi(x) + -\i\sin\varphi(x))$. \\ -Èíslu $\varphi(x)\in\left[ 0,2\pi \right)$ øíkáme {\I argument} -èísla~$x$, nìkdy znaèíme $\mathop{\rm arg} x$. - -\:Navíc $\varphi({\overline{x}}) = -\varphi(x)$. -\endlist - -\s{Exponenciální tvar} - -\itemize\ibull -\:Eulerova formule: $e^{i\varphi} = \cos\varphi + \i\sin\varphi$. - -\:Ka¾dé $x\in{\bb C}$ lze tedy zapsat jako $\vert x\vert \cdot e^{\i\cdot -\varphi(x)}$. - -\:Násobení: $xy = \left(\vert x\vert\cdot e^{\i\cdot\varphi(x)}\right) \cdot - \left(\vert y\vert\cdot e^{\i\cdot\varphi(y)}\right) = \vert x\vert - \cdot \vert y\vert \cdot e^{\i\cdot(\varphi(x) + \varphi(y))}$. \\ -(absolutní hodnoty se násobí, argumenty sèítají) - -\:Umocòování: $x^\alpha = \left(\vert x\vert\cdot e^{\i\cdot\varphi(x)}\right)^ - \alpha = {\vert x\vert}^\alpha\cdot e^{\i \alpha \varphi(x)}$. - -\:Odmocòování: $\root n\of x = {\vert x\vert}^{1/n} \cdot e^{\i\cdot -\varphi(x)/n}$. \\ -Pozor -- odmocnina není jednoznaèná: $1^4=(-1)^4=\i^4=(-\i)^4=1$. -\endlist - -\s{Odmocniny z~jednièky} - -\itemize\ibull -\:Je-li nìjaké $x\in{\bb C}$ $n$-tou odmocninou z~jednièky, musí platit: -$\vert x \vert = 1$, tak¾e $x=e^{\i\varphi}$ pro nìjaké~$\varphi$. -Proto $x^n = e^{\i\varphi n} = \cos{\varphi n} + \i\sin\varphi n = 1$. -Platí tedy $\varphi n = 2k\pi$ pro nìjaké $k\in{\bb Z}$. - -\:Z~toho plyne: $\varphi = 2k\pi/n$ \\ -(pro $k=0,\ldots,n-1$ dostáváme rùzné $n$-té odmocniny). - -\:Obecné odmocòování: $\root n \of x = {\vert x\vert}^{1/n} \cdot e^{\i\varphi - (x)/n} \cdot u$, kde $u=\root n\of 1$. - -\:Je-li $x$ odmocninou z 1, pak $\overline{x} = x^{-1}$ -- je toti¾ $1 = \vert -x\cdot \overline{x}\vert = x\cdot \overline{x}$. -\endlist - -\s{Primitivní odmocniny} - -\s{Definice:} $x$ je {\I primitivní} $k$-tá odmocnina z 1 $\equiv x^k=1 \land \forall j: 0Tuto definici splòují napøíklad èísla $\omega = e^{2\pi \i / k}$ a $\overline\omega = e^{-2\pi\i/k}$. -Platí toti¾, ¾e $\omega^j = e^{2\pi\i j/k}$, co¾ je rovno~1 právì tehdy, -je-li $j$ násobkem~$k$ (jednotlivé mocniny èísla~$\omega$ postupnì obíhají -jednotkovou kru¾nici). Analogicky pro~$\overline\omega$. - -\>Uka¾me si nìkolik pozorování fungujících pro libovolné èíslo~$\omega$, -které je primitivní $k$-tou odmocninou z~jednièky (nìkdy budeme potøebovat, -aby navíc $k$ bylo sudé): - -\itemize\ibull -\:Pro $0\leq jFFT($P$, $ \omega$) - -\>{\sl Vstup:} $p_{0}, \ldots , p_{n-1}$, koeficienty polynomu $P$ stupnì -nejvý¹e $n-1$, 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 èleny se sudými a lichými exponenty (jako v pùvodní -my¹lence) 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$, $\omega^2$ je -$n/2$-tá primitivní odmocnina, a mocniny $\omega^2$ jsou stále po dvou -spárované ($n$ je mocnina dvojky, a tedy i $n/2$ je sudé; popø. $n=2$ a je to -zøejmé). -\: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 Diskrétní 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}$$ -(DFT si lze mimo jiné pøedstavit jako funkci vyhodnocující polynom -s~koeficienty $x_k$ v~bodech $\omega^j$). Takovéto zobrazení je lineární -a~tedy popsatelné maticí $\Omega$ s~prvky $\Omega_{jk} = \omega^{jk}$. -Chceme-li umìt pøevádìt z~hodnot polynomu na koeficienty, zajímá nás inverze -této matice. - - -\ss{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¾ijme 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)}\hbox{.}$$ -\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é øady -s kvocientem $\omega ^{(j-k) }$ a~dostaneme ${{\omega^{(j-k)n} -1} \over -{\omega^{(j-k)} -1}} ={1-1 \over \neq 0} = 0$ ( -$\omega^{j-k} - 1$ je jistì $\neq 0$, nebo» $\omega$ je $n$-tá primitivní -odmoncina a $j-kNa¹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}$). - -Vyhodnocení polynomu lze provést vynásobením $\Omega$, pøevod do pùvodní -reprezentace vynásobením $\Omega^{-1}$. My jsme si ale v¹imli chytrého -spárování, a vyhodnocujeme polynom rychleji ne¾ kvadraticky (proto FFT, -jako¾e {\it fast}, ne jako {\it fuj}). Uvìdomíme-li si, ¾e $\overline \omega = -\omega^{-1}$ je také $n$-tá primitivní odmocnina z 1 (má akorát -úhel s opaèným znaménkem), tak mù¾eme stejným trikem vyhodnotit i~zpìtný -pøevod -- nejprve vyhodnotíme $A$ a $B$ v $\omega^j$, poté pronásobíme -hodnoty a~dostaneme tak hodnoty polynomu $C=A\cdot B$, a pustíme na nì -stejný algoritmus s~$\omega^{-1}$ (hodnoty $C$ -vlastnì budou v~algoritmu \uv{koeficienty polynomu}). Nakonec jen získané -hodnoty vydìlíme $n$ a~máme chtìné koeficienty. - -\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í FFT:} - -\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{Paralelní implementace FFT} - -\figure{img.eps}{Pøíklad prùbìhu algoritmu na~vstupu velikosti 8}{3in} - -Zkusme si prùbìh algoritmu FFT znázornit graficky (podobnì, jako jsme kreslili -hradlové sítì). Na~levé stranì obrázku se nachází vstupní vektor $x_0,\ldots,x_{n-1}$ -(v~nìjakém poøadí), na~pravé stranì pak výstupní vektor $y_0,\ldots,y_{n-1}$. -Sledujme chod algoritmu pozpátku: Výstup spoèítáme z~výsledkù \uv{polovièních} -transformací vektorù $x_0,x_2,\ldots,x_{n-2}$ a $x_1,x_3,\ldots,x_{n-1}$. -Èerné krou¾ky pøitom odpovídají výpoètu lineární kombinace $a+\omega^kb$, -kde $a,b$ jsou vstupy krou¾ku a $k$~nìjaké pøirozené èíslo závislé na poloze -krou¾ku. Ka¾dá z~polovièních transformací se poèítá analogicky z~výsledkù -transformace velikosti $n/4$ atd. Celkovì výpoèet probíhá v~$\log_2 n$ vrstvách -po~$\Theta(n)$ operacích. - -Jeliko¾ operace na~ka¾dé vrstvì probíhají na sobì nezávisle, mù¾eme je poèítat -paralelnì. Ná¹ diagram tedy popisuje hradlovou sí» pro paralelní výpoèet FFT -v~èase $\Theta(\log n\cdot T)$ a prostoru $\O(n\cdot S)$, kde $T$ a~$S$ znaèí -èasovou a prostorovou slo¾itost výpoètu lineární kombinace dvou komplexních èísel. - -\s{Cvièení:} Doka¾te, ¾e permutace vektoru $x_0,\ldots,x_{n-1}$ odpovídá bitovému -zrcadlení, tedy ¾e na pozici~$b$ shora se vyskytuje prvek $x_d$, kde~$d$ je -èíslo~$b$ zapsané ve~dvojkové soustavì pozpátku. - -\s{Nerekurzivní FFT} - -Obvod z~pøedchozího obrázku také mù¾eme vyhodnocovat po~hladinách zleva doprava, -èím¾ získáme elegantní nerekurzivní algoritmus pro výpoèet FFT v~èase $\Theta(n\log n)$ -a prostoru $\Theta(n)$: - -\algo -\algin $x_0,\ldots,x_{n-1}$ -\:Pro $j=0,\ldots,n-1$ polo¾íme $y_k\= x_{r(k)}$, kde $r$ je funkce bitového zrcadlení. -\:Pøedpoèítáme tabulku $\omega^0,\omega^1,\ldots,\omega^{n-1}$. -\:$b\= 1$ \cmt{velikost bloku} -\:Dokud $b