Month: June 2026

In my previous behemoth post I mentioned the struggle of implementing a comb + notch filter that looked passable for analog composite video. The blurring effect of adding the notch with FIR filtering was far too great, and there’s no way TVs at the time did it that way.

Nagging issues like these haunt me until I solve it, and I started experimenting with IIR filtering since it’s one of those topics I never really explored on GPUs, due to it being very GPU unfriendly. A notch IIR biquad filter is extremely effective on paper, and there has to be a way to do it well in parallel.

A biquad IIR

This is the building block of a ton of audio filters out there, since they can easily model classic analog filters and are easy to tweak parametrically. In the Z plane they have two zeroes and two poles:

// Z is a complex number representing where frequency response hits 0.
// P is a complex number representing a pole
// where the frequency response hits infinity (or pure resonance).
// Frequency response at a given frequency is measured by setting
// z = exp(2 * j * w)

// I probably screwed this formula up somehow,
// don't read too much into it.
H(z) = (z - Z) * (z - conj(Z)) / ((z - P) * (z - conj(P));

Having conjugated poles and zeroes ensures the filter coefficients are real, so the actual filter looks more like this in code:

out[x] =
  in[x] * b0 + in[x - 1] * b1 + in[x - 2] * b2 + // FIR portion
  out[x - 1] * a1 + out[x - 2] * a2 // IIR portion

The Z plane is often a little abstract, but in this case, we will design the notch directly in the Z domain. For example, we place zeroes where we want to remove frequency response, e.g. at the chroma subcarrier, then place poles close to it at the same frequency. The intuition is that close to the carrier, the response will fall to 0, and further away the pole and zero mostly negate each other, leaving a unity frequency response. The response can be made as sharp as desired, but sharper notches lead to more ringing.

E.g. for a pole set at 0.9 times the zero. Quite sharp. Getting this kind of response out of a normal convolution requires a comical number of taps. The response at DC here isn’t exactly 0 dB, but the FIR portion can be rescaled to ensure unity gain.

struct
{
  float b0, b1, b2, a1, a2; // a0 is implicitly 1.0.
} push = {};

// Compute a biquad IIR notch filter.
double w = 2.0 * muglm::pi<double>() / 27.0;
constexpr double fir_q = 0.99;
constexpr double iir_q = 0.85;

// Don't do a complete notch since that leads to really bad ringing.
// We just need a fairly narrow band-stop though ...

if (options.system == System::NTSC)
  w *= 315.0 / 88.0;
else
  w *= 4.43361875;

// Preconvolves two zeroes at exp(i * w) and exp(-i * w).
double b0 = 1.0;
double b1 = -2.0 * fir_q * cos(w);
double b2 = fir_q * fir_q;

// Preconvolves two poles at exp(i * w * Q) and exp(-i * w * Q).
// Flip sign here since we're designing the filter as B(z) / A(z) and
// when evaluating the filter we flip the signs.
double a1 = 2.0 * iir_q * cos(w);
double a2 = -iir_q * iir_q;

// Normalize the FIR gain to obtain a 0 dB gain at DC.
double fir_gain = b0 + b1 + b2;
double target_fir_gain = 1.0 - a1 - a2;
double fir_scale = target_fir_gain / fir_gain;

b0 *= fir_scale;
b1 *= fir_scale;
b2 *= fir_scale;

push.b0 = float(b0);
push.b1 = float(b1);
push.b2 = float(b2);
push.a1 = float(a1);
push.a2 = float(a2);

cmd.push_constants(&push, 0, sizeof(push));

The trick – a spicy parallel scan

While IIRs have infinite impulse response, we’re only really interested in the output of the active signal area, which is ~1440 pixels. A 1440-tap convolution is fine if we want to use FFT convolution, but that’s quite overkill for what’s literally a 5-tap filter. Surely there are ways to do IIRs effectively without falling back to serial algorithms.

Many parallel GPU algorithms boil down to scan operations in some way, or have scans as a primary ingredient of their madness. I found an older paper from 1990 (Prefix Sums and Their Applications, Guy E. Blelloch, section 1.4) which goes through this algorithm in a way that was easy enough to digest. It’s not the state of the art I think, but it gets the job more than adequately for my needs. For the IIR filter, the FIR portion of it is handled trivially with a filter kernel. We only need to concern ourselves with the feedback portion of the filter.

Usually we think of scans as just simple binary operators.

for i in range(len(data)):
  data[i] += data[i - 1]

This is common enough that we have special subgroup operations for them, like subgroupInclusiveAdd. This is not enough to implement IIRs since it cannot capture the “decay” of an impulse response. The two intuitions from the paper is that the binary operator can be as complex as we want it to be as long as it follows some properties that are expected of binary operators, and the data doesn’t have to be a scalar either. For an IIR with two histories, we need a 2×2 matrix and 2 element vector, then define a binary operator which is effectively:

struct Biquad
{
  mat2 A;
  vec2 B;
};

Biquad init(float v)
{
  // Note that GLSL is column major here.
  return Biquad(mat2(vec2(a1, a2), vec2(1, 0)), vec2(v, 0));
}

Biquad scan(Biquad prev, Biquad current)
{
  return Biquad(prev.A * current.A, prev.B * current.A + current.B);
}

The intuition is that the A matrix keeps track of impulse response decay over time, while the B vector serves as the history buffer. Once the scan is done, the result is extracted from B.x and written out to memory. The performance is excellent and processes my frame at like sub 10 microseconds. Can’t complain about that.

For the full shader example, see https://github.com/Arntzen-Software/parallel-gs/blob/main/analog/shaders/iir.comp, but the general gist of it is:

256 thread workgroup – work on N pixels per thread

The horizontal resolution is expected to be reasonable, so N is about 6 in my use case. Not having to split the scan over multiple kernels is good.

This means we can load 8 inputs and compute the FIR portion of it for 6 inputs all in registers. The initial scan can also be done in registers:

float xs[PER_THREAD + 2];
Biquad quad[PER_THREAD];

[[unroll]]
for (int i = -2; i < PER_THREAD; i++)
  xs[i + 2] = LoadInputInRow(PER_THREAD * int(index) + i);

// Do the FIR portion and produce 6 values.
[[unroll]]
for (int i = 0; i < PER_THREAD; i++)
  quad[i] = init(xs[i + 2] * b0 + xs[i + 1] * b1 + xs[i + 0] * b2);

[[unroll]]
for (int i = 1; i < PER_THREAD; i++)
  quad[i] = scan(quad[i - 1], quad[i]);

Custom scan with magic biquad operator

Then we do the custom scan over the subgroup using ShuffleUp.

for (uint step = 1; step < gl_SubgroupSize; step *= 2)
{
  Biquad prev;
  prev = subgroupShuffleUp(quad[PER_THREAD - 1], step);
  if (gl_SubgroupInvocationID >= step)
  {
    [[unroll]]
    for (int i = 0; i < PER_THREAD; i++)
      quad[i] = scan(prev, quad[i]);
  }
}

… and then complete the scan via groupshared memory with same strategy. Groupshared memory is kept at a minimum since we only need to store the last value per subgroup. That’s about it, and with that I have a decoded composite signal I’m happy with.

For very large IIR filters, this approach gets out of hand quickly since we need a big N x N matrix to hold the recurring state per element, but chaining together biquads is done trivially by just scanning multiple times. Since we can do the horizontal filter in one compute workgroup, the intermediate steps can easily live in shared memory if desired to avoid excessive memory bandwidth.