1. No category

Fouriertransformation och spektralmetoder
Michael Hanke
Skolan för teknikvetenskap
SF1538 Projekt i simuleringsteknik
1 (28)
Introduktion
Varför Fouriertransformation?
• Många processer är periodiska.
• Den snabba Fouriertransformationen är ett kraftfullt vertyg för att
undersöka sådana processer, både teoretiskt och numeriskt
I detta kapitel gör vi tre saker:
• Vi undersöker hur man kan approximera periodiska funktioner med
hjälp av sinus-/kosinusfunktioner och komplexa
exponentialfunktioner.
• Vi introducerar en effektiv algoritm för att göra det numeriskt (FFT).
• Vi använder Fouriertransformationens egenskaper för att konstruera
en ny numerisk metod för att lösa partiella differentialekvationer.
2 (28)
Periodiska funktioner
• Vi undersöker periodiska funktioner f definerad på R med perioden
2π, dvs f (x + 2π) = f (x) för alla x ∈ R. Därför räcker det med att
anta att funktionen är definerat på intervallet [−π, π].
• Funktionen kallas för jämn om f (x) = f (−x) för alla x ∈ R.
• Funktionen kallas för udd om f (x) = −f (−x) för alla x ∈ R.
3 (28)
Fourier sinusserie
• Prototyp av udda 2π-periodiska funktioner: sin(nx).
• Fourier sinusserie
S(x) =
∞
X
bn sin(nx)
n=1
Q: Vilka funktioner f : [−π, π] → R kan representeras bra
genom en sinusserie?
A: Väldigt många (udda) funktioner om man mäter felet i
minsta kvadratmening (minsta kvadratanpassning),
Z π
1
(SN (x) − f (x))2 dx.
kSN − f k2 =
2π −π
Här betecknar SN de första N termer i S.
4 (28)
Hur beräknar man koefficienterna?
• MVi har att:
(
0, om n 6= k,
sin(nx) sin(kx)dx =
π, om n = k.
−π
Z
π
• Antag att (f ojämn)
∞
X
f (x) =
bn sin(nx).
n=1
• Multiplicera ekvationen med sin(kx) och integrera:
Z
π
f (x) sin(kx)dx =
−π
∞
X
Z
bn =
1
π
Z
sin(nx) sin(kx)dx = bk π.
−π
n=1
• Därmed får vi
π
bn
π
f (x) sin(nx)dx.
−π
5 (28)
Exempel: Rektangulär våg
Bild fattas!!!!


−1, om x ∈ (−π, 0),
SW (x) = 1,
om x ∈ (0, π),


0
om x = −π, 0, π.
Fourier sinusserie:
sin 3x
sin 5x
sin 7x
4 sin x
SW (x) =
+
+
+
+ ···
π
1
3
5
7
6 (28)
Exempel (cont)
2
1.5
1
0.5
0
−0.5
−1
−1.5
−2
−3
−2
−1
0
1
2
3
7 (28)
Gibbs fenomen
2
1
1.5
1
0.8
0.5
0.6
0
−0.5
0.4
−1
−1.5
SW
N=9
N=15
N=21
0.2
−2
0
−3
−2
−1
0
1
2
3
2.2
2.4
2.6
2.8
3
3.2
3.4
Gibbs fenomen: Partialsummor håller konstant avstånd nära
diskontinuiteter!
(men konvergerar i minsta kvadratmening)
8 (28)
Fourier kosinusserie
• För jämna funktioner är det lämpligt att använda cosinusfunktioner
cos(nx).
• Fourier cosinusserie:
C (x) = a0 +
∞
X
an cos(nx).
n=1
• Vi har


om n 6= k,
0,
cos(nx) cos(kx)dx = 2π, om n = k = 0,

−π

π, om n = k > 0.
Z
π
• På samma sätt som ovan får vi (f jämn)
(
an =
Rπ
1
2πR −π f (x)dx,
1 π
π −π f (x) cos(nx)dx,
om k = 0,
om k > 0.
9 (28)
Exempel: Ramp
Funktionen fås genom att integrera SW :
RR(x) = |x|.
4
3.5
3
2.5
2
1.5
1
0.5
0
−0.5
−3
−2
−1
0
1
2
3
Fourier kosinusserie:
π π cos x
cos 3x
cos 5x
cos 7x
RR(x) = −
+
+
+
+
·
·
·
2
4
12
32
52
72
Obs: Koefficienterna fås genom integration av de för SW .
10 (28)
Fourierserie: Allmäna periodiska funktioner
Låt f : [−π, π] → R vara en periodisk funktion. Sedan gäller:
• f = fjämn + fojämn
• fjämn = 1/2(f (x) + f (−x))
• fojämn = 1/2(f (x) − f (−x))
Därför kan vi skriva f som en summa av sinus- och kosinusserier:
f (x) = fjämn + fojämn
∞
∞
X
X
= a0 +
an cos(nx) +
bn sin(nx).
n=1
n=1
11 (28)
Fourierserie: Den komplexa versionen
• Moivres teorem: e iα = cos α + i sin α.
• Låt ck = (ak − ibk )/2,
c−k = (ak + ibk )/2 = c̄k .
ck e ikx + c−k e −ikx = ck (cos kx + i sin kx) + c−k (cos kx − i sin kx)
= (ck + c−k ) cos kx + i(ck − c−k ) sin kx
= ak cos kx + bk sin kx.
• Därmed får vi
f (x) =
∞
X
cn e inx .
n=−∞
• Från och med nu använder vi alltid den komplexa versionen!
• Obs: Även om f är en reell funktion så är koefficienterna cn
vanligtvis komplexa.
12 (28)
Fourierserie (cont)
• Vanligt beteckning:
fˆn ≡ cn .
• Det är enkelt att verifiera: n = 0, ±1, ±2, . . .
1
fˆn =
2π
Z
π
f (x)e −inx dx
−π
13 (28)
Intressanta egenskaper
• Parsevals identitet
Z
π
|f (x)|2 dx = kf k22 = 2π
−π
∞
X
|fˆn |2
n=−∞
Därför gäller
∞
X
f ∈ L2 (−π, π) ⇔
|fˆn |2 < ∞.
n=−∞
• Derivator:
f (p) (x) =
∞
X
(in)p fˆn e inx .
n=−∞
• Låt
p
Hper
vara mängden av all 2π periodiska funktioner som är p
gånger (generaliserad) deriverbara. Då gäller
f ∈
p
Hper
⇔
∞
X
k 2p |fˆn |2 < ∞.
n=−∞
Koefficienterna avtar snabbare ju oftare f är deriverbar!
14 (28)
Intervalltransformationer
Q: Antag att vi använder “basintervallet” (0, 2π) istället för
(−π, π). Vad händer med koefficienterna?
A: Ingenting eftersom alla inblandade funktioner har samma
period 2π.
Q: Antag att funktionen f har perioden T istället för 2π. Hur
kommer vi fram till en lämplig Fourierserie?
A: Härledningen görs i två steg:
• Perioden T måste vara en “multiple” av perioden av
de komplexa exponentialfunktionerna. Låt ω = 2π/T .
Sedan görs ansatsen
∞
X
f (x) =
fˆn e inωx
n=−∞
• Genom ansatsen ovan fås
1
fˆn =
T
Z
T
f (x)e
0
−inωx
1
dx =
T
Z
T /2
f (x)e −inωx dx.
−T /2
15 (28)
Diskret Fouriertransformation (DFT)
• Standardversionen av DFT bygger på intervallet (0, 2π).
• Låt (0, 2π) vara indelat i N ekvidistanta delinterval,
2π
, xj = jh.
N
• Integralen approximeras genom trapetsregeln. Med
f (0) = f (2π) = f (xN ) fås
Z 2π
N−1
N−1
1
h X −inxj
1 X
f (x)e −inx dx ≈
fj e
=
fj (e −ih )jkn
2π 0
2π
N
h=
j=0
j=0
=
1
N
N−1
X
fj w̄ jn
j=0
=: f˜n
Här: fj = f (xj ) och w = e ih .
• Sats: (invers DFT)
fj =
N−1
X
n=0
f˜n w nj
16 (28)
Matlabs DFT
• Matlab har två funktioner för DFT och den inversa DFT: fft och
ifft.
• Båda använder en annan normering: Faktorn 1/N finns inbyggd i
ifft, medan den fattas i fft.
• Ibland delas faktorn så att båda transformationer ser mer
√
“symmetriska” ut: Istället för 1/N används faktorn 1/ N i DFT och
dess invers.
17 (28)
Snabb Fouriertransformation (FFT)
• Antalet räkneoperationer för DFT: N 2 komplexa multiplikationer och
additioner (vi antar att de f -oberoende koefficienterna w jk redan har
beräknats i förväg.) Väldigt dyrt!
• Om N är en potens av 2 så kan man härleda en
(divide-and-conquer) algoritm som bara har komplexiteten
O(N log N): Cooley och Tukeys FFT
• Exempel: N = 212
• naiv algoritm: 224 ≈ 1.7 · 107 operationer
• FFT: 6 × 1012 ≈ 2.4 · 104 operationer (nästan 1000 gånger snabbare)
18 (28)
Förskjuten FFT
• Att använda intervallet (0, 2π) leder till standardalgoritmen.
• Vad händer om man använder intervallet (−π, π) istället?
• Lite algebra leder till (N jämn!)
1
f˜n =
N
N/2−1
X
fl e −ihnl
l=−N/2
N/2−1
fl =
X
f˜n e ihnl
n=−N/2
• De diskreta koefficienterna kommer att ha olika ordningar (men är
identiska!):
DFT
förskjuten DFT
(f˜0 , f˜1 , . . . , f˜N−1 )
(f˜N/2 , f˜N/2+1 , . . . , f˜N−1 , f˜0 , . . . , f˜N/2−1 )
• Det är just vad matlabs fftshift åstadkommer!
19 (28)
Spektralinterpolation
• DFT och dess invers kan redan beräknas om bara de diskreta värden fj är
kända.
• Exponentialuttryck kan däremot enkelt fortsättas på hela intervallet
(0, 2π) resp (−π, π):
e ihnl = e inxl
• För givna fj kan vi beräkna deras DFT f˜l och sedan definera
ΠN f (x) :=
N−1
X
f˜n e inx
n=0
• Enligt Sats är Πf en interpolation till funktionen f !
Interpolationsfunktionen kallas för trigonometriskt polynom.
• På samma sätt kan vi använda den “symmetriska” versionen:
N/2−1
Πc,N f (x) :=
X
f˜n e inx
n=−N/2
• Obs: ΠN f 6= Πc,N f !!
• Mellan interpolationsnoderna xl har Πc,N vanligtvis bättre noggrannhet än
ΠN .
20 (28)
Spektralinterpolation (cont)
• Vi såg tidigare att koefficienterna i Fourierserien avtar väldigt snabbt
om funktionen f är många gånger deriverbar.
• Det antyder att Fourierserien konvergerar ju snabbare ju oftare
funktionen är deriverbar.
• Man kan visa: Om funktionen f tillhör C ∞ [−π, π], så är
kf − Πc,N k2 ≤ C1 e −C2 N
• Konvergensen är snabbare än vilken polynom som helst. Det kallas
för exponentiell konvergens.
21 (28)
Spektralinterpolation (cont)
Q: Kan en funktion entydigt återfås genom att använda
diskreta värden?
A: Nej! Frågan går ut på att kunna konstruera en entydig
interpolation.
Däremot går det bra om vi använder a-priori information
(t ex funktionen är ett polynom av ett visst grad).
Q: Vilken klass kan vi använder om vi bara tittar på
periodiska funktioner och trigonometriska polynom?
A: Tag den enkla periodiska funktionen f (t) = ae i(ωx+φ) . Vi
behöver (minst) två värden i intervallet [0, ω/(2π)) för att
kunna identifiera de två parametrarna entydigt.
22 (28)
Nyquist sampling theorem
• Låt T (= h) vara steglängden för “sampling” av signalen.
• Definition: Nyquist sampling rate: T = π/ω.
Sats:
Låt sampling rate T vara given. Frekvenserna som är högre än Nyquist
frequensen ωN = π/T kan inte bestämmas. En signal med högre frekvens
än ωN mappas till en signal med lägre frekvens. Denna effekt kallas för
aliasing.
23 (28)
Spektralmetoder: Derivering
Ide: Om vi har en funktion given vid diskreta noder så kan vi få
en approximation till derivatan genom att interpolera
funktionen med en deriverbar funktion p(x) och sedan
sätter u 0 (xj ) ≈ p 0 (xj ).
Exempel
uj +1 −uj
h
u −u
u 0 (xj ) ≈ j +12h j −1
• Stegvis linjär interpolation. u 0 (xj ) ≈
• Stegvis kvadratiskt interpolation:
24 (28)
Derivering (cont)
Använd spektralinterpolation,
N/2−1
X
p(x) = Πc,N u =
ũn e inx
n=−N/2
Derivatan:
N/2−1
p 0 (x) =
X
ũn ine inx
n=−N/2
Approximation:
N/2−1
p 0 (xj ) =
X
ũn ine injh
n=−N/2
25 (28)
Derivering: Algoritm
1
Låt vektorn u vara given.
2
Beräkna FFT ũ av u.
3
Använd fftshift för att förskjuta intervallet till (−π, π).
4
Beräkna vektorn (. . . , inũn , . . .).
5
Använd fftshift en gång till.
6
Till sist, beräkna p 0 genom ifft.
Obs: Man kan omordna algoritmen lite så att det bara behövs en
fftshift (jf anteckningarna).
26 (28)
Spektralmetoder för PDE
• Genom att använda deriveringsalgoritmen är det enkelt att
konstruera en linjemetod för PDE.
• Metoden konvergerar väldigt snabbt om lösningen är glatt och
periodiskt.
• Komplexiteten med FFT är O(N log N) jämfört med komplexiteten
O(N) för polynomiala derivator.
27 (28)
Sammanfattning
• Approximation av periodiska funktioner med hjälp av “enkla”
periodiska funktioner.
• Analys av periodiska fenomen med hjälp av DFT och
trigonometriska polynom.
• Spektralmetoder för partiella differentialekvationer.
28 (28)