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 "hloupý" algoritmus - jednodu¹e vezmu 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 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ím v¾dy vynásobit "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 \s{ Trochu algebry na zaèátek: }
   8
   9 Libovolný polynom P øádu n mù¾eme mít reprezentován 2 rùznými zpùsoby:
  10
  11 \itemize\ibull
  12 \: jeho koeficienty, èili èísly $a_{0}, a_{1}, \ldots ,a_{n}$
  13 \: jeho hodnotami v n + 1 rùzných bodech $x_{0}, x_{1}, \ldots , x_{n+1}$ èili èísly $P(x_{0} ), P(x_{1} ), \ldots ,P( x_{n+1} )$
  14 \:
  15 \endlist
  16
  17 Pov¹imnìme si jedné skuteènosti - máme-li 2 polynomy $A$ a $B$ øádu $n$ a $x_{0}\ldots, x_{k}$ body, pak platí $C(x_{k}) = A(x_{k}) * 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 (= rychlej¹í ne¾ kvadratický, jinak bychom si nepomohli oproti hloupému algoritmu) na pøevod polynomu z jedné reprezentace do druhé a zase zpìt.
  18
  19 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.
  20
  21 \s{ Idea, jak by mìl algorimus pracovat: }
  22 \itemize\ibull
  23 \: vybereme 2n + 2 bodù $x_{0}, x_{1}, \ldots , x_{2n+2}$
  24 \: v tìchto bodech vyhodnotíme polynomy $A$ a $B$
  25 \: nyní ji¾ v lineárním èase získám polynom C (viz vý±e)
  26 \: invernznì pøevedu hodnoty polynomu C v 2n+2 bodech na jeho koeficienty
  27 \endlist
  28
  29 Je asi vidìt, µe klíèové jsou kroky 2 a 4. Vybrání bodù jistì stihneneme pohodlnì v lineárním èas a vynásobení samotných hodnot té¾ (máme 2n+2 bodù, a $C(x_{k}) = A(x_{k}) * B(x_{k}), k = 0,1,2 \ldots 2n+2$, tak¾e na to nepotøebujeme více ne¾ $2n+2$ násobení).
  30
  31 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.
  32
  33 \s{ Vyhodnocení polynomu metodou rozdìl a panuj (algoritmus FFT): }
  34 Mìjme polynom P øádu $n$ a chceme jej vyhonotit 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}$.
  35
  36 Polynom P rozlo¾íme na dvì èásti, první obsahuje èleny se sudými koeficienty, druhá s lichými:
  37
  38 $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}$
  39
  40 $S(x^{2}) = p_{0}x^{0} + p_{2}x^{2} + ... + p_{n - 2}x^{n - 2}$
  41 $L(x^{2}) = p_{1}x^{1} + p_{3}x^{3} + ... + p_{n - 1}x^{n - 1}$
  42
  43 Tak¾e obecnì $P(x) = S(x^{2}) + xL(x^{2})$ a $P(-x) = S(x^{2}) - xL(x^{2})$
  44 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}$).
  45
  46 Pøíklad:
  47 $3 + 4x + 6x^{2} + 2x^{3} + x^{4} + 10x^{5} = (3 + 6x^{2} + x^{4}) + x(4 + 2x^{2} + 10x^{4}):$
  48
  49
  50 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 v u¾ v druhém levlu 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 1. Máme n n-tých odmocnin z 1, 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}$. 2 poznámky:
  51 \itemize\ibull
  52 \: n-té odmocniny z 1 jsou spárované, èili $\omega^{j} = -\omega^{n/2 + j} $
  53 \: umocníme-li v¹echny na druhou, vznikne nám n/2 n-pùltých odmocnin z 1, které jsou i nadále spárované
  54 \endlist
  55
  56 \s{Tak a teï koneènì ten slavný algoritmus: }
  57 FFT(P, $ \omega$)
  58 Vstup: $p_{0} \ldots p_{n-1}$ to jest vektor koeficientù polynomu P a $\omega$ co¾ je n-tá odmocina 1
  59 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})$
  60 \algo
  61 \: pokud n = 1, vra» $P_{0}$ a konec
  62 \: jinak rozdìl P na sudé a liché koeficienty a zarekurzi se FFT(S, $\omega^{2}$) a FFT(L, $\omega^{2}$)
  63 \: pro $j = 0 \ldots n - 1$ spoèítej: $P(\omega^{j}) = S(\omega^{2j}) + \omega^{j} * L(\omega^{2j})$
  64
  65 \endalgo
  66
  67
  68 \s{ Èasová slo¾itost:}
  69  $T(n)=2T({n \over 2} ) + O(n)  \Rightarrow$ slouitost  $O(n)$ ... stejnì jako MergeSort.
  70
  71
  72
  73
  74 Máme tedy algoritmus, který "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.
  75 Definuje me si DFT algoritmus, která vyu¾ívá maticovou reprezentaci a s jeho¾ pomocí získáme hledaný algoritmus.
  76
  77 \s{Definice:}
  78
  79  {\sc Diskretní Fourierova transformace} $(DFT)$
  80 je funkce $f:  { C ^n} -> { C ^n}$ , ze  $y=f(x) \equiv \forall j \  y_{j} = \sum \limits ^{n-1}_{k=0} x_{k} . \omega ^{k}  $
  81
  82 \s{Poznámka:}
  83
  84 Vezmeme polynom, který má $x_{kj}$ jako koeficienty a vyhodnotíme ho v bodì $\omega ^{j} [y_{j} = x(\omega^{j})] =>  {f}$  je linearní $=>$ mù¾eme  napsat $f(x) = \Omega_{x} ,\  \Omega _{jk} =\omega ^{jk}$, kde $\Omega$ je matice
  85
  86
  87 \s{Jak najít inverzní matici?} Víme ¾e $\Omega =\Omega ^{T}$ protoue $\omega ^{jk} = \omega ^{kj}$ .
  88
  89 \s{Jak vypadají øádky této matice?} Vyu¾ijeme následující lemma, které si ale napøed doká¾eme :)
  90
  91 \s{Lemma:}
  92
  93 \proof souèin
  94 $$\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} $$
  95 Proto¾e $ \overline{\omega^{kl}} = \overline{\omega} ^{kl} = {({1 \over \omega} )}^{kl} = \omega ^{-kl}$.
  96
  97 \itemize\ibull
  98 \: 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$
  99
 100 \: Pokud $j=k \sum \limits ^{n-1}_{l=0} (\omega ^{0}) ^{l} = n$.
 101 \endlist
 102 \qed
 103
 104 a nyní slibované a u¾ i dokázané lemma
 105
 106 \s{Lemma:} $\Omega _{j} . \Omega _{k} =$
 107
 108 \itemize\ibull
 109 \: $0$ pro $j\neq k$
 110
 111 \: $1$ pro $j=k$
 112 \endlist
 113
 114 \s{Dùsledek:}$$\Omega*\overline{\Omega} = nE $$
 115
 116
 117 \>Jedná se o skalární souèin (jako pøedtím, èile prvek na pozici $ij$ je $0$ nebo $n$) $=>\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*\overline{\Omega} = E$
 123
 124 $\Omega^{-1}_{jk} = {1 \over n}\overline{\omega^{jk}} = {1 \over n}\omega^{-jk} = {1 \over n}  {(\omega^{-1})}^{jk}  $
 125
 126 kde $\omega^{-1}$ je $\overline{\omega}$
 127
 128
 129 \>Ná¹ algoritmus poèítá tedy i inverzní transformaci, akorát 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ì po¾ít jinou matici a vydìlit $n$!
 130
 131 \s{Vìta:} pro $n= 2^k$ lze $DFT$ na $C^n$ spoèítat v èase $O(n \log n)$ a $DFT^{-1}$ takté¾
 132
 133 \s{Dùsledek:}
 134
 135 Polynomy stupnì $n$ lze násobit v èase $O(n \log n)$.
 136
 137 $O(n \log n)$ pro vyhodnocení, $O(n)$ pro vynásobení a $O(n \log n)$ pro pøevedení zpìt
 138
 139 \s{Pou¾ití:}
 140
 141 \itemize\ibull
 142
 143 \: zpracování signálu - rozklad na $\sin$ a $\cos$ o rùzných frekvencích $=>$ spektrální rozklad
 144 \: JPEG
 145 \: násobení dlouhých èísel v èase $O(n \log n)$
 146 \endlist
 147
 148
 149 \figure{img.eps}{Pøíklad prùbìhu algoritmu na vstupu velikosti 8}{3in}
 150
 151
 152 Je to schéma zapojení kombinaèního obvodu (tzv. "motýlek")
 153
 154 \s{Z toho:}
 155
 156 \itemize\ibull
 157 \: kombinaèní obvod pro DFT
 158         s $O(\log n)$ hladinami
 159         a $O(n)$ hradly na hladinì
 160 \: nerekurzivní algoritmus (postupujeme z leva) v èase $O(n \log n)$
 161
 162 èísla vstupu jsou èísla v binárním tvaru 0-7 pøeètená pozpátku
 163 \endlist
 164
 165
 166 \bye