# Same Answer, Two Centuries Apart *Task 5 of 6 · [The FFT Butterfly](https://gpu.rocks/learn/fft-butterfly-d4375da7.md) · GPU.js Learn* A permutation and a stack of index arithmetic produced *something*. The only way to know it is the Fourier transform is to compute the Fourier transform the way the definition says and compare, bin by bin. That is this task: write the naive DFT, run both at n = 8,192 — thirteen passes this time, since log₂(8,192) = 13 — and check they agree. Then count. The DFT visits every sample for every bin: n² = **67,108,864** complex multiplies. The FFT does log₂(n) passes of n/2 butterflies: **53,248**. A factor of 1,260 — and it is not a constant, it is a ratio that grows without bound. At n = 65,536, a window of about 1.5 seconds of CD audio, it is **8,192×**, which is the difference between a spectrogram that renders live and one that does not render at all. Gauss found the trick in 1805 and left it unpublished; Cooley and Tukey published it in 1965 and it went on to underwrite digital audio, JPEG's cousin the DCT, radio astronomy, MRI, and the modem you are reading this through. Now the stopwatch, and the reason this task is the size it is. **Kernel launches are not free.** A launch plus the readback around it costs about half a millisecond here, and the FFT needs fourteen of them where the definition needs one. At n = 512 that overhead *is* the transform, and the algorithm doing 114× less arithmetic loses the race outright. That is a real and permanent GPU lesson — a small transform is not worth sending to a GPU at all — but it is no longer this task's punchline, because the fix for it is the one change the code below makes to the FFT you built in task 4: every pass hands the next one a *texture* (`pipeline: true`) instead of a JavaScript array, so the spectrum crosses back to the CPU once instead of fourteen times. Measured, that single change takes the FFT from 8.2 ms to 1.5 ms. With both things settled — enough arithmetic that it outweighs the launches, and no copies that buy nothing — the race is not close. On the machine this was written on, in gpu mode: the FFT in **1.5 ms**, the definition in **3.5 ms**. Run it and read your own two numbers off the console. Then switch **Mode** from Auto to CPU, where there is no launch overhead left and nothing but the operation count decides it, and watch the gap open past **35×** — against a naive transform that has been handed only a 1,024-bin slice of the spectrum instead of all 8,192, because the whole thing is two seconds of a single thread. A machine with no GPU at all, running WebGL in software, reports about 135×. One last thing in the console, which is not about speed. The two answers agree to about six parts in ten thousand, and nearly all of that gap is the *naive* one being wrong: 8,192 float32 terms summed into one running total, against thirteen roundings for the FFT. The fast algorithm is also the accurate one, and that is not a coincidence. Measuring Speed Honestly is the module that owns the benchmarking argument; the ⏱ **Benchmark** chip answers a different question again, timing the whole file at once rather than either transform alone. ## Figures - **one pass over a square, or thirteen passes over a line** ## Goal **Goal:** write the naive DFT kernel, confirm it agrees with the FFT bin by bin, and log the two operation counts and both timings. ## Requirements - The DFT kernel is one thread per output cell, looping over all `this.constants.n` samples - The output is `bins` wide but every bin sums over all `n` samples — the output width and the loop bound are different numbers - Angles are `-2π·k·t / n` with `k = this.thread.x` — no scaling by `n` - Compare the two spectra bin by bin and `console.log` the verdict - `console.log` both operation counts and both elapsed times (already wired up) ## Hint 1 — which index is the loop, which is the thread? `this.thread.x` is the bin `k`, fixed for the whole thread. The loop variable `t` walks the samples. So the array lookup is `x[t]` — if `this.thread.x` appears inside the brackets, the two have been swapped. ## Hint 2 — the accumulation ```js const angle = (-2 * Math.PI * k * t) / this.constants.n; re += x[t] * Math.cos(angle); im += x[t] * Math.sin(angle); ``` — the input is real, so there is no fourth term. Task 6 needs the general version. ## Hint 3 — what "agree" should mean Not equality. The two computations do wildly different numbers of float32 operations and will not land on the same bits — the DFT accumulates 8,192 terms into one running total, the FFT thirteen. Compare the largest difference against the size of the spectrum, not against zero: `worst / peak < 5e-3` is a generous margin and a truthful one. Expect about `6e-4`, and expect nearly all of it to belong to the DFT: its phase `-2π·k·t / n` runs `k·t` up past 67 million, which float32 cannot even hold to the nearest whole number. ## Same idea elsewhere Checking a fast implementation against the slow definition on a small input is how every FFT library is tested — cuFFT, FFTW and rocFFT all ship exactly this comparison in their accuracy suites, usually reported as relative error against a reference computed at higher precision. It is also the honest answer to "is my GPU port correct": not "the pictures look the same", but the definition, at a size where you can afford it. The second habit here travels just as far: when a chain of kernels loses to a brute-force one, count the round trips before you touch the arithmetic. Keeping intermediates on the device is what `pipeline: true` is for in gpu.js, what CUDA graphs and Metal command buffers are for elsewhere, and it is very often the whole difference. ## Starter code ```js // The definition versus the algorithm. Same answer, different bill. const gpu = new GPU({ mode }); const n = 8192; // log2(8192) = 13 butterfly passes // The naive transform reads the WHOLE signal once for every bin it produces, so // a full spectrum is n * n work. A GPU eats that; one CPU thread needs about two // seconds of it. So on the CPU backend the definition is asked for a 1,024-bin // slice instead — an 8x head start, and it loses anyway. const bins = mode === 'gpu' ? n : 1024; // ---- the fast one, complete (tasks 3 and 4), with ONE change // Every pass hands the next one a TEXTURE instead of a JavaScript array. Copying // the spectrum off the card after each of thirteen passes costs far more than the // passes do, and that — not the arithmetic — is why a small FFT loses this race. // Two butterfly kernels take turns because a pipelined kernel writes into its own // texture and so cannot also be the one reading it. const stage = { output: [n, 2], pipeline: true }; const scramble = gpu.createKernel(function (x) { let v = this.thread.x; let reversed = 0; for (let b = 0; b < this.constants.bits; b++) { reversed = reversed * 2 + (v % 2); v = Math.floor(v / 2); } if (this.thread.y === 0) return x[reversed]; return 0; }, { ...stage, constants: { bits: 13 } }); const onePass = function (spectrum, half) { const i = this.thread.x; const j = Math.floor(i / half) % 2; const base = i - j * half; const r = i % half; const angle = (-Math.PI * r) / half; const wr = Math.cos(angle); const wi = Math.sin(angle); const ar = spectrum[0][base]; const ai = spectrum[1][base]; const br = spectrum[0][base + half]; const bi = spectrum[1][base + half]; const tr = wr * br - wi * bi; const ti = wr * bi + wi * br; if (this.thread.y === 0) { if (j === 0) return ar + tr; return ar - tr; } if (j === 0) return ai + ti; return ai - ti; }; const evenPass = gpu.createKernel(onePass, stage); const oddPass = gpu.createKernel(onePass, stage); function fft(samples) { let buffer = scramble(samples); let pass = 0; for (let half = 1; half < n; half *= 2) { buffer = (pass++ % 2 === 0 ? evenPass : oddPass)(buffer, half); } // one read back, at the very end. The CPU backend has no textures and has // handed back a plain array already, so there is nothing to convert there. return buffer.toArray ? buffer.toArray() : buffer; } // ---- the slow one const dft = gpu.createKernel(function (x) { const k = this.thread.x; let re = 0; let im = 0; for (let t = 0; t < this.constants.n; t++) { // TODO: accumulate x[t] · e^(-2πi·k·t / n) into re and im. } if (this.thread.y === 0) return re; return im; }, { output: [bins, 2], constants: { n } }); // Warm both up before timing: the first call compiles a shader, and a GPU handed // its first work of the day is still spinning its clocks up. fft(signal); dft(signal); fft(signal); dft(signal); const t0 = performance.now(); const fast = fft(signal); const t1 = performance.now(); const slow = dft(signal); const t2 = performance.now(); let worst = 0; let peak = 0; for (let k = 0; k < bins; k++) { worst = Math.max(worst, Math.abs(fast[0][k] - slow[0][k]), Math.abs(fast[1][k] - slow[1][k])); peak = Math.max(peak, Math.hypot(slow[0][k], slow[1][k])); } console.log('bins the definition was asked for:', bins, 'of', n); console.log('DFT complex multiplies:', bins * n); console.log('FFT complex multiplies:', (n / 2) * Math.log2(n)); console.log('worst bin difference:', worst, 'against a peak of', peak); console.log('agree:', worst / peak < 5e-3); console.log('fft ms:', (t1 - t0).toFixed(3)); console.log('dft ms:', (t2 - t1).toFixed(3)); console.log('the definition took', ((t2 - t1) / (t1 - t0)).toFixed(1), 'times as long'); ``` --- Interactive version: https://gpu.rocks/learn/fft-butterfly-d4375da7/5 [Previous task](https://gpu.rocks/learn/fft-butterfly-d4375da7/4.md) · [Next task](https://gpu.rocks/learn/fft-butterfly-d4375da7/6.md)