Fourier-Reihen, Teil 9 – komplexe Signale und Kurven in der Ebene

In den bisherigen Teilen haben wir uns mit der Fourier-Analyse reeller Signale beschäftigt. Dabei haben wir rotierende Zeiger unterschiedlicher Frequenzen addiert und die Projektion des Summenzeigers ergab unser zeitabhängiges Signal (s. Teil 1).

Der Summenzeiger hat dabei recht komplizierte Kurven in der komplexen Ebene beschrieben (s. speziell Teil 2). In diesem Teil stellen wir nun die Frage, wie wir geschlossene, ebene Kurven in eine Summe von rotierenden Zeigern verwandeln können.

Einfache Beispiele für solche Kurven sind Lissajous-Figuren wie in Abb. 1 gezeigt. Wir betrachten dabei die Bahnkurve eines Punktes, dessen x– und y-Koordinaten allgemeine Sinus-Funktionen der Zeit t sind. Wenn der Quotient der beiden Frequenzen rational ist, sind die Bahnen geschlossen – und damit periodisch.

Abb. 1: Bahn eines Punktes, dessen x– und y-Koordinaten allgemeine Sinus-Funktionen sind. Speziell ist $f'=0.2\,\text{Hz}$ und $f''=0.4\,\text{Hz}$ , was eine Periodendauer von $T=5\,\text{s}$ bedeutet.

Abb. 1: Bahn eines Punktes, dessen x– und y-Koordinaten allgemeine Sinus-Funktionen sind. Speziell ist $f'=0.2\,\text{Hz}$ und $f''=0.4\,\text{Hz}$ , was eine Periodendauer von $T=5\,\text{s}$ bedeutet.

Lissajous-Figuren

Interpretieren wir die x– bzw. y-Koordinate eines bewegten Punktes als Real- bzw. Imaginärteil eines komplexen Zeigers, erhalten wir das komplexe Signal

$\underline{s}(t)=x(t)+\underline{i}\,y(t)$ ,

das für geschlossene Kurven periodisch ist. Nachdem Real- und Imaginärteil rein reelle und periodische Funktionen sind, können wir deren Fourier-Koeffizienten

$\displaystyle\begin{aligned}\underline{X}_k&=\frac{1}{T}\int_0^Tx(t)\cdot e^{-\underline{i}k\omega_1t}\,\mathrm{d}t\\\underline{Y}_k&=\frac{1}{T}\int_0^Ty(t)\cdot e^{-\underline{i}k\omega_1t}\,\mathrm{d}t\end{aligned}$

berechnen. Aufgrund der Linearität des Integrals gilt weiter

$\displaystyle\begin{aligned}\underline{X}_k+\underline{i}\,\underline{Y}_k&=\frac{1}{T}\int_0^Tx(t)\cdot e^{-\underline{i}k\omega_1t}\,\mathrm{d}t+\underline{i}\,\frac{1}{T}\int_0^Ty(t)\cdot e^{-\underline{i}k\omega_1t}\,\mathrm{d}t\\&=\frac{1}{T}\int_0^T\left(x(t)+\underline{i}\,y(t)\right)\cdot e^{-\underline{i}k\omega_1t}\,\mathrm{d}t\\&=\frac{1}{T}\int_0^T\underline{s}(t)\cdot e^{-\underline{i}k\omega_1t}\,\mathrm{d}t=\underline{S}_k\,.\end{aligned}$

Wir erhalten also direkt die Fourier-Koeffizienten unseres komplexen Signals.

Im Gegensatz zu reellen Signalen gilt jetzt aber im Allgemeinen

$\displaystyle\begin{aligned}\underline{S}_{-k}&=\underline{X}_{-k}+\underline{i}\,\underline{Y}_{-k}\\&=\underline{X}_k^*+\underline{i}\,\underline{Y}_k^*\\&=\left(\underline{X}_k-\underline{i}\,\underline{Y}_k\right)^*\\&\neq\underline{S}_k^*\,.\end{aligned}$

Wir können uns also nicht mehr auf das einseitige Spektrum – sprich komplexe Amplituden – beschränken, sondern müssen negative Frequenzen mitberücksichtigen. Negative Frequenzen bedeuten aber einfach, dass der Zeiger im Uhrzeigersinn rotiert. Mit einer positiven Frequenz rotiert er gegen den Uhrzeigersinn.

In unserem Signal aus Abb. 1 ist das Frequenzverhältnis $f''/f'=2$ ganzzahlig. Deshalb legt die kleinere Frequenz die Grundfrequenz $f_1=f'$ fest und damit auch die Periodendauer $T=1/f_1=5\,\text{s}$ .

Unser Signal ist die Summe zweier allgemeiner Sinus-Funktionen:

$\underline{s}(t) = 2\sin(\tau f'\cdot t) + 2\underline{i}\sin(\tau f''\cdot t-\tfrac{\tau}{8})$ .

Im Spektrum können daher nur die Frequenzen $f'=f_1$ und $f''=2f_1$ vorkommen. D.h., außer für k = ±1 und k = ±2 müssen alle Fourier-Koeffizienten gleich 0 sein. Die exakte Rechnung ergibt

$\displaystyle\begin{aligned}\underline{S}_1&=-\underline{i}\quad &\underline{S}_{-1}&=\underline{i}\\\underline{S}_2&=\frac{1}{\sqrt{2}}-\frac{\underline{i}}{\sqrt{2}}\quad&\underline{S}_{-2}&=-\frac{1}{\sqrt{2}}-\frac{\underline{i}}{\sqrt{2}}\,.\end{aligned}$

Ich lasse die schneller rotierenden Zeiger gerne um die langsameren rotieren. Deshalb sind alle Fourier-Koeffizienten so sortiert wie oben gezeigt. Mit diesen Koeffizienten ergibt sich für unser Signal die Animation in Abb. 2.

Abb. 2: Die Lissajous-Figur aus Abb. 1 gezeichnet als Summe rotierender Zeiger. Während das Signal symmetrisch zur imaginären Achse ist, gilt das für die rotierenden Zeiger nicht.

Abb. 3 zeigt das zweiseitige Spektrum der Fourier-Koeffizienten der Lissajous-Figur aus Abb. 2. In diesem Fall merkt man $\underline{S}_{-k}\neq\underline{S}_k^*$ nur in den Phasen (und da auch nur für k = 2).

Abb. 3: Das zweiseitige Spektrum der Fourier-Koeffizienten der Kurve aus Abb. 1 bzw. Abb. 2.

Die Wikipedia-Seite über Lissajous-Figuren listet noch viele interessante Frequenz-Verhältnisse auf. Für nicht-ganzzahlige, aber rationale Frequenzverhältnisse können wir die Grundfrequenz und damit die Periodendauer wie bei den Schwebungen in Teil 5 bestimmen.

Ptolemäisches Weltbild und Kepler’s Ellipsen

Eine weitere, recht einfache Lissajous-Figur ist eine Ellipse. Ihre Standard-Parameterform ist

$\displaystyle\begin{aligned}\underline{s}(t)&=a\cos(\omega t)+\underline{i}\,b\sin(\omega t)\\&=a\sin(\omega t+\tau/4)+\underline{i}\,b\sin(\omega t)\,.\end{aligned}$

Real- und Imaginärteil haben also dieselbe Frequenz, sind aber $\tau/4=90^\circ$ phasenverschoben. Die Amplituden a und b sind die Längen der Ellipsen-Halbachsen. Die Exzentrizität e der Ellipse ist

$e=\sqrt{a^2-b^2}$ .

Dividieren wir das durch die große Halbachse (bei uns a), erhalten wir die numerische Exzentrizität

$\displaystyle\varepsilon=\sqrt{1-\frac{b^2}{a^2}}$ .

Die Ellipse in Abb. 4 hat die Halbachsen $a=2$ und $b=1$ , und damit die Exzentrizität $e =\sqrt{2^2-1^2}=\sqrt{3}\approx1.732$ bzw. $\varepsilon=\sqrt{1-1^2/2^2}=\sqrt{3}/2\approx0.866$ .

Abb. 4: Eine Ellipse mit großer Halbachse $a=2$ und kleiner Halbachse $b=1$ (in 1. Hauptlage). Die Periodendauer beträgt $T=4\,\text{s}$ , die dunklen Punkte auf der Ellipse haben einen Zeitabstand von $0.2\,\text{s}$ . Der orangene Punkt ist ein Brennpunkt der Ellipse, der den Abstand e vom Ursprung hat.

Abb. 4: Eine Ellipse mit großer Halbachse $a=2$ und kleiner Halbachse $b=1$ (in 1. Hauptlage). Die Periodendauer beträgt $T=4\,\text{s}$ , die dunklen Punkte auf der Ellipse haben einen Zeitabstand von $0.2\,\text{s}$ . Der orangene Punkt ist ein Brennpunkt der Ellipse, der den Abstand e vom Ursprung hat.

Bei einer Periodendauer von $T=4\,\text{s}$ beträgt die Grundfrequenz $f_1=0.25\,\text{Hz}$ . Mit demselben Argument wie oben können nur bei $k=\pm1$ Fourier-Koeffizienten herauskommen, also bei $f=\pm f_1=\pm1/T$ . Die konkrete Rechnung ergibt

$\displaystyle\underline{S}_1=\frac{a+b}{2}=\frac{3}{2}\qquad \underline{S}_{-1}=\frac{a-b}{2}=\frac{1}{2}$ .

Interessanterweise ergeben zwei aufeinander abrollende Kreise eine Ellipse, wenn ihre Radien das Verhältnis

$\displaystyle\frac{r_1}{r_2}=\frac{a+b}{a-b}$

haben.

Das entsprechende Spektrum zeigt Abb. 5. Würde sich der Planet im statt gegen den Uhrzeigersinn drehen, wären $\underline{S}_1$ und $\underline{S}_{-1}$ vertauscht.

Abb 5: Zweiseitiges Spektrum der Fourier-Koeffizienten der Ellipse aus Abb. 4. Die Amplituden bei allen Frequenzen außer bei plus/minus der Grundfrequenz sind gleich 0.

Im ptolemäischen (geozentrischen) Weltbild wurden alle Planetenbahnen mithilfe von Epizykeln (und weiteren Tricks) beschrieben. Für eine Ellipse braucht man nur einen Hauptkreis (Deferent) und einen Epizykel. Wunderbar, wir sind glücklich – oder?

Ein Planet bewegt sich zwar näherungsweise in einer Ellipse um die Sonne, aber nicht so, wie es in Abb. 4 gezeigt ist. Die dunklen Punkte in Abb. 4 zeigen jeweils gleiche Zeitabstände an. In der Nähe der Brennpunkte bewegt sich der Zeiger etwas langsamer als am Ende der kleinen Halbachsen.

Tatsächlich bewegt sich ein Planet in der Nähe der Sonne (eines Brennpunkts) schneller als weit weg von der Sonne (anderer Brennpunkt). Das ist das berühmte 3. Keplersche Gesetz, das auf der Drehimpulserhaltung beruht. Es gilt:

Bei Bahnkurven geht es nicht nur um die Form der Kurve, sondern es ist auch wichtig, wann das Objekt wo ist.

(Das sollte jeder beim Überqueren von Bahngleisen beherzigen!)

Ich habe keine Parametrisierung der Ellipse gefunden, die den Drehimpuls erhalten würde. Deshalb habe ich die entsprechenden Bewegungsgleichungen numerisch gelöst (Details im Anhang unten) und dann die Fourier-Transformation durchgeführt. Das Ergebnis ist die Animation in Abb. 6.

Abb. 6: Eine Kepler-Ellipse ist eine elliptische Bewegung, die den Drehimpuls erhält. Die dunklen Punkte auf der Ellipse haben wieder einen Zeitabstand von $0.2\,\text{s}$ . Der orangene Punkt ist die »Sonne«.

Abb. 6: Eine Kepler-Ellipse ist eine elliptische Bewegung, die den Drehimpuls erhält. Die dunklen Punkte auf der Ellipse haben wieder einen Zeitabstand von $0.2\,\text{s}$ . Der orangene Punkt ist die »Sonne«.

Man sieht an der Bewegung des roten Summenzeigers sehr schön, wie der Planet in der Nähe der Sonne schneller wird und nach seinem »Vorbeiflug« wieder langsamer. Auch die dunklen Punkte, die alle wieder im Zeitabstand 0.2 s gezeigt sind, verdeutlichen das sehr gut.

Wie die Abb. 6 zeigt braucht es jetzt mehr als zwei Kreise. Obwohl die Kreisradien relativ schnell gegen 0 gehen (s. Abb. 7), musste ich bis k = ±55 (±13.75 Hz) summieren, damit der ganz rechte Punkte der Ellipse erreicht wurde. Die Asymmetrie zwischen positiven und negativen Frequenzen ergibt sich wieder durch den gewählten Umlaufsinn.

Abb. 7: Zweiseitiges Spektrum der Fourier-Koeffizienten der Kepler-Ellipse aus Abb. 6.

Unsere Ellipse hat schon eine relativ große Exzentrizität von $\varepsilon\approx0.866$ . Die zu Keplers Zeiten bekannten Planeten haben viel kleinere Exzentrizitäten. Für die Erde hat man z.B. $\varepsilon\approx0.017$ gemessen und für den Mars $\varepsilon\approx0.094$ . Beide Bahnen könnte man mit freiem Auge fast nicht von einem Kreis unterscheiden.

Entsprechend braucht man im ptolemäischen Weltbild nicht ganz so viele Epizykel wie in Abb. 6 gezeigt. Allerdings ist Abb. 6 von einem ruhenden Beobachter außerhalb des Sonnensystems gesehen. Aus Sicht der Erde muss man die Epizykel der Erdbahn und der Planetenbahn kombinieren.

Beliebige geschlossene Kurven

Wie sieht es mit anderen geschlossenen Kurven, also periodischen komplexen Signalen $\underline{s}(t)$ aus?

Eines der bekanntesten Beispiele ist das Video Ptolemy and Homer (Simpson), in dem eine beliebte Zeichentrickfigur mit rotierenden Zeigern gezeichnet wird. Wer dieses Video gesehen hat weiß, dass praktisch keine Einschränkungen existieren – zumindest für zusammenhängende und geschlossene Kurven. Man kann den Stift zwar nicht absetzen, aber problemlos einen Teil der Kurve hin- und wieder zurückfahren.

Als künstlerisch Minderbegabter habe ich erst gar nicht versucht Homer Simpson zu zeichnen; meine Zeichenkünste kamen mit der Blume in Abb. 8 an ihre Grenzen. Gezeichnet (eher programmiert) habe ich sie mit Bézierkurven in Asymptote. Weil die Blume nicht durch eine bekannte Funktion beschrieben werden kann, müssen wir für die Berechnung der Zeiger auf die Diskrete Fourier-Transformation aus Teil 6 zurückgreifen.

Abb. 8: Eine Blume gezeichnet mit einer durchgehenden, geschlossenen Linie (links). 128 Punkte mit gleichem Abstand auf dieser Linie (rechts). Der zeitliche Verlauf der Kurve soll durch die Farben der Punkte angedeutet werden.

Dazu müssen wir die Kurve in eine Menge von Punkten umwandeln, wie es im rechten Teil von Abb. 8 gezeigt ist. Das kann man händisch machen, oder dem Computer überlassen. Wie wir gleich sehen werden, braucht es sehr viele Punkte – die zweite Alternative ist also empfehlenswert.

Nach der Wahl einer Periodendauer T haben die N Punkte jeweils den Zeitabstand $\Delta t=T/N$ . Wenn sie nicht auch denselben Längenabstand haben, läuft der Summenzeiger mal schneller, mal langsamer. Das ist kein prinzipielles Problem und könnte sogar gezielt eingesetzt werden. Die Periodendauer T selber sollte so gewählt werden, dass sich eine angenehme Animation der Zeiger ergibt.

Mit dem folgenden MATLAB/Octave-Code können nun die Zeiger berechnet werden. Hier kommt die Periodendauer T nirgendwo vor; die benötigen wir erst beim tatsächlichen Zeichnen der Zeiger.

% das komplexe Signal
sn = [0.11-1.87i,0.055217-1.70265i, ...
      % 124 Punkte ausgelassen
      -0.055793-1.697763i,-0.000949-1.867389i];

N = length(sn);

% Fourier-Koeffizienten berechnen
Sk = fft(sn)/N;

% nach dem Schema k = 0, 1, -1, 2, -2, 3, -3, ...
% sortieren
kmax = 2*63+1;
sortedSk = complex(zeros(1, kmax));
% MATLAB beginnt bei 1 zu zählen
sortedSk(1) = Sk(1);
for k = 2:2:kmax
    sortedSk(k) = Sk(k/2 + 1);
    sortedSk(k + 1) = Sk(N - k/2 + 1);
end

% für Asymptote aufbereitet speichern
fout = fopen('Sk.txt', 'w');
fprintf(fout, 'pair[] Sk = {(%.8e,%.8e),\n', ...
        real(sortedSk(1)), imag(sortedSk(1)));
for k = 2:2:(kmax - 1)
    fprintf(fout, '             (%.8e,%.8e), ', ...
            real(sortedSk(k)), imag(sortedSk(k)));
    fprintf(fout, '(%.8e,%.8e)', ...
            real(sortedSk(k + 1)), ...
            imag(sortedSk(k + 1)));
    if k ~= kmax - 1
        fprintf(fout, ',\n');
    else
        fprintf(fout, '};\n');
    end
end
fclose(fout);

Wie schon oben erwähnt, hätte ich die Zeiger gerne in der Reihenfolge k = 0, ±1, ±2, … Der fft-Befehl liefert sie aber in der Reihenfolge k = 0, 1, 2, 3, …, -3, -2, -1. Deshalb sortiere ich sie um, bevor sie in die Datei »Sk.txt« geschrieben werden. Sie sind so formatiert, dass Asymptote direkt damit arbeiten kann.

Aus N Messwerten erzeugt die DFT auch N Fourier-Koeffizienten. Abzüglich des konstanten Mittelwerts $\underline{S}_0$ bleiben von den 128 Punkten aus Abb. 8 noch 127 rotierende Zeiger übrig. Ich hätte jedoch gerne immer Paare mit plus/minus derselben Frequenz, wodurch sich die Zahl auf 63 Paare reduziert. Das Ergebnis der Rechnung zeigt die Animation in Abb. 9.

Abb. 9: Berechnet man die Fourier-Transformation der 128 Punkte aus Abb. 8, und zeichnet die Summe bis ±63-mal die Grundfrequenz, ergibt sich die rote Kurve in dieser Animation. Spitzen und Ecken kommen noch nicht gut heraus. Die Periodendauer wurde willkürlich auf 15 s festgesetzt.

Das Ergebnis ist schon ganz anschaulich, speziell die Blattspitzen und andere Ecken kommen jedoch noch nicht ganz heraus. Mehr können unsere 128 Punkte aber nicht liefern.

Um auf Nummer Sicher zu gehen, habe ich den Computer die Koordinaten von 4096 äquidistanten Punkten ermitteln lassen und daraus die DFT berechnet. Mit den ersten ±300 Fourier-Koeffizienten ergibt sich die Animation in Abb. 10. Die rote Kurve zeichnet die grüne fast perfekt nach.

Abb. 10: Berechnet man die Koordinaten von 4096 Punkten und summiert die Fourier-Koeffizienten bis ±300-mal die Grundfrequenz ergibt sich eine fast perfekte Nachzeichnung der grünen Kurve.

Abb. 11 zeigt das zweiseitige Spektrum der Fourier-Koeffizienten der Animation aus Abb. 10. Obwohl die Amplituden ab etwa ±3 Hz (k = ±45) fast 0 sind, wurde in Abb. 10 die Summe bis ±300 gebildet, damit die Spitzen und Ecken halbwegs vernünftig herauskommen.

Abb. 11: Das zweiseitige Spektrum der Fourier-Koeffizienten für die Animation in Abb. 10 bis ±4 Hz. Für eine Periodendauer von 15 s beträgt die Grundfrequenz 1/15 Hz. Es wird also das Spektrum bis k = ±60 gezeigt.

Auch wenn es vielleicht ungewohnt ist, die Blumenkurve und Abb. 11 enthalten exakt dieselbe Information.

Was die automatische Umwandlung von gezeichneten Kurven in Punktmengen betrifft, hat der Nutzer anderstood auf Mathematica StackExchange ein paar Ideen gezeigt. Perfekt ist das Ergebnis aber noch nicht.

Diskussion

Kurven sind 1-dimensional, weil wir entlang ihnen nur vor/zurück gehen können. Dennoch füllen sie eine 2-dimensionale Ebene. Wir werden uns bald mit 2- und mehr-dimensionalen Fourier-Transformationen beschäftigen. (Als Teaser ein Beitrag zu Schwingungen und Wellen im Raum.) Vorher wird es aber voraussichtlich ein paar Posts zu 3D-Computergraphik geben.

Anhang

Um die zeitabhängige Position $\underline{s}(t)=x(t)+\underline{i}\,y(t)$ eines Planeten auf einer Kepler-Ellipse zu erhalten, müssen wir wegen der Gravitation die Bewegungsgleichungen

$\displaystyle\begin{aligned}\ddot{x}&=-GM\cdot x/r^3\\\ddot{y}&=-GM\cdot y/r^3\end{aligned}$

mit $r=\sqrt{x^2+y^2}$ lösen. Der Stern steht dabei im Koordinatenursprung, ist also im Vergleich zu Abb. 4 um die Exzentrizität e nach links geschoben.

Wie man das numerisch recht einfach machen kann, hat Richard Feynman schon relativ zu Beginn seiner Lectures gezeigt. Für eine bessere Genauigkeit habe ich die Integration mit Mathematica gemacht.

Aber wie bekommen wir eine Ellipse mit vorgegebenen Halbachsen a und b und vorgegebener Periodendauer T? Zunächst gilt nach dem 3. Keplerschen Gesetz, dass das Produkt aus Gravitationskonstante G und Sternmasse M gleich

$GM = \left(\frac{\tau}{T}\right)^2\cdot a^3$

ist (wenn der Stern viel schwerer als der Planet ist). Im Prinzip legen wir mit der großen Halbachse und der Periodendauer die Sternmasse fest. Zusätzlich brauchen wir noch passende Anfangsbedingungen, z.B. starten wir in sonnennächster Position (Perihel) senkrecht nach oben. Die Bahngeschwindigkeit in diesem Punkt ist

$\sqrt{GM\cdot\left(\frac{2}{a-e}-\frac{1}{a}\right)}$ ,

weil der Abstand zu Sonne $r=a-e$ beträgt. Insgesamt sind die Anfangsbedingungen:

$\displaystyle\begin{aligned}x(0)&=a-e\quad&\dot{x}(0)&=0\\y(0)&=0&\dot{y}(0)&=\sqrt{GM\cdot\left(\frac{2}{a-e}-\frac{1}{a}\right)}\, .\end{aligned}$

Mathematica ermittelt zunächst für x und y je eine Interpolationsfunktion. Daraus habe ich dann 4096 Punkte in gleichen Zeitabständen berechnen lassen und ihre x-Koordinaten wieder um e nach rechts geschoben.

Mit dem Befehl Fourier[Punkte, FourierParameters -> {-1, 1}] berechnet Mathematica dann die diskrete Fourier-Transformation dieser Punkte. Die optionalen FourierParameters musste ich verwenden, weil Mathematica sonst Hin- und Rücktransformation durch $\sqrt{N}$ dividiert. Wir haben aber immer nur bei der Hintransformation durch N dividiert.

Natürlich sortiert Mathematica die Fourier-Koeffizienten wieder anders als MATLAB/Octave oder wie ich sie gerne hätte. Mathematica hat die Reihenfolge k = 0, -1, -2, -3, …, 3, 2, 1.

Eine Animation des Logos meiner Schule mit als Fourier-Reihe findet sich in Teil 9b.