From f7c2a9b76085fe8df71e19e19cb657cf055de3e6 Mon Sep 17 00:00:00 2001 From: Martin Mares Date: Tue, 1 Dec 2009 22:07:33 +0100 Subject: [PATCH] Korektury. --- 9-fft/9-fft.tex | 310 ++++++++++++++++++++++++++++++++---------------- 1 file changed, 210 insertions(+), 100 deletions(-) diff --git a/9-fft/9-fft.tex b/9-fft/9-fft.tex index c2249ce..53549b3 100644 --- a/9-fft/9-fft.tex +++ b/9-fft/9-fft.tex @@ -1,6 +1,7 @@ +\def\scharfs{\char"19} \input lecnotes.tex \def\imply{\Rightarrow} -\prednaska{8}{Fourierova transformace}{ \vbox{\hbox{(K.Jakubec, M.Polák +\prednaska{9}{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ý @@ -15,82 +16,73 @@ m Fourier Transform}. -\ss{Trochu algebry na~zaèátek:} -\>Libovolný polynom $P$ stupnì $n$ lze reprezentovat dvìma rùznými zpùsoby: +\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}$). -\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 +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: -\>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 +\>{\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ì $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$. +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 -\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 +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ý 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}$. +\: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 - (viz vý¹e). +\: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. 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í). - +\>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. 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. +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$ øá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}$. +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_{n-2}x^{n-2}) + (p_{1}x^{1} + - p_{3}x^{3} + \ldots + p_{n-1}x^{n-1})$$ +$$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(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}$$ +$$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}$$ -\>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}$). +\>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 + @@ -99,38 +91,143 @@ $3 + 4x + 6x^{2} + 2x^{3} + x^{4} + 10x^{5} = (3 + 6x^{2} + x^{4}) + x(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{\bf Komplexní intermezzo} \itemize\ibull -\:Pro $k\neq j,\ 0\leq k it +{\parskip12pt +\hfill{\bf Základní operace} +\: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$. +\vglue-8pt +\qquad 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$. +\vglue-8pt +\qquad $\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}$. + +\vglue-8pt +\qquad Také $\vert \alpha x \vert = \vert \alpha\vert \cdot \vert x \vert$. + +\:Dìlení: $x/y = (x\cdot \overline{y}) / (y \cdot \overline{y})$. + +\hfill{\bf Gau{\scharfs}ova rovina a goniometrický tvar} + +\: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 ({\it komplexní +jednotky\/}). + +\vglue-8pt +\qquad 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))$. + +\vglue-8pt +\qquad Èíslu $\varphi(x)\in\left[ 0,2\pi \right)$ øíkáme {\it argument\/} +èísla~$x$, nìkdy znaèíme $\mathop{\rm arg} x$. + +\:Navíc $\varphi({\overline{x}}) = -\varphi(x)$. + +\hfill{\bf Komplexní èísla: Exponenciální tvar} + +\: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))}$. + +\vglue-8pt +\qquad (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}$. + +\vglue-8pt +\qquad pozor -- odmocnina není jednoznaèná: $1^4=(-1)^4=\i^4=(-\i)^4=1$. + +\hfill{\bf Komplexní èísla: Odmocniny z~jednièky} + +\:Je-li nìjaké $x\in{\bb C}$ $n$-tou odmocninou z~jednièky, musí platit: +\vglue-8pt +\qquad $\vert x \vert = 1$, tak¾e $x=e^{\i\varphi}$ pro nìjaké~$\varphi$, +\vglue-8pt +\qquad $e^{\i\varphi n} = \cos{\varphi n} + \i\sin\varphi n = 1$, proèe¾ +$\varphi n = 2k\pi$ pro nìjaké $k\in{\bb Z}$. + +\:Z~toho plyne: $\varphi = 2k\pi/n$ +\vglue-8pt +\qquad (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}$. + +\vglue0pt minus4pt % hazel underfull vboxy +\: $x$ je $k$-tá {\it primitivní} odmocnina z 1 $\equiv (\forall j{\bf Konec intermezza} +% end komplex +\bigskip +Vra»me se nyní k algoritmu. Z poslední èásti komplexního intermezza se zdá, +¾e by nemusel být ¹patný nápad zkusit vyhodnocovat polynom v mocninách +$n$-té primitivní odmocniny z~jedné (tedy za $x_0,x_1,\ldots,x_{n-1}$ +z~pùvodního algoritmu zvolíme $\omega^0,\omega^1,\ldots,\omega^{n-1}$). +Aby nám v¹e vycházelo pìknì, zvolíme $n$ jako mocninu dvojky. \ss{Celý algoritmus bude vypadat takto:} \>FFT($P$, $ \omega$) -\>{\sl Vstup:} $p_{0}, \ldots , p_{n-1}$, koeficienty polynomu $P$, a~$\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 , @@ -139,9 +236,12 @@ 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. +\:Jinak rozdìlíme $P$ è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})$. @@ -166,27 +266,31 @@ K~tomu n \>{\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}$ +(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. -\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í +\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$ +\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)}$$ +\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é 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). +\: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-k