Task 5 of 5

Payoff: The Same Curve, via the FFT

Every lag of the correlation costs a pass over the signal, so the whole curve costs n passes over n samples: O(n²), with a large constant. There is a way out, and it is one of the prettiest results in the subject: the Wiener–Khinchin theorem. The autocorrelation of a signal is the inverse transform of its power spectrum. Transform, square the magnitudes, transform back — O(n log n), and the whole correlation comes out at once.

Two details make it true rather than nearly true. First, the transform believes your signal repeats forever, so a shifted copy wraps around from the end onto the beginning; the cure is to zero-pad to twice the length, which is why 16,384 samples go into a 32,768-point transform. Second, the inverse transform of this library — like most — leaves a factor of 32,768 behind, so divide by it at the end.

This task is bigger than the other four, and that is the whole point. Tasks 1–4 analysed one 512-sample window, because that is what pitch detection actually does. Here the buffer is 16,384 samples — two seconds of that same plucked string — and the question is the whole correlation, every one of its 16,384 lags. That is 16,384 threads each sweeping all 16,384 samples — 268,435,456 loop iterations, of which the 134,225,920 inside the shrinking overlap do a multiply-add — against 31 kernel launches over 32,768 points. At this size the asymptotics have stopped being a promise about large n and started being the thing you are waiting for.

One thing this route is not is fragile. The power spectrum throws phase away, so nothing here ever calls Math.atan2 — and a bin whose imaginary part is near zero has a phase that flips a whole turn on float32 noise. Squaring magnitudes has no such seam, which is why a 32,768-point round trip through 31 kernel launches still tracks the brute-force curve to better than one part in 10,000 — measured at one part in 86,000 on the hardware these notes were written on.

The FFT itself is written for you below — The FFT Butterfly derives that pass properly; here it is a given — and its shape should look familiar even so: a JavaScript loop calling one kernel 15 times with a doubling stride, which is the halving ladder from Reductions read backwards. Read its settings, though, because one of them is the difference between a GPU transform and 31 round trips to the GPU and back.

Then press the point home with a stopwatch — warm up first, and read the numbers the way Measuring Speed Honestly insists on. The direct route is not a strawman: it is one launch, it is perfectly parallel, and below a few thousand samples it wins, because 31 launches cost more than one until there is enough work to pay for them. Find out what happens when there is.

Two planes, one complex array. gpu.js has no complex number, so this track carries one as two planes of floats: a kernel declares output: [n, 2] and the result reads result[0][i] for the real part and result[1][i] for the imaginary part. Going the other way, a complex argument is a nested [2][n] array read as data[0][i] and data[1][i]. It is the same trick Optical Flow uses to return two flow components from one kernel — output: [w, h] is indexed [y][x], so [n, 2] is indexed [plane][i].

Goal: write the power-spectrum kernel, wire up the transform → power → inverse transform route, show its answer matches brute force, and time both.

Requirements

Hint 1 — the power kernel

The power spectrum is real, so plane 1 goes home empty:

const i = this.thread.x;
if (this.thread.y === 0) {
  const re = spec[0][i];
  const im = spec[1][i];
  return re * re + im * im;
}
return 0;

Note re * re + im * im and not Math.sqrt(...). The theorem wants the magnitude squared; a square root here produces a curve that is smooth, plausible and wrong.

Hint 2 — four lines of wiring
const spec = transform(padded(), -1);
const p = power(spec);
const back = readBack(transform(p, +1));
for (let lag = 0; lag < 16384; lag++) {
  out[lag] = back[0][lag] / 32768;
}

readBack is there because the ladder is pipelined: every pass hands the next one a texture that never leaves the GPU, and only the last result is downloaded. Drop it and you are indexing a Texture object, which has no [0]. If instead every value comes out exactly 32768 times too big, the division is the piece that is missing.

Hint 3 — timing it honestly

The first call to anything here compiles a shader, so throw it away:

correlate(signal); viaFFT();          // warm up, discard

let t0 = Date.now();
for (let i = 0; i < 5; i++) correlate(signal);
console.log('brute force:', (Date.now() - t0) / 5, 'ms');

Both routes end in exactly one download — correlate returns an array, viaFFT calls readBack once — so neither is being timed with its homework left on the GPU. The Benchmark button will still disagree with your own numbers, and it is not wrong: it times every kernel separately with a forced readback after each, which for a 15-pass ladder measures the ladder taken apart.

Same idea elsewhere

Convolution and correlation via the transform is the reason cuFFT, vkFFT and Apple's vDSP exist, and the reason every deep-learning framework once shipped an FFT-based convolution path. Both halves of what you just measured are load-bearing: the crossover is real, so small filters stay direct (below a few hundred taps the direct form wins on hardware that likes dense arithmetic), and so is the rule that got the transform route across it — a multi-pass GPU algorithm keeps its intermediates in device memory. Every one of those libraries batches its passes for the same reason, and a CUDA programmer calling cudaMemcpy between kernels is making the mistake this task's pipeline flag exists to avoid.

All tasks in Autocorrelation & Pitch

  1. Compare a Signal With Itself
  2. From Lag to Pitch
  3. Three Ways to Get It Wrong
  4. The Note That Isn't There
  5. Payoff: The Same Curve, via the FFT

This page is an interactive exercise — the editor, the GPU runner and your saved progress need JavaScript. The text above is the full brief.