021: Fourier’s Window

2026-10-18

Last time a wave went through a gap in a wall and made a pattern on the far wall. A single slit made a broad hump. Two slits made fringes. It looked as if each gap had its own pattern and each one needed working out.

They are all one rule. The pattern a long way behind a gap is the Fourier transform of the gap. A Fourier transform takes a picture and says how much of each ripple it is made of: a few broad ripples and a few fine ones, in every direction. A gap is a picture too, open where the wave gets through and shut where it does not. The wave leaves the gap in all directions and each direction sees a different ripple of the gap.

The picture has changed from last time. Then we looked down on the tank from above and the wall was a line with gaps in it. Now we look at the wall head-on, as a wave coming towards us would see it. The wall is a flat sheet and the gap is a shape cut out of it. The slits are tall bars and the pattern falls on a flat screen a long way behind, so it spreads in two directions and not one.

On the left below is the gap in the wall, white where it is open. On the right is the pattern on the screen, brightest in the middle. Paint on the gap with the mouse and the pattern follows.

Each point in the pattern stands for one ripple. Where it is says how many crests of the ripple fit across the gap and which way they run. How bright it is says how well the gap lines up with the ripple. Move the mouse over the pattern and the ripple for that point appears on the gap, red for the crests and blue for the troughs. It is not a wave arriving at the wall. It is a test pattern, the thing the gap is being compared with, and it pulses in place so that the crests and troughs swap. At a bright point the crests reach every open part of the gap together, so the open parts flash red and blue in step and the sum swings far. At a dark point some open parts are on crests while others are on troughs, they cancel and the sum stays near zero.

Try these:

• Two slits. The same fringes as last time, now with the rest of the pattern. The slits are tall and narrow, so the pattern is spread sideways and squeezed up and down. Move the mouse to a bright fringe and both slits flash red and blue together. Move it to the dark band between and one slit is red while the other is blue.
• Paint over half of one slit, or make a slit fatter. Whatever gets bigger in the gap gets smaller in the pattern. A wide gap makes a tight pattern and a narrow one a broad pattern.
• Circle. A round gap makes rings. It is how a telescope blurs a star.
• Grating. Five slits instead of two. The fringes get sharper.
• Clear, then draw anything you like.

viewof shape = {
  const ui = await import(new URL("../lib/ui.js", document.baseURI).href);
  return ui.segmented({
    label: "Gap",
    options: [
      ["one", "One slit"],
      ["two", "Two slits"],
      ["grating", "Grating"],
      ["circle", "Circle"],
    ],
    value: "two",
  });
}

{
  // Shared WebGPU boilerplate: setup, buffers and the animation loop (lib/gpu.js), and
  // the fast Fourier transform (lib/fft.js)
  const gpu = await import(new URL("../lib/gpu.js", document.baseURI).href);
  const fft = await import(new URL("../lib/fft.js", document.baseURI).href);
  const aperture = await import(new URL("../lib/aperture.js", document.baseURI).href);
  const { slits } = aperture;
  const ui = await import(new URL("../lib/ui.js", document.baseURI).href);

  // Controls are read live from their elements, so moving them does not restart
  // this cell.
  const shapeInput = viewof shape;
  const brush = 6; // radius of the brush, in cells

  // The gap on the left and the pattern on the right, each 512 by 512 cells
  const N = 512;
  const { canvas, device, context, format } = await gpu.setup({
    width: 1024,
    height: 512,
    resolution: "auto",
  });
  canvas.style.touchAction = "none";

  // The gap: 1 where it is open and 0 where it is wall, kept here and copied across
  const mask = new Float32Array(N * N);
  let area = 0; // how much of it is open
  let painted = false; // true once the shape has been changed by hand
  let dirty = true; // the gap has changed since the pattern was made
  let needRipple = true; // the ripple shown on the gap has to be worked out again

  function shapeMask(kind) {
    aperture.fill(mask, N, kind);
    painted = false;
    dirty = true;
  }

  // ------------------------------------------------------------------- the transform
  const transform = fft.create(device, N);
  const maskBuffer = gpu.storageBuffer(device, mask);

  // Copy the gap into the transform's input: real part from the gap, imaginary zero
  const loadModule = device.createShaderModule({
    code: `
      @group(0) @binding(0) var<storage, read> gap: array<f32>;
      @group(0) @binding(1) var<storage, read_write> data: array<vec2f>;

      @compute @workgroup_size(16, 16)
      fn load(@builtin(global_invocation_id) id: vec3u) {
        if (id.x >= ${N}u || id.y >= ${N}u) { return; }
        let i = id.y * ${N}u + id.x;
        data[i] = vec2f(gap[i], 0.0);
      }
    `,
  });
  const loadPipeline = device.createComputePipeline({
    layout: "auto",
    compute: { module: loadModule, entryPoint: "load" },
  });
  const loadGroup = device.createBindGroup({
    layout: loadPipeline.getBindGroupLayout(0),
    entries: [
      { binding: 0, resource: { buffer: maskBuffer } },
      { binding: 1, resource: { buffer: transform.data } },
    ],
  });

  // ---------------------------------------------------------------------- drawing
  const paramBuffer = device.createBuffer({
    size: 32,
    usage: GPUBufferUsage.UNIFORM | GPUBufferUsage.COPY_DST,
  });
  const cornerBuffer = gpu.unitSquare(device);
  const showModule = device.createShaderModule({
    code: `
      struct Params {
        area: f32,    // how much of the gap is open: the brightest the pattern can be is area squared
        kx: f32,      // the ripple picked out: how many fit across the gap sideways
        ky: f32,      // and up and down
        phase: f32,   // where its crests sit: on the first open cell of the gap
        pickX: f32,   // where it is in the pattern, in cells
        pickY: f32,
        pulse: f32,   // swings between 1 and -1 with time, so that crests and troughs swap
        pad1: f32,
      }

      @group(0) @binding(0) var<uniform> params: Params;
      @group(0) @binding(1) var<storage, read> gap: array<f32>;
      @group(0) @binding(2) var<storage, read> spectrum: array<vec2f>;

      struct Out {
        @builtin(position) position: vec4f,
        @location(0) uv: vec2f,
      }

      @vertex
      fn showVertex(@location(0) corner: vec2f) -> Out {
        return Out(
          vec4f(corner, 0.0, 1.0),
          vec2f(corner.x * 0.5 + 0.5, 0.5 - corner.y * 0.5),
        );
      }

      @fragment
      fn showFragment(in: Out) -> @location(0) vec4f {
        let n = ${N}u;
        let y = min(u32(in.uv.y * f32(n)), n - 1u);
        if (abs(in.uv.x - 0.5) < 0.0015) { return vec4f(0.5, 0.5, 0.5, 1.0); }

        let crest = vec3f(1.0, 0.282, 0.204);   // #FF4834
        let trough = vec3f(0.22, 0.60, 0.85);

        // The gap, on the left, with the ripple picked out in the pattern laid over it:
        // red where its crests fall and blue where its troughs do. Where the gap is open
        // the ripple shows in full and elsewhere only faintly
        if (in.uv.x < 0.5) {
          let x = min(u32(in.uv.x * 2.0 * f32(n)), n - 1u);
          let open = gap[y * n + x];
          let ripple = params.pulse * cos(6.283185 * (params.kx * f32(x) + params.ky * f32(y)) / f32(n) + params.phase);
          let tint = select(trough, crest, ripple > 0.0);
          let shut = mix(vec3f(0.10), tint, 0.12 * abs(ripple));
          let through = mix(vec3f(0.92), tint, abs(ripple));
          return vec4f(mix(shut, through, open), 1.0);
        }

        // The pattern, on the right: the zero frequency is the first element of the
        // transform, so shift by half to put it in the middle
        let x = min(u32((in.uv.x - 0.5) * 2.0 * f32(n)), n - 1u);
        let f = spectrum[((y + n / 2u) % n) * n + (x + n / 2u) % n];
        let intensity = (f.x * f.x + f.y * f.y) / max(params.area * params.area, 1.0);
        // Four decades of brightness, on a log scale
        let t = clamp((log(intensity + 1e-9) / 2.302585 + 4.0) / 4.0, 0.0, 1.0);
        let dark = vec3f(0.10, 0.10, 0.10);
        let white = vec3f(1.0, 0.95, 0.9);
        var colour = select(mix(dark, crest, t / 0.6), mix(crest, white, (t - 0.6) / 0.4), t > 0.6);
        // A ring round the point whose ripple is shown on the left
        let away = length(vec2f(f32(x), f32(y)) + 0.5 - vec2f(params.pickX, params.pickY));
        if (abs(away - 8.0) < 1.5) { colour = vec3f(0.95); }
        return vec4f(colour, 1.0);
      }
    `,
  });
  const showPipeline = device.createRenderPipeline({
    layout: "auto",
    vertex: { module: showModule, entryPoint: "showVertex", buffers: [gpu.squareLayout] },
    fragment: { module: showModule, entryPoint: "showFragment", targets: [{ format }] },
    primitive: { topology: "triangle-list" },
  });
  const showGroup = device.createBindGroup({
    layout: showPipeline.getBindGroupLayout(0),
    entries: [
      { binding: 0, resource: { buffer: paramBuffer } },
      { binding: 1, resource: { buffer: maskBuffer } },
      { binding: 2, resource: { buffer: transform.data } },
    ],
  });

  // ------------------------------------------------------------- picking a ripple
  // A point in the pattern stands for one ripple: kx crests across the gap sideways and
  // ky up and down. How much of it the gap contains is the sum, over the open cells, of
  // exp(-2 pi i (kx x + ky y) / N). The transform on the GPU does this for every ripple at
  // once. Here it is done for just the one picked, to show it on the gap
  let pick = { kx: 0, ky: 0 }; // the middle of the pattern, where the ripple is flat
  let brightness = 0; // how much of the ripple the gap holds, relative to the most it could
  let anchor = 0; // puts a crest of the ripple on the first open cell of the gap
  let open = new Uint32Array(0); // the open cells, in order, found again whenever the gap changes
  // cos and sin of 2 pi j / N for each whole j, so that the sum below needs no trigonometry
  const cosine = Float64Array.from({ length: N }, (_, j) => Math.cos((2 * Math.PI * j) / N));
  const sine = Float64Array.from({ length: N }, (_, j) => Math.sin((2 * Math.PI * j) / N));
  function updateRipple() {
    let re = 0;
    let im = 0;
    for (const cell of open) {
      const x = cell % N;
      const y = (cell - x) / N;
      // The ripple's phase at this cell, in steps of 2 pi / N
      const step = (((pick.kx * x + pick.ky * y) % N) + N) % N;
      re += cosine[step];
      im -= sine[step];
    }
    if (open.length) {
      const x = open[0] % N;
      const y = (open[0] - x) / N;
      anchor = (-2 * Math.PI * (pick.kx * x + pick.ky * y)) / N;
    } else {
      anchor = 0;
    }
    // The size of the sum is how bright the point is. As the ripple pulses on the gap
    // the sum swings up and down by that much
    brightness = (re * re + im * im) / Math.max(area * area, 1);
    rippleReadout.textContent =
      `ripple: ${Math.abs(pick.kx)} across, ${Math.abs(pick.ky)} up · brightness ${brightness.toFixed(3)}`;
  }

  // ------------------------------------------------------------------ painting
  // Brush strokes go into the gap, on the left half of the canvas
  let painting = false;
  let last = null;
  const place = (event) => {
    const box = canvas.getBoundingClientRect();
    const fx = (event.clientX - box.left) / box.width;
    const fy = (event.clientY - box.top) / box.height;
    return [fx * 2 * N, fy * N];
  };
  function stamp(x, y) {
    const r = brush;
    for (let j = Math.max(0, Math.floor(y - r)); j <= Math.min(N - 1, Math.ceil(y + r)); j++) {
      for (let i = Math.max(0, Math.floor(x - r)); i <= Math.min(N - 1, Math.ceil(x + r)); i++) {
        if ((i - x) ** 2 + (j - y) ** 2 <= r * r) mask[j * N + i] = 1;
      }
    }
    painted = true;
    dirty = true;
  }
  function stroke(from, to) {
    // Stamp along the way so that a fast drag leaves no gaps
    const steps = Math.max(1, Math.ceil(Math.hypot(to[0] - from[0], to[1] - from[1]) / 2));
    for (let s = 1; s <= steps; s++) {
      stamp(from[0] + ((to[0] - from[0]) * s) / steps, from[1] + ((to[1] - from[1]) * s) / steps);
    }
  }
  // Over the pattern: pick the ripple under the mouse
  function choose(at) {
    const kx = Math.min(N - 1, Math.floor(at[0] - N)) - N / 2;
    const ky = Math.min(N - 1, Math.max(0, Math.floor(at[1]))) - N / 2;
    if (kx !== pick.kx || ky !== pick.ky) {
      pick = { kx, ky };
      needRipple = true;
    }
  }
  function back() {
    if (pick.kx !== 0 || pick.ky !== 0) {
      pick = { kx: 0, ky: 0 };
      needRipple = true;
    }
  }
  canvas.addEventListener("pointerdown", (event) => {
    const at = place(event);
    if (at[0] >= N) return choose(at); // the pattern side is not paintable
    canvas.setPointerCapture(event.pointerId);
    painting = true;
    last = at;
    stamp(...at);
  });
  canvas.addEventListener("pointermove", (event) => {
    const at = place(event);
    if (painting) {
      stroke(last, at);
      last = at;
    } else if (at[0] >= N) {
      choose(at);
    } else if (event.pointerType !== "touch") {
      back();
    }
  });
  // With the mouse off the pattern the pick goes back to the middle
  canvas.addEventListener("pointerleave", (event) => {
    if (event.pointerType !== "touch") back();
  });
  const stop = () => (painting = false);
  canvas.addEventListener("pointerup", stop);
  canvas.addEventListener("pointercancel", stop);

  // -------------------------------------------------------------------- the charts
  // The middle row of the pattern (the zero frequency vertically) is the first row of
  // the transform
  const readback = gpu.readback(device, [{ buffer: transform.data, bytes: N * 8 }]);
  let row = []; // {k, I}: brightness along the middle row, relative to the centre
  let needRead = true;

  const [rowHost] = charts.children;
  const draw = ui.plotter(Plot);

  // What last time's formula gives for these slits: a hump from the width of each slit
  // times fringes from the spacing between them, with k the number of ripples that fit
  // across the 512 cells
  function formula() {
    const kind = shapeInput.value;
    if (painted || !(kind === "one" || kind === "two")) return [];
    const a = slits[kind].width;
    const d = kind === "two" ? slits[kind].centres[1] - slits[kind].centres[0] : 0;
    const rows = [];
    for (let k = -100; k <= 100; k += 0.5) {
      const alpha = (Math.PI * a * k) / N;
      const hump = alpha === 0 ? 1 : (Math.sin(alpha) / alpha) ** 2;
      rows.push({ k, I: hump * (d === 0 ? 1 : Math.cos((Math.PI * d * k) / N) ** 2) });
    }
    return rows;
  }

  function drawRow() {
    draw(rowHost, {
      height: 220,
      x: { label: "ripples across the gap (sideways in the pattern)", grid: true, domain: [-100, 100] },
      y: { label: "brightness (relative to the centre)", grid: true, domain: [0, 1.02] },
      marks: [
        Plot.line(row, { x: "k", y: "I", stroke: "#FF4834" }),
        Plot.line(formula(), {
          x: "k",
          y: "I",
          stroke: "currentColor",
          strokeOpacity: 0.7,
          strokeDasharray: "4 3",
        }),
      ],
    });
  }

  // ------------------------------------------------------------------- the picture
  const readout = ui.readout(); // how fast the GPU is
  const rippleReadout = ui.readout(); // which ripple is picked

  // How long a 512 by 512 transform takes on this GPU: a batch of them, timed, with the
  // first one thrown away. It overwrites the transform's data, so the pattern is made again
  // afterwards (the gap starts out marked as changed)
  async function benchmark() {
    const warm = device.createCommandEncoder();
    transform.encode(warm);
    device.queue.submit([warm.finish()]);
    await device.queue.onSubmittedWorkDone();
    const runs = 100;
    const encoder = device.createCommandEncoder();
    for (let r = 0; r < runs; r++) transform.encode(encoder);
    const start = performance.now();
    device.queue.submit([encoder.finish()]);
    await device.queue.onSubmittedWorkDone();
    return (performance.now() - start) / runs;
  }
  const ms = await benchmark();
  readout.textContent = `${N} by ${N} transform: ${ms.toFixed(2)} ms`;

  shapeMask(shapeInput.value);
  // A press on the shape selector builds the shape. That includes pressing the one that
  // is already chosen, which puts it back as it was
  const again = () => shapeMask(shapeInput.value);
  shapeInput.addEventListener("click", again);
  invalidation.then(() => shapeInput.removeEventListener("click", again));
  const clearButton = ui.button({
    label: "Clear",
    onClick: () => {
      mask.fill(0);
      painted = true;
      dirty = true;
    },
  });
  drawRow();

  function frame(ms) {
    const encoder = device.createCommandEncoder();
    let sample = false;
    if (dirty) {
      dirty = false;
      const cells = [];
      for (let i = 0; i < mask.length; i++) if (mask[i]) cells.push(i);
      open = Uint32Array.from(cells);
      area = open.length;
      device.queue.writeBuffer(maskBuffer, 0, mask);
      const load = encoder.beginComputePass();
      load.setPipeline(loadPipeline);
      load.setBindGroup(0, loadGroup);
      load.dispatchWorkgroups(N / 16, N / 16);
      load.end();
      transform.encode(encoder);
      needRead = true;
      needRipple = true;
    }
    if (needRipple) {
      needRipple = false;
      updateRipple();
    }
    // The ripple stays where it is and pulses: crests turn to troughs and back
    const centre = N / 2 + 0.5;
    device.queue.writeBuffer(
      paramBuffer,
      0,
      new Float32Array([area, pick.kx, pick.ky, anchor, pick.kx + centre, pick.ky + centre, Math.cos(ms / 350), 0]),
    );
    if (needRead && !readback.pending) {
      readback.copy(encoder);
      needRead = false;
      sample = true;
    }

    const pass = gpu.beginRenderPass(encoder, context);
    pass.setPipeline(showPipeline);
    pass.setBindGroup(0, showGroup);
    pass.setVertexBuffer(0, cornerBuffer);
    pass.draw(6);
    pass.end();
    device.queue.submit([encoder.finish()]);

    if (sample) {
      readback.read((buffer) => {
        const f = new Float32Array(buffer);
        const centre = f[0] ** 2 + f[1] ** 2 || 1;
        row = [];
        for (let k = -100; k <= 100; k++) {
          const i = (k + N) % N;
          row.push({ k, I: (f[2 * i] ** 2 + f[2 * i + 1] ** 2) / centre });
        }
        drawRow();
      });
    }
  }

  gpu.run(invalidation, device, frame);

  // The canvas with the controls laid over its corner
  yield ui.stage(canvas, [shapeInput, clearButton, rippleReadout, readout]);
}

Analysis

charts = {
  // An empty host; the simulation cell above draws into it as data arrives
  const ui = await import(new URL("../lib/ui.js", document.baseURI).href);
  return ui.chartHosts(["row"]);
}

The brightness along the middle row of the pattern, for the One slit and Two slits gaps. The solid line is what the transform gives and the dashed line is last time’s formula for the pattern behind a slit: a hump from the width of the slit and, for two, fringes from the space between them. They are the same thing. The formula was the Fourier transform of a slit worked out by hand and here the computer has done it for every point at once. It stops matching as soon as you paint on the gap.

A sum of ripples

In 1822 Joseph Fourier published his book on how heat flows. In it was the claim that any pattern along a line can be built by adding sine waves of different sizes. A picture is a pattern along two lines at once, so it can be built from ripples: waves that run across it in some direction with some number of crests.

To find how much of one ripple a gap contains you multiply the gap by the ripple, cell by cell and add up. That sum is one point of the pattern:

F(k_x, k_y) = \sum_{x, y} h(x, y) \, e^{-2\pi i (k_x x + k_y y) / N}.

Here h is 1 where the gap is open and 0 where it is wall. The exponential is the ripple, with k_x crests across the gap sideways and k_y up and down and N is 512. The brightness of the point is |F|^2. It is the sum you watch when you point at the pattern. Where the ripple puts a crest on every open part of the gap the terms all add and the sum is large. Where it puts crests on some and troughs on others they cancel.

How the computer does it

The page works out that sum by hand for the one ripple you point at, a loop over the open cells of the gap. For the whole pattern it would need a sum for every ripple. The GPU does that for all 262,144 of them at once by a shortcut. We’ll find out how later. The readout under the canvas is how long it takes on this machine, measured when the page loaded.

The shader that draws the pattern squares each result to get a brightness and puts it on a log scale, four factors of ten from the brightest to the dimmest, since the middle is very much brighter than the edges. The transform puts the zero frequency in the corner of the array and the drawing shifts it by half to the middle.

Sources

• Fourier, J. (1822), Théorie analytique de la chaleur (Paris: Firmin Didot).
• Adams, C. S. and Hughes, I. G. (2018), Optics f2f: From Fourier to Fresnel (Oxford University Press), for the pattern behind a gap as a Fourier transform.

Moving Figures