Praktikum II “Numerik am Computer” FS 2021
Prof. Dr. H. Harbrecht Universität Basel
A. Bot, J. Bopp, M. Fallahpour, H. Gubler, F. Rios, C. Santos
Beilage zur Serie 6
Trigonometrische Interpolation und schnelle Fourier-Transformation Trigonometrische Interpolation
Wir möchten eine Funktion interpolieren, von der wir wissen, dass sie 2π-periodisch ist, d.h.
f(x) =f(x+ 2πk)für alle x∈[0,2π) und k ∈Z,
Dazu seien eine gerade Anzahl von n = 2m Stützstellen(x0, y0),(x1, y1), . . . ,(xn−1, yn−1) gegeben, wobei die Knoten {xk}n−1k=0 äquidistant verteilt seien mit
xk:= 2πk
n , und yk :=f(xk) für k= 0, . . . , n−1.
Eine weitverbreitete Methode ist, ein Interpolationspolynom der Form P(x) = a0
2 +
m−1
X
`=1
[a`cos(`x) +b`sin(`x)] + am
2 cos(mx) (1)
durch die Stützstellen aufzustellen, das die Periodizität der Funktionensin(`x)undcos(`x) ausnutzt. Durch einen Wechsel in das Komplexe können wir die Koeffizienten a`, b` ∈ R mit Hilfe des trigonometrischen Polynoms
Π(x) = β0 +β1eix +β2e2ix+· · ·+βn−1e(n−1)ix mit β` ∈C und
Π(xk)=! yk für alle k = 0,1, . . . , n−1 bestimmen. Mit Hilfe der Beziehung
eix = cos(x) +isin(x)
findet man (siehe Vorlesung “Einführung in die Numerik”), dass
a0 = 2β0, am = 2βm, a` =β`+βn−`, b` =i(β`−βn−`), `= 1, . . . , m−1.
Beachte, dass zwar P(xk) = Π(xk) = yk für k = 0, . . . , n−1 gilt, im Allgemeinen aber P(x)6= Π(x)ist, wennx6=xk. Durch Nachrechnen kann man zeigen, dass die Koeffizienten β` gegeben sind durch
β` = 1 n
n−1
X
m=0
ymωn−m`, ` = 0,1, . . . , n−1, (2) wobei ωn:= e2πi/n die n-te komplexe Einheitswurzel bezeichnet.
1
Schnelle Fourier-Transformation
Alle Koeffizienten β`gemäss (2) zu berechnen hat einen Aufwand vonO(n2), was für viele Anwendungen (z.B. in der digitalen Signalverarbeitung) zu teuer ist. Solche Anwendungen benutzen stattdessen die schnelle Fourier-Transformation, die nur einen Aufwand von O(nlog2n)hat.
Schreibt man
Tn =
ωn0 ωn0 ωn0 . . . ωn0 ωn0 ωn1 ωn2 . . . ω(n−1)n
ωn0 ωn2 ωn4 . . . ω2(n−1)n
... ... ... ... ωn0 ω(n−1)n ωn2(n−1) . . . ω(n−1)(n−1)n
, β=
β0 β1 ... βn−1
und y=
y0 y1 ... yn−1
,
so ergibt sich aus (2), dass
β = 1 nT?ny.
Zur Erinnerung: T?n bedeutet, dass die Matrix Tn transponiert und komplex konjugiert wird. Für eine komplexe Zahl der Form a+ib ist die komplexe Konjugation definiert als a+ib :=a−ib, d.h. das Vorzeichen des Imaginärteils wird umgekehrt. In MATLAB für eine Matrix B: conj(B)zur komplexen Konjugation.
Die Abbildung y7→β wird als diskrete Fourier-Transformation bezeichnet, während ihre Umkehrung β 7→ y = Tnβ Fourier-Synthese genannt wird. Aus der Symmetrie von Tn kann man folgern, dassβ= 1nTnygilt, weshalb man die Fourier-Transformation mit Hilfe der Fourier-Synthese und vice versa berechnen kann. Die Fourier-Synthese entspricht der Auswertung des trigonimetrischen Polynoms p(x) an den äquidistanten Stellenxk. Das Matrix-Vektor-Produkt Tny lässt sich besonders effizient mithilfe einer Divide-and- Conquer-Strategie berechnen, vorausgesetzt die Anzahl der Stützstellen ist eine Zweier- potenz, d.h. n = 2j, j ∈ N (siehe Vorlesung “Enführung in die Numerik”). Dabei geht man nach folgendem Algorithmus vor:
Algorithmus (schnelle Fourier-Transformation).
Input:Funktionsauswertungen y∈Rn zu Stützstellen {2πk/n}n−1k=0, n= 2j, j ∈N Output: Koeffizientenvektorγ ∈Cn
Ist n= 2, so bilde
γ =T2y=
1 1 1 −1
y.
Andernfalls setze n:=n/2und
2
1. berechne die Matrix-Vektor-Produkte
a=Tn
y0 y2 ... y2n−2
und b=Tn
y1 y3 ... y2n−1
wiederum mit der schnellen Fourier-Transformation,
2. berechne c∈Rn für jedes k = 0,1, . . . , n−1als ck :=eπik/nbk, 3. setze
γ0 γ1 ... γn−1
=a+c und
γn γn+1
... γ2n−1
=a−c.
Sobald man den Koeffizientenvektor γ bestimmt hat, ist β = n1γ. Beachte, dass wenn y∈Cn,yzuerst komplex konjugiert werden muss, bevor man die Fourier-Transformation berechnet.
Rekursive Funktionen in MATLAB
Der Algorithmus für die schnelle Fourier-Transformation ruft sich im ersten Schritt (Punkt 1 im Algorithmus oben) zweimal selbst auf, d.h. er ist rekursiv. In diesem Schritt benutzt man einfach wie gewohnt die Funktion, als ob man sie von einem Skript her aufrufen würde. Als Beispiel für eine rekursive Funktion sei die folgende gegeben, welche rekursiv n! berechnet:
function n_fac = factorial(n) if n == 1
n_fac = 1;
else
n_fac = n*factorial(n-1);
end end
Allgemeine Informationen zum Praktikum befinden sich auf der Webseite
http://cm.dmi.unibas.ch/teaching/praktikumII/praktikumII.html
3
