Creating a Fourier Transform Animation

Alexander Schoch,Tue Jul 16 2024•math python english

If you're reading this, you might have already seen a few examples of these very pretty animations:

In this blog post, we're going to look at how they work and how to write such an animation. Specifically, we're going to

look at what a fourier series and a fourier transform are (go there),
use Inkscape to convert an SVG file into coordinate points we can use (go there),
use python to generate frames for the animation (go there), and
use ffmpeg to combine all frames into one video (go there).

We're not going to look at

how the fourier transform does its magic,
how to use python, and
the math basics (specifically: complex numbers, summation, exponential function, trigonometry).

Therefore, I expect you to

have some basic experience in python, and
know what complex numbers are and how to use them (we are gonna look at complex exponentials, though).

ℹ️

You can find the source code at the end of this page or on my GitLab (for all figures) (opens in a new tab).

The Fourier Transform

💡

If you don't care about the theory, you can jump directly to the practical part. The theory is not important for the code, but it's always nice to understand what you're doing :)

Let's look at an arbitrary function, a simple step function in this case:

f(t) = \begin{cases} 1 &\text{if } t \leq \pi \\ -1 &\text{otherwise} \end{cases}

Step function

The Fourier Series states that any periodic function can be expressed as an infinite sum of specifc sine waves (periodic --> we just loop this function after $2\pi$ - problem solved). Cosine waves are just shifted sine waves: $\cos x = \sin(x + \pi/2)$ . The more functions we add, the better our approximation will generally be.

The following three figures show three fourier series of the step function using 1, 3 and 30 summands. The thin lines depict the waves a fourier series is composed of. It is visible that the fourier series is able to construct the original function really well after a few iterations.(source code (opens in a new tab))

Fourier Series with one summand

Fourier Series with three summands

Fourier Series with thirty summands

The job of the Fourier Transform is to figure out the specific wave functions needed to construct any function as well as possible. This includes three quantities:

The frequency $f$ : How fast does a function oscillate?
The amplitude $a$ : How "large" is a function?
The phase $\phi$ : How much to the left or right is a function shifted?

Sine Wave: Frequency, Amplitude and Phase

Discretization

Contrary to the examples shown above, the input image will not be defined as smooth curves but as a set of points on a 2D plane. Because of this, we'll use a version of the fourier transfrom specifically made for this kind of data, the Discrete Fourier Transform.

The DFT (source: Wikipedia (opens in a new tab)) is defined as

X_k = \sum_{n = 0}^{N-1} x_n \cdot \exp\left\{ -i2\pi\frac{k}{N}n \right\},

where $x_n$ is a point on the input image, $X_k$ is the fourier transform at index $k$ and $N$ is the total number of input points.

This formula shows us a different issue: Our input points will be two-dimensional and contain a $x$ and $y$ coordinate.

From 1D to 2D

For the fourier transfrom, the easiest way to deal with two-dimensional numbers is to use complex numbers (it's not gonna get too difficult, I promise!).

Euler's Formula (opens in a new tab) states that any imaginary exponential function can be split into oscillating real and imaginarpy parts via

e^{it} = \cos(t) + i \sin(t).

Using this property, it is easy to see that the imaginary exponential describes an oscillating point in two dimensions. (source code (opens in a new tab))

Animation of Euler's Formula

If you look at the animation at the top of this page, you'll see many arrows that rotate in a circle with constant speed. This is the same thing as the animation above, but we look at the arrow in the $t$ direction. Imagine standing on the $\text{Real}$ axis label and looking at the arrow from the front. This rotating arrow is our two-dimensional wave function. According to the fourier series, the sum of many of those rotating arrows or wave functions equals our input function almost infinitely well. So, in our final animation, we will omit the time dimension to get these rotating arrows.

Sticking it all together

Similarly to the fourier series examples at the start of this post, we can calculate the frequencies (how fast an arrow spins), its amplitude (how large an arrow is) and its phase (what angle the arrow starts at) by applying the DFT to an input set. The "sum of waves" we saw earlier is extended to 2D by adding the individual arrows (= drawing them after each other, tip to end) together. If we do this for every point in time, we'll get an animation of a "drawing machine"!

Input Data

In this example, I'll use the logo of my university (opens in a new tab) as an input image:

If you want to follow along using my input data, you can download my already discretified input file and continue with Data Preprocessing.

The first step is to convert this SVG file into a series of about 1000 points we can use for further processing (anything between 100 and 5000 will do). A very handy program for this kind of task is Inkscape.

Install Inkscape

Inkscape can be downloaded from their website (opens in a new tab) or via your app store or package manager.

Install the ExportXY Extension

Download ExportXY.inx and ExportXY.py and store them in your inkscape extensions folder. I am writig this post on Fedora Linux, so I can't verify the paths for Windows and Mac.

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

Open input file in Inkscape

Open your file and remove everything that should not be used later on. Your .svg file should now only contain one object that can be drawn as one path.

Convert Object to a path

Select your object and convert it to a path by pressing Shift + Ctrl + C. You know that it worked if you click on the "Node Tool" on the left and seeing the nodes of your path as gray diamonds.

Inkscape Node Tool

Add more Nodes

Now, click on "Extensions" > "Modity Path" > "Add Nodes" and add points to each segment, such that your result contains about 1000 points.

Inkscape Add Nodes Inkscape Add Nodes Dialog

Your path should now look like this (Node Tool):

Inkscape with added nodes

Export

Click on "Extensions" > "Export" > "Export XY" and copy the contents to a new file input.dat.

Inkscape XY Export

Data Preprocessing

First, let's setup the python environment. For that, create a virtual environment, activate it and install the packages numpy, pandas, matplotlib, pillow and alive_progress.

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

Now, create a new python script and open it in any text editor or IDE of your choosing. Import the packages installed before.

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

Next, load the data that we generated in Inkscape. As this data is TAB-separated, we tell pandas to use \t as delimiter and not use the first line as a header for the columns. Afterwards, all input data points are converted to complex numbers (x value = real part, y value = imaginary part). Converting it to a numpy array allows us to grab real and imaginary parts of the whole array with e.g. input_complex.real.

# 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))])

The problem with this raw data is that the coordinate system is still arbitrary to us. The coordinates of each point are defined by the original SVG file, and we have to transform (as in "move and scale", not fourier transform) the input data to fit our future canvas.

Let's first define a canvas size as a complex number (in pixels). My image is is very landscape-y, so I use a wide canvas. Define a top left corner and a bottom right corner, between which the image will be drawn. Here, I define that 75% of my canvas width should be used for the image. Note that the top left corner is negative (here: $-300-150i$ ). The reason for this is that the fourier animation will start at (0, 0), and it looks better if this point is at the center of the image rather than in the top left. center, target_width, and target_height are useful constants we will use later.

# 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

Next, we calculate the extent our original image has by taking minimum and maximum values of the real and imaginary parts. Vertical values are converted back to complex numbers, such that adding e.g. topmost_point will only add a value to the $y$ axis.

Because the input data does not necessarily have the same aspect ratio as the bounding box enclosed with TOP_LEFT and BOTTOM_RIGHT, and we'll fit the width of the image to this bounding box, we calculate a parameter margin_top in order to vertically center the image without stretching along the $y$ axis.

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

Now, we transform the original image to the correct dimensions. First, we subtract leftmost_point and topmost_point, such that the top left corner of the image is at (0,0). Then, we divide by the width of the image and multiply by target_width in order for our image to have a width of, in this case, 600px. We then add TOP_LEFT to horizontally center the image and then add margin_top to also vertically center the image.

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

To test if this worked correctly, let's plot the transformed input. The reason for the inverted minus signs in plt.ylim is that matplotlib's coordinate system has an inverted $y$ axis (positive = up) compared to our input data (positive = down).

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()

Transformed Input Data Plot

Fourier Transform

In order to store the results of our fourier transform, pandas DataFrames come in handy. For that, create a new empty DataFrame with the columns number, frequency, amplitude and phase.

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

Now, we implement the formula stated above,

X_k = \sum_{n = 0}^{N-1} x_n \cdot \exp\left\{ -i2\pi\frac{k}{N}n \right\},

by initializing the sum to zero (= 0 + 0j) and adding $x_n \cdot \exp \{ \ldots \}$ to it. We then divide the sum by N in order to normalize the results by the input data size. This operation is allowed, as it only affects the amplitude, and all amplitudes are multiplied by the same factor. Last, we calculate the frequency as $k$ (remember: the fourier transform assigns an amplitude and a phase to each frequency), the amplitude as the absolute value of $X_k$ and the phase as the rotation of $X_k$ and store it in the DataFrame created above.

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]

Last, for our animation to look nice, we sort the DataFrame by amplitude. The reason for this is that our script will later draw about 1000 circles and lines, and most of them are really, really tiny. It thus makes sense that we sort all arrows by their size (= importance), such that larger arrows are drawn first and closer to the center.

This operation is, of course, allowed, as the sum of multiple functions does not care about the order in which they were added.

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

Drawing

Now, we want to generate a frame for each time step. Because this can take a while, I opted to include a progress bar to display the progress of the rendering. To use it, we have to wrap the whole time loop in

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

We then define a few variables that we need later. path will contain the sums of all previous time steps, such that we can show the path our drawing machine has generated so far. interval is the time step that is added to time after each iteration.

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

⚠️

The time variable needs to range from 0 to exactly $2\pi$ with exactly $N$ steps. If you do it any way differently (e.g. with $N-1$ steps, or not until $2\pi$ ), the animation will break, I promise you! I've been there, and I've spent hours debugging a simple for time in np.linspace(0, 2 * np.pi, N) where N+1 would've been correct.

Rendering out the images happens in a few steps:

Create new Empty image

We create an empty image with the dimensions CANVAS_SIZE and some background color.

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

Calculate/Draw Sum

We draw the sum of all functions at a point in time time. For that, we start at the center of the canvas and then loop over all points in the fourier transform. The end point of each line is calculated as the current sum value plus the function at index $k$ . This is maybe better understood by splitting the exponential into trigonometric functions:

a\exp\left\{ ift + i\phi \right\} = a\cos(ft + \phi) + ia\sin(ft + \phi)

where $a$ is the amplitude, $f$ is the frequency and $\phi$ is the phase. This exponential function is just a wave function, as seen at the beginning of this post, extended to two dimensions.

I then calculate a line thickness as a function of its length, as longer lines should be thicker than shorter ones. We then draw the line and the circle it rotates in, and update the start value for the next line. At the end, we append start (= the total sum) to path.

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)

It is also possible to loop $k$ until some other value ( $\leq N$ ), resulting in a less detailed animation. Here, $k$ goes from 0 to 50, which also looks really cool!

Fourier Animation Frame with lower resolution

Draw Path

Now, draw the path (all calculated sums so far). We stop at the second-to-last element, as each path also uses path[i+1] and thus we only have $M-1$ segments for $M$ points.

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')

Save Image

First, we generate a filename to ensure our operating system sorts the frames correctly. For that, we calculate the number digits the highest frame number will have via

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

This calculates the number of digits we need in order to store all frames using the same filename convention and length, e.g. frames/0000.png, frames/0001.png, ..., frames/1439.png.

We then zero-pad the frame index with this number to generate the filename. Make sure that the frames/ folder exists.

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

Increase Time Step

Increase the time step and push the progress bar to the next iteration.

time += interval
bar()

All in all, the drawing part looks like so:

# 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()

Creating the Video

Last, we can use ffmpeg to combine all frames into one video. -r 60 specifies the frame rate, -c:v libx264 is the video codec. You can, of course, also use the video editor of your choice for this.

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

And we're done!

You should now have a working fourier series animation!

Please feel free to share your finished animation with me! You can send it via email (schochal[at]ethz.ch) or post it on Mastodon and tag me (@alexander_schoch@linuxrocks.online). Also, I'd love to hear your feedback about this post!

Complete Code

You can either copy the code here, or download it directly from GitLab (opens in a new tab).

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()

Comments

Subscribe via RSS.