`synthio.Biquad` is inaccurate below ~300 Hz and above ~16 kHz (two independent causes)

[bdbarnett]

A biquad asked for a corner near either end of the audio band is not the
filter that was asked for, and it fails silently — a LOW_PASS at 50 Hz
returns silence, a HIGH_PASS at 30 Hz boosts, and a low shelf does
nothing.

There are two independent causes and neither can be fixed alone, which is
why they are in one report.

Cause 1: Q15 coefficients in an int32 accumulator

shared-module/synthio/Biquad.h:

#define BIQUAD_SHIFT (15)

shared-module/synthio/Biquad.c:

static int32_t biquad_scale_arg_float(mp_float_t arg) {
    return (int32_t)MICROPY_FLOAT_C_FUN(round)(MICROPY_FLOAT_C_FUN(ldexp)(arg, BIQUAD_SHIFT));
}
...
static inline int32_t biquad_filter_sample(int32_t input, int32_t a1, int32_t a2, int32_t b0, int32_t b1, int32_t b2, int32_t x0, int32_t x1, int32_t y0, int32_t y1) {
    return synthio_sat16((b0 * input + b1 * x0 + b2 * x1 - a1 * y0 - a2 * y1 + (1 << (BIQUAD_SHIFT - 1))), BIQUAD_SHIFT);
}

Both halves are the problem, for different reasons.

The coefficients are what quantize. A low-pass has b0 = (1 - cos W0)/2,
which goes to zero with W0. At 48 kHz:

corner	b0 (normalised)	in Q15	stored
1000 Hz	3.916e-3	128.32	128
400 Hz	6.608e-4	21.65	22
200 Hz	1.682e-4	5.51	6
100 Hz	4.244e-5	1.39	1
50 Hz	1.066e-5	0.35	0

By 100 Hz the numerator is carrying one and a half significant bits, and by
50 Hz it has rounded away entirely.

The accumulator is what overflows. a1 approaches −2 as W0 goes to
zero: at 100 Hz it is −1.9815, which is −64929 in Q15. Multiplied by a
full-scale sample that is 2.13e9, against an INT32_MAX of 2.147e9 — so one
of the five terms already fills the accumulator, with four still to add.

And the feedback state is rounded to whole samples. y0/y1 hold the
saturated int16 output, and a biquad low down has both poles close to the
unit circle: 1/A(z) at DC is about 4000 for a 100 Hz low-pass and 43000 for
a 30 Hz high-pass. That is the gain the loop applies to whatever error is fed
back into it, so half an LSB of rounding is not half an LSB by the time it
comes round again. This is what takes a 50 Hz low-pass from "wrong" to
"silent", and what makes a 30 Hz high-pass gain rather than merely mis-sit.

Cause 2: fast_sincos() fits one polynomial to two functions

#define FOUR_OVER_PI (4 / M_PI)
static void fast_sincos(mp_float_t theta, sincos_result_t *result) {
    mp_float_t x = (theta * FOUR_OVER_PI) - 1;
    ...
    mp_float_t evens = c4x4 + c2x2 + c0, odds = c5x5 + c3x3 + c1x;
    result->c = evens + odds;
    result->s = evens - odds;
}

One 5th-order fit to both sine and cosine over [0, pi/2], worth about
2.3e-5 of absolute error at worst. That is a fine approximation in
general, and wrong at both ends of the audio band for two different reasons.

At the bottom, every low-pass, notch and shelf coefficient is built from
1 - cos W0, which is a difference of two nearly-equal numbers:

corner (48 kHz)	1 - cos W0	cosine error	error as a fraction
1000 Hz	8.56e-3	2.3e-5	0.3 %
300 Hz	7.71e-4	1.4e-5	1.8 %
100 Hz	8.57e-5	5.5e-6	6.4 %
50 Hz	2.14e-5	2.8e-6	13.2 %
20 Hz	3.43e-6	1.2e-6	33.7 %

At the top, pi/2 in W0 is SR/4 — only 12 kHz at 48 kHz — and above
that the fit is extrapolating:

tone (48 kHz)	sine	error	1 + cos W0	error
16 kHz	+0.86914 vs +0.86603	+0.4 %	0.50221 vs 0.50000	+0.4 %
20 kHz	+0.53305 vs +0.50000	+6.6 %	0.15369 vs 0.13397	+14.7 %
22 kHz	+0.33865 vs +0.25882	+30.8 %	0.07766 vs 0.03407	+128 %

This half appears not to have been reported before. It is also the half that
is cheap to fix.

Worth noting for a board rather than a desktop: mp_float_t is single
precision by default (py/circuitpy_mpconfig.h,
MICROPY_FLOAT_IMPL_FLOAT), so the cancellation in 1 - cos W0 is worse
there than in the double-precision figures above — float epsilon alone costs
about 3.5 % at 100 Hz before the polynomial's error is counted.

What it measures

One biquad, 48 kHz, Q 0.707 unless stated, on a build of main. A corner
should read −3.01 dB and a passband 0.00:

asked for	reads	should be
LOW_PASS 50 Hz, at cutoff	silence	−3.01
LOW_PASS 100 Hz, at cutoff	−3.94	−3.01
LOW_PASS 200 Hz, at 50 Hz	−0.03	−0.02
HIGH_PASS 30 Hz, at cutoff	+9.03	−3.01
HIGH_PASS 100 Hz, at cutoff	−4.36	−3.01
LOW_SHELF 80 Hz +1.5 dB, at 20 Hz	−0.28	+1.49
LOW_PASS 16 kHz, at cutoff	−3.06	−3.01
LOW_PASS 20 kHz, at cutoff	−3.71	−3.01
LOW_PASS 22 kHz, at cutoff	−7.30	−3.01
HIGH_PASS 22 kHz, at 1 kHz (stopband)	−52.41	−82.56

The usable range is roughly 300 Hz to 16 kHz, and about 400 Hz to 12 kHz for
anything under half a decibel of error.

Why both causes have to go together

Turning each cause off in the closed form, at 48 kHz:

	as shipped	Q15 only	polynomial only	ideal
LOW_PASS 50 Hz @ cutoff	−6.09	−6.09	−1.96	−3.01
LOW_PASS 100 Hz @ cutoff	−3.95	−3.95	−2.48	−3.01
LOW_SHELF 80 Hz +1.5 @ 20 Hz	+0.08	+0.08	+1.50	+1.49
HIGH_PASS 30 Hz @ cutoff	−4.71	−4.71	−3.12	−3.01
LOW_PASS 20 kHz @ cutoff	−3.71	−3.01	−3.71	−3.01
LOW_PASS 22 kHz @ cutoff	−7.29	−3.01	−7.29	−3.01

Read the columns:

• Below 300 Hz, quantization swamps everything — "as shipped" and "Q15
  only" are identical, so the polynomial's error is currently invisible
  there. Widen the coefficients and it appears: about 1 dB at 50 Hz.
• Above 16 kHz, the polynomial is the whole story — "Q15 only" is exact
  and "polynomial only" reproduces the shipped error to the last digit.

So a fix to either half alone leaves half the band still wrong. (The closed
form is also optimistic at the very bottom: it says −6.09 dB at 50 Hz where
the build gives silence, and −4.71 at 30 Hz where the build gives +9.03. The
difference is the whole-sample feedback state described above, which a
transfer function cannot express.)

A fix, and what it costs

This is what we did in a downstream port of these files. Offering it as a
starting point rather than a patch to take, because the cost below is real
and it is your call whether to pay it.

1. Per-filter coefficient format instead of a fixed Q15. Take the largest
   of the five normalised coefficients and give them all as many fractional
   bits as that one has room for in an int32_t, capped at 30. A plain
   low-pass tops out near 2 (that is a1) and gets 29 bits; a 20 dB shelf
   reaches ~200 and gets 23. One fixed format for every filter means giving
   them all the shelf's. The chosen shift has to travel with the coefficients
   — they are meaningless apart — so it is cached in synthio_biquad_t
   alongside a1..b2.
2. int64_t accumulator, and a feedback state carrying ~12 bits below the
   sample grid. The second half is not optional, for the 1/A(z) reason
   above: rounding the feedback to whole samples hands a gain of thousands
   half an LSB to amplify.
3. Replace fast_sincos(): reflect into [0, pi/2]
   (sin(pi - t) = sin t, cos(pi - t) = -cos t) and evaluate the two Taylor
   series separately, seven terms each. Worst error across 20 Hz – 24 kHz
   falls from 2.3e-5 to 6.3e-9, and it stays sane past Nyquist, so a frequency
   asked for above Nyquist is merely wrong rather than absurd. About ten extra
   multiplies, paid once per coefficient update at block rate, not per sample.

With all three, every mode lands within 0.03 dB of the double-precision
closed form from 50 Hz to 22 kHz.

The cost, measured by cross-compiling the filter kernel at -Os and
counting instructions in the per-sample function:

core	before	after	note
Cortex-M4 / M7	37	76	no library calls — SMULL/SMLAL
Cortex-M0+	55	154	5 __aeabi_lmul + 1 __aeabi_lasr

Roughly 2x on anything with a long multiply, closer to 4x without one. One
stereo biquad at 22050 Hz goes from about 5 % of a 48 MHz M0+ to about 17 %.
That is the trade Q15 was making, and on an RP2040 it is not a small one. If
it matters, the products are 32x32 into 64 and could be hand-written for M0
rather than going through __aeabi_lmul.

Step 3 on its own is nearly free on every core, and it is the entire fix above
16 kHz — so it is worth considering separately even if steps 1 and 2 are not
wanted.

Repro

# What a biquad does at the ends of the audio band.
import array, math
import audiocore, audiofilters, synthio

RATE = 48000
FRAMES = 8192


def tone(hz, level=6000, frames=FRAMES):
    values = array.array("h")
    for frame in range(frames):
        values.append(int(level * math.sin(2 * math.pi * hz * frame / RATE)))
    return values


def through(mode, corner, hz, Q=0.7071067811865475, A=None, level=6000):
    values = tone(hz, level)
    source = audiocore.RawSample(values, sample_rate=RATE, channel_count=1)
    biquad = (synthio.Biquad(mode, corner, Q=Q) if A is None
              else synthio.Biquad(mode, corner, Q=Q, A=A))
    node = audiofilters.Filter(filter=biquad, sample_rate=RATE,
                               channel_count=1, bits_per_sample=16,
                               samples_signed=True, buffer_size=FRAMES * 2)
    node.play(source, loop=True)
    total = 0
    count = 0
    for block in range(4):
        data = bytes(audiocore.get_buffer(node)[1])
        if block < 2:                      # let the filter settle
            continue
        for index in range(0, len(data) - 1, 2):
            sample = data[index] | (data[index + 1] << 8)
            if sample >= 32768:
                sample -= 65536
            total += sample * sample
            count += 1
    out = math.sqrt(total / count)
    reference = level / math.sqrt(2)
    return 20 * math.log10(out / reference) if out > 0 else float("-inf")


FM = synthio.FilterMode
print("A corner should read -3.01 dB, and a passband 0.00.")
print()
for label, mode, corner, hz, A in (
        ("LOW_PASS   50 Hz, at cutoff", FM.LOW_PASS, 50.0, 50.0, None),
        ("LOW_PASS  100 Hz, at cutoff", FM.LOW_PASS, 100.0, 100.0, None),
        ("LOW_PASS  200 Hz, at 50 Hz ", FM.LOW_PASS, 200.0, 50.0, None),
        ("HIGH_PASS  30 Hz, at cutoff", FM.HIGH_PASS, 30.0, 30.0, None),
        ("HIGH_PASS 100 Hz, at cutoff", FM.HIGH_PASS, 100.0, 100.0, None),
        ("LOW_SHELF  80 Hz +1.5 dB, at 20 Hz",
         FM.LOW_SHELF, 80.0, 20.0, 10 ** (1.5 / 40)),
        ("LOW_PASS   16 kHz, at cutoff", FM.LOW_PASS, 16000.0, 16000.0, None),
        ("LOW_PASS   20 kHz, at cutoff", FM.LOW_PASS, 20000.0, 20000.0, None),
        ("LOW_PASS   22 kHz, at cutoff", FM.LOW_PASS, 22000.0, 22000.0, None),
        ("HIGH_PASS  22 kHz, at 1 kHz (stopband)",
         FM.HIGH_PASS, 22000.0, 1000.0, None)):
    print("%-40s %+8.2f dB" % (label, through(mode, corner, hz, A=A)))

A corner should read -3.01 dB, and a passband 0.00.

LOW_PASS   50 Hz, at cutoff                  -inf dB
LOW_PASS  100 Hz, at cutoff                 -3.94 dB
LOW_PASS  200 Hz, at 50 Hz                  -0.03 dB
HIGH_PASS  30 Hz, at cutoff                 +9.03 dB
HIGH_PASS 100 Hz, at cutoff                 -4.36 dB
LOW_SHELF  80 Hz +1.5 dB, at 20 Hz          -0.28 dB
LOW_PASS   16 kHz, at cutoff                -3.06 dB
LOW_PASS   20 kHz, at cutoff                -3.71 dB
LOW_PASS   22 kHz, at cutoff                -7.30 dB
HIGH_PASS  22 kHz, at 1 kHz (stopband)     -52.41 dB

audiocore.get_buffer is not in upstream; on a board, play each filter into
an AudioOut and capture instead. The effects run from 1 dB to total silence,
so they survive any reasonable measurement.

---

Verified still present in 10.3.0-rc.0 (2026-08-28) by inspection of the tagged source.

[dhalbert]

The LLM output here is hard to understand. Please rewrite it for a human to read, and say what motivated this investigation. Were you dissatisfied with the filter results. What is your use case?

What do you want to be done to fix this? Do you have a fix? As the text mentions, there are tradeoffs of accuracy vs code size and computation, and one might say the current implementation is good enough.

[bdbarnett]

Hey there — taking another run at this, shorter.

We hit this building a synthesizer/effects library on top of synthio: musical filter sweeps routinely close a low-pass below 200 Hz, and there the filter isn't the one asked for. At 48 kHz, a LOW_PASS at 50 Hz reads silence at its own cutoff instead of −3 dB, and a HIGH_PASS at 30 Hz boosts +9 dB. Above ~16 kHz there's an independent problem: a 22 kHz HIGH_PASS leaks ~30 dB into its stopband. Usable range today is roughly 300 Hz–16 kHz.

Two independent causes: Q15 coefficient quantization (by 100 Hz the smallest numerator coefficient is down to about one bit), and the shared sine/cosine polynomial extrapolating above SR/4.

We have a working fix in our downstream port and can send it as a PR — per-filter coefficient scaling with a wider accumulator, plus a corrected trig fit. The cost is real on small chips: per stereo biquad, roughly 2× instructions on M4/M7 and ~3× on M0+ (five 64-bit multiply helper calls, ~5% → ~17% of a 48 MHz M0+).

Your call on appetite: (a) the full fix, (b) the full fix gated per-port or config so M0+ keeps today's cost, or (c) just the cheap half — the trig fix, which fully repairs the >16 kHz end and leaves low frequencies as they are. Happy to open whichever as a PR; if the answer is "current is good enough," no hard feelings — we'll carry it downstream.

Full measurements available if useful.

[dhalbert]

@relic-se: Interested in your thoughts.

[relic-se]

Although the implementation of synthio.Biquad has come a long way from where it started, it is definitely not perfect. These issues that you have identified are not a surprise, and I think it is a good idea to continue to enhance biquad functionality which is critical to many components of the synthio and audio systems. A lot of the math here is over my head (and I do agree that the LLM language is a bit wordy), but thank you for providing detailed analysis for us to work off of.

1. Per-filter coefficient format instead of a fixed Q15. Take the largest of the five normalised coefficients and give them all as many fractional bits as that one has room for in an int32_t, capped at 30. A plain low-pass tops out near 2 (that is a1) and gets 29 bits; a 20 dB shelf reaches ~200 and gets 23. One fixed format for every filter means giving them all the shelf's. The chosen shift has to travel with the coefficients — they are meaningless apart — so it is cached in synthio_biquad_t alongside a1..b2.

This could definitely be a solution, but I can't say that it is exactly elegant. I'd be very interested to see what kind of performance we'd get by just using single-precision floating point on platforms with an FPU (ie: RP2350 / M33) which I would assume would also solve these problems. Otherwise, I can see us moving forward with this implementation.

2. int64_t accumulator, and a feedback state carrying ~12 bits below the sample grid. The second half is not optional, for the 1/A(z) reason above: rounding the feedback to whole samples hands a gain of thousands half an LSB to amplify.

I don't think I fully understand what you're proposing here. Are you saying to change the properties of biquad_filter_state to int64_t? I'd be interested to see what your code looks like here in a PR and if there'd be some other alternative here which wouldn't require 64-bit calculations which aren't ideal on 32-bit ALUs.

3. Replace fast_sincos(): reflect into [0, pi/2] (sin(pi - t) = sin t, cos(pi - t) = -cos t) and evaluate the two Taylor series separately, seven terms each. Worst error across 20 Hz – 24 kHz falls from 2.3e-5 to 6.3e-9, and it stays sane past Nyquist, so a frequency asked for above Nyquist is merely wrong rather than absurd. About ten extra multiplies, paid once per coefficient update at block rate, not per sample.

Yes, this does seem like the solution here, but I wish there'd be a way to optimize it further. Once again, seeing the code in practice here would help that evaluation.

Step 3 on its own is nearly free on every core, and it is the entire fix above 16 kHz — so it is worth considering separately even if steps 1 and 2 are not wanted.

This should definitely be submitted as a separate PR for further evaluation.

audiocore.get_buffer is not in upstream; on a board, play each filter into an AudioOut and capture instead. The effects run from 1 dB to total silence, so they survive any reasonable measurement.

Rather than doing a specific build with the audiocore debug stuff enabled, audiofilewriter.AudioFileWriter is available to write to a WAV file and then read it back for analysis (https://docs.circuitpython.org/en/latest/shared-bindings/audiofilewriter/index.html). But the example provided is definitely good though.

---

TL;DR

Go ahead and prepare these updates (steps 1-3) and push them out as individual PRs if possible for individual code analysis and testing.

---

Out of curiosity, what platform are you targeting? That might help target some of these efforts. You've included M0 in a lot of your analysis, but that platform is extremely limited in flash availability. It may be necessary to enhance synthio.Biquad only on some platforms depending on the effect on flash, ram and performance.

[dhalbert]

Which LLM proposed these fixes? We might point a frontier model at this and see what it says?

[bdbarnett]

It was Claude Fable 5 ultracode. This is over my head too, but, we've completely ported all CircuitPython audio to micropython and CPython. We are testing on ESP32-S3, ESP32-P4, micropython on Linux, micropyton.exe on Windows, CPython on everything - Windows, Linux (WSL), WASM, and android. I don't want to burden you. We will keep the code either way. Just wanted to offer what we can.

[tannewt]

we've completely ported all CircuitPython audio to micropython and CPython.

Why have you ported it to micropython and CPython? It'd be great as part of Blinka.