Skip to content

FIx: Use full-precision multiply for biquad feedback coefficients#950

Closed
rt-rtos wants to merge 1 commit into
shorepine:mainfrom
rt-rtos:fix/filter-fixedpoint-precision
Closed

FIx: Use full-precision multiply for biquad feedback coefficients#950
rt-rtos wants to merge 1 commit into
shorepine:mainfrom
rt-rtos:fix/filter-fixedpoint-precision

Conversation

@rt-rtos

@rt-rtos rt-rtos commented Jul 20, 2026

Copy link
Copy Markdown
Contributor

Use full-precision multiply for biquad feedback coefficients

Discovered while investigating the root cause of intermittent clipping caused by resonant
HPF in sub the 100hz range.

Symptom

An HPF with cutoff below ~100 Hz and Q > 1 produces massive sustained clipping,
worse at higher Q. Also reproducible on the LPF when a high-Q resonance overlaps
the note fundamental, though it's less audible there because the low cutoff
attenuates most of the band. The repro is erratic - some cutoff/Q combinations
ring up violently, others nearby lose their resonance instead.

Root cause

FILT_MUL_SS is SMULR6, which truncates both operands to 12 fractional bits
(>>11 on s8.23) before multiplying.

The split-feedback kernel dsps_biquad_f32_ansi_split_fb exists precisely to
protect pole precision - it computes feedback as
y0 = w0 + 2*y1 - y2 - e*y1 + f*y2 with the corrections e = 2 + a1 and
f = 1 - a2 deliberately small. But small is exactly what 12-bit truncation
destroys. At fc = 65 Hz, Q = 4, e ~ f ~ 0.0011, so after quantization they
retain only 2-3 significant bits.

So the mechanism meant to preserve the pole is the one that loses it. With the
corrections down to a couple of bits, the pole's frequency and damping land
quasi-randomly around their targets: overshoot and it rings up ~+12 dB into
sustained clipping - full-scale, since an HPF passes the whole spectrum -
overdamp and the resonance vanishes instead. Which side the rounding falls on
decides which you get, which is why the repro felt erratic.

Independently, the low-cutoff LPF numerator (1-c)/2 ~ 3e-5 quantizes to
zero, which is why LPF below ~100 Hz goes silent.

Ruled out: block-floating-point state overflow (instrumented - zero
near-overflow events in all failing static-cutoff runs).

Evidence

Verified on a host-side replica of the fixed-point path: s8.23 SAMPLE,
SMULR6, the split_fb kernel, and filter_process's normbits BFP logic, all
driven by a retriggered 110 Hz saw at 0.4 FS and run against a double-precision
reference on the same coefficient stream. 4-second runs, static cutoff.

Filter fc Q peak (SMULR6) peak (reference) clipped samples peak (SMUL64R)
HPF 65 4 4.128 0.664 166394 0.664
HPF 80 4 3.418 0.871 129961 0.871
HPF 80 8 1.593 0.912 51449 0.913
HPF 65 8 0.450 (resonance lost) 0.664 0 0.664
LPF 65 4 0.000 (silent) 0.445 0 0.444
LPF 200 8 1.085 0.918 974 0.918
HPF 500 4 0.850 0.824 0 0.824

With the exact multiply, every row matches the double-precision reference to 3
decimals and all clipping disappears. Above ~200 Hz the old path was already
fine, which is consistent with the diagnosis - e and f only get small at low
cutoff.

The one case that still clips

HPF fc=10Hz Q=8.0 sweep=6oct   SMULR6:   peak=7.337  ref=1.562  clipped=34005  near-wrap=464
HPF fc=10Hz Q=8.0 sweep=6oct   SMUL64R:  peak=1.562  ref=1.562  clipped=12075  near-wrap=0

The double-precision reference clips
identically at the same peak. That's legitimate resonance gain from a 6-octave
sweep on a Q=8 resonator.

Recreating it

The harness is a single self-contained C file replicating the fixed-point path
(no AMY build required):

Happy to attach the harness as a gist or paste it in-thread if that's useful for
independent verification - it's ~190 lines.

The fix

In the fixed-point build only, FILT_MUL_SS becomes SMUL64R, an exact
32x32->64 rounding multiply:

SMUL64R(a, b) = ((int64_t)a * b + 2^22) >> 23

The float build keeps SMULR6 untouched.

Cost is roughly 2 extra Xtensa instructions per multiply across the biquad
kernels. Note this includes the parametric EQ path, which was a deliberate
SMULR6 choice for CPU reasons - if EQ cost measurably regresses on target, the
EQ path can be specialized back to SMULR6, since its fixed mid-band
coefficients are not precision-sensitive. Happy to do that pre-emptively if
preferred.

Hardware-verified on ESP32-S3: the original repro is gone and
resonant low-frequency filters sound better subjectively.

Related, not addressed here

  • dsps_biquad_gen_hpf_f32 / _bpf_ lack the w0 floor (~70 Hz) and Q floor
    (0.51) that _lpf_ has. With the precise multiply they're stable down to the
    4.8 Hz LOWEST_RATIO clamp, so those floors are now a voicing question rather
    than a stability one.

FILT_MUL_SS was SMULR6, which truncates both operands to 12 fractional
bits before multiplying. For low-cutoff biquads the split-feedback pole
corrections e = 2 + a1 and f = 1 - a2 are ~2^-10, so they retained only
2-3 significant bits and the pole landed almost randomly around its
target:

- sub-150 Hz HPF/BPF with Q >= 2 rang up ~+12 dB into sustained
  full-scale clipping, or lost its resonance entirely depending on
  where the quantized pole fell;
- the sub-100 Hz LPF numerator (1-c)/2 ~ 3e-5 rounded to zero,
  producing silent output.

Replace it, in the fixed-point build only, with SMUL64R: an exact
32x32->64 rounding multiply. Costs a couple of extra instructions per
multiply in the biquad kernels.

reduced comment verbosity
dpwe added a commit that referenced this pull request Jul 20, 2026
filters.c uses SMUL64R (direct 64x64->64) on desktop and ESP32S3.

Incorporates and supersedes #950.
@dpwe

dpwe commented Jul 20, 2026

Copy link
Copy Markdown
Collaborator

The int64_t mul speedup is incredible on the EPS32-S3! I incorporated into the "main" LPF24 filter too, and made it my own submit #951. Sorry to appropriate your massive improvement!

@dpwe dpwe closed this Jul 20, 2026
@rt-rtos

rt-rtos commented Jul 20, 2026

Copy link
Copy Markdown
Contributor Author

No worries. I am glad you could claw back more performance than was lost for correctness.

Appending my A/B bench. Delayed because it required new scenes.

Scene Wire string (first osc) G type F cutoff R Q
filt_hpf65 v30w2f110G3F65R4l0.4Z 3 = HPF 65 Hz 4
filt_bpf90 v36w2f110G2F90R8l0.4Z 2 = BPF 90 Hz 8
filt_lpf60 v42w2f110G1F60R2l0.4Z 1 = LPF 60 Hz 2
filt_eq x6,-6,6Z + v48w2f110l0.4Z none - -

Which code each filter reaches

flowchart TD
    G1["G1 LPF"] --> SFB["dsps_biquad_f32_ansi_split_fb"]
    G2["G2 BPF"] --> SFB
    G3["G3 HPF"] --> SFB
    G4["G4 LPF24"] --> TW["split_fb_twice (BFP)"]
    EQ["parametric EQ"] --> PEQ["parametric_eq_process"]

    SFB --> FMS["FILT_MUL_SS"]
    TW --> T16["top16SMUL + scan_max"]
    PEQ --> TEQ["FILT_MUL_SS_EQ = top16SMUL"]

    FMS --> B950["#950: SMULR6 to SMUL64R<br/>slower, more accurate"]
    T16 --> B951["#951: whole BFP path replaced<br/>by nobfp kernel on SMUL64R"]
    TEQ --> NC["untouched by both"]

    B950 --> R1["+4.6 to +6.5%"]
    B951 --> R2["-23.8%"]
    NC --> R3["-0.12%, CRC same"]
Loading

Stale comments discovered in the process

  • The EQ has its own macro (FILT_MUL_SS_EQ, filters.c:541). The comment
    above FILT_MUL_SS warning that "parametric_eq was chewing major CPU" is
    stale - that was already solved by giving the EQ a separate macro. It reads
    as a live constraint on this change and is not one. Worth fixing (I took at at face value).

Also dead, so irrelevant to both PRs despite using FILT_MUL_SS: hpf_buf
(defined, zero callers repo-wide) and dsps_biquad_f32_ansi_commuted (its only
call site, filters.c:773, is commented out).


Measurements

ESP32-S3 240 MHz, fixed-point, 48 kHz, block 256, free pacing, 5 captures/side,
3 passes, block budget 5333 us. Harness: amy-bench-public. Every run's
baseline is 2291ff8, so all columns are directly comparable.

#951 as merged (2291ff8 -> d427432)

Scene Path A cyc B cyc delta noise headroom output
juno6 LPF24 1,150,097 1,043,114 -9.30% ±0.00% 10.2% -> 18.5% CHANGED
saw_lpf6 LPF24 450,022 342,971 -23.79% ±0.00% 64.9% -> 73.2% CHANGED
filt_lpf60 G1 320,411 341,352 +6.54% ±0.01% 75.0% -> 73.4% CHANGED
filt_bpf90 G2 322,343 340,888 +5.75% ±0.02% 74.8% -> 73.4% CHANGED
filt_hpf65 G3 326,180 341,210 +4.61% ±0.01% 74.5% -> 73.4% CHANGED
filt_eq EQ 348,365 347,952 -0.12% ±0.01% 72.8% -> 72.8% same

juno6 headroom 10.2% -> 18.5% is the most consequential single number
here: it nearly doubles the margin on the tightest scene in the set.

#950 alone (2291ff8 -> 042b5be)

Scene delta output
filt_lpf60 +6.51% CHANGED
filt_bpf90 +5.54% CHANGED
filt_hpf65 +4.60% CHANGED
filt_eq +0.00% same
saw_lpf6 -0.00% same
juno6, dx76, sine8, fx_sine8, idle +0.00 to +0.24% same

The 12 dB regressions are identical across #950 and #951 (within 0.2 pp), which
confirms they come from the multiply swap alone and not from anything else in
#951. The LPF24 gain is entirely #951's

Conclusion:

  • Huge LPF24 perf win
  • Minor LPF,BPF,HPF perf loss. Improved correctness.

@rt-rtos

rt-rtos commented Jul 20, 2026

Copy link
Copy Markdown
Contributor Author

I don´t have both RPI Picos to test on device, but codegen seems like a win there too.
Epecially given the reduced complexity and correctness:

LPF24 path M0+ (RP2040) insns M0+ libcalls / sample M33 (RP2350) insns M33 libcalls / sample
BFP (top16SMUL, block floating-point) 473 __clzsi2 322 none
SMUL64R (PR #951) 209 __aeabi_lmul 112 none

arm-none-eabi-gcc 14.3, -O2, -mthumb. insns = whole-function static count; the per-sample loop dominates over a 64-sample block. Libcalls/sample are loop-body call sites (BFP also has 4 more __clzsi2 in once-per-block setup: 9 total). M33 has hardware SMULL + CLZ, so both paths are libcall-free there.

@dpwe

dpwe commented Jul 20, 2026

Copy link
Copy Markdown
Collaborator

I tried it out on my Pimoroni 2 (RP2350) board. I expected it not to help, and I sort-of got that confirmed (shutdown with fewer notes under SMUL64) but I have no faith that the Arduino UI is actually recompiling when I want it to, or that it's using the optimization flags I requested.

I should adapt the "make speedtest" target to run for my RP2350 setup, then I could measure for sure.

@dpwe

dpwe commented Jul 20, 2026

Copy link
Copy Markdown
Collaborator

From what I read, the RP2040 multiply is 32x32->32, but on the RP2350 it's 32x32->64, which would make all the difference. I must go deeper!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants