Yesterday I wrote a very brief introduction to z-transforms and finite-impulse-response filters. Today I’ll try to do infinite-impulse-response (IIR) filters. Obviously, one blog post can’t do justice to the topic, but I’ll try to get far enough that my son can implement some simple filters for his project. If he needs more detail, the tutorial at CCRMA by Julius Smith is probably worth reading.
The basic idea of an IIR filter is that the transfer function is not a polynomial but a rational function. For example, the standard “biquad” digital filter unit is the ratio of two quadratic formulas: 
There are a lot of different ways to implement the biquad section, trading off various properties, but one that has particularly nice numeric properties is the “direct form I” implementation.
Direct form I for implementing a biquad filter section. Taken from https://ccrma.stanford.edu/~jos/filters/Direct_Form_I.html
The transfer function can be checked by looking at the summation in the center and taking z-transforms: 
The left half of the implementation computes the numerator of the transfer function, which may have either two real zeros or a conjugate pair of complex ones (or a single real zero, if b2=0). It is often useful to look at these zeros in polar form, 
The right half of the implementation computes the denominator of the transfer function producing a pair of poles, and it is again useful to look at the poles in polar form. For stability of the filter we need to have the poles within the unit circle— even getting them too close to the unit circle can cause the output to get very large and cause numeric overflow.
Note that any rational function with real coefficients can be factored into biquad units with real coefficients, just by factoring the polynomials into a product of roots, and pairing the complex roots in their conjugate pairs. Because the direct implementation of a high degree polynomial is very sensitive to rounding errors in the parameters, it is common to factor the transfer function into biquad units.
About a year ago, I was playing with a biquad filter for an optical pulse monitor, where I wanted a roughly 0.4—2.5 Hz bandpass filter, with a sampling rate of 30Hz, using only small integer coefficients. I put the zeros at 1 and —1, and used a filter-design site on the web to get real-valued coefficients for the denominator. I then approximated the real-valued coefficients as integers divided by powers of two, and tried tweaking the coefficients up and down to get close to the response I wanted, using gnuplot to plot the gain:
Biquad design for a bandpass filter with peak around ω=π/15 using only small integer coefficients.
The code I used for implementing the biquad bandpass filter on a microcontroller was fairly simple:
y_0 = b0*(x_0-x_2) - a1*y_1 -a2*y_2;
Output(y_0);
x_2 = x_1;
x_1 = x_0;
y_2 = y_1;
y_1 = y_0/a0;
Note that I’m outputting the high-precision output from the summation node, which has an extra gain of a0. The memory elements (x_1, x_2, y_1, y_2) have the same precision as the input—only the summation is done at higher precision. The division in the last line is not really a division but a shift, since I restricted a0 to be a power of 2. This shifting could be done after the multiplications of y_1 and y_2 in the first line, if y_1 and y_2 are carried to higher precision, but then the intermediate results in the summation could get too large—a lot depends on how big the word size is compared to the input values and the a0 scaling.
My son is probably interested in two digital filters—one to remove the DC bias (so a high-pass with cutoff frequency around 10–20Hz) and one to pick out the beat from a non-linearly modified signal (absolute value or squared signal). The beat he is interested in picking out is probably around 60–180bpm or 1–3Hz—similar to what I did for picking out heartbeats (I used 0.4—2.5Hz as a rough design goal on the pulse monitor). Unfortunately, as you go to higher sampling frequencies, you need to have higher precision in the filter coefficients.
The biquad coefficient calculator at http://www.earlevel.com/main/2013/10/13/biquad-calculator-v2/ seems pretty good, though it is a bit tedious to rescale and round the coefficients, then check the result. So I wrote a Python script to use the scipy.signal.iirfilter function to design filters, then scaled and rounded the coefficients. The scaling factor had to get much larger as the sampling frequency went up to get a good bandpass filter near 1Hz, otherwise rounding errors resulted in a low-pass filter rather than a bandpass filter (perhaps one of the poles ended up at 1 rather than as a complex pair?). To make a 0.33–3Hz bandpass filter, I needed to scale by 230 at 40kHz, 221 at 10kHz, 216 at 3kHz, 215 at 1kHz, 211 at 300Hz, and 28 at 100Hz. The scaling factors needed for 40kHz sampling would exceed the 32-bit word size, so this approach did not look very promising.
It may be more effective to use two separate biquad sections: a low-pass filter with a fairly high cutoff to downsample the signal, then a bandpass filter at the lower sampling rate. This approach allows using much lower-precision computations. I wrote a little Python program to test this approach also, aiming for a 0.33–3Hz bandpass filter.
Here is the filter response of a low-pass filter, followed by down-sampling and a bandpass filter. Note that the scaling factor is only 214, but the filter response is quite nice.
So the overall design of the loudness sensor will probably be a series of filters:
- high-pass filter at 20Hz to remove DC bias, leaving sound
- nonlinear operation (squaring or absolute value)
- low-pass filter at 200Hz
- down-sample to 500Hz
- bandpass filter at 0.33–3Hz
One possible problem with this approach—look at the phase change when the beat is not 1Hz. At 0.33Hz, the phase change is about 37º, which is about 0.3s—a larger shift in the beat than we’d want to see. We may have to look at more complicated filter design that has smaller phase shifts. (The -49º phase shift at 3Hz is only about 45msec, and so not a perceptual problem.)
Here is a possible filter design for the initial high-pass filter:
Getting a 20Hz cutoff frequency required a somewhat larger scaling factor than for the other filters, but still quite reasonable for 14-bit values on a 32-bit machine.
Here is a version of the two-stage filter design program:
#!/usr/bin/env python
from __future__ import division,print_function
import numpy as np
import scipy
from scipy.signal import iirfilter,freqz,lfilter,tf2zpk
import matplotlib.pyplot as plt
pi=3.1415926535897932384627
fs = 40.e3 # sampling frequency
low_fs = 1000. # sampling frequency after downsampling
cutoffs = (0.33,3.0) # cutoff frequencies for bandpass filter in Hz
scale_factor = 2.**14
def scaled(b,a):
"""scaled b and a by scale_factor and round to integers.
Temporarily increase scale factor so that b[0] remains positive.
"""
extra_scale=1.
b_scaled = [int(np.round(scale_factor*c*extra_scale)) for c in b]
while b_scaled[0]==0:
b_scaled = [int(np.round(scale_factor*c*extra_scale)) for c in b]
extra_scale *= 2
a_scaled = [int(round(scale_factor*c*extra_scale)) for c in a]
print (" b=",b, "a=",a)
z,p,k = tf2zpk(b,a)
print ( "zeros=",z, "poles=",p, "gain=",k)
print ("scaled: b=",b_scaled, "a=",a_scaled)
z,p,k = tf2zpk(b_scaled,a_scaled)
print ( "zeros=",z, "poles=",p, "gain=",k)
return b_scaled,a_scaled
def plot(b1,a1, b2,a2, frequencies=None):
"""Plot the gain (in dB) and the phase change of the
concatentation of filters sepecified by b1,a1 and b2,a2.
The b1,a1 filter is designed to run at sampling rate fs
The b2,a2 filter is designed to run at sampling rate low_fs
Both are designed with the filter type specified in global filter.
"""
if frequencies is None:
worN =[pi*10.**(k/200.) for k in range(-1000,0)]
else:
worN = [2*pi*f/fs for f in frequencies]
freq,response1 = freqz(b1,a1, worN=worN)
freq2,response2 = freqz(b2,a2, worN=[f*fs/low_fs for f in worN])
freq *= fs/(2*pi)
response = response1*response2
fig=plt.figure()
plt.title ('{}: b1={} a1={} fs={}\n b2={} a2={} fs={}'.format(filter,b1,a1,fs,b2,a2,low_fs))
ax1 = fig.add_subplot(111)
plt.semilogx(freq,20*np.log10(np.abs(response)), 'b')
plt.ylabel('Amplitude [dB]', color='b')
plt.grid()
plt.legend()
plt.xlabel('Frequency [Hz]')
ax2 = ax1.twinx()
angles = [ 180/pi*ang for ang in np.unwrap(np.angle(response)) ]
plt.plot(freq,angles, 'g')
plt.ylabel('Phase change [degrees]', color='g')
plt.show()
for low_fs in [100., 200., 400., 500., 800., 1000., 2000.]:
filter='bessel'
low_cut = 0.4*low_fs
b1,a1 = iirfilter(2, low_cut/fs, btype='lowpass', ftype=filter, rp=0.1, rs=0.1)
print (filter, "lowpass for fs=",fs, "cutoff=", low_cut)
b1_scaled,a1_scaled = scaled(b1,a1)
b2,a2 = iirfilter(1, [2*f/low_fs for f in cutoffs], btype='bandpass', ftype=filter, rp=0.1, rs=0.1)
print(filter, "bandpass for fs=", low_fs, "cutoffs=",cutoffs)
b2_scaled,a2_scaled = scaled(b2,a2)
plot(b1_scaled, a1_scaled, b2_scaled, a2_scaled, [10.**(k/200.) for k in range(-400,min(401,int(200.*np.log10(low_fs/2))))])
, where 
. Note that the z values that correspond to sinusoids are along the unit circle, from DC at 
up to the Nyquist frequency 
at 
.
Gas station without pumps 