Mitigate floating-point numerical errors for very-low corner frequency low-pass filter with DSP

Question

I'm designing a low-pass filter for a digital signal processing application that ideally just passes a very small bandwidth above DC. I'm using an IIR biquad filter for this, where the coefficients are derived using the instructions here. A smaller bandwidth leads to a longer filtering time (larger time constant) but yields a more accurate result whereas a larger bandwidth can be filtered faster but is less accurate. Both of these are valid use cases.

Here's the code I've got

#!/usr/bin/env python
import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import freqz

# calculates filter coefficients using link above
# fc is corner frequency, fs is sample freq
def iir_lp_coeffs(fc, fs):
    w0 = 2 * np.pi * fc / fs
    q = 1 / np.sqrt(2)
    alpha = np.sin(w0) / (2 * q)
    b0 = (1 - np.cos(w0)) / 2
    b1 = 1 - np.cos(w0)
    b2 = b0
    a0 = 1 + alpha
    a1 = -2 * np.cos(w0)
    a2 = 1 - alpha
    b0 /= a0
    b1 /= a0
    b2 /= a0
    a1 /= a0
    a2 /= a0
    a0 /= a0
    return (
        np.array([b0, b1, b2], dtype=np.float64),
        np.array([a0, a1, a2], dtype=np.float64),
    )


fc = 2  # low pass corner frequency (Hz)
fsample = 500e3
b, a = iir_lp_coeffs(fc, fsample)

w, h = freqz(b, a, worN=int(1e6), fs=fsample)
fig, ax = plt.subplots()
ax.plot(w, 20 * np.log10(abs(h)))
ax.set_ylim(-40, 10)
ax.set_xscale("log")
plt.show()
print(w[0:10])
print(abs(h[0:10]))

The current settings use 64-bit floating point with a cutoff frequency of $2\,\text{Hz}$. This all works fine, and I can even decrease the corner frequency substantially as long as I increase the granularity of freqz (with worN=).

For instance here's a plot of the gain response with the above code (note that I've cut the x axis off at the higher frequencies):

However, my actual application requires 32-bit floating point. When I do this (set dtype of iir_lp_coeffs to np.float32), I get non-unity gain in the passband. For instance, here's a gain response with fc=10 using 32-bit:

If I set the corner frequency higher, the gain response looks correct again (e.g., fc=100 looks fine).

Am I running up against the limit of what's possible with 32-bit FP? Or, is there another strategy that would allow me to get away with the lower precision of 32-bit? Have I correctly diagnosed this issue as a floating-point issue?

Answer 1

I think your issue might be coefficient quantization and filter topology. A direct form biquad has poor quantization effects about 0 and π radians. It’s easier to analyze such effects in fixed point, but even though floating point has a much larger range, it still has shortcomings. In particular, if you add a very small number to a very large number, the small number disappears because it can’t be aligned for the operation in the available number of mantissa bits. This can cause the order of operations to affect the result. For instance, where S is a small number and L is large, L - L + S = S, but L + S - L = 0.

Udo Zolzer covers the differences between several filter structures in his book Digital Audio Signal Processing. I borrow direct form quantization effects on pole locations from the book:

See how precision is lost near 0 and π. Other filter topologies might have higher precision near 0, while being much worse near π, which may be a good tradeoff for uses such as yours. The Gold and Rader form has a very even distribution, it looks like a perfect grid.

Another simple and popular filter that has good quantization characteristics at low frequencies is the "Chamberlin" state variable filter. There are improved versions of this filter, as it has problems at higher frequencies (from about one-sixth of the sample rate and up), but the plain Chamberlin is very good at low frequencies, where you need it.

See my article on the Chamberlin state variable filter here:

The digital state variable filter

Zolzer presents modified Chamberlin structures here:

The Modified Chamberlin and Zölzer Filter Structures

In particular, see the graph of quantization effect near zero for the Chamberlin structure—very dense near zero, at the expense of poorer performance at high frequencies, compared to the direct form graph:

Answer 2

For the cookbook LPF, i would make use of this trig identity:

$$ \cos(\omega_0) = 1 - 2\sin^2 \left( \frac{\omega_0}{2} \right) $$

what you're doing is subtracting a number very close to one from one and all of the information is in the difference.