8-fft/8-fft.tex

   1 \input lecnotes.tex
   2 \prednaska{8}{Fourierova transformace}{(K.Jakubec, M.Polák a G.Ocsovszky)}
   3
   4 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 -- jednodu¹e vezmeme koeficienty prvního polynomu a ka¾dým z nich pøenásobím v¹echny koeficienty toho druhého. Pokud øád prvního polynomu je $n$ a druhého $m$, tak èasová slo¾itost nám vyjde nìco jako $\O(mn)$. To není a¾ tak ¹patné, v nejhor¹ím pøípadì se dostaneme na $\O(n^{2})$ (pokud $m = n$). Na první pohled se mù¾e zdát, ¾e rychleji to prostì nejde (pøeci musíme v¾dy vynásobit \uv{ka¾dej 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}.
   5
   6
   7 \ss{Trochu algebry na zaèátek:}
   8 \>Libovolný polynom $P$ øádu $n$ mù¾eme být reprezentován dvìma rùznými zpùsoby:
   9
  10 \itemize\ibull
  11 \:svými koeficienty, èili èísly $a_{0}, a_{1}, \ldots ,a_{n}$, nebo
  12 \:svými hodnotami v $n + 1$ rùzných bodech $x_{0}, x_{1}, \ldots , x_{n}$, èili èísly $P(x_{0}),$ $P(x_{1}),$ $\ldots , P(x_{n})$.
  13 \endlist
  14
  15 \>Pov¹imnìme si jedné skuteènosti -- máme-li dva polynomy $A$ a $B$ øádu $n$ a body $x_{0}, \ldots, x_{k}$, pak platí $C(x_{k}) = A(x_{k}) \cdot B(x_{k}), k = 0,1,2, \ldots, n+1.$ 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.
  16
  17 Dále bychom si mìli uvìdomit, ¾e stupeò na¹eho výsledného polynomu $C$ bude $\leq 2n+1$ (kde $n$ je stupnìm výchozích polynomù). To snad netøeba nijak vysvìtlovat, ka¾dý si to snadno ovìøí, jen dodáme, ¾e pokud chceme polynom $C$ reprezentovat pomocí jeho hodnot v bodech, musíme vzít $2n + 2$ bodù. Tímto konèí malá algebraická vsuvka.
  18
  19 \s{Idea, jak by mìl algoritmus pracovat:}
  20 \algo
  21 \:Vybereme $2n + 2$ bodù $x_{0}, x_{1}, \ldots , x_{2n+1}$.
  22 \:V~tìchto bodech vyhodnotíme polynomy $A$ a $B$.
  23 \:Nyní ji¾ v lineárním èase získáme polynom $C$ (viz vý¹e).
  24 \:Inverznì pøevedeme hodnoty polynomu $C$ v $2n+2$ bodech na jeho koeficienty.
  25 \endalgo
  26
  27 \>Je asi 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+2$ 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+2$ násobení).
  28
  29 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 stránce Matrina Mar¹e jsou k dispozici slajdy, zde to bude zapsáno o~trochu struènìji.
  30
  31 \ss{Vyhodnocení polynomu metodou rozdìl a panuj (algoritmus FFT):}
  32 Mìjme polynom $P$ øádu $n$ a chceme 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_{i}$ shodují s~druhými mocninami $-x_{i}$.
  33
  34 Polynom $P$ rozlo¾íme na dvì èásti, první obsahuje èleny se sudými exponenty, druhá s~lichými:
  35
  36 $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}$
  37
  38 $S(x^{2}) = p_{0}x^{0} + p_{2}x^{2} + \ldots + p_{n - 2}x^{n - 2}$,
  39 $L(x^{2}) = p_{1}x^{1} + p_{3}x^{3} + \ldots + p_{n - 1}x^{n - 1}$
  40
  41 \>Tak¾e obecnì $P(x) = S(x^{2}) + xL(x^{2})$ a $P(-x) = S(x^{2}) - xL(x^{2})$.
  42 Jinak øeèeno, vyhodnocování $P(x)$ v $n$ bodech se nám smrskne na vyhodnocení $S(x)$ a $L(x)$ (oba mají polovièní stupeò ne¾ $P(x)$) v $n/2$ bodech (proto¾e $(x_{i})^{2} = (-x_{i})^{2}$).
  43
  44 \s{Pøíklad:}
  45 $3 + 4x + 6x^{2} + 2x^{3} + x^{4} + 10x^{5} = (3 + 6x^{2} + x^{4}) + x(4 + 2x^{2} + 10x^{4})$.
  46
  47
  48 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 $n$-tou komplexní odmocninu z~jedné. Máme $n$ $n$-tý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é odmociny vypadají takto: $1, \omega, \omega^{2}, \ldots , \omega^{n - 1} $, kde $\omega = e^{2 \pi i/ n}$.
  49
  50 \s{Dvì poznámky:}
  51 \itemize\ibull
  52 \:$n$-té odmocniny z~jednièky jsou spárované, èili $\omega^{j} = -\omega^{n/2 + j}$,
  53 \:umocníme-li v¹echny na druhou, vznikne nám $n/2$ $n/2$-tých odmocnin z~jedné, které jsou i nadále spárované.
  54 \endlist
  55
  56 \ss{Tak a teï koneènì ten slavný algoritmus:}
  57 \>FFT($P$, $ \omega$)
  58
  59 \>{\sl Vstup:} $p_{0}, \ldots , p_{n-1}$, koeficienty polynomu $P$, a $\omega$, $n-$tá odmocina z~jedné.
  60
  61 \>{\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})$.
  62
  63 \algo
  64 \:Pokud $n = 1$, vra» $P_{0}$ a konec.
  65 \:Jinak rozdìl $P$ na sudé a liché koeficienty a zarekurzi se do FFT($S$, $\omega^{2}$) a FFT($L$, $\omega^{2}$).
  66 \:Pro $j = 0, \ldots , n - 1$ spoèítej: $P(\omega^{j}) = S(\omega^{2j}) + \omega^{j} \cdot L(\omega^{2j})$.
  67
  68 \endalgo
  69
  70
  71 \s{Èasová slo¾itost:}
  72 \>$T(n)=2T({n \over 2} ) + \O(n)  \Rightarrow$ slo¾itost  $\O(n \log n)$, stejnì jako MergeSort.
  73
  74
  75
  76
  77 Máme tedy algoritmus, který \uv{pøevede} koeficienty polynomu na hodnoty tohoto polynomu v~námi zadaných bodech. Ale potøebujeme také algoritmus, který doká¾e reprezentaci polynomu pomocí hodnot pøevést zpìt na koeficienty polynomu. Tedy nìjaký inverzní algoitmus.
  78 Definuje me si algoritmus DFT, která vyu¾ívá maticovou reprezentaci a s~jeho¾ pomocí získáme hledaný algoritmus.
  79
  80 \s{Definice:}
  81 \>{\I Diskretní Fourierova transformace} $(DFT)$
  82 je funkce $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} . \omega ^{jk}$.
  83
  84 \s{Poznámka:}
  85 Vezmeme polynom, který má $x_{kj}$ jako koeficienty a vyhodnotíme ho v~bodì $\omega ^{j} [y_{j} = x(\omega^{j})] \Rightarrow  {f}$  je linearní $\Rightarrow$ mù¾eme  napsat $f(x) = \Omega_{x} ,\  \Omega _{jk} =\omega ^{jk}$, kde $\Omega$ je matice.
  86
  87
  88 \s{Jak najít inverzní matici?} Víme, ¾e $\Omega =\Omega ^{T}$ proto¾e $\omega ^{jk} = \omega ^{kj}$.
  89
  90 \ss{Jak vypadají øádky této matice?}
  91 Vyu¾ijeme následující lemma, které si ale napøed doká¾eme :)
  92
  93 \ss{Lemma:}
  94
  95 \proof Souèin
  96 $$\Omega _{j}  \Omega _{k} = \sum \limits ^{n-1}_{l=0} \Omega _{jl}  \overline{\Omega _{kl}} = \sum \limits _{l}  \omega ^{jl}   \overline{\omega ^{kl}}  = \sum \limits _{l} \omega ^{jl}  \omega ^{-kl} = \sum \limits _{l} \omega ^{(j-k)l } = \sum \limits ^{n-1}_{l} (\omega^{j-k}) ^{l}, $$
  97
  98 proto¾e $ \overline{\omega^{kl}} = \overline{\omega} ^{kl} = {({1 \over \omega} )}^{kl} = \omega ^{-kl}$.
  99
 100 \itemize\ibull
 101 \: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 @ -1} = {0 \over \neq 0} = 0$.
 102
 103 \:Pokud $j=k \sum \limits ^{n-1}_{l=0} (\omega ^{0}) ^{l} = n$.
 104 \endlist
 105 \qed
 106
 107 \>A nyní slibované a u¾ i dokázané lemma:
 108
 109 \s{Lemma:} \quad $\Omega _{j}  \cdot \Omega _{k} = \left\{
 110 {\displaystyle 0 \ldots j\neq k}\atop
 111 {\displaystyle 1 \ldots j=k}
 112 \right.$.
 113
 114 \s{Dùsledek:} \quad $\Omega \cdot \overline{\Omega} = nE$.
 115
 116
 117 \>Jedná se o skalární souèin (jako pøedtím, èili prvek na pozici $ij$ je $0$ nebo $n$) $\Rightarrow\Omega^{-1} = {1 \over n} \overline{\Omega}$.
 118
 119
 120 \>Na¹li jsme inverzi:
 121
 122 $\Omega({1 \over n}  \overline{\Omega}) = {1 \over n}\Omega \cdot \overline{\Omega} = E$, \quad
 123 $\Omega^{-1}_{jk} = {1 \over n}\overline{\omega^{jk}} = {1 \over n}\omega^{-jk} = {1 \over n}  {(\omega^{-1})}^{jk}$, \quad
 124 kde $\omega^{-1}$ je $\overline{\omega}$.
 125
 126
 127 \>Ná¹ algoritmus poèítá tedy i inverzní transformaci, pouze místo $\omega_n$  pou¾ijeme $\overline{\omega_n}$ a vydìlíme $n$. Co¾ je skvìlé -- staèí znát pouze jeden algoritmus u~kterého staèí v~jednom pøípadì pou¾ít jinou matici a vydìlit $n$.
 128
 129 \s{Vìta:} Pro $n= 2^k$ lze DFT na ${\bb C}^n$ spoèítat v~èase $\O(n \log n)$ a DFT$^{-1}$ takté¾.
 130
 131 \s{Dùsledek:}
 132
 133 \>Polynomy stupnì $n$ lze násobit v èase $\O(n \log n)$:
 134 $\O(n \log n)$ pro vyhodnocení, $\O(n)$ pro vynásobení a $\O(n \log n)$ pro pøevedení zpìt.
 135
 136 \s{Pou¾ití:}
 137
 138 \itemize\ibull
 139
 140 \:Zpracování signálu -- rozklad na siny a cosiny o~rùzných frekvencích $\Rightarrow$ spektrální rozklad.
 141 \:JPEG.
 142 \:Násobení dlouhých èísel v èase $\O(n \log n)$.
 143 \endlist
 144
 145
 146 \figure{img.eps}{Pøíklad prùbìhu algoritmu na vstupu velikosti 8}{3in}
 147
 148
 149 \>To je schéma zapojení kombinaèního obvodu (tzv. \uv{motýlek}).
 150
 151 \s{Z toho:}
 152
 153 \itemize\ibull
 154 \:Kombinaèní obvod pro DFT
 155         s~$\O(\log n)$ hladinami
 156         a $\O(n)$ hradly na hladinì.
 157 \:Nerekurzivní algoritmus (postupujeme zleva) v~èase $\O(n \log n)$.
 158    Èísla vstupu jsou èísla v~binárním tvaru pøeètená pozpátku.
 159 \endlist
 160
 161
 162 \bye