How Fourier works

To recap, we take our inputs and throw those into Fourier. We get a bunch of complex numbers out, of which we throw away the second half (because the second half is the same as the first half, but then inverted). Using (im,re) as a vector we calculate the direction and the strength. We do still multiply the strength by two, because we threw away half our results.

Now let's look into what Fourier actually does. Let's first look at how the implementation is laid out. It's laid out in two nested for loops. These are:

• For each frequency
  • Initialize something
  • For each sample
    • Do something
  • Add result to output array

So for every single possible outcome frequency we do something for every single sample. And then after doing this, we throw away half the samples. Now that sounds useless, and unless you're doing really complex stuff with Fourier (in which case you already know more than I could ever tell you) it's safe to ignore. So let's put a note here to ourselves when we get to making Fourier faster in a next post, we shouldn't calculate it for more than half the frequencies.

Within the outer for loop we initialize a complex number, with im=0 and re=0. After the inner loop we add the result to the output array. Nothing interesting here.

Now the inner loop is where the magic happens. We grab the amplitude at the current index, or the current measurement, and put it into a variable. After that we take what we call the current rotation angle for this point in time for the frequency we're checking. The best way to understand this is with an example.

With 20 measurement points, at a frequency of 2, we'll "rotate" a full circle up and down twice. Let's assume we're currently calculate at measurement 5.
-1 * (2 * Math.PI) * frequency * (timer / N) becomes
-2PI * 2 * 5/20, or -2PI * (2 * 0.25), which becomes -PI.
This means that over the course of 20 measurements, our rotational angle goes from 0 to -PI to -2PI to -3PI to -4PI, or 2 full circles, which correlates to a frequency of two again.

So now we've got the amplitude (which was the current measurement) and the rotation angle. Next we throw the rotation angle into a complex number, using cos(angle) as the real part and sin(angle) as the imaginary part. We then multiply it with a complex number with an imaginary part of 0, and a real part of the amplitude. The new real part becomes the real parts multiplied minus the imaginary parts multiplied. As the imaginary part is 0, this means the new real part simply becomes the real parts multiplied. The new imaginary part becomes it's own real part multiplied by the other imaginary part plus it's own imaginary part multiplied by the other real part. Striping away anything multiplied by 0 turns this into it's own imaginary part times the other real part. To summarize, this multiplication by the current amplitude simply multiplies both the real and the imaginary part, so both cos(angle) and sin(angle), by the current amplitude. As such, this code could also have been written as such:

const dataPointContribution = new ComplexNumber({
    re: Math.cos(rotationAngle) * currentAmplitude,
    im: Math.sin(rotationAngle) * currentAmplitude
});

We then add this contribution to the sum for this given frequency. Adding two complex numbers is simply adding their real parts and their imaginary parts.

After this we still divide it by the amount of inputs, basically calculating an average. Division with complex numbers looks a lot more complicated in the code we started off with, so let's see what's actually going on. First I'll repeat the code as it was when we started off.

  divide(divider) {
    const complexDivider = this.toComplexNumber(divider);
    const dividerConjugate = this.conjugate(complexDivider);
    const finalDivident = this.multiply(dividerConjugate);
    const finalDivider = (complexDivider.re ** 2) + (complexDivider.im ** 2);
    return new ComplexNumber({
      re: finalDivident.re / finalDivider,
      im: finalDivident.im / finalDivider,
    });
  }
  conjugate(number) {
    const complexNumber = this.toComplexNumber(number);
    return new ComplexNumber({
      re: complexNumber.re,
      im: -1 * complexNumber.im,
    });
  }

Now again, we'll be dividing by a simple number. That means dividing by a complex number that only has a real part. Using this knowledge we can simplify the code a lot for this specific use-case, to this:

  divide(divider) {
    const complexDivider = {im:0,re:divider};
    const dividerConjugate = complexDivider;
    const finalDivident = this.multiply(dividerConjugate);//So {im:this.im*divider, re:this.re*divider}
    const finalDivider = (complexDivider.re ** 2) + (0);
    return new ComplexNumber({
      re: finalDivident.re / finalDivider,
      im: finalDivident.im / finalDivider,
    });
  }

Now lets simplify this a bit further.

  divide(divider) {
    return new ComplexNumber({
      re: this.re * divider / (divider ** 2),
      im: this.im * divider / (divider ** 2),
    });
  }

Okay, one more step to go.

  divide(divider) {
    return new ComplexNumber({
      re: this.re / divider,
      im: this.im / divider,
    });
  }

What's interesting to note is that although Fourier is based on logic found in complex numbers, when implementing you have absolutely nothing to do with any complex number logic.

Using all our findings, let's rewrite our code a bit to something that's easier to play around with. And while we're at it, let's also make it so that we get out a frequency and offset radian right away.

function dft(inputAmplitudes, iCloseToZeroTreshold) {
  const N = inputAmplitudes.length;
  const signals = [];
  for (let frequency = 0; frequency < N; frequency += 1) {
    let frequencySignal = {cos:0,sin:0};
    
    for (let timer = 0; timer < N; timer += 1) {
      const currentAmplitude = inputAmplitudes[timer];
      const rotationAngle = -1 * (2 * Math.PI) * frequency * (timer / N);
      frequencySignal.cos += Math.cos(rotationAngle) * currentAmplitude;
      frequencySignal.sin += Math.sin(rotationAngle) * currentAmplitude;
    }
    frequencySignal.cos /= N;
    frequencySignal.sin /= N;
    var iRad = Math.atan2(e.cos, e.sin);
    var iStrength = e.cos / Math.sin(iRad) * 2;
    if (iStrength >= iCloseToZeroTreshold){
      signals.push({iRad:iRad,iStrength:iStrength});
    }
  }
  return signals;
}

Now this allows for another optimization, which is moving the division to the strength calculation.

function dft(inputAmplitudes, iCloseToZeroTreshold) {
  const N = inputAmplitudes.length;
  const signals = [];
  for (let frequency = 0; frequency < N; frequency += 1) {
    let frequencySignal = {cos:0,sin:0};
    
    for (let timer = 0; timer < N; timer += 1) {
      const currentAmplitude = inputAmplitudes[timer];
      const rotationAngle = -1 * (2 * Math.PI) * frequency * (timer / N);
      frequencySignal.cos += Math.cos(rotationAngle) * currentAmplitude;
      frequencySignal.sin += Math.sin(rotationAngle) * currentAmplitude;
    }
    var iRad = Math.atan2(frequencySignal.cos, frequencySignal.sin);
    var iStrength = frequencySignal.cos / Math.sin(iRad) * 2 / N;
    if (iStrength >= iCloseToZeroTreshold){
      signals.push({iRad:iRad,iStrength:iStrength});
    }
  }
  return signals;
}

Okay, that's a whole lot more readable. But do we understand it yet? I don't. Let's look at what actually causes this to work by plotting it out on a canvas on a small number of samples with a simple formula, sin(x), for only the correct frequency.

What we see here is that the average of this is an arrow pointing straight up. Now let's see what happens if we'd try to measure frequency 2, which is not in the signal.

Now we've got arrows pointing in all directions, with the same strength, which ends up being a vector of (0,0). So this example shows that attempting to detect the wrong signal ends up canceling out, while when detecting the correct signal the vectors amplify each other.

Now let's give this another shot, but with a different signal.

Now we see that all vectors cancel each other out again. Let's detect the correct signal.

This amplifies again. I'm starting to understand why Fourier outputs the correct frequencies, but not the wrong ones, and I hope you are too. But lets see what happens with a more complicated signal.

Looking at these plots I'm noticing something. It seems as if these plots are just the sums of the plots with the individual simpler signals. Let me plot them together. Red will be Fourier with the signal Math.sin(i), blue will be Math.sin(2*i), and black will be Math.sin(i)+Math.sin(i*2).

Okay, apparently it's not the sum, it's the average. But it still explains exactly how Fourier works. Let's still see what happens if we try to measure a signal that is't in there.

And here we see that some signals already cancel themselves out, and the total of all of them resolves to the vector (0,0) again.

From this, we can draw a general conclusion of how Fourier works, and why it works.

Fourier works by generating the output signal to be tested, and multiplying this by the measurements at the same position. Given no start offset (e.g. sin(5x)*3) these two signals will multiply in a way that they strengthen each other, in a sense resonating with each other. When measuring with an input signal of sin(5x) this will give us the output vector (0,3).
Given a start offset (e.g. sin(x+0.5*PI)*3, shifting the signal 90 degrees) the aforementioned will still apply, but instead we'll get out an output vector (3,0). The direction of the vector will indicate the "start offset", while the length of the vector will indicate the strength of the signal.