Erstellen einer Animation mittels Fouriertransformation (German)

Wenn du dies liest, hast du vielleicht schon ein paar Beispiele für diese sehr hübschen Animationen gesehen:

In diesem Blogpost werden wir uns ansehen, wie sie funktionieren und wie man eine solche Animation schreibt. Konkret werden wir

• uns ansehen, was eine Fourier-Reihe und eine Fourier-Transformation sind (Link),
• Inkscape verwenden, um eine SVG-Datei in Koordinatenpunkte umzuwandeln, die wir verwenden können (Link),
• Python verwenden, um Bilder für die Animation zu erzeugen (Link), und
• ffmpeg verwenden, um alle Frames zu einem Video zu kombinieren (Link).

Wir werden uns nicht ansehen,

• wie die Fourier-Transformation ihre Magie macht,
• wie man Python benutzt, und
• wie die benutzten mathematischen Grundlagen funktionieren (spezifisch: komplexe Zahlen, Summation, Exponentialfunktion, Trigonometrie).

Daher gehe ich davon aus, dass du

• einige Grundkenntnisse in Python hast, und
• weisst, was komplexe Zahlen sind und wie man sie benutzt (wir werden uns allerdings komplexe Exponentialfunktionen anschauen).

Die Fouriertransformation

Betrachten wir eine beliebige Funktion, in diesem Fall eine einfache step function:

Die Fourier-Reihe besagt, dass jede periodische Funktion als eine unendliche Summe bestimmter Sinuswellen ausgedrückt werden kann (periodisch —> wir loopen diese Funktion einfach nach $2\pi$. Problem gelöst). Kosinuswellen sind einfach verschobene Sinuswellen: $\cos x = \sin(x + \pi/2)$. Je mehr Funktionen wir hinzufügen, desto besser wird unsere Annäherung im Allgemeinen sein.

Die folgenden drei Abbildungen zeigen drei Fourierreihen der step function mit 1, 3 und 30 Summanden. Die dünnen Linien stellen die Wellen dar, aus denen sich eine Fourierreihe zusammensetzt. Es ist zu erkennen, dass die Fourierreihe die ursprüngliche Funktion nach einigen Iterationen sehr gut beschreiben kann.(Quellcode)

Die Aufgabe der Fourier-Transformation besteht nun darin, die spezifischen Wellenfunktionen zu ermitteln, die benötigt werden, um eine beliebige Funktion so gut wie möglich zu konstruieren. Dazu gehören drei Grössen:

• Die Frequenz $f$: Wie schnell schwingt eine Funktion?
• Die Amplitude $a$: Wie “gross” ist eine Funktion?
• Die Phase $\phi$: Um wie viel nach links oder rechts ist eine Funktion verschoben?

Diskretisierung

Im Gegensatz zu den oben gezeigten Beispielen wird das Input-Bild nicht als glatte Kurven, sondern als eine Reihe von Punkten auf einer 2D-Ebene definiert. Aus diesem Grund werden wir eine Version der Fourier-Transformation verwenden, die speziell für diese Art von Daten entwickelt wurde, die Diskrete Fourier-Transformation.

Die DFT (Quelle: Wikipedia) ist definiert als

wobei $x_n$ ein punkt auf dem Input-Bild, $X_k$ die Fouriertransformation am index $k$ und $N$ die Anzahl Input-Punkte darstellt.

Diese Formel zeigt uns aber bereits ein weiteres Problem: Unser Input wird zweidimensional sein und eine $x$- und eine $y$-Koordinate enthalten.

Von 1D zu 2D

Für die Fourier-Transformation ist der einfachste Weg, mit zweidimensionalen Zahlen umzugehen, die Verwendung komplexer Zahlen (es wird nicht zu schwierig, versprochen!).

Die Eulersche Formel besagt, dass jede imaginäre Exponentialfunktion in oszillierende reelle und imaginäre Teile aufgeteilt werden kann:

Mit dieser Eigenschaft ist sichtbar, dass die komplexe Exponentialfunktion eine Oszillation in zwei Dimensionen beschreibt.(source code)

Wenn du dir die Animation ganz oben auf dieser Seite ansiehst, siehst du viele Pfeile, die sich mit konstanter Geschwindigkeit im Kreis drehen. Dies ist dasselbe wie die Animation oben, aber wir betrachten den Pfeil in $t$-Richtung. Stell dir vor, du stehst auf der $\text{Real}$-Achsenbeschriftung und betrachtest den Pfeil von vorne. Dieser rotierende Pfeil ist unsere zweidimensionale Wellenfunktion. Gemäss der Fourier-Reihe entspricht die Summe vieler dieser rotierenden Pfeile oder Wellenfunktionen fast unendlich gut unserer Input-Funktion. In unserer endgültigen Animation lassen wir also die Zeitdimension weg, um diese rotierenden Pfeile zu erhalten.

Un nun, alles zusammen

Ähnlich wie bei den Fourierreihen-Beispielen zu Beginn dieses Beitrags können wir die Frequenzen (wie schnell sich ein Pfeil dreht), seine Amplitude (wie gross ein Pfeil ist) und seine Phase (bei welchem Winkel der Pfeil beginnt) berechnen, indem wir die DFT auf eine Input-Menge anwenden. Die “Summe der Wellen”, die wir vorhin gesehen hatten, wird auf 2D ausgedehnt, indem die einzelnen Pfeile addiert werden (d. h. sie werden hintereinander gezeichnet, von Spitze zu Ende). Wenn wir dies für jeden Punkt in der Zeit tun, erhalten wir eine Animation einer “Zeichenmaschine”!

Input-Daten

In diesem Beispiel verwende ich das Logo von meiner Universität als Input:

Wenn du meine Eingabedaten verwenden möchtest, kannst du meine bereits diskretisierte Input-Datei herunterladen und mit Vorbereitung der Daten fortfahren.

Der erste Schritt besteht darin, diese SVG-Datei in eine Serie von etwa 1000 Punkten zu konvertieren, die wir für die weitere Verarbeitung verwenden können (alles zwischen 100 und 5000 ist tiptop). Ein sehr praktisches Programm für diese Art von Aufgabe ist Inkscape.

C:\Users\username\AppData\Roaming\Inkscape\Extensions

~/.config/inkscape/extensions/

~/.config/inkscape/extensions/

Vorbereitung der Daten

Zuerst richten wir die Python-Environment ein. Dazu erstellen wir eine virtuelle Umgebung, aktivieren sie und installieren die Packages numpy, pandas, matplotlib, pillow und alive_progress.

python -m venv venv
source venv/bin/activate
 
pip install numpy pandas matplotlib pillow alive_progress

Erstelle nun ein neues Python-Skript und öffne es in einem beliebigen Texteditor oder einer IDE deiner Wahl. Importiere die zuvor installierten Packages.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
 
from PIL import Image, ImageDraw
from alive_progress import alive_bar

Als nächstes laden wir die Daten, die wir in Inkscape erzeugt haben. Da diese Daten TAB-separiert sind, weisen wir Pandas an, \t als Trennzeichen zu verwenden und die erste Zeile nicht als header für die Spalten zu verwenden. Anschliessend werden alle eingegebenen Datenpunkte in komplexe Zahlen umgewandelt (x-Wert = Realteil, y-Wert = Imaginärteil). Die Konvertierung in ein numpy-Array ermöglicht es uns, Real- und Imaginärteil des gesamten Arrays mit z.B. input_complex.real zu erhalten.

# Input Data
input = pd.read_csv('input.dat', delimiter='\t', header=None)
input_complex = np.array([complex(input[0][i], input[1][i]) for i in range(len(input))])

Das Problem mit diesen Rohdaten ist, dass das Koordinatensystem für uns immer noch arbiträr ist. Die Koordinaten der einzelnen Punkte sind durch die ursprüngliche SVG-Datei definiert, und wir müssen die Eingabedaten transformieren (im Sinne von “verschieben und skalieren”, nicht Fourier-Transformation), damit sie in unseren zukünftigen canvas passen.

Definieren wir zunächst die Grösse des canvas als komplexe Zahl (in Pixeln). Mein Bild ist sehr querformatig, daher verwende ich einen breiten canvas. Definiere eine linke obere Ecke und eine rechte untere Ecke, zwischen denen das Bild gezeichnet werden soll. Hier lege ich fest, dass 75% der Canvasbreite für das Bild verwendet werden sollen. Beachte, dass die obere linke Ecke negativ ist (hier: $-300-150i$). Der Grund dafür ist, dass die Fourier-Animation bei (0, 0) beginnt, und es sieht besser aus, wenn dieser Punkt in der Mitte des Bildes liegt als oben links. center, target_width und target_height sind nützliche Konstanten, die wir später verwenden werden.

# Constants
CANVAS_SIZE = 800 + 400j
 
TOP_LEFT = -CANVAS_SIZE / 8 * 3
BOTTOM_RIGHT = CANVAS_SIZE / 8 * 3
 
center = (TOP_LEFT + BOTTOM_RIGHT) / 2
target_width = BOTTOM_RIGHT.real - TOP_LEFT.real
target_height = BOTTOM_RIGHT.imag - TOP_LEFT.imag

Als Nächstes berechnen wir die Ausdehnung unseres Originalbildes, indem wir die Minimal- und Maximalwerte des Real- und Imaginärteils nehmen. Vertikale Werte werden wieder in komplexe Zahlen umgewandelt, sodass das Hinzufügen von z.B. topmost_point nur einen Wert auf der $y$-Achse ergibt.

Da die Eingabedaten nicht notwendigerweise das gleiche Seitenverhältnis haben wie die mit TOP_LEFT und BOTTOM_RIGHT eingeschlossene Bounding Box, und wir die Breite des Bildes an diese Bounding Box anpassen werden, berechnen wir einen Parameter margin_top, um das Bild vertikal zu zentrieren, ohne es entlang der $y$-Achse zu strecken.

leftmost_point = min(input_complex.real)
rightmost_point = max(input_complex.real)
topmost_point = complex(0, min(input_complex.imag))
bottommost_point = complex(0, max(input_complex.imag))
 
width = rightmost_point - leftmost_point
height = bottommost_point - topmost_point
margin_top = (target_height - height.imag / width * target_width) / 2

Jetzt transformieren wir das Originalbild in die richtigen Dimensionen. Zuerst subtrahieren wir leftmost_point und topmost_point, sodass die obere linke Ecke des Bildes bei (0,0) liegt. Dann dividieren wir durch die width des Bildes und multiplizieren mit der target_width, sodass unser Bild eine Breite von 600 Pixeln hat. Dann fügen wir TOP_LEFT hinzu, um das Bild horizontal zu zentrieren, und addieren dann margin_top, um das Bild auch vertikal zu zentrieren.

input_transformed = (input_complex - leftmost_point - topmost_point) / width * target_width + TOP_LEFT + complex(0, margin_top)

Um zu testen, ob dies korrekt funktioniert hat, plotten wir den transformierten Input. Der Grund für die umgekehrten Minuszeichen in plt.ylim ist, dass das Koordinatensystem von matplotlib eine umgekehrte $y$-Achse (positiv = oben) im Vergleich zu unseren Eingabedaten (positiv = unten) hat.

plt.plot(input_transformed.real, input_transformed.imag, color='black')
plt.xlim(-CANVAS_SIZE.real / 2, CANVAS_SIZE.real / 2)
plt.ylim(CANVAS_SIZE.imag / 2, -CANVAS_SIZE.imag / 2)
plt.show()

Fourier-Transformation

Um die Resultate unserer Fourier-Transformation zu speichern, bieten sich pandas DataFrames an. Erstelle dazu einen neuen leeren DataFrame mit den Spalten number, frequency, amplitude und phase.

N = len(input_transformed)
fourier_transform = pd.DataFrame(columns = ['number', 'frequency', 'amplitude', 'phase'], index = np.arange(N))

Jetzt implementieren wir die Formel von oben,

indem wir die Summe mit Null initialisieren (= 0 + 0j) und $x_n \cdot \exp \{ \ldots \}$ dazuaddieren. Dann teilen wir die Summe durch N, um die Ergebnisse mit der Grösse der Input-Daten zu normalisieren. Dieser Vorgang ist erlaubt, da er sich nur auf die Amplitude auswirkt und alle Amplituden mit demselben Faktor multipliziert werden. Zuletzt berechnen wir die Frequenz als $k$ (zur Erinnerung: die Fourier-Transformation ordnet jeder Frequenz eine Amplitude und eine Phase zu), die Amplitude als Absolutwert von $X_k$ und die Phase als Drehung von $X_k$ und speichern sie in dem oben erstellten DataFrame.

for k in range(N):
    X_k = 0 + 0j
    for n in range(N):
        X_k += input_transformed[n] * np.exp(-1j * 2 * np.pi * k * n / N)
    X_k /= N
    frequency = k
    amplitude = abs(X_k)
    phase = np.angle(X_k)
    fourier_transform.iloc[k] = [X_k, frequency, amplitude, phase]

Damit unsere Animation schön aussieht, sortieren wir den DataFrame nach Amplitude. Der Grund dafür ist, dass unser Skript später etwa 1000 Kreise und Linien zeichnen wird, und die meisten davon sind sehr, sehr klein. Es macht also Sinn, dass wir alle Pfeile nach ihrer Grösse (= Wichtigkeit) sortieren, sodass grössere Pfeile zuerst gezeichnet werden und näher an der Mitte liegen.

Diese Operation ist natürlich erlaubt, da die Summe mehrerer Funktionen unabhängig von der Reihenfolge der Summanden ist.

fourier_transform.sort_values(by = 'amplitude', inplace=True, ascending = False)

Zeichnung

Nun wollen wir für jeden Zeitschritt ein Bild erzeugen. Da dies eine Weile dauern kann, habe ich mich dafür entschieden, einen Fortschrittsbalken einzubauen, der den Fortschritt des Render-Prozesses anzeigt. Um diesen zu verwenden, müssen wir die gesamte Zeitschleife folgenden Code packen:

with alive_bar(N) as bar:
    [...]

Dann definieren wir ein paar Variablen, die wir später brauchen. path wird die Summe aller vorherigen Zeitschritte enthalten, sodass wir den Pfad anzeigen können, den unsere Zeichenmaschine bis jetzt erzeugt hat. interval ist der Zeitschritt, der nach jeder Iteration zu time hinzugefügt wird.

path = [] # stores the points of the already calculated sums
 
interval = 2 * np.pi / N
time = 0

Das Rendern der Bilder erfolgt in mehreren Schritten:

image = Image.new("RGB", (int(CANVAS_SIZE.real), int(CANVAS_SIZE.imag)), 'white')
draw = ImageDraw.Draw(image)

start = CANVAS_SIZE / 2
for k in range(N):
    point = fourier_transform.iloc[k]
    end = start + point.amplitude * np.exp(1j * point.frequency * time + 1j * point.phase)
 
    thickness = round(5 * point.amplitude / max(fourier_transform.amplitude))
    draw.line([start.real, start.imag, end.real, end.imag], width=thickness, fill='black')
    draw.circle((start.real, start.imag), point.amplitude, outline='gray', width=1)
 
    start = end
path.append(start)

for i in range(len(path)-1):
    draw.line([path[i].real, path[i].imag, path[i+1].real, path[i+1].imag], width=5, fill='black')

num_digits = int(np.log10(N)) + 1 # For filename generation

filename = "frames/%s.png" % str(i).zfill(num_digits)
image.save(filename)

time += interval
bar()

Alles in allem sieht der Zeichnungsteil so aus:

# Animation
num_digits = int(np.log10(N)) + 1 # For filename generation
with alive_bar(N) as bar:
    for i in range(N):
        image = Image.new("RGB", (int(CANVAS_SIZE.real), int(CANVAS_SIZE.imag)), 'white')
        draw = ImageDraw.Draw(image)
 
        # path.append(draw_epicycles(time, draw))
        start = CANVAS_SIZE / 2
        for k in range(N):
            point = fourier_transform.iloc[k]
            end = start + point.amplitude * np.exp(1j * point.frequency * time + 1j * point.phase)
            thickness = round(5 * point.amplitude / max(fourier_transform.amplitude))
            draw.line([start.real, start.imag, end.real, end.imag], width=thickness, fill='black')
            draw.circle((start.real, start.imag), point.amplitude, outline='gray', width=1)
            start = end
        path.append(start)
 
        for i in range(len(path)-1):
            draw.line([path[i].real, path[i].imag, path[i+1].real, path[i+1].imag], width=5, fill='black')
 
        time += interval
 
        filename = "frames/%s.png" % str(i).zfill(num_digits)
        image.save(filename)
        bar()

Erstellen des Videos

Zuletzt können wir ffmpeg verwenden, um alle Bilder zu einem Video zu kombinieren. -r 60 gibt die Framerate an, -c:v libx264 ist der Videocodec. Du kannst dafür natürlich auch den Videoeditor deiner Wahl verwenden.

ffmpeg -r 60 -pattern_type glob -i 'frames/*.png' -c:v libx264 fourier_transform_animation.mp4

Und… Fertig!

Damit solltest du nun eine funktionierende Fourier-Serien-Animation haben!

Schick mir gerne deine fertige Animation! Du kannst mir dafür eine Mail zusenden (schochal[at]ethz.ch) oder sie auf Mastodon posten und mich markieren (@alexander_schoch@linuxrocks.online). Ich würde mich ausserdem sehr über Feedback zu diesem Post freuen!

Kompletter Code

Du kannst den Code entweder hier kopieren oder ihn direkt von GitLab herunterladen.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
 
from PIL import Image, ImageDraw
from alive_progress import alive_bar
 
# Input Data
input = pd.read_csv('input.dat', delimiter='\t', header=None)
input_complex = np.array([complex(input[0][i], input[1][i]) for i in range(len(input))])
 
# Constants
CANVAS_SIZE = 800 + 400j
 
TOP_LEFT = -CANVAS_SIZE / 8 * 3
BOTTOM_RIGHT = CANVAS_SIZE / 8 * 3
 
center = (TOP_LEFT + BOTTOM_RIGHT) / 2
target_width = BOTTOM_RIGHT.real - TOP_LEFT.real
target_height = BOTTOM_RIGHT.imag - TOP_LEFT.imag
 
# Transforming (moving / scaling) input to desired dimensions
leftmost_point = min(input_complex.real)
rightmost_point = max(input_complex.real)
topmost_point = complex(0, min(input_complex.imag))
bottommost_point = complex(0, max(input_complex.imag))
 
width = rightmost_point - leftmost_point
height = bottommost_point - topmost_point
margin_top = (target_height - height.imag / width * target_width) / 2
 
input_transformed = (input_complex - leftmost_point - topmost_point) / width * target_width + TOP_LEFT + complex(0, margin_top)
 
# Discrete Fourier Transform
N = len(input_transformed)
fourier_transform = pd.DataFrame(columns = ['number', 'frequency', 'amplitude', 'phase'], index = np.arange(N))
 
for k in range(N):
    X_k = 0 + 0j
    for n in range(N):
        X_k += input_transformed[n] * np.exp(-1j * 2 * np.pi * k * n / N)
    X_k /= N
    frequency = k
    amplitude = abs(X_k)
    phase = np.angle(X_k)
    fourier_transform.iloc[k] = [X_k, frequency, amplitude, phase]
 
fourier_transform.sort_values(by = 'amplitude', inplace=True, ascending = False)
 
path = [] # stores the points of the already calculated sums
 
num_digits = int(np.log10(N)) + 1 # For filename generation
 
interval = 2 * np.pi / N
time = 0
 
# Animation
with alive_bar(N) as bar:
    for i in range(N):
        image = Image.new("RGB", (int(CANVAS_SIZE.real), int(CANVAS_SIZE.imag)), 'white')
        draw = ImageDraw.Draw(image)
 
        start = CANVAS_SIZE / 2
        for k in range(N):
            point = fourier_transform.iloc[k]
            end = start + point.amplitude * np.exp(1j * point.frequency * time + 1j * point.phase)
            thickness = round(5 * point.amplitude / max(fourier_transform.amplitude))
            draw.line([start.real, start.imag, end.real, end.imag], width=thickness, fill='black')
            draw.circle((start.real, start.imag), point.amplitude, outline='gray', width=1)
            start = end
        path.append(start)
 
        for i in range(len(path)-1):
            draw.line([path[i].real, path[i].imag, path[i+1].real, path[i+1].imag], width=5, fill='black')
 
        time += interval
 
        filename = "frames/%s.png" % str(i).zfill(num_digits)
        image.save(filename)
        bar()

Subscribe via RSS.

MIT 2024 © Alexander Schoch.