Ever wondered how computers process and analyze images and shapes? You might have seen images stored as .png, .jpg, .bmp, or .svg while working with graphics or image processing. Broadly, image files can be represented using two fundamental formats: raster graphics and vector graphics.

• Raster graphics store images as a grid of pixels, where each pixel contains color or intensity information.
• Vector graphics represent shapes using mathematical equations — lines, curves, and polygons — allowing them to scale infinitely without losing quality.

In mathematical terms, an image in a two-dimensional space can be represented in two ways: spatial domain and frequency domain. In the spatial domain, the image is represented as a grid of pixels, where each pixel has a specific intensity value. In the frequency domain, the image is represented as a combination of different frequencies, which can be analyzed using Fourier analysis.

Fourier analysis is a mathematical technique used to decompose periodic functions into their constituent frequencies and analyze the amplitude and phase of each frequency component. Some of the key applications of Fourier analysis include audio processing, audio filtering, noise reduction, signal reconstruction, image compression, feature extraction, and recognition, object detection, and speech recognition among others. This means AirPods in your ears use some form of Fourier analysis to filter out the background noise when you activate noise cancellation in a subway or on a plane.

At its core, Fourier analysis is about decomposing a time-domain signal into its constituent frequencies, and inferring the amplitude and phase of each frequency component.

Draw a shape on the left canvas and watch Fourier epicycles trace it on the right! ^[1]

Left: Draw your shape | Right: Fourier epicycles tracing it

Table of Contents

1. Recap
   • Sinusoidal Waves
   • A Harmonic Wave
   • Functions of Waves and Signals
2. Fourier Series
   • Periodic Signals
   • Approximating Complex Waves with Fourier Series
3. Drawing with Fourier Series: Epicycles
   • Sketches as Functions
   • Fourier Approximation of the Bat Silhouette
   • Epicycles: Rotating Vectors
   • Complex Exponential Notation
   • Computing Epicycles with the DFT
   • Retracing our drawing with Epicycles
   • The Mathematics Behind Epicycles
   • Why This Works
   • Generalizing: From Periodic to Non-Periodic Signals
4. Fast Fourier Transform (FFT)
   • The Performance Problem
   • Enter the FFT
   • How FFT Works (Intuition)
   • Why This Matters
   • Libraries and Implementation

Fourier Analysis

Instead of jumping straight into the math, I’ll try to build an intuitive understanding of how waves and signals work which would make it easier to grasp the need for Fourier analysis.

---

Here’s a quick overview recap on how waves and signals work, followed by an introduction to Fourier series, Fourier transform, and discrete Fourier transform.

Recap

A wave is a disturbance that travels through a medium, and a signal is a wave that carries information.

Sinusoidal Waves

The most fundamental wave pattern in nature is the sinusoidal wave (or sine wave). A sine wave is a smooth, continuous oscillation that repeats periodically. It’s described mathematically by the sine function, which produces a characteristic S-shaped curve when plotted.

The basic sine wave equation is:

$y = \sin(x)$

where $y$ oscillates smoothly between -1 and +1 as $x$ increases.

This simple mathematical function appears everywhere in nature—from the oscillation of a pendulum to the vibration of guitar strings, from alternating current in electrical circuits to the propagation of light and sound waves.

Sinusoidal waves are special because they are the building blocks of all other periodic waves. No matter how complex a wave pattern looks, it can be broken down into a sum of simple sine (and cosine) waves—this is the fundamental insight of Fourier analysis.

A Harmonic Wave:

The equation $f(t) = A \sin(\omega t)$ represents a simple harmonic wave or oscillation with the following parameters:

• Amplitude (A): The maximum displacement from the equilibrium position. It determines how “tall” the wave peaks are. In our example, $A = 40$ units.
• Angular Frequency (ω): The rate of oscillation in radians per second. It is related to the frequency $f$ (in Hz) by the formula $\omega = 2\pi f$. In our example, if $f = 2$ Hz, then $\omega = 4\pi$ radians per second.
• Time (t): The independent variable, usually representing time in seconds.

Suppose we have a wave with frequency $f = 2$ Hz; the wave oscillates 2 times per second,

1. $T = \frac{1}{f} = \frac{1}{2} = 0.5$ seconds. (One cycle is completed in 0.5 seconds)
2. $\omega = 2\pi f = 2\pi \times 2 = 4\pi$ radians per second. (2 cycles are completed in 1 second)

The example equation $f(t) = 40 \sin(4\pi t)$ describes a wave with an amplitude of 40 units and a frequency of 2 Hz. This means the wave oscillates between -40 and +40 units, completing 2 full cycles every second.

This equation describes many periodic phenomena, such as a vibrating string, sound waves, light waves (for small angles in the wave equation), etc.

Functions of Waves and Signals

In mathematics, when we add two functions (say, $f(x)$ and $g(x)$) together), we simply add their output values at each input point. This is expressed as:

$(f + g)(x) = f(x) + g(x)$

This simple property of functions has profound implications in physics. When two waves occupy the same space simultaneously, they don’t collide or bounce off each other like physical objects. Instead, they superpose—their effects combine according to the principle of superposition.

The total displacement at any point is simply the sum of the individual wave displacements:

$y_{total}(t) = y_1(t) + y_2(t) + y_3(t) + ...$

Wave Interference

And this leads us to the fascinating phenomenon of wave interference.

Constructive Interference occurs when two waves are in phase (peaks align with peaks). Their amplitudes add together, creating a stronger wave:

$\sin(\omega t) + \sin(\omega t) = 2\sin(\omega t)$

For example, when two sound waves constructively interfere, they produce louder noise. In light waves, constructive interference creates bright spots on reflective surfaces.

Destructive Interference occurs when waves are out of phase (peaks align with troughs). They cancel each other out:

$\sin(\omega t) + \sin(\omega t + \pi) = 0$

This is the principle behind noise-cancelling headphones—they detect incoming sound waves and generate an exact opposite wave to cancel them out. When light waves destructively interfere, they appear as dim patches.

These wave behaviors demonstrate that mathematical function addition directly corresponds to physical wave superposition, bridging abstract mathematics with observable phenomena in our universe.

---

Now that we have a basic understanding of waves and signals, imagine being given a complex wave pattern, like a sound wave or an image signal, and being asked to figure out what simple sine and cosine waves make up that complex pattern. This is exactly where Fourier analysis comes in. Fourier analysis provides the mathematical tools to decompose complex signals into their constituent sine and cosine waves, allowing us to analyze the frequency content of the signal.

This means when you record a sound using your phone, Fourier analysis can identify which low amplitude high-frequency signals are present in the complex sound wave, apply a 180-degree phase shift to them, and add them back to the original sound wave to cancel out the noise.

---

Fourier Series

The Fourier series is a powerful tool that represents any reasonable^[2] periodic function as a sum of sine and cosine waves with different frequencies, amplitudes, and phases.

Periodic Signals

For periodic signals (functions that repeat over time, with period T), a signal 𝑓(𝑡) can be written as a sum of sine and cosine terms:

$f(t) = a_0 + \sum_{n=1}^{\infty} a_n \cos(2\pi n f_0 t) + b_n \sin(2\pi n f_0 t)$

where $f_0$ is the fundamental frequency (the lowest frequency component), and $a_0$, $a_n$, and $b_n$ are the Fourier coefficients that determine the amplitude of each frequency component:

• $a_n$ represents the amplitude of the cosine wave at frequency $n f_0$
• $b_n$ represents the amplitude of the sine wave at frequency $n f_0$

Alternatively, this can be written in amplitude-phase form:

$f(t) = a_0 + \sum_{n=1}^{\infty} A_n \cos(2\pi n f_0 t + \phi_n)$

where $A_n = \sqrt{a_n^2 + b_n^2}$ is the amplitude and $\phi_n = \arctan(-b_n/a_n)$ is the phase of the $n$-th harmonic component.

Sound waves are periodic signals, and the Fourier series can be used to decompose them into their constituent frequencies.

Approximating Complex Waves with Fourier Series

Fourier series can approximate any periodic function by summing up sine and cosine waves. The more terms we include in the series, the better our approximation becomes. With an infinite number of terms, we can perfectly reconstruct the original function. (In practice, we use a finite number of terms to get a good approximation.)

Let’s see this in action with a square wave - a function that abruptly alternates between +1 and -1:

As you increase the number of terms, notice how the red Fourier approximation gets closer to the gray target square wave. With just a few terms, we can already capture the general shape. With more terms, the approximation becomes remarkably accurate, though you’ll always see some overshoot at the discontinuities (known as the Gibbs phenomenon)—this happens because finite sums of smooth waves can’t perfectly match a jump discontinuity. Even functions with sharp corners and discontinuities can be represented as sums of smooth, continuous sine waves!

Drawing with Fourier Series: Epicycles

Coming back to recreating images using mechanical arms above, we can use Fourier series to break down complex 2D paths and use the same principles of wave decomposition to reconstruct them.

Sketches as Functions

When you draw a picture, your pen traces a path through 2D space. You go in two directions simultaneously: horizontally (x-axis) and vertically (y-axis). At any moment in time $t$, your pen has:

• An x-coordinate: $x(t)$
• A y-coordinate: $y(t)$

This means any drawing can be represented as two separate functions of time - one tracking horizontal movement, the other tracking vertical movement.

Just like we decomposed sound waves, we can analyze the path by separating it into two time-domain signals.

Fourier Approximation of the Bat Silhouette

We have two functions now: $x(t)$ and $y(t)$. Using Fourier series, we can approximate both functions as sums of sine and cosine waves. Which means, two sets of Fourier coefficients - one for $x(t)$ and one for $y(t)$ can be used to retrace the path of the drawing.

With just a few terms, we get a rough outline. As we increase the number of terms, the Fourier reconstruction gets closer and closer to the original bat silhouette.

Epicycles: Rotating Vectors

There are primarily two ways to represent coordinates in 2D space: Cartesian coordinates $(x, y)$ and Polar coordinates $(r, \theta)$. In Cartesian coordinates, a point is defined by its horizontal (x) and vertical (y) distances from the origin; P(x, y). In Polar coordinates, a point is defined by its distance from the origin (radius $r$) and the angle ($\theta$) it makes with the positive x-axis; P(r, θ).

Polar coordinates use trigonometric functions to relate to Cartesian coordinates: $x = r \cos(\theta)$, $y = r \sin(\theta)$. In the polar coordinates system, a point can be represented as a vector originating from the origin, with length $r$ and angle $\theta$.

Look at this example of $\sin(t)$ represented as a rotating vector in the complex plane:

Notice how the tip of the rotating vector traces out the sine wave over time. Similarly, a periodic function can be represented as a sum of rotating vectors (epicycles) in the complex plane. Attach an epicycle on the tip of another epicycle, and so on.

Complex Exponential Notation

Earlier, we wrote the Fourier Transform and DFT using separate sine and cosine terms. There’s a more elegant way to write this using Euler’s formula, which establishes a beautiful connection between exponentials and trigonometry:

$e^{i\theta} = \cos(\theta) + i\sin(\theta)$

This remarkable equation means a complex exponential represents a rotating vector in the complex plane. From this, we can derive:

$\cos(\theta) = \frac{e^{i\theta} + e^{-i\theta}}{2}$

$\sin(\theta) = \frac{e^{i\theta} - e^{-i\theta}}{2i}$

Why this matters: Instead of writing separate $A(\omega)\cos(\omega t) + B(\omega)\sin(\omega t)$ terms, we can combine them into a single complex exponential:

$F(\omega) = \int_{-\infty}^{\infty} f(t) e^{-i\omega t} dt$

This is the compact form of the Fourier Transform you’ll see in textbooks. It’s the same math as before, just written more concisely.

Similarly, the DFT can be written as:

$X[k] = \sum_{n=0}^{N-1} x[n] e^{-i2\pi nk/N}$

Now, coming back to the Fourier series representation of our drawing, we combine x and y into one complex signal:

$z(t) = x(t) + i y(t)$

Fourier coefficients can be expressed as:

$z(t) = \sum_{n=-\infty}^{\infty} c_n e^{i n \omega t}$

Each term $e^{i n \omega t}$ represents a rotating vector (an epicycle) at frequency $n$. The magnitude $|c_n|$ gives the radius, and $\arg(c_n)$ gives the starting phase.

Each term in the Fourier series can be visualized as a circle rotating at a specific frequency. When you connect these circles end-to-end (placing the center of each circle at the tip of the previous one), their combined motion traces out your original drawing.

Computing Epicycles with the DFT

Your drawing is sampled at discrete pixel points—not a continuous mathematical curve. For such sampled periodic signals, we use the Discrete Fourier Transform (DFT) to extract the frequency components.

The DFT takes $N$ sampled points and computes $N$ frequency components, each with amplitude and phase:

$X[k] = \sum_{n=0}^{N-1} x[n] \, e^{-i2\pi nk/N}$

Each complex number $X[k]$ tells us:

• Magnitude $|X[k]|$: radius of epicycle $k$
• Phase $\arg(X[k])$: starting angle of epicycle $k$

Since our drawing is 2D, we combine $x$ and $y$ coordinates into one complex signal: $z[n] = x[n] + iy[n]$. The DFT treats this closed path as one period of a repeating pattern—that’s why the epicycles loop perfectly!

Retracing our drawing with Epicycles

Summing up, we have:

• A sketch
• Decomposed into two time-domain signals: $x(t)$ and $y(t)$, representing horizontal and vertical movement over time
• Each horizontal and vertical signal approximated using Fourier series

Since your drawing is sampled at discrete points (not a continuous function), we use the Discrete Fourier Transform (DFT)—the sampled version of the Fourier transform.

Step 1: Analyzing the Drawing (DFT)

We compute complex coefficients $c_n$ by projecting the sampled path $(x_k + i y_k)$ onto each rotating basis $e^{-i2\pi nk/N}$:

for (let n = -numFrequencies; n <= numFrequencies; n++) {
    let realSum = 0, imagSum = 0;

    for (let k = 0; k < N; k++) {
        const angle = (-2 * Math.PI * n * k) / N;
        realSum += x[k] * Math.cos(angle) - y[k] * Math.sin(angle);
        imagSum += x[k] * Math.sin(angle) + y[k] * Math.cos(angle);
    }

    // Convert to polar form: amplitude and phase
    amplitude = sqrt(realSum² + imagSum²) / N;
    phase = atan2(imagSum, realSum);
}

This gives us:

• Amplitude: Radius of each rotating circle
• Phase: Starting angle of each circle
• Frequency: How fast each circle rotates

Step 2: Reconstructing the Path (Inverse DFT)

To redraw at time $t$ (from 0 to 1):

x(t) = Σ amplitude[n] * cos(2π * freq[n] * t + phase[n])
y(t) = Σ amplitude[n] * sin(2π * freq[n] * t + phase[n])

Each term is one rotating epicycle. Sum them all up to get the final pen position.

Step 3: Animation

for (t = 0; t < 1; t += 0.01) {
    position = evaluateEpicycles(coefficients, t);
    drawPoint(position.x, position.y);
}

As $t$ goes from 0 → 1, we complete one full cycle and trace the entire drawing!

Now, for each Fourier coefficient, we create an epicycle (rotating circle) in the complex plane. By connecting these epicycles tip-to-tail, the final tip traces out the original drawing.

If you’ve followed this far, you’d have a fair understanding of how we used Fourier series to decompose a 2D drawing into two time-domain signals, approximated those signals using sine and cosine waves, and represented those waves as rotating vectors (epicycles) in the complex plane.

The Mathematics Behind Epicycles

For a closed curve sampled at $N$ points $(x_k, y_k)$, we calculate the Fourier coefficients:

$c_n = \frac{1}{N} \sum_{k=0}^{N-1} (x_k + iy_k) e^{-2\pi i nk/N}$

Each coefficient $c_n$ can be converted to polar form:

• Radius: $|c_n|$ (the size of the rotating circle)
• Frequency: $n$ (how fast it rotates)
• Phase: $\arg(c_n)$ (where it starts)

At any time $t$, the position is:

$z(t) = \sum_{n=-M}^{M} c_n e^{2\pi int}$

where $M$ is the number of terms used in the approximation, and $z(t) = x(t) + iy(t)$ represents the pen’s position in the complex plane.

Why This Works

This technique works because:

1. Any periodic path can be approximated by summing sinusoids
2. Rotating vectors in the complex plane naturally create sinusoidal motion
3. Adding vectors tip-to-tail creates the epicycle visualization
4. The more terms we include, the closer we get to the original drawing

Generalizing: From Periodic to Non-Periodic Signals

So far, we’ve treated your drawing as a periodic signal—one complete cycle that the DFT assumes repeats infinitely. This is why closed curves work perfectly with epicycles.

But what about signals that truly don’t repeat? A recording of your voice, a single burst of noise, or any one-time event?

Fourier Transform for Non-Periodic Signals

The Fourier Transform generalizes the Fourier series to non-periodic signals. Instead of discrete harmonics (1×, 2×, 3× the fundamental frequency), we now have a continuous spectrum of all possible frequencies.

Using complex exponential notation, the Fourier Transform is:

Forward Transform (Signal → Frequency):

$F(\omega) = \int_{-\infty}^{\infty} f(t) e^{-i\omega t} \, dt$

Inverse Transform (Frequency → Signal):

$f(t) = \frac{1}{2\pi} \int_{-\infty}^{\infty} F(\omega) e^{i\omega t} \, d\omega$

where:

• $f(t)$ is the time-domain signal
• $F(\omega)$ is the frequency-domain representation (a complex number with amplitude and phase)
• $\omega$ is the angular frequency

The Fourier Transform is essential for analyzing non-repeating signals in audio processing, image filtering, and countless other applications.

The DFT: Sampled and Periodic

The Discrete Fourier Transform (DFT) we used for drawings combines two properties:

1. Discrete: Works with sampled points (not continuous functions)
2. Periodic: Assumes the sample sequence repeats infinitely

For a sequence of $N$ samples $x[0], x[1], \ldots, x[N-1]$:

$X[k] = \sum_{n=0}^{N-1} x[n] \, e^{-i2\pi nk/N}$

Inverse DFT (reconstructing the signal):

$x[n] = \frac{1}{N} \sum_{k=0}^{N-1} X[k] \, e^{i2\pi nk/N}$

Each $X[k]$ decomposes into:

• Amplitude: $|X[k]| = \sqrt{A[k]^2 + B[k]^2}$
• Phase: $\arg(X[k]) = \arctan(-B[k]/A[k])$

where $A[k]$ and $B[k]$ are the cosine and sine components.

The DFT is the foundation of all digital signal processing—analyzing digital audio, compressing images (JPEG), processing video, and more.

Fast Fourier Transform (FFT)

While our visualizations use the standard DFT algorithm shown above, real-world applications use a much faster algorithm called the Fast Fourier Transform (FFT).

The Performance Problem

The DFT we implemented has a critical limitation: computational complexity.

For a signal with $N$ points:

• DFT complexity: $O(N^2)$ operations
  • We loop through $N$ frequencies
  • For each frequency, we sum over $N$ points
  • Total: $N \times N = N^2$ operations

Example: For a 1-second audio clip at 44,100 Hz (CD quality):

• $N = 44,100$ samples
• DFT requires: $(44,100)^2 = 1.9$ billion operations
• On a modern CPU: several seconds to compute

This is too slow for real-time applications like audio processing, video streaming, or live music visualizers.

Enter the FFT

In 1965, Cooley and Tukey rediscovered an algorithm (originally found by Gauss in 1805) that computes the exact same result as DFT, but much faster:

FFT complexity: $O(N \log N)$ operations

Same example (44,100 samples):

• FFT requires: $44,100 \times \log_2(44,100) \approx 700,000$ operations
• Speedup: 2,700× faster than DFT!
• On a modern CPU: milliseconds instead of seconds

How FFT Works (Intuition)

The FFT uses a divide-and-conquer strategy:

1. Divide: Split the signal into even and odd indexed samples
2. Conquer: Recursively compute FFT on each half
3. Combine: Merge the results using clever math

The key insight: many of the calculations in DFT are redundant. The FFT eliminates this redundancy by reusing intermediate results.

// DFT: Computes everything from scratch
for (n in frequencies) {
    for (k in points) {
        // N² operations total
    }
}

// FFT: Divides work recursively
function FFT(points) {
    if (points.length === 1) return points;

    even = FFT(points[0, 2, 4, ...]);  // Recursion
    odd  = FFT(points[1, 3, 5, ...]);  // Recursion

    return combine(even, odd);  // N log N operations total
}

Performance deep dive:

For production systems, three details matter:

1. Power-of-two sizes: FFT is fastest when $N = 2^k$ (e.g., 512, 1024, 2048). This is because the divide-and-conquer splits evenly at each level. Non-power-of-two sizes work but require different algorithms (Bluestein’s, mixed-radix) that are slower.

2. Iterative vs recursive: The code above shows recursive FFT for clarity, but production libraries use iterative FFT (bit-reversal + in-place butterfly operations). Why? Recursive FFT has function call overhead and poor cache locality. Iterative FFT processes data in cache-friendly sequential passes.

3. Cache friendliness: Modern FFT implementations (like FFTW) use clever memory layouts and loop tiling to keep data in L1/L2 cache. A well-tuned FFT is often an order of magnitude faster in practice than a naive implementation, even with the same $O(N \log N)$ complexity.

This is why audio buffer sizes are always powers of two (256, 512, 1024)—it’s not arbitrary, it’s for FFT performance!

Why This Matters

DFT is great for:

• Learning and understanding Fourier analysis
• Small signals (< 1,000 points)
• Our interactive visualizations

FFT is essential for:

• Audio processing: Real-time effects, compression (MP3, AAC)
• Image processing: JPEG compression, filters, Instagram effects
• Video streaming: Netflix, YouTube use FFT for compression
• Telecommunications: 4G/5G networks, WiFi
• Medical imaging: MRI and CT scans
• Scientific computing: Weather prediction, astronomy

Libraries and Implementation

Most programming languages provide optimized FFT libraries:

// JavaScript (using Web Audio API)
const fft = new FFT(bufferSize);
fft.forward(audioData);  // O(N log N) instead of O(N²)

// Python (NumPy)
import numpy as np
frequencies = np.fft.fft(signal)  // Uses FFTW library

// MATLAB
X = fft(x);  // Built-in FFT

These libraries use highly optimized implementations like FFTW (Fastest Fourier Transform in the West), which can be 100× faster than a naive FFT implementation.

Equipped with this you can now go on to explore the vast west of signal processing. This new weaponry can help you build something cool, like:

• Implement noise reduction in sounds.
• Sharpen images with denoising.
• Create music visualizers that respond to audio frequencies.

---

1. This article including the image retracing is inspired from injuly’s fourier series blog. ↩︎

2. “Reasonable” functions are those that are piecewise continuous and have a finite number of discontinuities within a period. Not every periodic function is representable in a strong sense. There are periodic functions that are so pathological that: 1.the Fourier series fails to converge at some points, or even 2. diverges almost everywhere. These exist (classic results by du Bois-Reymond, later strengthened by others). You won’t meet them in engineering, but mathematically they kill the fully-general “any periodic function” claim. ↩︎

---