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

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.

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.
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.
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.

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.

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.

@relic-se

@relic-se: Interested in your thoughts.

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.

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.

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.

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.

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

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.

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.

Why? Because your audio code was the right fit for my project. Your audio implementation is fantastic for constrained boards but also holds up extremely well on cpython.

You are welcome to the port - https://github.com/pydevices/audioif
On testpypi as pydevices-audioif.

I think the additions may be more valuable to you:

audiodynamics (compression, limiting, expansion, gating, transient shaping)
audioroute (fan one stream out to parallel branches)
audiomath (multiply one stream by another — ring and amplitude modulation)
audioecho (a delay with a filter, a soft-clip and a cross-feed inside its feedback loop)
audioconvolve (apply a measured or synthesized impulse response, by partitioned FFT; not for constrained boards!)

If you'd rather not discuss on this issue, let me know where. I'd be happy to talk with you about it. github:PyDevices/pydevices has a discussions page if you'd rather use it.

If you'd rather not discuss on this issue, let me know where.

Those are really interesting. How about another "enhancement" issue here?

Hi I was using Claude to help me with some audio FX testing. When playing some stuff with the whole FX chain, I kept hearing a low rumbling, farty sound. We narrowed it down to the filters causing the sound, and Claude found a bug in how synthio.Biquad rounds numbers.

I believe that this issue I am facing might be related to this and decided not to open a new GitHub issue. I asked Claude to simplify the issue 4 times and I hope it is understandable: -

The very short version
A rounding error makes negative numbers come out 1 unit (between -32,768 and 32,768) too high on every sample. For a Biquad with a low frequency setting, the filter’s feedback magnifies that tiny error hundreds of times, turning something far below the noise floor into a steady offset and low-frequency wobble that is loud enough to hear.

What happened
On an RP2350 (CircuitPython 10.3.1) I added a master EQ of three biquads (200 Hz, 1 kHz, 5 kHz), all set to gain 1.0 (A=1.0), so it should do nothing. It gave my kick drum low thumping pops; bypassing it fixed them. With the EQ on, the output also never went back to silence after a hit: the DAC data line kept toggling for 3+ seconds (if EQ was bypassed, it was completely quiet).

Commenter's note: Claude used a shelf filter as an example as it said that it's a very clear example -- a filter that should not be doing anything is introducing noise due to the aforementioned rounding issue.

Why (from the source and a simulation, not a patched build)
biquad_filter_sample adds half a unit and then calls synthio_sat16, which already rounds toward zero, so negative results end up about one unit too high. The filter's feedback magnifies that small error, and the lower the filter's frequency, the more it is magnified: about 4 times bigger each time the frequency is halved. As a share of the loudest possible signal, the leftover error is about 0.003% at 4 kHz, 0.05% at 1 kHz, 0.2% at 500 Hz, 1.1% at 200 Hz (what I heard as pops) and 2.2% at 150 Hz.

Possible fix
Remove the extra half step (or round evenly). In my simulation that makes the flat filter exactly neutral, and it costs nothing at run time.

Claude also spat the following check code that it used to validate the bug, which shows that after feeding noise through the formula, the signal does not reset to 0: -

def sat16(n, r=15):                    # same rounding as synthio_sat16
    if n < 0:
        n += (1 << r) - 1
    n >>= r
    return max(-32768, min(32767, n))

import random
rnd = random.Random(1)
a1, a2 = -63921, 31193                 # 200 Hz low shelf, A=1.0, 36 kHz (numerator = denominator)
b0, b1, b2 = 32768, a1, a2
for extra in (1 << 14, 0):             # as in 10.3.1, then with the extra half step removed
    x0 = x1 = y0 = y1 = 0
    tail = []
    for n in range(36000):
        inp = int(rnd.gauss(0, 100)) if n < 18000 else 0     # half a second of quiet noise, then silence
        out = sat16(b0*inp + b1*x0 + b2*x1 - a1*y0 - a2*y1 + extra)
        x1, x0, y1, y0 = x0, inp, y0, out
        if n >= 34000:
            tail.append(out - inp)
    print(extra, sum(tail) / len(tail))    # prints about 372 for the first, 0 for the second

Navigation Menu

Repository navigation

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

Cause 1: Q15 coefficients in an int32 accumulator

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

What it measures

Why both causes have to go together

A fix, and what it costs

Repro

TL;DR

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

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

Description

Cause 1: Q15 coefficients in an int32 accumulator

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

What it measures

Why both causes have to go together

A fix, and what it costs

Repro

Activity

dhalbert commented on Aug 31, 2026

bdbarnett commented on Aug 31, 2026

dhalbert commented on Aug 31, 2026

relic-se commented on Sep 1, 2026

TL;DR

dhalbert commented on Sep 1, 2026

bdbarnett commented on Sep 2, 2026

tannewt commented on Sep 2, 2026

bdbarnett commented on Sep 4, 2026

tannewt commented on Sep 4, 2026

GlasierXplor commented on Oct 10, 2026

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions

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

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