# Multiplying Spectra Convolves in a Circle

*Task 5 of 5 · [Filtering in the Frequency Domain](https://gpu.rocks/learn/frequency-filtering-8c225e10.md) · GPU.js Learn*

Task 1 had an accident built into it: the signal occupied 20 of 64 samples and the
filter 5, so the 24-sample answer had room to spare. Take the room away and the theorem
bites.

Multiplying two length-`n` spectra does not produce the convolution you want.
It produces the **circular** convolution — the one that treats both buffers as
loops, so whatever runs off the end comes back on at the beginning:

```js
circular[i] = linear[i] + linear[i + n]
```

Here `sig` is 32 samples and `filt` carries 9 taps, so the honest
answer is `32 + 9 − 1 = 40` samples long. Multiply their 32-bin spectra and the
last 8 samples fold back onto the first 8, corrupting them by as much as
**0.37** against a signal whose largest value is 0.66. The signature is
unmistakable once you have seen it: *the tail of the answer appears at its head*.

The fix is to give the answer somewhere to live. Zero-pad both buffers to at least
`n + m − 1` before transforming — 40 here, and since every real FFT wants a
power of two, **64**. The overhang then lands on padding instead of on your
data. The starter transforms at 32 and prints the damage; fix `PAD` and the
guard inside `pad`, and the error should collapse to float32 rounding — a few
millionths on a GPU, a few hundred-thousandths on a software renderer.

This is also where the cost story gets honest. Padding to a power of two means
transforming 64 points to convolve 32 with 9 — and 9 taps is short enough that the sliding
window wins outright, on any hardware. Frequency-domain convolution starts paying somewhere
around a few dozen taps, and exactly where depends on the machine. The Benchmark button will
price this pipeline for you on both backends; Measuring Speed Honestly owns the rest of the
methodology, including why a single timing is a rumour.

## Figures

- **eight samples with nowhere to go, and they do not go quietly**

## Goal

**Goal:** pick a `PAD` long enough to hold the whole answer,
and finish `pad` so it copies the first 32 samples and zeroes the rest.

## Requirements

- Set `PAD` to at least `32 + 9 − 1 = 40`; round up to `64`, the next power of two
- Inside `pad`, return `src[this.thread.x]` only while `this.thread.x` is below `this.constants.src`, and `0` beyond it
- Guard that read: past the end of the source there is nothing to read, and gpu.js answers `NaN` on the CPU backend and garbage on WebGL
- Leave the rest of the pipeline alone — the transforms and the complex multiply are task 1's, unchanged

## Hint 1 — how long is the answer?

Convolving `n` samples with `m` taps produces
`n + m − 1` samples: the filter's last tap is still hanging over the
signal's last sample `m − 1` steps after the signal has ended. With 32 and
9 that is 40. Any transform length of 40 or more works; 64 is what you would use in
practice, because that is what an FFT wants.

## Hint 2 — the guard

Two lines, and the `if` matters — a ternary would still describe an
out-of-range read to the GLSL compiler:

```js
if (this.thread.x < this.constants.src) {
  return src[this.thread.x];
}
return 0;
```

`filt` already arrives in a 32-sample buffer with 9 taps and 23 zeros, so
the same kernel pads it correctly with no special case.

## Hint 3 — what wrapping looks like

Before you fix it, compare the first eight samples the starter prints against
the last eight of `direct`. They are the same numbers. That is not a
coincidence and it is not rounding — it is the tail of the convolution, delivered to
the front of the buffer.

## Same idea elsewhere

Every FFT convolution in production carries this bookkeeping. cuFFT, FFTW and
`scipy.signal.fftconvolve` all pad internally and hand back the
`n + m − 1` answer; `numpy.fft` does not, which is why the "why is my
filtered audio clicking at the block boundaries" question outlives every generation of
programmers. Streaming filters solve the same problem with overlap-add or overlap-save:
pad each block, convolve, and add the overhang into the next block instead of throwing it
at the start of this one.

## Starter code

```js
// Zero-padding, or: where does the answer live?
const SRC = 32;   // both source buffers are 32 samples long; filt has 9 taps

// TODO: SRC is exactly the source length, which leaves the answer nowhere to
// go. How long is the convolution of 32 samples with 9 taps?
const PAD = SRC;

const gpu = new GPU({ mode });

// GIVEN — the forward transform (see task 1).
const dft = gpu.createKernel(function (sig) {
  const k = this.thread.x;
  let re = 0;
  let im = 0;
  for (let t = 0; t < this.constants.n; t++) {
    const angle = (-2 * Math.PI * k * t) / this.constants.n;
    re += sig[t] * Math.cos(angle);
    im += sig[t] * Math.sin(angle);
  }
  if (this.thread.y === 0) return re;
  return im;
}, { output: [PAD, 2], constants: { n: PAD } });

// GIVEN — the inverse.
const idft = gpu.createKernel(function (spec) {
  const t = this.thread.x;
  let re = 0;
  let im = 0;
  for (let k = 0; k < this.constants.n; k++) {
    const angle = (2 * Math.PI * k * t) / this.constants.n;
    const c = Math.cos(angle);
    const s = Math.sin(angle);
    re += spec[0][k] * c - spec[1][k] * s;
    im += spec[0][k] * s + spec[1][k] * c;
  }
  if (this.thread.y === 0) return re / this.constants.n;
  return im / this.constants.n;
}, { output: [PAD, 2], constants: { n: PAD } });

// GIVEN — task 1's complex multiply, unchanged.
const mulSpectra = gpu.createKernel(function (a, b) {
  const i = this.thread.x;
  const ar = a[0][i];
  const ai = a[1][i];
  const br = b[0][i];
  const bi = b[1][i];
  if (this.thread.y === 0) return ar * br - ai * bi;
  return ar * bi + ai * br;
}, { output: [PAD, 2] });

// YOUR JOB — copy the source into a longer buffer and zero the rest.
const pad = gpu.createKernel(function (src) {
  // TODO: only the first this.constants.src samples exist. Past that there
  // is nothing to read — return 0 instead of reading off the end.
  return src[this.thread.x];
}, { output: [PAD], constants: { src: SRC } });

const result = idft(mulSpectra(dft(pad(sig)), dft(pad(filt))));

let worst = 0;
const compared = Math.min(40, result[0].length);
for (let i = 0; i < compared; i++) worst = Math.max(worst, Math.abs(result[0][i] - direct[i]));
console.log('worst error vs the direct convolution:', worst.toFixed(9));
console.log('first 8 of result:', Array.from({ length: 8 },
  (unused, i) => Number(result[0][i].toFixed(3))));
console.log('last 8 of direct: ', Array.from({ length: 8 },
  (unused, i) => Number(direct[32 + i].toFixed(3))));
```

---

Interactive version: https://gpu.rocks/learn/frequency-filtering-8c225e10/5

[Previous task](https://gpu.rocks/learn/frequency-filtering-8c225e10/4.md)