21
2 FOURIER ANALYSE EN DE “FAST FOURIER TRANSFORM”
2
Fourier analyse en de “Fast Fourier Transform”
Zij f een continue 2π-periodieke funktie op IR (eventueel met complexe waarden), dan kunnen we f ontwikkelen in een Fourier-reeks, f (x) =
+∞ X
fk eikx ,
(2.1)
k=−∞
fk := Zonder bewijs vermelden we:
( f , eikx ) 1 = 2 ikx 2π ke k2
Z
2π
f (x) e−ikx dx .
(2.2)
0
Stelling. Als f ∈ C(IR), en f (x + 2π) = f (x), ∀ x ∈ IR, dan geldt: (i) k f −
Pn
fk eikx k2 → 0 als n → ∞,
k=−n
(ii) Een parti¨ele som Σnk=−n fk eikx is de beste benadering van f in L2 -zin in de deelruimte van trigonometrische polynomen van graad ≤ n (dit is de ruimte opgespannen door { eikx | k ∈ ZZ , |k| ≤ n} ). (iii) De rij van CESARO-sommen N → ∞.
1 N
PN
n=0
Pn
k=−n
fk eikx convergeert uniform naar f voor
(iv) Als f ∈ C m (IR), dan dalen de Fourier-co¨effici¨enten van f minstens met de m-de macht van 1/k, | fk | ≤ k −m k f (m) k∞ . (2.3) We kunnen de transformatie (2.2) discretiseren met behulp van de trapeziumregel, Z
b
f (x) dx =
a
(b − a)3 ′′ b−a (f (a) + f (b)) + f (ξ) . 2 12
Door het interval [0, 2π] te verdelen in n deelintervallen van lengte 2π/n en de trapeziumregel op ieder deelinterval toe te passen, vinden we 1 2π
Z
0
2π
f (x) e−ikx
n−1 X 1 f (0) f (2jπ/n) e− + dx = n 2 j=1
=
X 1 n−1 f (2jπ/n) e− n j=0
2ikjπ/n
2ikjπ/n
f (2π) + Rest + 2
(2.4)
+ Rest .
Daar f periodiek is met het integratie-interval als periode, kunnen we bewijzen, dat de rest van orde O(n−m k f (m) k∞ ) is als f ∈ C m (IR), zie het hoofdstuk over integratie. Als we de funktiewaarden f (2πj/n) vervangen door xj , de j-de komponent van de vektor x ∈ IRn , en de discrete Fourierco¨effici¨enten yk noemen, de k-de komponent van de vektor y ∈ IRn , dan kunnen we de discrete Fouriertransformatie Fn op n punten neerschrijven als (Fn y)j := xj :=
n−1 X
yk e2ikjπ/n .
k=0
Fn is een transformatie van C n naar C n , waarvoor geldt:
(2.5)
2 FOURIER ANALYSE EN DE “FAST FOURIER TRANSFORM”
22
1
Stelling. n− 2 Fn is een unitaire transformatie van C n en voor de inverse van Fn geldt: yk = (Fn∗ x)k =
X 1 n−1 xj e−2ikjπ/n . n j=0
(2.6)
Bewijs. De bewering volgt onmiddellijk uit de orthogonaliteit van de collectie vektoren {ek | k = 0, · · · , n − 1} met componenten ekj := e2ikjπ/n t.o.v. het discrete inprodukt (x, y) :=
n−1 X
xj yj
(2.7)
j=0
en uit het feit dat deze vektoren lengte
√
n hebben.
Fourier-transformatie wordt zeer veel gebruikt voor spektrale analyse van (bijna) periodieke verschijnselen en wordt in vele vakgebieden op grote schaal toegepast. Een van de redenen hiervoor ligt in het bestaan van een zeer snel algoritme voor het berekenen van zo’n transformatie voor het geval dat n een macht van twee is, namelijk de Fast Fourier Transform (FFT) van Cooley & Tuckey. We zullen hier een afleiding geven van deze algoritme (voor een standaard sequenti¨ele processor). Laat het aantal punten n van Fn even zijn, zeg n = 2m , dan kunnen we de discrete Fourier-som (2.6) xj =
2m−1 X
yk eikjπ/m ,
(2.8a)
k=0
verdelen in een som over de termen met even index en een som over de termen met oneven index, xj =
m−1 X
y2k e
2ikjπ/m
+ e
ijπ/m
m−1 X
y2k+1 e2ikjπ/m .
(2.8b)
k=0
k=0
Hiermee is het probleem (2.7a) met 2m punten gesplitst in twee problemen met m punten. D.w.z. als u de vektor van komponenten met even index is en v de vektor van komponenten met oneven index (u, v ∈ C m en y ∈ C 2m ), dan geldt (ga na!) voor j = 0, · · · , m − 1 : = (Fm u)j
+ (Fm v)j eijπ/m
xj+m = (F2m y)j+m = (Fm u)j
− (Fm v)j eijπ/m
xj
=
F2m y)j
(2.9)
Het aantal arithmetische operaties, dat nodig is om F2m y uit te rekenen, is dus gelijk aan tweemaal het aantal operaties voor Fm u plus m maal (een vermenigvuldiging plus een evaluatie van een emacht plus twee optellingen). Als n een tweemacht is, n = 2t , dan is het totale aantal dus gelijk aan (2t + 2 · 2t−1 + · · · ) maal de basiseenheid (1 V + 1 E + 2 O) (= n 2 log n maal de basiseenheid). Om met de reduktieformules (2.9) een nette algoritme neer te schrijven, herschrijven we de opdracht “maak een 2n-punts Fouriertransformatie (n is een macht van 2) van de het (complexe) array y[p] · · · y[p + 2n − 1]” in de recursieve procedure
2 FOURIER ANALYSE EN DE “FAST FOURIER TRANSFORM”
23
procedure F2T( var x : vector; p, m : integer ); (∗ overschrijf het deelarray x[p · · · p + 2 ∗ m − 1] met zijn 2 ∗ m–punts Fourier-getransformeerde ∗) var y : vector; j : integer; begin (∗ splits array in twee stukken; elementen met even index vooraan en oneven achteraan ∗) for j := 0 to m − 1 do (∗ Verwisselingsstap ∗) y[p + j] := x[p + 2 ∗ j]; y[p + m + j] := x[p + 2 ∗ j + 1]; end; if m > 2 then (∗ FT op beide delen apart ∗) F2T(y, p, m div 2); F2T(y, p + m, m div 2); end for j := 0 to m − 1 do (∗ Assemblagestap ∗) x[p + j] := y[p + j] + exp(ijπ/m) ∗ y[p + m + j]; x[p + m + j] := y[p + j] − exp(ijπ/m) ∗ y[p + m + j]; end; end;
De recursieve algoritme bestaat dus uit een verwisselingsstap, gevolgd door een recursieve oproep van FT op beide helften. Na terugkeer uit deze recursie volgt een assemblagestap. Op het diepere recursieniveau is dit hetzelfde. In het totaal worden dus eerst alle verwisselingen gedaan op steeds meer steeds kortere deelarrays, totdat de recursie is aangekomen op het niveau n = 2. Daarna volgen alle assemblagestappen, te beginnen met deelarrays van lengte 2. De verwijdering van de recursiviteit leidt dus tot de volgende algoritme: procedure FFT( var x : vector; n : integer ); (∗ overschrijf het array x[0 .. n] met zijn n–punts Fourier-getransformeerde. ∗) var y : vector; j, m, p : integer; begin m := n ÷ 2; while m > 1 do (∗ splits alle deelarrays van lengte 2 ∗ m ∗) for p := 0 to n − 1 by 2 ∗ m do for j := 0 to m − 1 do y[p + j] := x[p + 2 ∗ j]; y[p + m + j] := x[p + 2 ∗ j + 1]; end; end; x := y ; m := m ÷ 2; end; while m < n do (∗ assembleer deelarrays van lengte m ∗) for p := 0 to n − 1 by 2 ∗ m do for j := 0 to m − 1 do y[p + j] := x[p + j] + exp(ijπ/m) ∗ x[p + m + j]; y[p + m + j] := x[p + j] − exp(ijπ/m) ∗ x[p + m + j]; end; end; x := y ; m := m ∗ 2; end; end;
Bij een verwisselingsstap worden de elementen met even index naar de eerste helft verplaatst en die met oneven index naar de tweede helft. De binaire representatie van een even getal (van p binaire cijfers) eindigt met een “nul” en de binaire representatie van een element in de eerste helft (van 0 · · · 2p − 1) begint met een nul; evenzo eindigt de binaire representatie van een oneven getal met een “een” en die van een element van de tweede helft begint met een “een”. Een verwisselingsstap op een array van 2p elementen (met indices 0 · · · 2p − 1) betekent dus, dat ieder element een nieuwe plaats krijgt met een index, die we eenvoudig uit de binaire representatie van de oude kunnen afleiden door het laaste bit vooraan te plaatsen. In de volgende verwisselingsstap op het diepere recursieniveau werken we met array van halve lengte en gebruiken we dus alleen de
24
2 FOURIER ANALYSE EN DE “FAST FOURIER TRANSFORM”
laatste p − 1 bits. Het p-de bit wordt op het diepere niveau dus niet meer veranderd. De uitkomst van alle verwisselingen tesamen is dus, dat ieder element verplaatst wordt naar een positie met een index, waarvan de binaire representatie het gespiegelde is van de oorspronkelijke index. We hoeven het array dus maar ´e´en keer te doorlopen om alle verwisselingen te doen. In het geval, dat ons uitgangsarray lengte 8 heeft, kunnen we alle verwisselingsstappen alsvolgt tekenen: 7
111
6
110 Z
- 7
1
: »» »» » Z» » 101 * © Z Z©© 100 H ©©Z ½ HH Z> © ½ ~ © H½ Z 011 ½HH j H 010 XX½ ½ XXX XX z 001 ½
0
000
5 4 3 2
- 7
111 110
5
: 6 110 XXX »»» XX »»» XX z 5 » 101
4
100
100
3
011
- 4 - 3
2
: 2 010 XXX »»» XX »»» XX » z 1 001
010
- 0
000
6
1
- 0
111
000
101 011 001
Figure 1: Verwisselingspatroon voor een array van lengte 8. De indices zijn binair uitgeschreven. We zien, dat het bitpatroon van alle indices precies wordt omgedraaid. Bij het uitvoeren van alle assemblagestappen kunnen we tenslotte een hoop werk uitsparen op het berekenen van de exponenten door alle assemblages, waarin exp(ijπ/n) gebruikt wordt, tegelijk te doen en door de beide for-lussen te verwisselen en exp(ijπ/n) in de j-de slag te berekenen als de exponent van de vorige slag maal de basisexponent exp(iπ/n) . Dit geeft de effici¨entere algoritme: procedure FFT( var x : vector; n : integer ); (∗ overschrijf het array x[0 .. n − 1] met zijn n–punts Fourier-getransformeerde. ∗) var j, m, p : integer; E, Z, W, : real begin for p := 0 to n − 1 do j := gespiegelde(p); if j > p then verwissel(x[ j ] , x[ p ]) end end; m:=1; while m < n do (∗ assembleer deelarrays van lengte m ∗) E := exp(iπ/n); Z := 1; for j := 0 to m − 1 do for p := 0 to n − 1 by 2 ∗ m do W := Z ∗ x[p + m + j]; x[p + m + j] := x[p + j] − W ; x[p + j] := x[p + j] + W ; end; Z = Z ∗ E; end; m := m ∗ 2; end; end;
Opgave 2.1: Analoog aan de hier afgeleide FFT op N punten met N = 2t een tweemacht, welke gebaseerd is op herhaalde tweedeling, kunnen we ook een snel algoritme afleiden op basis van driedeling als N een macht van 3 is. Leid hiervoor de algoritme af.
REFERENCES
112
References [1] M. Hestenes & E. Stiefel, Methods of conjugate gradients for solving linear systems, J. Research NBS, 49, pp. 409 – 436, 1952. [2] C. Lanczos, An iteration method for the solution of the eigenvalue problem of linear differential and integral operators, J. Research NBS, 45, pp. 255 – 282, 1950. [3] J.K. Reid, On the method of conjugate gradients for the solution of large sparse systems of linear equations, Proc. Conf. on Large Sparse Sets of Linear Equations, Academic Press, New York, 1971. [4] J.A. Meijerink and H.A. van der Vorst, An iterative solution method for linear systems of which the coefficient matrix is a symmetric M-matrix, Math.of Comp., 31, pp. 148 – 162, 1977. [5] G.H. Golub & C.F. Van Loan, Matrix Computations, The Johns Hopkins University Press, Baltimore, Maryland, USA, 1ste druk, 1983, 2de druk, 1988, 3de druk, 1995. [6] R. Bulirsch & J. Stoer, Introduction to Numerical Analysis, Springer Verlag, Berlin, 1977. (Ook verkrijgbaar in een goedkope duitstalige pocketeditie). [7] D. Kincaid & W. Cheney, Numerical Analysis, Brooks & Cole Publishing Company, Pacific Grove, California, USA, 1991; 2de druk, 1996.
