torvalds-GuitarPedal/Documentation/biquad-precision

Biquads in single precision
===========================

Everything in this pedal is float32, because the FPU on an RP2354A is
single precision only and a double costs hundreds of cycles in software
helpers.  For most of the DSP that is not even close to a constraint.
For a biquad at the bottom of the audio range it is the whole story, and
this is an account of where exactly it hurts, taken from measurements
rather than from reasoning about ulps.

The short version: **three separate errors, of wildly different sizes,
in a part of the code small enough that they all look like one.**  Two
are fixed; the third is a floor.

	angle recovered from a cosine    53 %      at 20 Hz   - fixed
	float32 coefficient storage      0.4 %     at 20 Hz   - floor
	1-cos cancellation in the LPF    0.87 %    at 20 Hz   - fixed

They are not the same kind of thing.  The first was a defect and was
audible.  The second cannot be removed without changing what a
'struct biquad_coeff' is.  The third is a classic and turns out to be
worth almost nothing here, which is itself worth knowing.

A fourth thing sits underneath the first, and is section 4: the sine
table's own interpolation, which is where the 53 % actually came from
and which has since been fixed at source.  That does not make the half
angle in section 1 redundant, but it does change what it is for - see
the end of section 3.


1. The angle, which was the real bug
------------------------------------

Every cookbook constructor places its poles and zeros with -2*cos(w0).
So what it needs out of the sine table is an *angle*, and it asks for it
by way of a cosine.

fastsincos() is 256 entries a quarter, linearly interpolated.  Near zero
phase the sine is straight, so the chord and the arc agree - the error at
30 Hz measures 1.4e-08.  The cosine at the same place is at its maximum,
where the chord sits furthest under the arc, and the error is 3.8e-06.

That would be a rounding error if cos were being used as a number.  It is
being used as an angle, and

	d(cos)/dw = -sin(w0)

vanishes exactly where the error is worst.  So the angle that comes back
is wrong by err/sin(w0), which grows without limit as the frequency
falls:

	asked      built      error
	   20    30.68 Hz     +53 %
	   30    37.49 Hz     +25 %
	   60    63.58 Hz      +6 %
	  100   101.27 Hz    +1.3 %
	 4000  4000.06 Hz   0.001 %

**This is not a Q-dependent error, and that is the thing most easily got
wrong about it.**  The band lands in the wrong place by the same
percentage whatever its Q; a sweep of nine constructors at six Q values
found the same 53 % at 20 Hz in every column.  What Q changes is how much
a given misplacement costs in dB, because a narrow filter has a steeper
skirt to be wrong on.

Measured at real note frequencies, for filters the pedal actually builds,
the error before the fix was:

	                                     B0     E1     A1     E2
	                                     31     41     55     82 Hz
	boost Basscut HPF Q.707 @30Hz      0.48   0.70   0.92   0.56
	boost Basscut HPF Q.707 @10Hz      3.34   1.94   1.08   0.48
	tone Bass LOSHELF Q.707 @20Hz      4.47   0.74   1.34   1.15
	tone Mid PEAKING Q4 @50Hz          0.17   0.82   1.66   0.07

Note that the worst row is a Q of 0.707 and the Q of 4 is not
particularly bad.  Frequency is the variable, not Q.  A dedicated hum
notch at Q 20 is where it becomes total - the notch misses the hum
entirely and attenuates 60 Hz by 0.7 dB instead of 37 - but nothing in
the pedal asks for that yet, which is why this survived so long.

The fix is to take the *half* angle, whose sine is the well-behaved one,
and come back with the double-angle identities:

	cos(w) = 1 - 2 sin(w/2)^2
	sin(w) = 2 sin(w/2) cos(w/2)

Both fall out of the one table lookup that was happening anyway.
cos(w) is now built out of a squared sine, which is the entire point: a
squared sine carries the angle, and a cosine near 1.0 does not.


2. The floor, which is the coefficients themselves
--------------------------------------------------

With the angle fixed, what is left is float32.

A biquad's pole angle is carried in a1, which for a low-frequency
section is a number near -2.  One ulp there is 1.19e-07, and the angle
is recovered from it through the same vanishing derivative as before, so
the placement cannot be better than about ulp/(2 sin w).

Swept every semitone with a peaking section at Q 1, comparing against
the same filter built in double:

	band          f64 coeffs      f32 coeffs, exact w0    as built
	              max     rms      max        rms         max     rms
	  20- 80 Hz   0 %     0 %      0.43 %     0.15 %      0.86 %  0.23 %
	  80-320 Hz   0 %     0 %      0.025 %    0.008 %     0.020 % 0.008 %
	 320-1280     0 %     0 %      0.002 %    0.001 %     0.002 % 0.001 %

So below 80 Hz there is a floor around 0.4 % that belongs to the storage
format, and the path used to get there costs about a factor of two on
top of it, from rounding cos to float32 near 1.0 on the way past.

**Restructuring the arithmetic does not recover that factor of two.**
Handing the constructors (1-cos) directly, so that a value near 1.0 is
never formed, measures identically to the current code - because a1 has
to end up near -2 whatever route it took, and that final store is where
the information goes.  Getting past this needs a different coefficient
set or a different topology, not a cleverer way to compute the same five
floats.  It is 0.4 % at 20 Hz, a third of a semitone, and the frequency
pots step by a whole semitone, so nothing is waiting on it.

**Nor does fixing the sine table**, which is the more interesting of the
two, because it was predicted to.  Once section 4 is in, the plain
cos(w0) route and the half-angle route measure the same thing:

	half angle, corrected table    worst 20-80 Hz   0.855 %
	plain cos,  corrected table    worst 20-80 Hz   0.855 %

both through bench/coeff, both real code.  A model said the plain route
would come out at 0.43 % and it does not, for a reason the model left
out: 'struct sincos' holds floats, so the plain route stores cos(w0) as
a float32 near 1.0 exactly as the half-angle route stores 1 - 2s^2.  The
last rounding before the coefficient is the same rounding either way,
and it is the one that sets the floor.


3. The cancellation, which is famous and barely matters
-------------------------------------------------------

	static inline void _biquad_lpf(..., float freq, float Q)
	{
		...
		float b1 = (1 - w0.cos) * a0_inv;

cos(20 Hz) is 0.99999657, so (1 - w0.cos) is the textbook case of
subtracting two nearly equal numbers: the result is 3.4e-06 with about
two significant figures left of the eight that went in.

This is the one place in the file where that happens.  The high-pass
wants (1 + cos), which is near 2 at the bottom of the range and only
degenerates approaching Nyquist, where the loss is a factor of seven and
worth nothing.  The shelves and the peaking sections never form either
quantity.  **One line, one filter type, one end of the range.**

What makes it disappear is that cos is no longer an independent input.
With biquad_w0() inlined into the constructor, the compiler sees

	1 - (1 - 2*s*s)

and -ffast-math folds it to 2*s*s.  The 1.0 never appears, and the small
quantity is computed at full relative precision.  On the target that is
two instructions:

	vmul.f32  s15, s0, s0        ; s*s
	vadd.f32  s15, s15, s15      ; 2*s*s   = 1 - cos

and the cosine is materialised separately, only for the pole
coefficient, where a value near 1.0 is what is actually wanted:

	vsub.f32  s14, s12, s15      ; 1 - 2s^2  = cos
	vmul.f32  s14, s14, s9       ; * -2

**This only works if the constructor can see where cos came from**,
which is why biquad_w0() is called inside the constructors rather than
by the caller.  Measured across 1440 coefficient sets through
bench/coeff - which calls the constructors through a function pointer
and so cannot inline them - moving the computation inside changed 200 of
them, every one a low-pass numerator:

	20 Hz  0.87 %      440 Hz  0.0017 %
	30 Hz  0.28 %     1000 Hz  0.0002 %
	60 Hz  0.086 %

In the firmware itself this had already been happening, because there is
one translation unit and the constructors inlined into their callers.
The point of moving it is that the precision no longer depends on the
optimiser making that choice.

And the size of it: 0.87 % of a low-pass numerator is 0.075 dB of
passband gain.  It is worth having because it is free.  It was never the
audible problem, and anyone reading (1 - w0.cos) and reaching for the
smelling salts - which is the natural reaction, and was the first guess
here - should measure it before rewriting anything around it.

**This is what the half angle is for now.**  Section 4 removes the
reason it was introduced: with the table corrected, asking for cos(w0)
directly is no longer the bad question, and the two routes place a
filter identically.  What the half angle still buys is this fold, which
the plain route cannot have - cos straight out of the table is not
algebraically 1 - 2s^2, so there is nothing for the compiler to cancel:

	numerator b0     plain cos      half angle
	  20 Hz        1.695593e-06   1.708975e-06    +0.783 %
	  30 Hz        3.863572e-06   3.884833e-06    +0.547 %
	  60 Hz        1.532269e-05   1.535918e-05    +0.238 %

So it stays, having swapped jobs: it was protection against the table
and is now the thing that keeps one line of _biquad_lpf exact.


4. The chord, which was under all of it
---------------------------------------

Section 1 fixed the biquads by not asking the table the question it was
bad at.  This is the same problem attacked at the other end: making the
question safe to ask.

fastsincos() interpolates linearly, and **a chord under an arc is always
nearer zero than the function it approximates**.  So the error is not
noise - it has a sign, and the table reads uniformly small.  That signed
part is what matters, because it does not average out: it is what made
the test tone come out 3 ppm quiet, and it is what got divided by
sin(w0) and became 53 %.

Linear interpolation misses by (h^2/2)*t*(1-t)*f''.  For a sinusoid f''
is -f, so **the miss is proportional to the value the chord just
produced** - which makes the correction a pure relative scale, one
multiply-add each:

	float defect = QUARTER_SINE_DEFECT * phase * (1.0f - phase);
	x += defect * x;
	y += defect * y;

Being proportional is what makes it hold over the whole quarter rather
than only in the middle: at a zero crossing the correction vanishes
exactly as the error does.  Measured over 16384 points through the
function itself:

	worst error   4.74e-06 -> 1.32e-07     (about two float32 ulp)
	mean bias    -2.00e-06 -> -3.24e-09

Nine instructions, 41 to 50, of which eight are floating point.  The
constant is h^2/2 and is emitted by quarter_sine.py, which is what knows
h - written beside the user it would go stale the first time the table
size moved, with nothing to say so.

Cubic Hermite was measured and rejected.  Its slopes are the *other*
output's endpoints, so it needs no new table reads, and its mathematical
error is h^4/384 = 3.7e-12 - but it measures 3.4e-08, because at that
point it is measuring float32 and not interpolation.  It costs 26
instructions instead of nine to buy 1.8x of pure rounding noise.

**Nothing audible changed.**  The generated sine's worst spur goes from
-120 dB to about -169 dB, against a best-case board noise floor of
-102 dBFS.  What changed is that the plain, obvious way to ask this
table for a small angle now works: with this in, the original
formulation of biquad_w0() passes the whole check-biquad sweep on its
own.  The half angle went from being the fix to being an optimisation.

analyze.h's Hann window and lfo.h's sine walk the same table and do not
get the correction, deliberately.  Neither is ever asked for an angle,
so the chord is good enough - which is the same "good enough until it
isn't" as the biquads, in the two places where it stays good enough.


What to take away
-----------------

- A cosine is a bad way to carry a small angle, and the giveaway is that
  the error is divided by sin(w0) on the way out.  If a quantity is
  going to be used as an angle, keep it as one.

- **The code that hit this was not wrong.**  It asked for a cosine and
  expected the table's advertised precision, which is what anyone would
  do, and nothing in the interface said that cosine-near-zero is the one
  place that expectation fails.  A primitive with a quiet exception in
  it will be walked into again, by somebody who has no reason to
  suspect it - so an approximation is better fixed than documented.

- **Measure the real thing, not a model of it.**  Three times in one day
  a numpy model of this arithmetic disagreed with bench/coeff, and every
  time the model was the optimistic one.  Twice it was a sampling
  artifact and once - the 0.43 % in section 2 - it was leaving out a
  float32 store that the real struct performs.  The model is for
  deciding what to measure.

- Which of three co-located errors dominates is not obvious by reading,
  and here the famous one was the smallest by a factor of sixty.  The
  53 % was invisible in the source and the 0.87 % was written in a form
  that draws the eye.

- An interface that takes a derived quantity invites somebody to derive
  it differently.  parametric_eq.h did exactly that and built its filters
  the old way for a day after the fix landed, because it called the
  constructors directly instead of through the macros that knew about
  _w0().  Taking the frequency removes the choice.

- 'make check-biquad' sweeps the constructors and did not catch that,
  because bench/coeff.c called _w0() too.  A test that shares the
  harness's idea of the right way to call something tests the primitive
  and not its adoption; see the note on measurement in ISSUES.md 140.