Hitra Fourierjeva transformacija
Algoritem, ki je spremenil civilizacijo: od O(N²) k O(N log N) in zakaj je to razlika med nemogočim in instantnim.
Predpogoji: 3. stopnja. Potrebno: razumevanje DFT in kompleksnih eksponentov.
Leto je 1965. Računalniki so velikosti sobe, vsaka operacija stane dragocene sekunde strojnega časa. Jemata Cooley in Tukey objavita kratek članek v matematični reviji. Strokovna skupnost ga prebere z zmernim zanimanjem. Nihče ne prepozna, da gledajo enega najpomembnejših algoritmov v zgodovini računalništva — algoritma, ki bo desetletja kasneje stekel v vsakem mobilnem telefonu, vsakem MP3 predvajalniku, vsakem radarju in v vsaki medicinski napravi za slikanje možganov.
Cooley-Tukey algoritem za hitro Fourierjevo transformacijo (HFT, angleško FFT — Fast Fourier Transform) ni novo matematično orodje. Je mnogo hitrejši način izračuna istega rezultata. Razlika med počasnim in hitrim načinom je razlika med »ne bom dočakal konca izračuna« in »to traja manj kot milisekundo«.
Najprej: zakaj je navadna DFT počasna?
Spomnimo se: diskretna Fourierjeva transformacija (DFT) vzame N vzorcev signala in vrne N frekvenčnih komponent. Vsaka izhodna vrednost je vsota z vsemi vhodnimi vrednostmi. Za vsako od N izhodnih vrednosti opravimo N množenj in seštevanj. Skupaj: N × N = N² operacij.
za vsak k od 0 do N−1 → N izhodov × N operacij = N² skupaj
Za majhna N je to sprejemljivo. Toda signal dolžine tisoč vzorcev zahteva milijon operacij. Signal dolžine milijona vzorcev? Bilijon. Signal, ki se v realnem času obdeluje v telefonu? Nemogoče.
Premikajte drsnik. Za N = 256 je razlika že faktor 18. Za N = 1.048.576 (en milijon vzorcev) DFT zahteva bilijon operacij, FFT pa manj kot dvajset milijonov — razlika faktor 50.000. S sodobno frekvenco procesorja je DFT za tak signal počasna 15 minut, FFT pa 0,02 sekunde.
Ključna opažba: simetričnost rotatorjev
Zakaj je sploh možno pohitriti? Ker DFT opravlja ogromno nepotrebnega dela. Vsi tisti eksponenti e−j2πkn/N — imenujmo jih rotatorji ali twiddle faktorji — niso vsi različni. Ko jih narišemo na enotski krožnici, vidimo, da se ujemajo: rotator za k=1, n=1 je isti kot rotator za k=3, n=3 (modulo N). Polovica izračuna je podvajanje.
Opazite: pri N = 8 imamo le 8 možnih položajev rotatorjev. Toda DFT bi jih izračunala 64-krat (8×8). Cooley in Tukey sta opazila, da bi z malo pameti vsak rotator izračunali le enkrat in ga večkrat uporabili — in da se ta logika da gnezditi rekurzivno.
Ideja razpolavljanja: »razdeli in vladaj«
Osrednji trik FFT je preprost za opisati, genialen za opaziti: DFT dolžine N razdelimo na dve DFT dolžine N/2 — eno za sode in eno za lihe vzorce. Vsaka od teh se razdeli naprej. In naprej. Dokler ne pridemo do DFT dolžine 1 — ki je trivialna (ničkolikokrat je X[0] = x[0]).
Na vsakem nivoju drevesa opravimo skupaj N operacij. Drevo ima log₂N nivojev. Skupaj: N · log₂N operacij namesto N². To je celota triku.
Algoritem je elegantno kratek:
N = dolžina(x)
če N == 1: vrni x
sode = FFT(x[0, 2, 4, ..., N-2])
lihe = FFT(x[1, 3, 5, ..., N-1])
za k od 0 do N/2 - 1:
W = e−j2πk/N
X[k] = sode[k] + W · lihe[k]
X[k+N/2] = sode[k] − W · lihe[k]
vrni X
Korak »sode[k] + W · lihe[k]« in »sode[k] − W · lihe[k]« se imenuje metuljček (ang. butterfly). Vsak nivo drevesa je sestavljenec N/2 metuljčkov. V praksi se rekurzijo zamenja z iterativno različico, ki podatke razvršča z bitno obratnim zaporedjem (bit-reversal permutation) — da se izogne rekurzivnim klicem, ki bi zapravljali pomnilnik.
Zanimivost: algoritem zahteva, da je N potenca števila 2. Kadar signal ni take dolžine, se dopolni z ničlami (zero-padding) do naslednje potence 2. To je razlog, zakaj zvočni inženirji pogosto govorijo o oknih dolžin 512, 1024, 2048 — vse potence 2.
Metuljček: osnovna enota FFT
Vsak »metuljček« vzame dve vrednosti, jima pripiše rotator in vrne dve novi vrednosti. Vizualno to res spominja na metuljčka — dve črti, ki se prekrižata v sredini.
Metuljček je jedro vsega. Celoten FFT algoritem za N = 1024 vzorcev ni nič drugega kot 5120 metuljčkov, zaporedno urejenih v 10 nivojih. Vsak metuljček: dve množitvi, dve seštevalni operaciji. Nič magičnega — samo pametna organizacija enakega dela.
Vizualni dokaz: FFT v praksi
Poglejmo, kako FFT »vidi« signal. Spodaj ustvarite mešanico do treh sinusov, FFT pa bo v eni potezi prepoznal, kateri frekvenci sta prisotni — kot sopranica, ki iz hrupnega orkestra sliši posamezne instrumente.
Opazite: četudi dodate šum, ki naredi signal zgoraj videti kaotičen, spekter spodaj še vedno jasno pokaže vrhove pri pravih frekvencah. Prav to je moč Fourierjeve transformacije — iz hrupa izvleče red.
Cooley-Tukejev algoritem v radix-2 različici zahteva, da je N potenca 2, ker se N vsakič natanko prepolovi. Toda sveta dvojka ni edina možnost. Obstajajo različice za N, ki je potenca 3 (radix-3), potenca 4 (radix-4, ki je v praksi hitrejša od radix-2, ker izkoristi strukturo imaginarne enote), ali pa za N brez posebne strukture (mešani radix, Bluestein-ov algoritem, Rader-jev algoritem za praštevila).
Sodobne knjižnice kot je FFTW (Fastest Fourier Transform in the West, razvita na MIT) samodejno izberejo najboljšo strategijo glede na N in arhitekturo procesorja — in znajo vnaprej »načrtovati« optimalne vzorce dostopov do pomnilnika. FFTW je med najpogosteje citiranimi akademskimi programi v zgodovini.
Na grafičnih procesornikih (GPU) je situacija drugačna: masivna vzporednost GPE pomeni, da se splača razmeroma preproste DFT prilagoditi za vzporedno izvajanje namesto vlagati v rekurzivno strukturo. Toda osnovna ideja ostane: izogibati se nepotrebnemu delu z izkoriščanjem simetrij rotatorjev.
Ironija zgodbe je, da algoritem ni bil nov. Gauss je opisal v bistvu enako metodo že leta 1805 — toda v rokopisu, ki je bil objavljen šele posthumno in v latinščini. Različne posebne primere so odkrili še v vmesnih desetletjih.
Cooley in Tukey sta prišla do odkritja med projektom za Pentagon: ameriški seizmografi so zbirali podatke o podzemnih jedrskih preizkusih Sovietske zveze, toda analiza signalov je bila prepočasna. John Tukey je na sestanku Kennedyjevega svetovalnega odbora predlagal nov pristop. John Cooley, takrat mlad matematik pri IBM, ga je implementiral v tednih. Članek so objavili leta 1965 v reviji Mathematics of Computation.
V naslednjih mesecih je postal eden najpogosteje citiranih matematičnih del. Seveda: vsak inženir, ki se je kdajkoli ubadal z analizo signalov, je nenadoma imel orodje, ki je dotlej nedosegljive izračune naredilo trivialne. Pravijo, da je bil to eden redkih primerov, ko je matematični članek opazno spremenil industrijo v manj kot enem letu.
Časovna zahtevnost v številkah
Oglejmo si konkretne primere — koliko operacij dejansko stane DFT v primerjavi z FFT pri realnih vrednostih N:
| N (vzorci) | DFT operacije (N²) | FFT operacije (N·log₂N) | Pohitritev |
|---|---|---|---|
| 64 | 4.096 | 384 | 10,7× |
| 1.024 | 1.048.576 | 10.240 | 102× |
| 1.048.576 | 1.099.511.627.776 | 20.971.520 | 52.429× |
Milijon vzorcev zvoka pri 44.100 Hz (standardna kakovost CD) traja dobrih 22 sekund. DFT bi to obdelala v urah. FFT to stori v delčku sekunde na navadnem telefonu. Brez FFT ni predvajanja glasbe v realnem času, ni Zoom klica, ni radarskega sistema, ki bi deloval dovolj hitro, da bi bil koristen.
Zakaj N mora biti potenca 2?
Razpolavljanje deluje gladko le, kadar se N da deliti z 2 večkrat, vse do 1. To pomeni, da N mora biti potenca 2: 2, 4, 8, 16, …, 1024, 2048, … Kadar naš signal nima te dolžine, ga preprosto dopolnimo z ničlami (zero-padding) do naslednje potence 2. V praksi to niti ne škoduje — ničle ne dodajajo nobene informacije, le omogočijo algoritem.
Za prav ljubitelje teorije: obstajajo FFT algoritmi za poljuben N (Bluestein-ov algoritem), toda v praksi se jim izogibamo prav s tehniko ničel.
Kar zdaj veste
FFT ni nova matematika — je DFT izračunana z drastično manj dela, ker izkoriča simetričnost rotatorjev. Z razpolavljanjem problema v dve polovici (in rekurzivno naprej) dosežemo O(N log N) namesto O(N²). Na vsakem nivoju drevesa operirajo metuljčki — pari seštevanj s fazo. Rezultat je isti, pot pa 50.000-krat krajša. Na naslednji stopnji si ogledamo, kam vse ta algoritem seže: od MRI aparatov do umetne inteligence.