Task 4 of 5
Play a note through a telephone, or a small speaker, or a church organ's mixture stop, and something strange happens: the fundamental is gone — the hardware simply cannot move that slowly — and you hear it anyway. Your auditory system is not reading the loudest peak off a spectrum. It is finding the period the partials share.
signal below is built from harmonics 2, 3 and 4 of 128 Hz — that is
256, 384 and 512 Hz — and there is no energy at 128 Hz at all. Not
attenuated: absent, by construction. Take its spectrum and the loudest bin is
256 Hz, an octave above the note. Take its autocorrelation and the first peak is at
lag 64, which is 128 Hz, because 256, 384 and 512 all come back into step every 64
samples whether or not anything is oscillating at that rate.
So you need a spectrum to accuse. Write one — a direct transform, one thread per bin, which is the O(n²) way and perfectly fine at 512 samples. The DFT, Honestly has a great deal more to say about what a transform is for and what it costs; here it is only the witness for the prosecution.
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].
output: [512, 2] — one thread per (plane, bin); this.thread.x is the bin, this.thread.y is the planethis.constants.n samples with angle −2π · bin · i / n: plane 0 accumulates signal[i] · cos(angle), plane 1 signal[i] · sin(angle)Math.hypot(re, im) per bin and find the loudest one below bin 256; its frequency is bin * sampleRate / 512console.log('loudest partial:', hz, 'Hz') and console.log('autocorrelation:', hz, 'Hz')Each thread owns a single output cell, so it computes a single sum — its own plane's, for its own bin:
const bin = this.thread.x;
let re = 0;
let im = 0;
for (let i = 0; i < this.constants.n; i++) {
const angle = (-2 * Math.PI * bin * i) / this.constants.n;
re += signal[i] * Math.cos(angle);
im += signal[i] * Math.sin(angle);
}
if (this.thread.y === 0) return re;
return im;
Both sums get computed by both threads and one of them is thrown away. That is wasteful and it is also the shape gpu.js gives you — and at 512 bins the whole transform is one launch, so you will not notice. The next task runs the same idea over a buffer 32 times longer, which is where the O(n²) stops being a remark and becomes the thing you are waiting for.
The result is two rows, not two columns:
const spec = spectrum(signal);
const mag = Math.hypot(spec[0][bin], spec[1][bin]);
Plane first, bin second — output: [n, 2] is indexed [y][x]
like any other 2D output.
Watch the two conversions in this task, because they are reciprocals of each
other and mixing them up is the classic slip. A lag is a period, so
hz = sampleRate / lag. A bin is already a frequency count, so
hz = bin * sampleRate / n. With 8192 Hz over 512 samples each bin
is 16 Hz wide, and the answer you are looking for is bin 16.
float2 buffer in CUDA (which is exactly what cuFFT's
cufftComplex is), half2 pairs in Metal. The layout question —
interleaved [re, im, re, im…] or planar [re…][im…] — is the same
one everywhere, and planar is what coalesces.
This page is an interactive exercise — the editor, the GPU runner and your saved progress need JavaScript. The text above is the full brief.