5.6. Trigonometrische Interpolation

Wir betrachten in diesem Abschnitt noch einen Spezialfall der Polynominterpolation, nämlich die Interpolation auf äquidistanten Stützstellen auf dem komplexen Einheitskreis. Dieses sehr spezielle Interpolationsproblem führt uns zu einem der wichtigsten Hilfsmittel der angewandten Mathematik - der diskreten Fouriertransformation.

Da wir im Gegensatz zu den vorigen Abschnitten nun mit Polynomen über dem Körper \(\C\) arbeiten wollen, bezeichnen wir ab nun mit \(\Pi_n\) die Menge der Polynome mit komplexen Koeffizienten vom Grad \(n \in \N\). Ein Polynom \(P \in \Pi_n\) mit \(P \colon \C \rightarrow \C\) ist nun von der Gestalt

\[P(x) \ = \ \sum_{k=0}^n a_kx^k, \qquad a_k \in \C.\]

Es sei angemerkt, dass sich die Aussage zur Eindeutigkeit des Interpolationspolynoms aus Satz 5.1 direkt vom Körper \(\R\) auf \(\C\) verallgemeinern lässt.

Wir wählen im Folgenden \(n \in \N\) Stützstellen \(x_j \in \mathbb{S}^1, j=0,\ldots,n-1\) auf dem komplexen Einheitskreis, die gegeben sind durch

\[x_j \ = \ e^{it_j} \ = \ \cos t_j + i \sin t_j, \qquad t_j \ = \ \frac{2 \pi j}n, \qquad j=0,\ldots,n-1 .\]

Darüberhinaus nehmen wir an, dass wir zu den \(n\) Stützstellen oben ebenfalls \(n\) Stützwerte \(f_j, j=0,\ldots,n-1\) gegeben haben.

Da die Stützstellen \(x_j, j=0,\ldots,n-1\) paarweise unterschiedlich sind existiert also ein eindeutiges Interpolationspolynom \(P \in \Pi_{n-1}\) dessen unbekannte Koeffizienten \(a_k\) wir im Folgenden mit \(\hat{f}_k\) bezeichnen wollen und das die Form hat

\[P(x) \ = \ \sum_{k=0}^{n-1} \hat{f}_k x^k.\]

Da es sich um interpolierendes Polynom handelt muss also folgende Bedingung gelten

\[P(x_j) \ = \ \sum_{k=0}^{n-1} \hat f_k x_j^k \ = \ \sum_{k=0}^{n-1} \hat f_k e^{ikt_j} \ = \ \sum_{k=0}^{n-1} \hat f_k e^{2 \pi i k j/n} \ \overset{!}{=} \ f_j, \qquad j=0,\ldots,n-1.\]

In diesem Fall können wir eine sehr einfache Lösungsformel für die unbekannten Koeffizienten \(\hat{f}_k\) angeben, wie das folgende Theorem zeigt.

Satz 5.5 (Diskrete Fouriertransformation)

Es seien \(n\) Datenpunkte \((x_j, f_j) \in \C^2\) gegeben, wobei die Stützstellen \(x_j = e^{2\pi i j / n}\) äquidistant auf dem komplexen Einheitskreis liegen. Sei \(P(x) \ \coloneqq \ \sum_{k=0}^{n-1} \hat{f}_k x^k\) das eindeutig bestimmte Interpolationspolynom.

Dann lassen sich die unbekannten Koeffizienten \(\hat f_k, k=0,\ldots,n-1\) von \(P\) wie folgt ausrechnen

(5.14)\[ \hat f_k \ = \ \frac{1}n \sum_{j=0}^{n-1} f_j e^{-2 \pi i k j/n}.\]

Diese Darstellung der Koeffizienten nennt man die diskrete Fourier-Transformation der Länge \(n\) zu den Daten \(f_j, j=0,\ldots,n-1\).

Proof. Wir zeigen die Identität des Lemmas indem wir die diskrete Fourier-Transformation der Daten in das Interpolationspolynom einsetzen. Sei also \(x_j \in \mathbb{S}^1\) eine beliebige Stützstelle auf dem komplexen Einheitskreis. Dann gilt für das eindeutig bestimmte Interpolationspolynom \(P\)

\[\begin{split} P(x_j) \ &= \ \sum_{k=0}^{n-1} \hat f_k e^{2 \pi i k j/n} \ = \ \sum_{k=0}^{n-1} \frac{1}n \sum_{\ell=0}^{n-1} f_\ell e^{-2 \pi i k \ell/n} e^{2 \pi i k j/n} \\ & = \ \frac{1}n \sum_{k=0}^{n-1} \sum_{\ell=0}^{n-1} f_\ell e^{2 \pi i (j - \ell) k/n} \ = \ \frac{1}n \sum_{\ell=0}^{n-1} f_\ell \sum_{k=0}^{n-1} \left( e^{2 \pi i (j - \ell)/n}\right)^k \end{split}\]

Die innere Summe kann als Partialsumme einer geometrischen Reihe \(s_{n-1}(j, \ell) \coloneqq \sum_{k=0}^{n-1} q^k(j, \ell)\) für \(q(j, \ell) \coloneqq e^{2 \pi i (j - \ell)/n}\) interpretiert und berechnet werden, so dass wir erhalten

\[s_{n-1}(j, \ell) \ = \ \frac{q^n(j, \ell) - 1}{q(j, \ell) - 1} \ = \ \frac{e^{2 \pi i (j - \ell)} - 1}{e^{2 \pi i (j - \ell)/n} - 1} \ = \ \begin{cases} \ 0, \quad &\text{ falls } j \neq \ell, \\ \ n, \quad &\text{ falls } j = \ell. \end{cases}\]

Damit berechnen wir also insgesamt

\[P(x_j) \ = \ \frac{1}n \sum_{\ell=0}^{n-1} f_\ell s_{n-1}(j, \ell) \ = \ f_j.\]

Da dies für jede beliebige Stützstelle \(x_j, j=0,\ldots,n-1\) gilt, haben wir gezeigt, dass die Koeffizienten \(\hat f_k\) das eindeutige Interpolationspolynom ergeben. ◻

Korollar 5.4 (Orthogonalitätseigenschaft)

Im Beweis von Satz 5.5 können wir sehr schön eine Orthogonalitätseigenschaft der trigonometrischen Funktionen ablesen, denn es gilt

\[\frac{1}{n} \sum_{j=0}^{n-1} e^{2\pi jk / n} \ = \ \begin{cases} \ 1, \quad &\text{ falls } k \in n\mathbb{Z}, \\ \ 0, \quad &\text{ sonst }.\\ \end{cases}\]

Die unbekannten Koeffizienten \(\hat f_k\) in der diskreten Fouriertransformation werden auch Fourierkoeffizienten zu den Daten \(f_j, j=0,\ldots,n-1\) genannt. Aus dem Beweis von Satz 5.5 wird klar, dass man die Berechnungen auch umkehren kann, d.h., dass man aus einer Menge von Fourierkoeffizienten \(\hat f_k, k=0,\ldots,n-1\) die ursprünglichen Stützwerde \(f_j, j=0,\ldots,n-1\) ausrechnen kann. Dies lässt sich nämlich wie folgt realisieren

(5.15)\[f_j \ = \ \sum_{k=0}^{n-1} \hat f_k e^{2 \pi i j k/n}.\]

Diese Identität nennt man daher auch die inverse diskrete Fouriertransformation der Länge \(n\) zu den Fourierkoeffizienten \(\hat f_k, k=0,\ldots, n-1\). Man beachte, dass je nach wissenschaftlicher Disziplin die Definition der diskreten Fouriertransformation und ihrer Inversen leicht abweichen kann, wie z.B. die Vertauschung der Vorzeichens im Exponenten der Stützstellen oder die Normalisierung mit \(1/n\).

Bemerkung 5.10 (Rechenaufwand der diskreten Fourier-Transformation)

Der Rechenaufwand zur Berechnung der diskreten Fourier-Transformation in (5.14) und ihrer Umkehrabbildung in (5.15) benötigt bei Vorberechnung der \(n\) Koeffizienten \(t_j = \frac{2\pi j}{n}\) genau \(n^2\) FLOPS und somit liegt der Rechenaufwand in \(\mathcal{O}(n^2)\).

5.6.1. Schnelle Fouriertransformation

Die Berechnung der Fourierkoeffizienten für gegebene Daten erlaubt es diese in ein Spektrum aus Frequenzen zu zerlegen und darzustellen. Anders ausgedrückt kann die Fouriertransformation als Zerlegung eines Signals in seine verschiedenen Schwingungsanteile interpretiert werden und hat deshalb vielfache Anwendungen. Daher ist es naheliegend, dass die Fouriertransformation besonders in der Signalverarbeitung eine zentrale Rolle spielt. Ein Großteil der modernen Kompressionsalgorithmen, wie zum Beispiel MP3 oder JPEG verwenden die diskrete Fouriertransformation oder eine ihrer Varianten (z.B. die diskrete Kosinustransformation), um Daten effizient in ihrer Größe zu reduzieren ohne große Qualitätseinbußen hinnehmen zu müssen.

Durch die vielseitige Anwendbarkeit der diskreten Fouriertransformation ist es besonders wichtig, diese so effizient wie möglich zu berechnen. Eine formelle Beschreibung eines effizienten Algorithmus, der sogenannten schnellen Fouriertransformation (im Englischen: Fast Fourier Transform (FFT)) basiert auf der Pionierarbeit von Cooley-Tukey aus dem Jahr 1965 und basiert auf Ideen von Gauss. Die zentrale Idee des Algorithmus ist es die diskrete Fouriertransformation der Größe \(n \in \mathbb{N}\) in zwei diskrete Fouriertransformationen der Größe \(n/2\) aufzuteilen und diese Zerlegung rekursiv fortzusetzen.

Sei im Folgenden \(n=2m\) eine gerade natürliche Zahl. Dann können wir die diskrete Fouriertransformation der Länge \(n\)

\[\hat f_k \ = \ \frac{1}{n}\sum_{j=0}^{n-1} f_{j} q^{jk}, \qquad \text{ mit } \qquad q \ \coloneqq \ e^{-2\pi i / n},\]

zur Berechnung eines Fourierkoeffizienten \(f_k \in \C\) in gerade und ungerade Indizes wie folgt aufteilen:

\[\begin{aligned} n \hat f_k \ &= \ \sum_{j=0}^{m-1} f_{2j} q^{(2j)k} + \sum_{j=0}^{m-1} f_{2j+1} q^{(2j+1)k} \\ &= \ \sum_{j=0}^{m-1} f_{2j} (q^2)^{jk} + q^k \sum_{j=0}^{m-1} f_{2j+1} (q^2)^{jk}\\ &=: \ \hat g_k + q^k \hat u_k. \end{aligned}\]

Man sieht also, das wir die ursprüngliche diskrete Fouriertransformation \(\hat f_k, k=0\,ldots, n-1\) der Länge \(n\) durch zwei diskrete Fouriertransformationen \(\hat g_k, k=0,\ldots,m-1\) und \(\hat u_k, k=0,\ldots,m-1\) der Länge \(m = n/2\) berechnen können. Dies sieht man ein, da gilt

\[q^2 \ = \ e^{(-2\pi i / n)^2} \ = \ e^{-4\pi i / n } \ = \ e^{-2\pi i / (n / 2) } \ = \ e^{-2\pi i / m}.\]

Wie man nachrechnen kann gilt außerdem

\[\begin{split} \hat g_{k+m} \ &= \ \sum_{j=0}^{m-1} f_{2j} (q^2)^{j(k+m)} \ = \ \sum_{j=0}^{m-1} f_{2j} e^{-2\pi i j(k+m) / m} \\ &= \ \sum_{j=0}^{m-1} f_{2j} e^{-2\pi i jk / m}\cdot \underbrace{e^{-2\pi i j}}_{= \, 1} \ = \ \hat g_k \end{split}\]

und analog \(\hat u_{k+m} = \hat u_k\). Wegen \(m = n/2\) gilt außerdem

\[q^{k + m} \ = \ e^{(-2\pi i / n)^{(k+m)}} \ = \ e^{-2\pi i k / n} \cdot \underbrace{e^{-2\pi i (n/2) / n}}_{= \, -1} \ = \ -q^k.\]

Vernachlässigen wir den Normalisierungsfaktor \(n\) vor den Fourierkoeffizienten, so lassen sich alle ursprünglichen Fourierkoeffizienten \(\hat{f}_l=0,\ldots,n-1\) wie folgt über die Fouriertransformationen \(g_k\) und \(u_k\) für \(k=0,\ldots, m-1\) berechnen:

(5.16)\[\begin{split} \hat f_k \ = \ &\hat g_k + q^k \hat u_k,\\ \hat f_{k+m} \ = \hat g_{k+m} + q^{k+m} \hat u_{k+m} \ = \ &\hat g_k - q^k \hat u_k. \end{split}\]

Mit Hilfe dieser typischen Divide & Conquer Strategie können wir den Rechenaufwand für die schnelle Fouriertransformation bestimmen.

Lemma 5.1 (Rechenaufwand der schnellen Fouriertransformation)

Der Rechenaufwand der schnellen Fouriertransformation der Länge \(n = 2^p, p \in \N\) liegt in \(\mathcal{O}(n \log n)\).

Proof. Sei \(M_n = M_{2^p}\) der gesamte Rechenaufwand zur Berechnung der diskreten Fouriertransformation der Länge \(n = 2^p\). Wie bemerken zunächst, dass wir alle \(n\) Terme \(q^k, k=0,\ldots,n-1\) einmalig als fixe Einheitswurzeln vorberechnen können mit einem Rechenaufwand von \(\mathcal{O}(n)\).

Dann gilt mit der Zerlegung aus der schnellen Fouriertransformation in (5.16)

\[M_n \ = \ 2 M_{n/2} + n + n/2 \ = \ 2 M_m + 3m,\]

da wir \(n\) zusätzliche Additionen und \(m\) Multiplikationen durchführen müssen um alle Fourierkoeffizienten zu berechnen. Da \(n=2^p\) als eine Zweierpotenz gewählt ist, können wir mit \(m = 2^{p-1}\) rekursiv folgern

(5.17)\[\begin{split} M_{2^p} \ &= \ 2 M_{2^{p-1}} + 3 \cdot 2^{p-1} \ = \ 2( 2 M_{2^{p-2}} + 3\cdot 2^{p-2} )+3\cdot 2^{p-1} \\ &= \ 4 M_{2^{p-2}} + 2\cdot 3 \cdot 2^{p-1} \ \overset{\cdots}{=} \ 2^p M_1 + 3\cdot p \cdot 2^{p-1}. \end{split}\]

Da für die Fouriertransformation der Länge \(1\) gilt, dass man die Daten selbst erhält, d.h.

\[\hat{f_j} \ = \ \frac{1}{1} \sum_{j=0}^0 f_j q^{j\cdot 0} \ = \ f_j,\]

ist der Rechenaufwand also \(M_1=1\). Und somit erhalten wir aus (5.17), dass \(M_{2^p} = {\cal O}(p \cdot 2^p)\) bzw. \(M_n = {\cal O}(n \log n)\) gilt. ◻