Convolution is usually introduced with its formal definition:

Yikes. Let's start without calculus: Convolution is fancy multiplication.

Part 1: Hospital Analogy

Imagine you manage a hospital treating patients with a single disease. You have:

• A treatment plan: [3] Every patient gets 3 units of the cure on their first day.
• A list of patients: [1 2 3 4 5] Your patient count for the week (1 person Monday, 2 people on Tuesday, etc.).

Question: How much medicine do you use each day? Well, that's just a quick multiplication:

Plan  *  Patients       = Daily Usage
[3]   *  [1 2 3 4 5]    = [3 6 9 12 15]

Multiplying the plan by the patient list gives the usage for the upcoming days: [3 6 9 12 15]. Everyday multiplication (3 x 4) means using the plan with a single day of patients: [3] * [4] = [12].

Intuition For Convolution

Let's say the disease mutates and requires a multi-day treatment. You create a new plan: Plan: [3 2 1]

That means 3 units of the cure on the first day, 2 on the second, and 1 on the third. Ok. Given the same patient schedule [1 2 3 4 5], what's our medicine usage each day?

Uh... shoot. It's not a quick multiplication:

• On Monday, 1 patient comes in. It's her first day, so she gets 3 units.
• On Tuesday, the Monday gal gets 2 units (her second day), but two new patients arrive, who get 3 each (2 * 3 = 6). The total is 2 + (2 * 3) = 8 units.
• On Wednesday, it's trickier: The Monday gal finishes (1 unit, her last day), the Tuesday people get 2 units (2 * 2), and there are 3 new Wednesday people... argh.

The patients are overlapping and it's hard to track. How can we organize this calculation?

An idea: imagine flipping the patient list, so the first patient is on the right:

           Start of line
 5 4 3 2 1

Next, imagine we have 3 separate rooms where we apply the proper dose:

 Rooms     3 2 1              

On your first day, you walk into the first room and get 3 units of medicine. The next day, you walk into room #2 and get 2 units. On the last day, you walk into room #3 and get 1 unit. There's no rooms afterwards, and your treatment is done.

To calculate the total medicine usage, line up the patients and walk them through the rooms:

 Monday
 ----------------------------
 Rooms                  3 2 1                                      
 Patients       5 4 3 2 1

 Usage                  3

On Monday (our first day), we have a single patient in the first room. She gets 3 units, for a total usage of 3. Makes sense, right?

On Tuesday, everyone takes a step forward:

 Tuesday
 ----------------------------
 Rooms                  3 2 1              
 Patients ->      5 4 3 2 1

 Usage                  6 2      = 8

The first patient is now in the second room, and there's 2 new patients in the first room. We multiply each room's dose by the patient count, then combine.

Every day we just walk the list forward:

 Wednesday
 ----------------------------
 Rooms                  3 2 1              
 Patients ->        5 4 3 2 1

 Usage                  9 4 1    = 14


 Thursday
 -----------------------------
 Rooms                  3 2 1              
 Patients ->          5 4 3 2 1

 Usage                 12 6 2    = 20


 Friday
 -----------------------------
 Rooms                  3 2 1              
 Patients ->            5 4 3 2 1

 Usage                 15 8 3    = 26

Whoa! It's intricate, but we figured it out, right? We can find the usage for any day by reversing the list, sliding it to the desired day, and combining the doses.

The total day-by-day usage looks like this (don't forget Sat and Sun, since some patients began on Friday):

Plan      *  Patient List   = Total Daily Usage

[3 2 1]   *  [1 2 3 4 5]    = [3 8 14 20 26 14 5]
              M T W T F        M T W  T  F  S  S

This calculation is the convolution of the plan and patient list. It's a fancy multiplication between a list of input numbers and a "program".

Interactive Demo

Here's a live demo. Try changing F (the plan) or G (the patient list). The convolution $c(t)$ matches our manual calculation above.

(We define functions $f(x)$ and $g(x)$ to pad each list with zero, and adjust for the list index starting at 1.)

You can do a quick convolution with Wolfram Alpha:

ListConvolve[{3, 2, 1}, {1, 2, 3, 4, 5}, {1, -1}, 0]
{3, 8, 14, 20, 26, 14, 5}

(The extra {1, -1}, 0 aligns the lists and pads with zero.)

Application: COVID Ventilator Usage

I started this article 5 years ago (intuition takes a while...), but unfortunately the analogy is relevant today.

Let's use convolution to estimate ventilator usage for incoming patients.

• Set $f(x)$ as the percent of patients needing ventilators. For example, [.05 .03 .01] means 5% of patients need ventilators the first week, 3% the second week, and 1% the third week.
• Set $g(x)$ as the weekly incoming patients, in thousands.
• The convolution $c(t) = f * g$, shows how many ventilators are needed each week (in thousands). $c(5)$ is how many ventilators are needed 5 weeks from now.

Let's try it out:

• F = [.05, .03, .01] is the ventilator use percentage by week
• G = [10, 20, 30, 20, 10, 10, 10], is the incoming hospitalized patients. It starts at 10k per week, rises to 30k, then decays to 10k.

With these numbers, we expect a max ventilator use of 2.2k in 2 weeks:

The convolution drops to 0 after 9 weeks because the patient list has run out. In this example, we're interested in the peak value the convolution hits, not the long-term total.

Other plans to convolve may be drug doses, vaccine appointments (one today, another a month from now), reinfections, and other complex interactions.

The hospital analogy is the mental model I wish I had when learning. Now that we've tried it with actual numbers, let's add the Math Juice and convert the analogy to calculus.

Part 2: The Calculus Definition

So, what happened in our example? We had a list of patients and a plan. If the plan were simple (single day [3]), regular multiplication would have worked. Because the plan was complex, we had to "convolve" it.

Time for some Fun Facts™:

• Convolution is written $f * g$, with an asterisk. Yes, an asterisk usually indicates multiplciation, but in advanced calculus class, it indicates a convolution. Regular multiplication is just implied ($fg$).

• The result of a convolution is a new function that gives the total usage for any day ("What was the total usage on day $t=3$?"). We can graph the convolution over time to see the day-by-day totals.

Now the big aha: Convolution reverses one of the lists! Here's why.

Let's call our treatment plan $f(x)$. In our example, we used [3 2 1].

The list of patients (inputs) is $g(x)$. However, we need to reverse this list as we slide it, so the earliest patient (Monday) enters the hospital first (first in, first out). This means we need to use $g(-x)$, the horizontal reflection of $g(x)$. [1 2 3 4 5] becomes [5 4 3 2 1].

Now that we have the reversed list, pick a day to compute ($t = 0, 1, 2...$). To slide our patient list by this much, we use: $g(-x + t)$. That is, we reverse the list ($-x$) and jump to the correct day ($+t$).

We have our scenario:

• $f(x)$ is the plan to use
• $g(-x + t)$ is the list of inputs (flipped and slid to the right day).

To get the total usage on day $t$, we multiply each patient with the plan, and sum the results (an integral). To account for any possible length, we go from -infinity to +infinity.

Now we can describe convolution formally using calculus:

(Like colorized math? There's more.)

Phew! That's quite few symbols. Some notes:

• We use a dummy variable $\tau$ (tau) for the intermediate computation. Imagine $\tau$ as knocking on each room ($\tau={0, 1, 2, 3...}$), finding the dosage [$f(\tau)$], the number of patients [$g(t - \tau)$], multiplying them, and totaling things in the integral. Yowza. The so-called "dummy" variable $\tau$ is like i in a for loop: it's temporary, but does the work. (By analogy, $t$ is a global variable has a fixed value during the loop: it's the day we're calculating the usage for, such as t = Day 5).
• In the official definition, you'll see $g(t - \tau)$ instead of $g(- \tau+ t)$. The second version shows the flip ($-\tau$) and slide ($+t$). Writing $g(t - \tau)$ makes it seem like we're interested in the difference between the variables, which confused me.
• The treatment plan (program to run) is called the kernel: you convolve a kernel with an input.

Not too bad, right? The equation is a formal description of the analogy.

Part 3: Mathematical Properties of Convolution

We can't discover a new math operation without taking it for a spin. Let's see how it behaves.

Convolution is commutative: f * g = g * f

In our computation, we flipped the patient list and kept the plan the same. Could we have flipped the plan instead?

You bet. Imagine the patients are immobile, and stay in their rooms: [1 2 3 4 5]. To deliver the medicine, we have 3 medical carts that go to each room and deliver the dose. Each day, they slide forward one position.

  Carts ->
  1 2 3

      1 2 3 4 5
      Patients

As before, though our plan is written [3 2 1] (3 units on the first day), we flip the order of the carts to[1 2 3]. That way, a patient gets 3 units on their first day, as we expect. Checking with Wolfram Alpha, the calculation is the same.

ListConvolve[{1, 2, 3, 4, 5}, {3, 2, 1}, {1, -1}, 0]
{3, 8, 14, 20, 26, 14, 5}

Cool! It looks like convolution is commutative:

and we can decide to flip either $f$ or $g$ when calculating the integral. Surprising, right?

The integral of the convolution

When all treatments are finished, what was the total medicine usage? This is the integral of the convolution. (A few minutes ago, that phrase would have you jumping out of a window.)

But it's a simple calculation. Our plan gives each patient sum([3 2 1]) = 6 units of medicine. And we have sum([1 2 3 4 5]) = 15 patients. The total usage is just 6 x 15 = 90 units.

Wow, that was easy: the usage for the entire convolution is just the product of the subtotals!

I hope this clicks intuitively. Note that this trick works for convolution, but not integrals in general. For example:

If we separate $x \cdot x$ into two integrals we get:

• $ \int (x \cdot x) = \int x^2 = \frac{1}{3} x^3 $
• $\int x \cdot \int x = \frac{1}{2}x^2 \cdot \frac{1}{2}x^2 = \frac{1}{4}x^4$

and those aren't the same. (Calculus would be much easier if we could split integrals like this.) It's strange, but $\int (f * g)$ is probably easier to solve than $\int (fg)$.

Impulse Response

What happens if we sent a single patient through the hospital? The convolution would just be that day's plan.

Plan    * Patients = Convolution
[3 2 1] * [1]      = [3 2 1]

In other words, convolving with [1] gives us the original plan.

In calculus terms, a spike of [1] (and 0 otherwise) is the Dirac Delta Function. In terms of convolutions, this function acts like the number 1 and returns the original function:

We can delay the delta function by T, which delays the resulting convolution function too. Imagine our single patient shows up a week late ($\delta(t - T)$), so our medicine usage gets delayed for a week too:

Part 4: Convolution Theorem & The Fourier Transform

The Fourier Transform (written with a fancy $\mathscr{F}$) converts a function $f(t)$ into a list of cyclical ingredients $F(s)$:

As an operator, this can be written $\mathscr{F}\lbrace f \rbrace = F$.

In our analogy, we convolved the plan and patient list with a fancy multiplication. Since the Fourier Transform gives us lists of ingredients, could we get the same result by mixing the ingredient lists?

Yep, we can: Fancy multiplication in the regular world is regular multiplication in the fancy world.

In math terms, "Convolution in the time domain is multiplication in the frequency (Fourier) domain."

Mathematically, this is written:

or

where $f(x)$ and $g(x)$ are functions to convolve, with transforms $F(s)$ and $G(s)$.

We can prove this theorem with advanced calculus, that uses theorems I don't quite understand, but let's think through the meaning.

Because $F(s)$ is the Fourier Transform of $f(t)$, we can ask for a specific frequency ($s = 2\text{Hz}$) and get the combined interaction of every data point with that frequency. Let's suppose:

That means after every data point has been multiplied against the 2Hz cycle, the result is $3 + i$. But we could have kept each interaction separate:

Where $c_t$ is the contribution to the 2Hz frequency from datapoint $t$. Similarly, we can expand $G(s)$ into a list of interactions with the 2Hz ingredient. Let's suppose $G(2) = 7 - i$:

The Convolution Theorem is really saying:

Our convolution in the regular domain involves a lot of cross-multiplications. In the fancy frequency domain, we still have a bunch of interactions, but $F(s)$ and $G(s)$ have consolidated them. We can just multiply $F(2)G(2) = (3 + i)(7-i)$ to find the 2Hz ingredient in the convolved result.

By analogy, suppose you want to calculate:

It's a lot of work to cross-multiply every term: $(1 \cdot 5) + (1\cdot 6) + (1\cdot 7) + ...$

It's better to consolidate the groups into $(1 + 2 + 3 + 4) = 10$ and $(5 + 6 + 7 + 8) = 26$, and then multiply to get $10 \cdot 26 = 260$.

This nuance caused me a lot of confusion. It seems like $FG$ is a single multiplication, while $f * g$ involves a bunch of intermediate terms. I forgot that $F$ already did the work of merging a bunch of entries into a single one.

Now, we aren't quite done.

We can convert $f * g$ in the time domain into $FG$ in the frequency domain, but we probably need it back in the time domain for a usable result:

You have a riddle in English ($f * g$), you translate it to French ($FG$), get your smart French friend to work out that calculation, then convert it back to English ($\mathscr{F}^{-1}$).

And in reverse...

The convolution theorem works this way too:

Regular multiplication in the regular world is fancy multiplication in the fancy world.

Cool, eh? Instead of multiplying two functions like some cave dweller, put on your monocle, convolve the Fourier Transforms, and and convert to the time domain:

I'm not saying this is fun, just that it's possible. If your French friend has a gnarly calculation they're struggling with, it might look like arithmetic to you.

Mini proof

Remember how we said the integral of a convolution was a multiplication of the individual integrals?

Well, the Fourier Transform is just a very specific integral, right?

So (handwaving), it seems we could swap the general-purpose integral $\int$ for $\mathscr{F}$ and get

which is the convolution theorem. I need a deeper intuition for the proof, but this helps things click.

Part 5: Applications

The trick with convolution is finding a useful "program" (kernel) to apply to your input. Here's a few examples.

Moving averages

Let's say you want a moving average between neighboring items in a list. That is half of each element, added together:

This is a "multiplication program" of [0.5 0.5] convolved with our list:

ListConvolve[{1, 4, 9, 16, 25}, {0.5, 0.5}, {1, -1}, 0] 
{0.5, 2.5, 6.5, 12.5, 20.5, 12.5}

We can perform a moving average with a single operation. Neat!

A 3-element moving average would be [.33 .33 .33], a weighted average could be [.5 .25 .25].

Derivatives

The derivative finds the difference between neighboring values. Here's the plan: [1 -1]

ListConvolve[{1, 2, 3, 4, 5}, {1, -1}, {1, -1}, 0] 
{1, 1, 1, 1, 1, -5}                         // -5 since we ran out of entries

ListConvolve[{1, 4, 9, 16, 25}, {1, -1}, {1, -1}, 0] 
{1, 3, 5, 7, 9, -25}                        // discrete derivative is 2x + 1

With a simple kernel, we can find a useful math property on a discrete list. And to get a second derivative, just apply the derivative convolution twice:

F * [1 -1] * [1 -1]

As a shortcut, we can precompute the final convolutions ([1 -1] * [1 -1] ) and get:

ListConvolve[{1, -1}, {1,-1}, {1, -1}, 0] 
{1, -2, 1}

Now we have a single kernel [1, -2, 1] that gets the second derivative of a list:

ListConvolve[{1, 4, 9, 16, 25}, {1, -2, 1}, {1, -1}, 0] 
{1, 2, 2, 2, 2, -34, 25}

Excluding the boundary items, we get the expected second derivative:

Blurring / unblurring images

An image blur is essentially a convolution of your image with some "blurring kernel":

The blur of our 2D image requires a 2D average:

Can we undo the blur? Yep! With our friend the Convolution Theorem, we can do:

Whoa! We can recover the original image by dividing out the blur. Convolution is a simple multiplication in the frequency domain, and deconvolution is a simple division in the frequency domain.

A short while back, the concept of "deblurring by dividing Fourier Transforms" was gibberish to me. While it can be daunting mathematically, it's getting simpler conceptually.

More reading:

• Wikipedia on deconvolution, deblurring.
• Image Restoration Lecture, another one, and example for above
• Blind deconvolution is when you don't know the blur kernel (make a guess)

Algorithm Trick: Multiplication

What is a number? A list of digits:

1234 = 1000 + 200 + 30 + 4 = [1000 200 30 4]
5678 = 5000 + 600 + 70 + 8 = [5000 600 70 8]

And what is regular, grade-school multiplication? A digit-by-digit convolution! We sweep one list of digits by the other, multiplying and adding as we go:

Source

We can perform the calculation by convolving the lists of digits (wolfram alpha):

ListConvolve[{1000, 200, 30, 4}, {8, 70, 600, 5000}, {1, -1}, 0]
{8000, 71600, 614240, 5122132, 1018280, 152400, 20000}

sum {8000, 71600, 614240, 5122132, 1018280, 152400, 20000}
7006652

Note that we pre-flip one of the lists (it gets swapped in the convolution later), and the intermediate calculations are a bit different. But, combining the subtotals gives the expected result.

Faster Convolutions

Why convolve instead of doing a regular digit-by-digit multiplication? Well, the convolution theorem lets us substitute convolution with Fourier Transforms:

The convolution ($f * g$) has complexity $O(n^2)$. We have $n$ positions to process, with $n$ intermediate multiplications at each position.

The right side involves:

• Two Fourier Transforms, which are normally $O(n^2)$. However, the Fast Fourier Transform (a divide-and-conquer approach) makes them $O(n\log(n))$.
• Pointwise multiplication of the final result of the transforms ($\sum a_n \cdot b_n$), which is $O(n)$
• An inverse transform, which is $O(n\log(n))$

And the total complexity is: $O(n\log(n)) + O(n\log(n)) + O(n) + O(n\log(n)) = O(n\log(n))$

Regular multiplication in the fancy domain is faster than a fancy multiplication in the regular domain. Our French friend is no slouch. (More)

Convolutional Neural Nets (CNN)

Machine learning is about discovering the math functions that transform input data into a desired result (a prediction, classification, etc.).

Starting with an input signal, we could convolve it with a bunch of kernels:

Given that convolution can do complex math (moving averages, blurs, derivatives...), it seems some combination of kernels should turn our input into something useful, right?

Convolutional Neural Nets (CNNs) process an input with layers of kernels, optimizing their weights (plans) to reach a goal. Imagine tweaking the treatment plan to keep medicine usage below some threshold.

CNNs are often used with image classifiers, but 1D data sets work just fine.

• Nice writeup: https://ujjwalkarn.me/2016/08/11/intuitive-explanation-convnets/
• Digit classifier demo: https://cs.stanford.edu/people/karpathy/convnetjs/demo/mnist.html

LTI System Behavior

A linear, time-invariant system means:

• Linear: Scaling and combining inputs scales and combines outputs by the same amount
• Time invariant: Outputs depend on relative time, not absolute time. You get 3 units on your first day, and it doesn't matter if it's Wednesday or Thursday.

A fancy phrase is "A LTI system is characterized by its impulse response". Translation: If we send a single patient through the hospital [1], we'll discover the treatment plan. Then we can predict the usage for any sequence of patients by convolving it with the plan.

If the system isn't LTI, we can't extrapolate based on a single person's experience. Scaling the inputs may not scale the outputs, and the actual calendar day, not relative day, may impact the result (imagine fewer rooms available on weekends).

Engineering Analogies

From David Greenspan: "Suppose you have a special laser pointer that makes a star shape on the wall. You tape together a bunch of these laser pointers in the shape of a square. The pattern on the wall now is the convolution of a star with a square."

Regular multiplication gives you a single scaled copy of an input. Convolution creates multiple overlapping copies that follow a pattern you've specified.

Real-world systems have squishy, not instantaneous, behavior: they ramp up, peak, and drop down. The convolution lets us model systems that echo, reverb and overlap.

Now it's time for the famous sliding window example. Think of a pulse of inputs (red) sliding through a system (blue), and having a combined effect (yellow): the convolution.

(Source)

Summary

Convolution has an advanced technical definition, but the basics can be understood with the right analogy.

Quick rant: I study math for fun, yet it took years to find a satisfying intuition for:

• Why is one function reversed?
• Why is convolution commutative?
• Why does the integral of the convolution = product of integrals?
• Why are the Fourier Transforms multiplied point-by-point, and not overlapped?

Why'd it take so long? Imagine learning multiplication with $f \times g = z$ instead of $3 \times 5 = 15$. Without an example I can explore in my head, I could only memorize results, not intuit them. Hopefully this analogy can save you years of struggle.

Happy math.

Yikes. Rather than jumping into the symbols, let's experience the key idea firsthand. Here's a plain-English metaphor:

• What does the Fourier Transform do? Given a smoothie, it finds the recipe.
• How? Run the smoothie through filters to extract each ingredient.
• Why? Recipes are easier to analyze, compare, and modify than the smoothie itself.
• How do we get the smoothie back? Blend the ingredients.

Here's the "math English" version of the above:

• The Fourier Transform takes a time-based pattern, measures every possible cycle, and returns the overall "cycle recipe" (the amplitude, offset, & rotation speed for every cycle that was found).

Time for the equations? No! Let's get our hands dirty and experience how any pattern can be built with cycles, with live simulations.

If all goes well, we'll have an aha! moment and intuitively realize why the Fourier Transform is possible. We'll save the detailed math analysis for the follow-up.

This isn't a force-march through the equations, it's the casual stroll I wish I had. Onward!

From Smoothie to Recipe

A math transformation is a change of perspective. We change our notion of quantity from "single items" (lines in the sand, tally system) to "groups of 10" (decimal) depending on what we're counting. Scoring a game? Tally it up. Multiplying? Decimals, please.

The Fourier Transform changes our perspective from consumer to producer, turning What do I have? into How was it made?

In other words: given a smoothie, let's find the recipe.

Why? Well, recipes are great descriptions of drinks. You wouldn't share a drop-by-drop analysis, you'd say "I had an orange/banana smoothie". A recipe is more easily categorized, compared, and modified than the object itself.

So... given a smoothie, how do we find the recipe?

Well, imagine you had a few filters lying around:

• Pour through the "banana" filter. 1 oz of bananas are extracted.
• Pour through the "orange" filter. 2 oz of oranges.
• Pour through the "milk" filter. 3 oz of milk.
• Pour through the "water" filter. 3 oz of water.

We can reverse-engineer the recipe by filtering each ingredient. The catch?

• Filters must be independent. The banana filter needs to capture bananas, and nothing else. Adding more oranges should never affect the banana reading.

• Filters must be complete. We won't get the real recipe if we leave out a filter ("There were mangoes too!"). Our collection of filters must catch every possible ingredient.

• Ingredients must be combine-able. Smoothies can be separated and re-combined without issue (A cookie? Not so much. Who wants crumbs?). The ingredients, when separated and combined in any order, must make the same result.

See The World As Cycles

The Fourier Transform takes a specific viewpoint: What if any signal could be filtered into a bunch of circular paths?

Whoa. This concept is mind-blowing, and poor Joseph Fourier had his idea rejected at first. (Really Joe, even a staircase pattern can be made from circles?)

And despite decades of debate in the math community, we expect students to internalize the idea without issue. Ugh. Let's walk through the intuition.

The Fourier Transform finds the recipe for a signal, like our smoothie process:

• Start with a time-based signal
• Apply filters to measure each possible "circular ingredient"
• Collect the full recipe, listing the amount of each "circular ingredient"

Stop. Here's where most tutorials excitedly throw engineering applications at your face. Don't get scared; think of the examples as "Wow, we're finally seeing the source code (DNA) behind previously confusing ideas".

• If earthquake vibrations can be separated into "ingredients" (vibrations of different speeds & amplitudes), buildings can be designed to avoid interacting with the strongest ones.

• If sound waves can be separated into ingredients (bass and treble frequencies), we can boost the parts we care about, and hide the ones we don't. The crackle of random noise can be removed. Maybe similar "sound recipes" can be compared (music recognition services compare recipes, not the raw audio clips).

• If computer data can be represented with oscillating patterns, perhaps the least-important ones can be ignored. This "lossy compression" can drastically shrink file sizes (and why JPEG and MP3 files are much smaller than raw .bmp or .wav files).

• If a radio wave is our signal, we can use filters to listen to a particular channel. In the smoothie world, imagine each person paid attention to a different ingredient: Adam looks for apples, Bob looks for bananas, and Charlie gets cauliflower (sorry bud).

The Fourier Transform is useful in engineering, sure, but it's a metaphor about finding the root causes behind an observed effect.

Think With Circles, Not Just Sinusoids

One of my giant confusions was separating the definitions of "sinusoid" and "circle".

• A "sinusoid" is a specific back-and-forth pattern (a sine or cosine wave), and 99% of the time, it refers to motion in one dimension.
• A "circle" is a round, 2d pattern you probably know. If you enjoy using 10-dollar words to describe 10-cent ideas, you might call a circular path a "complex sinusoid".

Labeling a circular path as a "complex sinusoid" is like describing a word as a "multi-letter". You zoomed into the wrong level of detail. Words are about concepts, not the letters they can be split into!

The Fourier Transform is about circular paths (not 1-d sinusoids) and Euler's formula is a clever way to generate one:

Must we use imaginary exponents to move in a circle? Nope. But it's convenient and compact. And sure, we can describe our path as coordinated motion in two dimensions (real and imaginary), but don't forget the big picture: we're just moving in a circle.

Following Circular Paths

Let's say we're chatting on the phone and, like usual, I want us to draw the same circle simultaneously. (You promised!) What should I say?

• How big is the circle? (Amplitude, i.e. size of radius)
• How fast do we draw it? (Frequency. 1 circle/second is a frequency of 1 Hertz (Hz) or 2*pi radians/sec)
• Where do we start? (Phase angle, where 0 degrees is the x-axis)

I could say "2-inch radius, start at 45 degrees, 1 circle per second, go!". After half a second, we should each be pointing to: starting point + amount traveled = 45 + 180 = 225 degrees (on a 2-inch circle).

Every circular path needs a size, speed, and starting angle (amplitude/frequency/phase). We can even combine paths: imagine tiny motorcars, driving in circles at different speeds.

The combined position of all the cycles is our signal, just like the combined flavor of all the ingredients is our smoothie.

Here's a simulation of a basic circular path:

(Based on this animation, here's the source code. Modern browser required. Click the graph to pause/unpause.)

The magnitude of each cycle is listed in order, starting at 0Hz. Cycles [0 1] means

• 0 amplitude for the 0Hz cycle (0Hz = a constant cycle, stuck on the x-axis at zero degrees)
• 1 amplitude for the 1Hz cycle (completes 1 cycle per time interval)

Now the tricky part:

• The blue graph measures the real part of the cycle. Another lovely math confusion: the real axis of the circle, which is usually horizontal, has its magnitude shown on the vertical axis. You can mentally rotate the circle 90 degrees if you like.
• The time points are spaced at the fastest frequency. A 1Hz signal needs 2 time points for a start and stop (a single data point doesn't have a frequency). The time values [1 -1] shows the amplitude at these equally-spaced intervals.

With me? [0 1] is a pure 1Hz cycle.

Now let's add a 2Hz cycle to the mix. [0 1 1] means "Nothing at 0Hz, 1Hz of amplitude 1, 2Hz of amplitude 1":

Whoa. The little motorcars are getting wild: the green lines are the 1Hz and 2Hz cycles, and the blue line is the combined result. Try toggling the green checkbox to see the final result clearly. The combined "flavor" is a sway that starts at the max and dips low for the rest of the interval.

The yellow dots are when we actually measure the signal. With 3 cycles defined (0Hz, 1Hz, 2Hz), each dot is 1/3 of the way through the signal. In this case, cycles [0 1 1] generate the time values [2 -1 -1], which starts at the max (2) and dips low (-1).

Oh! We can't forget phase, the starting angle! Use magnitude:angle to set the phase. So [0 1:45] is a 1Hz cycle that starts at 45 degrees:

This is a shifted version of [0 1]. On the time side we get [.7 -.7] instead of [1 -1], because our cycle isn't exactly lined up with our measuring intervals, which are still at the halfway point (this could be desired!).

The Fourier Transform finds the set of cycle speeds, amplitudes and phases to match any time signal.

Our signal becomes an abstract notion that we consider as "observations in the time domain" or "ingredients in the frequency domain".

Enough talk: try it out! In the simulator, type any time or cycle pattern you'd like to see. If it's time points, you'll get a collection of cycles (that combine into a "wave") that matches your desired points.

But… doesn't the combined wave have strange values between the yellow time intervals? Sure. But who's to say whether a signal travels in straight lines, or curves, or zips into other dimensions when we aren't measuring it? It behaves exactly as we need at the equally-spaced moments we asked for.

Making A Spike In Time

Can we make a spike in time, like (4 0 0 0), using cycles? I'll use parentheses () for a sequence of time points, and brackets [] for a sequence of cycles.

Although the spike seems boring to us time-dwellers (one data point, that's it?), think about the complexity in the cycle world. Our cycle ingredients must start aligned (at the max value, 4) and then "explode outwards", each cycle with partners that cancel it in the future. Every remaining point is zero, which is a tricky balance with multiple cycles running around (we can't just "turn them off").

Let's walk through each time point:

• At time 0, the first instant, every cycle ingredient is at its max. Ignoring the other time points, (4 ? ? ?) can be made from 4 cycles (0Hz 1Hz 2Hz 3Hz), each with a magnitude of 1 and phase of 0 (i.e., 1 + 1 + 1 + 1 = 4).

• At every future point (t = 1, 2, 3), the sum of all cycles must cancel.

Here's the trick: when two cycles are on opposites sides of the circle (North & South, East & West, etc.) their combined position is zero (3 cycles can cancel if they're spread evenly at 0, 120, and 240 degrees).

Imagine a constellation of points moving around the circle. Here's the position of each cycle at every instant:

Time 0 1 2 3 
------------
0Hz: 0 0 0 0 
1Hz: 0 1 2 3
2Hz: 0 2 0 2
3Hz: 0 3 2 1

Notice how the the 3Hz cycle starts at 0, gets to position 3, then position "6" (with only 4 positions, 6 modulo 4 = 2), then position "9" (9 modulo 4 = 1).

When our cycle is 4 units long, cycle speeds a half-cycle apart (2 units) will either be lined up (difference of 0, 4, 8…) or on opposite sides (difference of 2, 6, 10…).

OK. Let's drill into each time point:

• Time 0: All cycles at their max (total of 4)
• Time 1: 1Hz and 3Hz cancel (positions 1 & 3 are opposites), 0Hz and 2Hz cancel as well. The net is 0.
• Time 2: 0Hz and 2Hz line up at position 0, while 1Hz and 3Hz line up at position 2 (the opposite side). The total is still 0.
• Time 3: 0Hz and 2Hz cancel. 1Hz and 3Hz cancel.
• Time 4 (repeat of t=0): All cycles line up.

The trick is having individual speeds cancel (0Hz vs 2Hz, 1Hz vs 3Hz), or having the lined-up pairs cancel (0Hz + 2Hz vs 1Hz + 3Hz).

When every cycle has equal power and 0 phase, we start aligned and cancel afterwards. (I don't have a nice proof yet -- any takers? -- but you can see it yourself. Try [1 1], [1 1 1], [1 1 1 1] and notice the signals we generate: (2 0), (3 0 0), (4 0 0 0)).

In my head, I label these signals as "time spikes": they have a value for a single instant, and are zero otherwise (the fancy name is a delta function.)

Here's how I visualize the initial alignment, followed by a net cancellation:

Moving The Time Spike

Not everything happens at t=0. Can we change our spike to (0 4 0 0)?

It seems the cycle ingredients should be similar to (4 0 0 0), but the cycles must align at t=1 (one second in the future). Here's where phase comes in.

Imagine a race with 4 runners. Normal races have everyone lined up at the starting line, the (4 0 0 0) time pattern. Boring.

What if we want everyone to finish at the same time? Easy. Just move people forward or backwards by the appropriate distance. Maybe granny can start 2 feet in front of the finish line, Usain Bolt can start 100m back, and they can cross the tape holding hands.

Phase shifts, the starting angle, are delays in the cycle universe. Here's how we adjust the starting position to delay every cycle 1 second:

• A 0Hz cycle doesn't move, so it's already aligned
• A 1Hz cycle goes 1 revolution in the entire 4 seconds, so a 1-second delay is a quarter-turn. Phase shift it 90 degrees backwards (-90) and it gets to phase=0, the max value, at t=1.
• A 2Hz cycle is twice as fast, so give it twice the angle to cover (-180 or 180 phase shift -- it's across the circle, either way).
• A 3Hz cycle is 3x as fast, so give it 3x the distance to move (-270 or +90 phase shift)

If time points (4 0 0 0) are made from cycles [1 1 1 1], then time points (0 4 0 0) are made from [1 1:-90 1:180 1:90]. (Note: I'm using "1Hz", but I mean "1 cycle over the entire time period").

Whoa -- we're working out the cycles in our head!

The interference visualization is similar, except the alignment is at t=1.

Test your intuition: Can you make (0 0 4 0), i.e. a 2-second delay? 0Hz has no phase. 1Hz has 180 degrees, 2Hz has 360 (aka 0), and 3Hz has 540 (aka 180), so it's [1 1:180 1 1:180].

Discovering The Full Transform

The big insight: our signal is just a bunch of time spikes! If we merge the recipes for each time spike, we should get the recipe for the full signal.

The Fourier Transform builds the recipe frequency-by-frequency:

• Separate the full signal (a b c d) into "time spikes": (a 0 0 0) (0 b 0 0) (0 0 c 0) (0 0 0 d)
• For any frequency (like 2Hz), the tentative recipe is "a/4 + b/4 + c/4 + d/4" (the amplitude of each spike is split among all frequencies)
• Wait! We need to offset each spike with a phase delay (the angle for a "1 second delay" depends on the frequency).
• Actual recipe for a frequency = a/4 (no offset) + b/4 (1 second offset) + c/4 (2 second offset) + d/4 (3 second offset).

We can then loop through every frequency to get the full transform.

Here's the conversion from "math English" to full math:

A few notes:

• N = number of time samples we have
• n = current sample we're considering (0 .. N-1)
• x_n = value of the signal at time n
• k = current frequency we're considering (0 Hertz up to N-1 Hertz)
• X_k = amount of frequency k in the signal (amplitude and phase, a complex number)
• The 1/N factor is usually moved to the reverse transform (going from frequencies back to time). This is allowed, though I prefer 1/N in the forward transform since it gives the actual sizes for the time spikes. You can get wild and even use $1/\sqrt{N}$ on both transforms (going forward and back creates the 1/N factor).
• n/N is the percent of the time we've gone through. 2 * pi * k is our speed in radians / sec. e^-ix is our backwards-moving circular path. The combination is how far we've moved, for this speed and time.
• The raw equations for the Fourier Transform just say "add the complex numbers". Many programming languages cannot handle complex numbers directly, so you convert everything to rectangular coordinates and add those.

Onward

This was my most challenging article yet. The Fourier Transform has several flavors (discrete/continuous/finite/infinite), covers deep math (Dirac delta functions), and it's easy to get lost in details. I was constantly bumping into the edge of my knowledge.

But there's always simple analogies out there -- I refuse to think otherwise. Whether it's a smoothie or Usain Bolt & Granny crossing the finish line, take a simple understanding and refine it. The analogy is flawed, and that's ok: it's a raft to use, and leave behind once we cross the river.

I realized how feeble my own understanding was when I couldn't work out the transform of (1 0 0 0) in my head. For me, it was like saying I knew addition but, gee whiz, I'm not sure what "1 + 1 + 1 + 1" would be. Why not? Shouldn't we have an intuition for the simplest of operations?

That discomfort led me around the web to build my intuition. In addition to the references in the article, I'd like to thank:

• Scott Young, for the initial impetus for this post
• Shaheen Gandhi, Roger Cheng, and Brit Cruise for kicking around ideas & refining the analogy
• Steve Lehar for great examples of the Fourier Transform on images
• Charan Langton for her detailed walkthrough
• Julius Smith for a fantastic walkthrough of the Discrete Fourier Transform (what we covered today)
• Bret Victor for his techniques on visualizing learning

Today's goal was to experience the Fourier Transform. We'll save the advanced analysis for next time.

Happy math.

Appendix: Projecting Onto Cycles

Stuart Riffle has a great interpretation of the Fourier Transform:

Imagine spinning your signal in a centrifuge and checking for a bias. I have a correction: we must spin backwards (the exponent in the equation above should be $e^{-i 2 \pi...}$). You already know why: we need a phase delay so spikes appear in the future.

Appendix: Another Awesome Visualization

Lucas Vieira, author of excellent Wikipedia animations, was inspired to make this interactive animation:

Fourier Toy - Click to download, requires flash

(Detailed list of control options)

The Fourier Transform is about cycles added to cycles added to cycles. Try making a "time spike" by setting a amplitude of 1 for every component (press Enter after inputting each number). Fun fact: with enough terms, you can draw any shape, even Homer Simpson.

Check out http://www.jezzamon.com/fourier/ for a great tool to draw any shape using epicycles.

Appendix: Article with R code samples

João Neto made a great writeup, with technical (R) code samples here:

http://www.di.fc.ul.pt/~jpn/r/fourier/fourier.html

Appendix: Using the code

All the code and examples are open source (MIT licensed, do what you like).

• Interactive example (view source)
• Github gist
• Reddit discussion on details of the computation, I'm pb_zeppelin