Discrete Fourier Transform

Warning- a whole lot of text and math speak ahead. Continue at your own risk.

Okay, let’s start with the Fourier Transform. The core idea is that a signal (= sequence of numbers over time) can be decomposed into a sum of sinusoids. To maybe sort of get an intuition for why, think about any point in a two-dimensional plane, such as (3,4), or (17,42). Both of these points, and every other point on the plane, can be represented as a sum of multiples of (1,0) and (0,1):

Notice that the reason you only need to know two things is that there are only two degrees of freedom: how far right or left, and how far up or down. In fact, I can uniquely determine any point in the plane given “any” two properties of that point as long as they are sufficiently not redundant (this collection of “properties” is referred to as a basis). For example, the point \((1,1)\) in Cartesian Coordinates, signifying a horizontal offset of positive one and a negative offset of positive one, can also be denoted as \((\sqrt{2},\frac{\pi}{4})\) in Polar Coordinates, signifying a radial offset of \(\sqrt{2}\) and an angular offset from the positive x-axis of \(\frac{\pi}{4}\) (radians). This reinterpretation of a unique point in some N-space is called a change of basis. One can think of the Fourier Transform as a type of change of basis.

In order to extend this understanding to signals, you must realize that there are infinitely many degrees of freedom: how big the signal is at the start, how big the signal is one infinitesimal moment later, how big the signal is the moment after that, and so on and so forth. Like we said before, since there are infinitely many degrees of freedom, we can uniquely classify any signal with infinitely many “properties”. It may seem like this doesn’t help at all, but defining the basis as the set of all harmonics of a sine wave has a nice real world analogue in the way we humans perceive sounds. And, as you will see, the cyclic nature of sinusoids also allows for a tangential interpretation.

Side note: if instead you defined the basis as the set of all powers of x, you get the Taylor series.

Congrats, you passed Fourier 101. Onto the main course:

The idea is that a parametric curve in 2-dimensions \(f(t) = (a(t),b(t))\) can be just as easily written as a complex function of one variable (ie. a signal) \(g(t) = a(t)+b(t)i\), after which I can take the Discrete Fourier Transform (DFT) to get a list of coefficients of the complex exponential basis vectors (by the way, in this context, a function is just a vector whose components are the values of the function at every point). The reason I’m using the DFT instead of the at least three other versions of the Fourier Transform is because it takes a discrete periodic signal in the time domain and returns a discrete periodic signal in the frequency domain. I am not really sure why, and most of the resources online either assume it or are really hand-wavey. What I do know is that a discreteness in one domain causes periodicity in the other, and vice versa. Therefore, since I want discreteness in my frequency domain (ie. a finite number of complex exponentials), I impose periodicity onto my time signal (ie. the parametric curve must be a closed loop). And since I can only get discrete time samples, my frequency domain will be periodic (I think this causes/is caused by/is related to aliasing).

Onto some implementation details. Since the DFT is pretty ubiquitous in many digital applications, it wasn’t too difficult to find a working implementation of it online. First of all, the input time-domain signals are obtained and formatted as a list of list of two elements representing the real and imaginary components. Once the input signals are obtained, they are passed into the dft function, which returns a new list. I adapted the code from the above website which was written in C++ to JavaScript. The site has a pretty easy to follow explanation of the implementation, so I won’t bother to go over it here.

function dft(time) {
  var ret = []
  var n = time.length
  for (var i = 0 ; i < n ; ++i) {
    var sreal = 0, simag = 0, angle = 0
    for (var j = 0 ; j < n ; ++j) {
      angle = tpi * i * j / n; // tpi = 2*Math.PI
      sreal +=  time[j][0] * Math.cos(angle) + time[j][1] * Math.sin(angle);
      simag += -time[j][0] * Math.sin(angle) + time[j][1] * Math.cos(angle);
    }
    ret.push([sreal/n,simag/n,i])
  }
  return ret
}

There is one implementation detail to note, however. Notice that my return type is a list of lists of three elements, not two. In addition to the real and imaginary components of the complex coefficient, I also push on the index, or the frequency value, because later on I sort the list in terms of radius, which mixes up the order.

Now comes the fun part. Here’s the formula for the Inverse DFT (IDFT):

(That’s actually not entirely correct; the IDFT is only defined at the sampled points, but here I’m allowing it to be defined continuously on \(0<t<2\pi\).)

There are a few other things to realize before moving on:

1. complex addition = vector addition
2. complex multiplication = rotation (ie. the product of two complex numbers is a complex number whose angle is the sum of the two angles and whose magnitude is the product of the two magnitudes)
3. complex exponentials can be interpreted as points around a unit circle.
4. From (2) and (3), it follows that multiplying a complex number by a complex exponential \(e^{i\theta}\) is equivalent to rotating the prior by \(\theta\)

First, take a look at a single term of the sum. Since \(X[k]\) is a complex number, the whole term can just be interpreted as a point at \(X[k]\) rotating at a constant speed in the counter clockwise direction. Extending this to multiple terms, it turns out that the IDFT is actually just the sum of rotating vectors.

Enter the Ptolemaic planetary system. Basically, this was way back when people still believed that the Earth was at center of the universe, and reasoned the crazy complicated paths of the other celestial bodies by saying that their motion can be described by picturing them as the very end of a long chain of satellites. For a clearer explanation, see this Mathematics StackExchange post (in fact, that post describes much, if not all, of what I try to accomplish in this project).

Now to combine the two ideas, the key is to understand that a satellite’s orbit around another celestial body is simply a rotating vector. If the discourse from the above seven-ish paragraphs is now combined, the conclusion is that any closed path can be drawn by an infinite number of planets orbiting the next in an infinitely long dizzying planetary conga line. Or, for practical purposes, any closed path can be approximated by a finite number of planets (perhaps you can see now why the Ptolemaic planetary system was an absurd argument for the geocentric system).

While the argument was absurd for explaining the motions of the planets in our sky, it’s perfect for drawing pictures with circles. After calculating the coefficients of the complex exponentials, the only thing left was to draw each circle with a radius of that circle’s coefficient, and rotate it with a speed of that circle’s frequency.

for (var i = 0 ; i < keyframes.length-options.compression; ++i ) {
  c = Complex.Polar(1,animt/360*tpi*keyframes[i][1])['*'](keyframes[i][0])['*'](Complex(w/2/scale))
  circle(ctx, 0, 0, c.abs, 'black', false)
  ctx.beginPath()
  ctx.moveTo(0,0)
  ctx.translate(c.real, c.imag)
  ctx.lineTo(0,0)
  ctx.stroke()
  pctx.translate(c.real, c.imag)
}

I use a complex number library here, but it isn’t strictly necessary, just a convenience. animt is the current time step, and c represents the vector of that given circle at time animt. Afterwards, it’s just a matter of drawing the circle, and translating the canvas so that the origin of the canvas is at the end of the vector of the circle that was just drawn.

After drawing some grid lines and adding some basic GUI, the pretty fun to play with Fourier drawing application is done. If you’re still here, I applaud you for your determination, perserverance, and/or insanity.