Gas station without pumps

2014 July 3

Digital filters: infinite-impulse-response filters

Filed under: Uncategorized — gasstationwithoutpumps @ 17:13
Tags: , , , , , ,

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: H(z) = \frac{b_{0} + b_{1}z^{-1} + b_{2}z^{-2}}{1 + a_{1}z^{-1} + a_{2}z^{-2}}.

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

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: Y(z) = X(z)(b_{0} + b_{1}z^{-1} + b_{2}z^{-2}) - Y(z)(a_{1}z^{-1} + a_{2}z^{-2}), followed by simple algebra to rearrange to get the transfer function above. One nice property of this particular implementation for fixed-point computation is that as long as the output stays within range, there are no problems with overflow in intermediate computations.

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, Ae^{i\theta}, since the resulting filter will have minima at \omega=\theta, and how sharp the minimum is depends on how close A gets to 1. One popular choice for bandpass filters is to put the zeros at 1 and –1, so that there are zeroes at DC (ω=0) and at the Nyquist frequency (ω=π), which gives 1-z^{-2} as the numerator. For a low-pass filter, both zeros are often put at the Nyquist frequency, giving (1+z^{-1})^2, while for a high-pass filter, both zeros are often put at DC, giving (1-z^{-1})^2.

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.

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 2 to the 14, but the  filter response is quite nice.

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.

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))))])

2013 July 27

Failed attempt at pulse oximeter

In Optical pulse monitor with little electronics  and Digital filters for pulse monitor, I developed an optical pulse monitor using an IR emitter, a phototransistor, 2 resistors, and an Arduino.  On Thursday, I decided to try to extend this to a pulse oximeter, by adding a red LED (and current-limiting resistor) as well.  Because excluding ambient light is so important, I decided to build a mount for everything out of a block of wood:

Short piece of 2x2 wood, with a 3/4" diameter hoe drilled with a Forstner bit part way through the block.  Two 1/8" holes drilled for 3mm LEDs on top, and one for a 3mm phototransistor on the bottom (lined up with the red LED).  Wiring channels were cut with the same 1/8" drill bit, and opened up a with a round riffler.  Electrical tape holds the LEDs and phototransistor in place (removed here to expose the diodes).

Short piece of 2×2 wood, with a 3/4″ diameter hole drilled with a Forstner bit partway through the block. Two 1/8″ holes drilled for 3mm LEDs on top, and one for a 3mm phototransistor on the bottom (lined up with the red LED). Wiring channels were cut with the same 1/8″ drill bit, and opened up a with a round riffler. Electrical tape holds the LEDs and phototransistor in place (removed here to expose the diodes).

My first test with the new setup was disappointing.  The signal from the IR LED swamped out the signal from the red LED, being at least 4 times as large. The RC discharge curves for the phototransistor for the IR signal was slow enough that I would have had to go to a very low sampling rate to see the red LED signal without interference from the discharge from the IR pulse.  I could reduce the signal for the IR LED to only twice the red output by increasing the IR current-limiting resistor to 1.5kΩ, and reduce the RC time constant of the phototransistor by reducing the pulldown resistor for it to 100kΩ The reduction in the output of the IR LED and decreased sensitivity of the phototransistor made about a 17-fold reduction in the amplitude of the IR signal, and the red signal was about a thirtieth of what I’d previously been getting for the IR signal.  Since the variation in amplitude that made up my real signal was about 10 counts before, it is substantially less than 1 count now, and is  too small to be detected even with the digital filters that I used.

I could probably solve this problem of a small signal by switching from the Arduino to the KL25Z, since going from a 10-bit ADC to a 16-bit ADC would allow a 64 times larger signal-to-noise ratio (that is, +36dB), getting me back to enough signal to be detectable even with the reductions..  I’ve ordered headers from Digi-Key for the KL25Z, so next week I’ll be able to test this.

I did do something very stupid yesterday, though in a misguided attempt to fix the problem.  I had another red LED (WP710A10ID) that was listed on the spec sheet as being much brighter than the one I’d been using (WP3A8HD), so I soldered it in.  The LED was clearly much brighter, but when I put my finger in the sensor, I got almost no red signal!  What went wrong?

A moment’s thought explained the problem to me (I just wish I had done that thinking BEFORE soldering in the LED).  Why was the new LED brighter for the same current?  It wasn’t that the LED was more efficient at generating photons, but that the wavelength of the light was shorter, and so the eye was more sensitive to it.

Spectrum of the WP3A8HD red LED that I first used.  It has a peak at 700nm and dominant wavelength at 660nm.  I believe that the "dominant wavelength" refers to the peak of the spectrum multiplied by the sensitivity of the human eye.

Spectrum of the WP3A8HD red LED that I first used. It has a peak at 700nm and dominant wavelength at 660nm. I believe that the “dominant wavelength” refers to the peak of the spectrum multiplied by the sensitivity of the human eye.  Spectrum copied from Kingbright preliminary specification for WP3A8HD.

Spectrum of the WP710A10ID brighter red LED.  The peak is at 627nm and the "dominant wavelength" is 617nm.  The extra brightness is coming from this shorter wavelength, where the human eye is more sensitive.

Spectrum of the WP710A10ID brighter red LED that didn’t work for me. The peak is at 627nm and the “dominant wavelength” is 617nm. The extra brightness is coming from this shorter wavelength, where the human eye is more sensitive. Image copied from the Kingbright spec sheet.

CIE_1931 luminosity curve, representing a stndardized sensitivity of the human eye under high-light conditions (scotopic vision).  Copied from Wikipedia: http://upload.wikimedia.org/wikipedia/commons/thumb/7/72/CIE_1931_Luminosity.png/640px-CIE_1931_Luminosity.png

1931 CIE luminosity curve, representing a standardized sensitivity of the human eye with bright lighting (photopic vision). The peak is at 555nm. Note that there are better estimates of human eye sensitivity now available (see the discussion of newer ones in the Wikipedia article on the Luminosity function).
Image copied from Wikipedia.

The new LED is brighter, because the human eye is more sensitive to its shorter wavelength, but the optimum sensitivity of the phototransistor is at longer wavelengths, so the phototransistor is less sensitive to the new LED than to the old one.

Typical spectral sensitivity of a silicon photodiode or phototransistor.  This curve does not take into account any absorption losses in the packaging of the part, which can substantially change the response.  Unfortunately, Kingbright does not publish a spectral sensitivity curve for their WP3DP3B phototransistor.  Image copied from https://upload.wikimedia.org/wikipedia/commons/4/41/Response_silicon_photodiode.svg

Typical spectral sensitivity of a silicon photodiode or phototransistor. This curve does not take into account any absorption losses in the packaging of the part, which can substantially change the response. Note that the peak sensitivity is in the infrared, around 950nm, not in the green around 555nm as with the human eye. Unfortunately, Kingbright does not publish a spectral sensitivity curve for their WP3DP3B phototransistor, so this image is a generic one copied from https://upload.wikimedia.org/wikipedia/commons/4/41/Response_silicon_photodiode.svg

This sensitivity is much better matched to the IR emitter (WP710A10F3C) than to either of the red LEDs:

Spectrum for the WP710A10F3C IR emitter, copied from the Kingbright spec sheet.  The peak is at 940nm with a 50nm bandwidth.  There is no "dominant wavelength", because essentially all the emissions are outside the range of the human eye.

Spectrum for the WP710A10F3C IR emitter, copied from the Kingbright spec sheet. The peak is at 940nm with a 50nm bandwidth. There is no “dominant wavelength”, because essentially all the emissions are outside the range of the human eye.

Furthermore, blood and flesh is more opaque at the shorter wavelength, so I had more light absorbed and less sensitivity in the detector, making for a much smaller signal.

Scott Prahl's estimate of oxyhemoglobin and deoxyhemoglobin molar extinction coefficients, copied from http://omlc.ogi.edu/spectra/hemoglobin/summary.gif The higher the curve here the less light is transmitted.  Note that 700nm has very low absorption, but 627nm has much higher absorption.

Scott Prahl’s estimate of oxyhemoglobin and deoxyhemoglobin molar extinction coefficients, copied from http://omlc.ogi.edu/spectra/hemoglobin/summary.gif
Tabulated values are available at http://omlc.ogi.edu/spectra/hemoglobin/summary.html and general discussion at http://omlc.ogi.edu/spectra/hemoglobin/
The higher the curve here the less light is transmitted. Note that 700nm has very low absorption (290), but 627nm has over twice as high an absorption (683).  Also notice that in the infrared

I had to go back to the red LED (WP3A8HD) that I started with. Here is an example of the waveform I get with that LED, dropping the sampling rate to 10Hz:

The green waveform is the voltage driving the red LED and through a 100Ω resistor.  The red LED is on for the 1/30th of second that the output is low, then the IR LED is on (through a 1.5kΩ resistor) for 1/30th of a second, then both are off.  THe yellow trace shows the voltage at the phototransistor emitter with a 680kΩ pulldown. This signal seems to have too little amplitude for the variation to be detected with the Arduino.

The green waveform is the voltage driving the red LED and through a 100Ω resistor. The red LED is on for the 1/30th of second that the output is low, then the IR LED is on (through a 1.5kΩ resistor) for 1/30th of a second, then both are off. THe yellow trace shows the voltage at the phototransistor emitter with a 680kΩ pulldown.
This signal seems to have too little amplitude for the variation to be detected with the Arduino (the scale is 1v/division with 0v at the bottom of the grid).

I can try increasing the signal by using 2 or more red LEDs (though the amount of current needed gets large), or I could turn down the IR signal to match the red signal and use an amplifier to get a big enough signal for the Arduino to read.  Sometimes it seems like a 4.7kΩ resistor on the IR emitter matches the output, and sometimes there is still much more IR signal received, depending on which finger I use and how I hold it in the device.

I was thinking of playing with some amplification, but I could only get a gain of about 8, and even then I’d be risking saturation of the amplifier.  I think I’ll wait until the headers come and I can try the KL25Z board—the gain of 64 from the higher resolution ADC is likely to be more useful.  If that isn’t enough, I can try adding gain also.  I could also eliminate the “off-state” and just amplify the difference between IR illumination and red illumination.  I wonder if that will let me detect the pulse, though.

2013 July 24

Digital filters for pulse monitor

In Optical pulse monitor with little electronics, I talked a bit about an optical pulse monitor using the Arduino and just 4 components (2 resistors, an IR emitter, and a phototransistor).  Yesterday, I had gotten as far as getting good values for resistors, doing synchronous decoding, and using a very simple low-pass IIR filter to clean up the noise.  The final result still had problems with the baseline shifting (probably due to slight movements of my finger in the sensor):

With digital low-pass filtering, the pulse signal is much cleaner, but the sharp downward transition at the start of each pulse has been rounded off by the filter.

(click to embiggen) Yesterday’s plot with digital low-pass filtering, using y(t) = (x(t) + 7 y(t-1) )/8.  There is not much noise, but the baseline wobbles up and down a lot, making the signal hard to process automatically.

Today I decided to brush off my digital filter knowledge, which I haven’t used much lately, and see if I could design a filter using only small integer arithmetic on the Arduino, to clean up the signal more. I decided to use a sampling rate fs = 30Hz on the Arduino, to avoid getting any beating due to 60Hz pickup (not that I’ve seen much with my current setup). The 30Hz choice was made because I do two measurements (IR on and IR off) for each sample, so my actual measurements are at 60Hz, and should be in the same place in any noise waveform that is picked up. (Europeans with 50Hz line frequency would want to use 25Hz as their sampling frequency.)

With the 680kΩ resistor that I selected yesterday, the 30Hz sampling leaves plenty of time for the signal to charge and discharge:

The grid line in the center is at 3v.  The green trace is the signal to on the positive side of the IR LED, so the LED is on when the trace is low (with 32mA current through the pullup resistor).  The yellow trace is the voltage at the Arduino input pin: high when light is visible, low when it is dark.  This recording was made with my middle finger between the LED and the phototransistor.

The grid line in the center is at 3v. The green trace is the signal to on the positive side of the IR LED, so the LED is on when the trace is low (with 32mA current through the pullup resistor). The yellow trace is the voltage at the Arduino input pin: high when light is visible, low when it is dark. This recording was made with my middle finger between the LED and the phototransistor.

I decided I wanted to replace the low-pass filter with a passband filter, centered near 1Hz (60 beats per minute), but with a range of about 0.4Hz (24 bpm) to 4Hz (240bpm). I don’t need the passband to be particularly flat, so I decided to go with a simple 2-pole, 2-zero filter (called a biquad filter). This filter has the transfer function

H(z) = \frac{b_{0} + b_{1}z^{-1} + b_{2}z^{-2}}{1+a_{1}z^{-1}+a_{2}z^{-2}}

To get the gain of the filter at a frequency f, you just compute \left| H( e^{i \omega} ) \right|, where \omega = 2 \pi f / f_{s}.  Note that the z values that correspond to sinusoids are along the unit circle, from DC at e^{0} = 1 up to the Nyquist frequency f_{s}/2 at e^{\pi} = -1.

The filter is implemented as a simple recurrence relation between the input x and the output y:

y(t) = b_{0} x(t) + b_{1}x(t-1) + b_{2}x(t-2) - a_{1}y(t-1) - a_{2}y(t-2)

This is known as the “direct” implementation.  It takes a bit more memory than the “canonical” implementation, but has some nice properties when used with small-word arithmetic—the intermediate values never get any further from 0 than the output and input values, so there is no overflow to worry about in intermediate computations.

I tried using an online web tool to design the filter http://www-users.cs.york.ac.uk/~fisher/mkfilter/, and I got some results but not everything on the page is working.  One can’t very well complain to Tony Fisher about the maintenance, since he died in 2000. I tried using the tool at http://digitalfilter.com/enindex.html to look at filter gain, but it has an awkward x-axis (linear instead of logarithmic frequency) and was a bit annoying to use.  So I looked at results from Tony Fisher’s program, then used my own gnuplot script to look at the response for filter parameters I was interested in.

The filter program gave me one obvious result (that I should not have needed a program to realize): the two zeros need to be at DC and the Nyquist frequency—that is at ±1.  That means that the numerator of the transfer function is just 1-z^{-2}, and b0=1, b1=0, and b2=–1.  The other two parameters it gave me were a2=0.4327386423 and a1=–1.3802466192.  Of course, I don’t want to use floating-point arithmetic, but small integer arithmetic, so that the only division I do is by powers of 2 (which the compiler turns into a quick shift operation).

I somewhat arbitrarily selected 32 as my power of 2 to divide by, so that my transfer function is now

H(z) = 32 \frac{1 - z^{-2}}{32+A_{1}z^{-1}+A_{2}z^{-2}}

and my recurrence relation is

y(t) = \left(32 \left( x(t) - x(t-2) \right) - A_{1} y(t-1) - A_{2} y(t-2) \right)/32

with A1 and A2 restricted to be integers.  Rounding the numbers from Fisher’s program suggested A1=-44 and A2=14, but that centered the filter at a bit higher frequency than I liked, so I tweaked the parameters and drew plots to see what the gain function looked like.  I made one serious mistake initially—I neglected to check that the two poles were both inside the unit circle (they were real-valued poles, so the check was just applying the quadratic formula).  My first design (not the one from Fisher’s program) had one pole outside the unit circle—it looked fine on the plot, but when I implemented it, the values grew until the word size was exceeded, then oscillated all over the place.  When I realized what was wrong, I checked the stability criterion and changed the A2 value to make the pole be inside the unit circle.

I eventually ended up with A1=-48 and A2=17, which centered the filter at 1, but did not have as high an upper frequency as I had originally thought I wanted:

The gain of the filter that I ended up implementing has -3dB points at about 0.43 and 2.15 Hz.

(click to embiggen) The gain of the filter that I ended up implementing has -3dB points at about 0.43 and 2.15 Hz.

Here is the gnuplot script I used to generate the plot—it is not fully automatic (the xtics, for example, are manually set). Click it to expand.

fs = 30	# sampling frequency
A0=32.  # multiplier (use power of 2)
b=16.

A1=-(A0+b)
A2=b+1

peak = fs/A0	# approx frequency of peak of filter

set title sprintf("Design of biquad filter, fs=%3g Hz",fs)

set key bottom center
set ylabel "gain [dB]"
unset logscale y
set yrange [-20:30]

set xlabel "frequency [Hz]"
set logscale x
set xrange [0.01:0.5*fs]

set xtics add (0.43, 2.15)
set grid xtics

j=sqrt(-1)
biquad(zinv,b0,b1,b2,a0,a1,a2) = (b0+zinv*(b1+zinv*b2))/(a0+zinv*(a1+zinv*a2))
gain(f,b0,b1,b2,a0,a1,a2) = abs( biquad(exp(j*2*pi*f/fs),b0,b1,b2,a0,a1,a2))
phase(f,b0,b1,b2,a0,a1,a2) = imag(log( biquad(exp(j*2*pi*f/fs),b0,b1,b2,a0,a1,a2)))

plot 20*log(gain(x,A0,0,-A0,  A0,A1,A2)) \
		title sprintf("%.0f (1-z^-2)/(%.0f+ %.0f z^-1 + %.0f z^-2)", \
			A0, A0, A1, A2), \
	20*log(gain(peak,A0,0,-A0,  A0,A1,A2))-3 title "approx -3dB"

I wrote a simple Arduino program to sample the phototransistor every 1/60th of a second, alternating between IR off and IR on. After each IR-on reading, I output the time, the difference between on and off readings, and the filtered difference. (click on the code box to view it)

#include "TimerOne.h"

#define rLED 3
#define irLED 5

// #define CANONICAL   // use canonical, rather than direct implementation of IIR filter
// Direct implementation seems to avoid overflow better.
// There is probably still a bug in the canonical implementation, as it is quite unstable.

#define fs (30) // sampling frequency in Hz
#define half_period (500000L/fs)  // half the period in usec

#define multiplier  32      // power of 2 near fs
#define a1  (-48)           // -(multiplier+k)
#define a2  (17)            // k+1

volatile uint8_t first_tick;    // Is this the first tick after setup?
void setup(void)
{
    Serial.begin(115200);
//    pinMode(rLED,OUTPUT);
    pinMode(irLED,OUTPUT);
//    digitalWrite(rLED,1);  // Turn RED LED off
    digitalWrite(irLED,1); // Turn IR LED off

    Serial.print("# bandpass IIR filter\n# fs=");
    Serial.print(fs);
    Serial.print(" Hz, period=");
    Serial.print(2*half_period);
    Serial.print(" usec\n#  H(z) = ");
    Serial.print(multiplier);
    Serial.print("(1-z^-2)/(");
    Serial.print(multiplier);
    Serial.print(" + ");
    Serial.print(a1);
    Serial.print("z^-1 + ");
    Serial.print(a2);
    Serial.println("z^-2)");
#ifdef CANONICAL
    Serial.println("# using canonical implementation");
#else
    Serial.println("# using direct implementation");
#endif
    Serial.println("#  microsec raw   filtered");

    first_tick=1;
    Timer1.initialize(half_period);
    Timer1.attachInterrupt(half_period_tick,half_period);
}

#ifdef CANONICAL
// for canonical implementation
 volatile int32_t w_0, w_1, w_2;
#else
// For direct implementation
 volatile int32_t x_1,x_2, y_0,y_1,y_2;
#endif

void loop()
{
}

volatile uint8_t IR_is_on=0;    // current state of IR LED
volatile uint16_t IR_off;       // reading when IR is off (stored until next tick)

void half_period_tick(void)
{
    uint32_t timestamp=micros();

    uint16_t IR_read;
    IR_read = analogRead(0);
    if (!IR_is_on)
    {   IR_off=IR_read;
        digitalWrite(irLED,0); // Turn IR LED on
        IR_is_on = 1;
        return;
    }

    digitalWrite(irLED,1); // Turn IR LED off
    IR_is_on = 0;

    Serial.print(timestamp);
    Serial.print(" ");

    int16_t x_0 = IR_read-IR_off;
    Serial.print(x_0);
    Serial.print(" ");

 #ifdef CANONICAL
    if (first_tick)
    {  // I'm not sure how to initialize w for the first tick
       w_2 = w_1 = multiplier*x_0/ (1+a1+a2);
       first_tick = 0;
    }
 #else
    if (first_tick)
    {   x_2 = x_1 = x_0;
        first_tick = 0;
    }
#endif

#ifdef CANONICAL
    w_0 = multiplier*x_0 - a1*w_1 -a2*w_2;
    int32_t y_0 = w_0 - w_2;
    Serial.println(y_0);
    w_2=w_1;
    w_1=w_0;
#else
     y_0 = multiplier*(x_0-x_2) - a1*y_1 -a2*y_2;
     Serial.println(y_0);
     y_0 /= multiplier;
     x_2 = x_1;
     x_1 = x_0;
     y_2 = y_1;
     y_1 = y_0;
#endif
}

Here are a couple of examples of the input and output of the filtering:

    The input signals here are fairly clean, but different runs often get quite different amounts of light through the finger, depending on which finger is used and the alignment with the phototransistor. Note that the DC offset shifts over the course of each run.

(click to embiggen) The input signals here are fairly clean, but different runs often get quite different amounts of light through the finger, depending on which finger is used and the alignment with the phototransistor. Note that the DC offset shifts over the course of each run.

IR-biquad-filtered

(click to embiggen) After filtering the DC offset and the baseline shift are gone. The two very different input sequences now have almost the same range. There is a large, clean downward spike at the beginning of each pulse.

Overall, I’m pretty happy with the results of doing digital filtering here. Even a crude 2-zero, 2-pole filter using just integer arithmetic does an excellent job of cleaning up the signal.

%d bloggers like this: