# Where the Saving Comes From

*Task 1 of 6 · [The FFT Butterfly](https://gpu.rocks/learn/fft-butterfly-d4375da7.md) · GPU.js Learn*

The discrete Fourier transform asks one question n times: *how much of the
signal is a wave that fits exactly k times into the window?* The answer is a sum over
every sample, so n answers cost n² terms:

```js
X[k] = Σ  x[t] · e^(-2πi·k·t / n)      t = 0 … n-1
```

At n = 512 that is 262,144 sine-cosine terms per transform, and the count grows with the
square. Now split the sum by the **parity of t** — the even-indexed samples in
one pile, the odd-indexed in the other. Pull the shared factor out of the odd pile and
something remarkable is left standing:

```js
X[k]       = E[k] + W·O[k]        W = e^(-2πi·k / n)
X[k + n/2] = E[k] - W·O[k]
```

where `E` and `O` are the (n/2)-point DFTs of the even and the odd
samples. Two half-length transforms and one multiplication give you the whole thing — and
*they* can be split the same way, and so on down. That recursion is the entire
algorithm: n² becomes n·log₂n. Do not take that on faith. Check it.

Complex numbers arrive here as **two planes of floats**, the convention every
module in this track shares: a kernel with `output: [n, 2]` is indexed
`result[p][i]`, plane `0` real and plane `1` imaginary.
`this.thread.y` tells a thread which plane it owns.

## Goal

**Goal:** compute the two 4-point half-transforms of `signal`
with one kernel, combine them by hand, and confirm the result is the 8-point DFT the
`dft8` kernel already computes.

## Requirements

- One `halfDft` kernel, called twice — `offset 0` for the even samples, `offset 1` for the odd
- It reads every SECOND sample: element `t` is `x[2 * t + offset]`
- Its angles divide by 4, not 8 — it is a 4-point transform
- Combine with `W = e^(-2πi·k/8)` and `console.log` the largest disagreement

## Hint 1 — which samples are mine?

The even pile is `x[0], x[2], x[4], x[6]`; the odd pile is
`x[1], x[3], x[5], x[7]`. Both are "every second sample", differing only in
where they start — so element `t` of the pile is
`x[2 * t + offset]`, and `offset` is an ordinary kernel
argument.

## Hint 2 — multiplying two complex numbers

There is no complex type, so you write it out. With
`w = wr + i·wi` and `o = or + i·oi`:

```js
// the minus lives in the REAL part
const tr = wr * or - wi * oi;
const ti = wr * oi + wi * or;
```

## Hint 3 — the check

For each `k` in 0…3, all four of these have to hold:

```js
E[0][k] + tr  →  X[0][k]
E[1][k] + ti  →  X[1][k]
E[0][k] - tr  →  X[0][k + 4]
E[1][k] - ti  →  X[1][k + 4]
```

Four bins of `E` and `O` reconstruct all eight bins of
`X`. That is the saving, in miniature.

## Same idea elsewhere

Every FFT in the world is this recursion made iterative — FFTW, cuFFT, rocFFT,
vDSP, Metal Performance Shaders. The radix varies (4, 8 and split-radix are common, and
cuFFT picks one per size), and Bluestein's algorithm rescues the sizes that will not
factor, but the move is always the same: turn one transform into several smaller ones plus
a twiddle.

## Starter code

```js
// Two 4-point transforms and one multiplication should rebuild all 8 bins.
const gpu = new GPU({ mode });

// The definition, written out: 8 bins × 8 samples = 64 terms.
const dft8 = gpu.createKernel(function (x) {
  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 += x[t] * Math.cos(angle);
    im += x[t] * Math.sin(angle);
  }
  if (this.thread.y === 0) return re;
  return im;
}, { output: [8, 2], constants: { n: 8 } });

const halfDft = gpu.createKernel(function (x, offset) {
  const k = this.thread.x;
  let re = 0;
  let im = 0;
  for (let t = 0; t < this.constants.m; t++) {
    const angle = (-2 * Math.PI * k * t) / this.constants.m;
    // TODO: read every SECOND sample, starting at `offset`.
    const v = x[t];
    re += v * Math.cos(angle);
    im += v * Math.sin(angle);
  }
  if (this.thread.y === 0) return re;
  return im;
}, { output: [4, 2], constants: { m: 4 } });

const X = dft8(signal);
const E = halfDft(signal, 0);
const O = halfDft(signal, 1);
console.log('E real:', E[0]);
console.log('O real:', O[0]);

let worst = 0;
for (let k = 0; k < 4; k++) {
  const angle = (-2 * Math.PI * k) / 8;
  const wr = Math.cos(angle);
  const wi = Math.sin(angle);

  // TODO: t = W · O[k] — a complex multiply, two lines.
  const tr = 0;
  const ti = 0;

  worst = Math.max(worst,
    Math.abs(E[0][k] + tr - X[0][k]), Math.abs(E[1][k] + ti - X[1][k]),
    Math.abs(E[0][k] - tr - X[0][k + 4]), Math.abs(E[1][k] - ti - X[1][k + 4]));
}

console.log('largest disagreement:', worst);
console.log('identity holds:', worst < 1e-3);
```

---

Interactive version: https://gpu.rocks/learn/fft-butterfly-d4375da7/1

[Next task](https://gpu.rocks/learn/fft-butterfly-d4375da7/2.md)