audiobiquad.Biquad: direct-form I with float state loses half a decibel of peak gain below 100 Hz at high Q

[bdbarnett]

audiobiquad.Biquad's recursion is direct-form I with float state
(src/shared/audioif_filter_f32.c:220-238):

y0 = b0 * x0 + b1 * x1 + b2 * x2 - a1 * y1 - a2 * y2;

At a low centre frequency the two feedback terms nearly cancel, and single
precision cannot hold the difference. A BAND_PASS section at 20 Hz / Q 32
on a 48 kHz graph has a1 = -1.999911308, a2 = +0.999918163, w0 = 0.002618
rad: 1 + a1 + a2 = 6.9e-6, so the wanted per-sample increment is about
2e-5 of the numbers being differenced, and the rounding error accumulates over
the Q*F_s/(pi*f0) samples the resonator rings for.

What it costs, measured. RBJ's constant-0 dB-peak band-pass has unity gain
at its own centre at every Q. Rendered through audiobiquad.Biquad, a 20 Hz
tone at that setting comes back -0.50 dB, and the sign of the error moves
with the probe level (+0.09 dB at -3 dBFS, -0.48 at -12, +0.44 at -20), which
is the tell that it is round-off and not a coefficient error:

$ PYTHONPATH=lib .venv/bin/python tools/phase2_probes/bandpass_lowcorner.py arithmetic
  RBJ closed form, double, exact                    +0.0000 dB
  the kernel's float32 coefficients, evaluated exactly -0.0031 dB
  those coefficients, float64 recursion              -0.0031 dB
  those coefficients, float32 recursion (the kernel) -0.4801 dB
  the class itself, the same 10 s render              -0.4801 dB

The four rows use the same five coefficients, read off the node with
Biquad.coefficients. Coefficient quantization costs 0.003 dB; the state's
precision costs the other 0.477. That rules out both the coefficient
arithmetic and audioif_sincos_reflect.

How far it reaches. Worst |gain at f0| over probe levels -3/-12/-20/-40
dBFS, one section, 48 kHz, audiocomponents BandPass over its own macro
grid (... bandpass_lowcorner.py map):

  f0 \ Q         0.5      0.707          2          4          8         16         32
    20.0    +0.0100    -0.0044    -0.0103    -0.0619*   +0.0434    -0.1679*   -0.4986*
    25.0    -0.0021    +0.0035    +0.0039    +0.0514*   -0.1284*   +0.1431*   -0.1869*
    31.5    -0.0018    +0.0096    -0.0140    +0.0148    -0.0470    -0.1746*   -0.3323*
    40.0    -0.0008    +0.0013    +0.0105    -0.0149    -0.0323    -0.0383    -0.1104*
    63.0    +0.0007    -0.0008    -0.0012    +0.0081    +0.0028    -0.0187    -0.0613*
   100.0    +0.0002    +0.0002    -0.0007    -0.0005    -0.0030    -0.0145    -0.0095
   200.0    +0.0014    +0.0009    +0.0016    +0.0013    +0.0021    +0.0015    -0.0067
  1000.0    +0.0013    +0.0027    +0.0027    +0.0027    +0.0027    +0.0027    +0.0027

* = outside 0.05 dB. Eleven of 56 cells, every one of them f0 <= 63 Hz with
Q >= 4. Above 100 Hz nothing is out at any width the surface reaches.

Why it matters beyond one class. This is the arithmetic every
audiobiquad-tier filter shares. audiocomponents' BandPass is parked on
it (PyDevices/audiocomponents#39), and the gate audit records LowPass T1 at
-0.444 dB and Notch at f0 20 Hz Q 32 reading -0.724 dB - the same corner, on
the same kernel. A sub-bass resonance is a setting players use; the shelf and
peaking modes carry the same cancellation at low w0.

What would fix it, in rough order of cost:

1. Transposed direct form II in audioif_biquad_f32_process_s16(). Same
   five coefficients, same two state words per channel, but the accumulation
   is not the difference of two nearly-equal large numbers - the standard
   answer to exactly this defect, and it costs nothing at run time.
2. double state for the four memory words (still float coefficients).
   Costs 16 bytes per channel and the ESP32-S3 has no double FPU, so this is
   the expensive one.
3. Leave it and document the corner in docs/upstream-diff.md, so a class
   that wants a 20 Hz Q 32 resonance knows the peak is not 0 dB there.

Reproduction is audiocomponents tools/phase2_probes/bandpass_lowcorner.py
on effects/p2b-bandpass; the class is unchanged by it (the file only reads
Biquad.coefficients and renders).