# Payoff: The Same Curve, via the FFT *Task 5 of 5 · [Autocorrelation & Pitch](https://gpu.rocks/learn/autocorrelation-b159433f.md) · GPU.js Learn* 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 **Goal:** write the power-spectrum kernel, wire up the transform → power → inverse transform route, show its answer matches brute force, and time both. ## Requirements - `power` takes the two-plane spectrum and returns `output: [32768, 2]`: plane 0 is `re² + im²`, plane 1 is `0` - Wire up `transform(padded(), −1)` → `power` → `transform(…, +1)`, read the result back, and divide the real plane by `32768` - Log the largest disagreement with brute force: `console.log('max difference:', d)` - Warm up, then time five runs of each and log `'brute force:'` and `'fft route:'` in milliseconds ## Hint 1 — the power kernel The power spectrum is real, so plane 1 goes home empty: ```js 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 ```js 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: ```js 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. ## Starter code ```js // Two routes to the same curve. One is O(n²), the other O(n log n). const gpu = new GPU({ mode }); // ---- route 1: brute force, exactly as in task 1 -------------------------- // 16,384 threads, each sweeping all 16,384 samples: // 268,435,456 loop iterations for the whole curve. const correlate = gpu.createKernel(function (signal) { const lag = this.thread.x; let sum = 0; for (let i = 0; i < this.constants.n; i++) { if (i + lag < this.constants.n) { sum += signal[i] * signal[i + lag]; } } return sum; }, { output: [16384], constants: { n: 16384 }, }); // ---- route 2: one radix-2 butterfly pass, run 15 times ------------ // GIVEN, and the same two-plane layout as before: data[0][i] real, // data[1][i] imaginary, output [32768, 2] read back as result[plane][i]. // // LADDER is what makes this a transform ON the GPU rather than 31 trips to // it and back: // pipeline + immutable — a pass hands the next one a texture that stays in // device memory, so the ladder pays for ONE download, not 31. // optimizeFloatMemory — four floats to a texel, which is also what keeps a // 32768-wide output inside the 16384-texel texture limit. // correlate() downloads its answer once too, so the race below is fair. const LADDER = { pipeline: true, immutable: true, optimizeFloatMemory: true }; const fftPass = gpu.createKernel(function (data, ns, sign) { const i = this.thread.x; const block = Math.floor(i / (2 * ns)); const r = i - block * 2 * ns; const k = r % ns; const t = block * ns + k; const partner = t + this.constants.half; const angle = (sign * Math.PI * k) / ns; const wr = Math.cos(angle); const wi = Math.sin(angle); const ar = data[0][t]; const ai = data[1][t]; const br = data[0][partner]; const bi = data[1][partner]; const tr = br * wr - bi * wi; const ti = br * wi + bi * wr; if (r < ns) { if (this.thread.y === 0) return ar + tr; return ai + ti; } if (this.thread.y === 0) return ar - tr; return ai - ti; }, { output: [32768, 2], constants: { half: 16384 }, ...LADDER, }); // TODO: the one kernel you have to write. const power = gpu.createKernel(function (spec) { // plane 0 → re * re + im * im, plane 1 → 0 return 0; }, { output: [32768, 2], ...LADDER }); // A pipelined kernel returns a gpu.js Texture; the CPU backend has no textures // and returns the array itself. This is the one download either route makes. function readBack(result) { return result.toArray ? result.toArray() : result; } // The ladder: 15 passes, the stride doubling each time. function transform(data, sign) { let cur = data; for (let ns = 1; ns < 32768; ns *= 2) cur = fftPass(cur, ns, sign); return cur; } // 16,384 samples zero-padded to 32,768, plus an all-zero // imaginary plane. Float32Array from the start: gpu.js locks an argument's // container type on the first call. function padded() { const re = new Float32Array(32768); const im = new Float32Array(32768); for (let i = 0; i < signal.length; i++) re[i] = signal[i]; return [re, im]; } function viaFFT() { const out = new Float32Array(16384); // TODO: transform → power → inverse transform → readBack → divide by 32768. return out; } const brute = correlate(signal); const fast = viaFFT(); let maxDiff = 0; for (let lag = 0; lag < 16384; lag++) { maxDiff = Math.max(maxDiff, Math.abs(fast[lag] - brute[lag])); } console.log('max difference:', maxDiff, 'on a peak of', brute[0]); // TODO: warm both routes up, then time five runs of each and log // console.log('brute force:', ms, 'ms'); // console.log('fft route:', ms, 'ms'); ``` --- Interactive version: https://gpu.rocks/learn/autocorrelation-b159433f/5 [Previous task](https://gpu.rocks/learn/autocorrelation-b159433f/4.md)