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

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?

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: • I'm getting excellent results with the Chamberlin filter in your article. Nice solution! Oct 19, 2020 at 15:55
• Great! I'm glad it worked for you. It was from his book, 40 years ago... Oct 19, 2020 at 17:09

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.

• This is clever. However, I'm not noticing an improvement and I think it's because even when I'm using 32-bit I'm first calculating the coefficients as 64-bit and then converting to 32-bit. Since I don't have issues with 64-bit anyway, getting a more precise 64-bit value before converting to 32-bit is unlikely to help much. Still, using this seems to only have upside. For anyone else coming here, the right way to apply it to the function I posted is b0=np.sin(w0/2)**2 and b1=2*b0. I think that's right, but correct me if I've misinterpreted. Oct 15, 2020 at 0:21
• you also have b2=b0, don't you? and what about a1? Oct 15, 2020 at 3:07
• and i think you want to check its behavior for fc set to something higher. Oct 15, 2020 at 3:12
• Yes, I left b2=b0. I've left a1=-2*np.cos(w0). That doesn't need to be changed does it? This will start to work anyway as I edge fc up. Are you saying you think the trig identity should get the gain response to improve for less increase of fc? Oct 15, 2020 at 3:15
• a1 will be very close to 2 when w0 is so small anyway.but so also will a0 and a2 be very close to 1. 2 Hz out of 500 Mhz is small. and there are only 25 bits of precision in the mantissa. that's your problem, it is too much for single precision without changing the definitions of the coefficients. 25 bits is only 16 or 32 million. much lower than 500 or 250 million. Oct 15, 2020 at 3:24