Task 5 of 6

Same Answer, Two Centuries Apart

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.

one pass over a square, or thirteen passes over a line

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

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.

All tasks in The FFT Butterfly

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