Gas station without pumps

2016 February 26

Digital filter lecture

Filed under: freshman design seminar — gasstationwithoutpumps @ 00:00
Tags: , , ,

On Wed 2016 Feb 24, I gave a lecture in the freshman design seminar on digital filters, covering about 3 weeks worth of material in an hour.  Needless to say, this was a very rushed and handwavey introduction to digital filtering, but it should be enough that the students can (with some additional scaffolding) implement bandpass filters for the pulse monitor project and the ultrasonic rangefinder project.

On Monday, I had covered the notion of digital signals as having discrete values and discrete times with a uniform sampling frequency, so I started with a signal x_0, x_1, \ldots , x_t and introduced the z-transform X(z) = x_0 + x_1 z^{-1} + x_1 z^{-2} +\cdots .  I explained that this was a linear transform, so that if we multiplied x by a constant, we would multiply X(z) by the same amount, and that the z-transform of the sum of two signals was the sum of their transforms.

I also showed that the signal 0, x_0, x_1, \ldots, which is x delayed by one tick, has the z-transform z^{-1} X(z).

I claimed (without proof) that linear filters consisted of delay elements, multiplication by constants, and addition, so that the z-transforms of the input (X(z)) and output (Y(z)) were related by a transfer function: H(z)= Y(z)/X(z), and that H(z) is a rational function for linear filters.

I first showed them a finite-impulse-response (FIR) filter:

Small finite-impulse response filter. (Block diagram drawn with http://www.draw.io/)

Small finite-impulse response filter. (Block diagram drawn with http://www.draw.io/)

I showed them how this filter had H(z) = b_0 + b_1 z^{-1} + b_2 z^{-2}, and could be implemented easily in pseudocode:

    x0 ← new value
    y ← b0 * x0 + b1*x1 + b2*x2
    x2 ← x1
    x1 ← x0

I explained, briefly, what an impulse response was (the values put out by a filter whose input is 1 at time 0 and 0 at all other times), and showed that the filter coefficients were the impulse response.

I also showed them a biquad filter element:

This biquad filter element uses the type 1 direct implementation, which has the advantage of not having any internal overflows, as long as the inputs and outputs remain within bounds. (Block diagram drawn with http://www.draw.io/)

This biquad filter element uses the type 1 direct implementation, which has the advantage of not having any internal overflows, as long as the inputs and outputs remain within bounds. (Block diagram drawn with http://www.draw.io/)

I gave them a simple implementation of the biquad element:

    x0 ← new value
    y ← (b0 * x0 + b1*x1 + b2*x2 - a1*y1 -a2*y2)/a0
    x2 ← x1
    x1 ← x0
    y2 ← y1
    y1 ← y0

I pointed out that multiplication and addition (of integers) was cheap on the Teensy boards, but division or floating-point arithmetic is expensive. Because of binary representation, division by powers of two is cheap, so we can keep the computation fast if we restrict a0 to powers of 2.

I showed how the recurrence relation in the code could be rearranged with simple algebra to get the transfer function H(z) = \frac{b_0 + b_1 z^{-1} + b_2 z^{-2}}{a_0 + a_1 z^{-1} + a_2 z^{-2}}. Somewhere in the lecture I mentioned that the poles of H(z) (that is, the zeros of $latex a_0 + a_1 z^{-1} + a_2 z^{-2}$) had to remain within the unit circle to keep the filter from oscillating, but I didn’t explain why.

I told the students that we would represent sinusoids with e^{j\omega t} = \cos(\omega t) + j \sim(\omega t), giving a brief explanation of the advantages of using exponentials rather than trig functions and reminding them of the popular abbreviation using angular frequency \omega = 2 \pi f instead of frequency f.

I claimed, without proof, that the response of a filter to a sinusoid with angular frequency \omega was just H(e^{j \omega}). I then shared with them a gnuplot script for plotting the response of a biquad element:

#band-pass for pulse monitor
fs = 60
A0= 256
A1= -388
A2= 141
B0 = A0
B1 = 0
B2 = -B0

A0_2 = 256
A1_2 = -449
A2_2 = 199

unset arrow
set arrow nohead from 40000,-20 to 40000,20

unset label
pole1a = ( -A1 + sqrt(A1*A1 - 4 * A0 * A2)) / (2*A0)
pole1b = ( -A1 - sqrt(A1*A1 - 4 * A0 * A2)) / (2*A0)
set label sprintf("pole magnitude %.3f", abs(pole1a)) at 1,-3
set label sprintf("pole1a at %.3f + %.3f j", real(pole1a), imag(pole1a)) at 1,-5
set label sprintf("pole1b at %.3f + %.3f j", real(pole1b), imag(pole1b)) at 1,-7


pole2a = ( -A1_2 + sqrt(A1_2*A1_2 - 4 * A0_2 * A2_2)) / (2*A0_2)
pole2b = ( -A1_2 - sqrt(A1_2*A1_2 - 4 * A0_2 * A2_2)) / (2*A0_2)
set label sprintf("pole2b at %.3f + %.3f j", real(pole2b), imag(pole2b)) at 1,-17
set label sprintf("pole2a at %.3f + %.3f j", real(pole2a), imag(pole2a)) at 1,-15
set label sprintf("pole magnitude %.3f", abs(pole2a)) at 1,-13

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

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

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

set grid xtics
set mxtics 10

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

set samples 5000

plot \
    20*log10(gain(x,B0,B1,B2, A0,A1,A2)) \
	title sprintf("%.0f (1 + %.0f z^-1 + %.0f z^-2)/(%.2f+ %.3f z^-1 + %.3f z^-2)", \
			B0, B1/B0, B2/B0, A0, A1, A2), \
    20*log10(gain(x,B0,B1,B2, A0_2,A1_2,A2_2)) \
	title sprintf("%.0f (1 + %.0f z^-1 + %.0f z^-2)/(%.2f+ %.3f z^-1 + %.3f z^-2)", \
			B0, B1/B0, B2/B0, A0_2, A1_2, A2_2) lt 3

I walked them through the code and showed them the result:

This is a pair of bandpass filters designed for a pulse monitor. The red curve may be a better choice, as it does not have such a sharp peak around 1.5Hz, but still reasonably suppresses the high frequencies and the DC drift.

This is a pair of bandpass filters designed for a pulse monitor. The red curve may be a better choice, as it does not have such a sharp peak around 1.5Hz, but still reasonably suppresses the high frequencies and the DC drift.

I showed them why any bandpass biquad has essentially the same numerator: we want to have a zero at DC (frequency 0, so at e^{j 2 \pi 0}=1) and at the Nyquist frequency $e^{j 2 \pi \frac{1}{2}}=-1$, so the numerator is always a multiple of (1-z^{1})(1+z^{-1})= 1-z^{-2}.

I also showed them the effect of having just the numerator (setting the denominator polynomial to 1), using both log and linear frequency scales to show the peak at half the Nyquist frequency (one quarter the sampling frequency):

A linear frequency scale make the symmetric frequency response of the FIR filter 256 (1–z^-2) clear.

A linear frequency scale make the symmetric frequency response of the FIR filter 256 (1–z^-2) clear.

I told them that I did not really know the details of how to choose specific filter parameters, and shared with them the code I used from the scipy.signal package for choosing the parameters:

#!/usr/bin/env python3

from __future__ import print_function, division

from scipy import signal
from cmath import sqrt  # complex square root

print("set parameters by giving name=value lines:")

sampling_freq = 30 # Hz
low_cutoff = 0.3        # Hz
high_cutoff = min(150, 0.49*sampling_freq) # Hz
fixed_point_scaling=64

while True:
    print("sample=",sampling_freq, "low=", low_cutoff, "high=", high_cutoff,
        "scale=", fixed_point_scaling)

    hi_over_Nyquist = high_cutoff/(0.5*sampling_freq)
    lo_over_Nyquist = low_cutoff/(0.5*sampling_freq)

    filter_float = signal.bessel(1,
                        [lo_over_Nyquist,hi_over_Nyquist], 
                        btype="bandpass")
    A= [int(x*fixed_point_scaling+0.5) for x in filter_float[1]]
    discr = A[1]*A[1] -4*A[0]*A[2]
    poles= [(-A[1]+sqrt(discr))/(2*A[0]), (-A[1] -sqrt(discr))/(2*A[0])]
    
    print("\tA=",A, "B=[1,0,-1]", "poles=",poles)
    
    line=input("param=value?")
    fields = line.split("=")
    if len(fields)==0: continue
    if len(fields)!=2:
        print("No parameter change found, exiting")
        break
    keyword=fields[0].strip().lower()
    value=float(fields[1])
    if keyword=="low" or keyword=="lo":
        low_cutoff = value
    elif keyword=="high" or keyword=="hi":
        high_cutoff = value
    elif keyword=="sample" or keyword=="freq":
        sampling_freq = value
    elif keyword=="scale":
        fixed_point_scaling=int(value)
    else:
        print("unrecognized keyword:", keyword)

I only walked them through a little of this code (mainly showing the scaling of frequencies to the Nyquist frequency, since that is what the signal package wants, and the call to the bessel function to get a Bessel filter). I explained that the fixed_point_scaling parameter was to get more resolution in the parameters when doing integer arithmetic, but didn’t really have time to explain what that meant. I did demonstrate setting the parameters to get one of the filters shown in the graph above.

On Friday, I plan to give them a little code snippet that I use:

static volatile     int32_t x_0, x_1, x_2;
static volatile     int32_t y_0, y_1, y_2;

// filter parameters for biquad bandpass filter

// selected for approx 0.66--6Hz with 60Hz sampling
#define SAMPLE_FREQ	(60)	// sampling frequency in Hz
#define a0  (256)
#define a1  (-388)  
#define a2  (141) 

#define gain (1)
  // b0= - b2= gain*a0
  // b1=0

#define DELAY_XY (x_2=x_1, x_1=x_0, y_2=y_1, y_1=y_0)
#define GENERAL_BANDPASS (y_0 =  ((gain*a0)*(x_0-x_2) -a1*y_1 -a2*y_2)/a0, DELAY_XY)

To use the filter, I include

    x_0=analogRead(MONITOR_PIN);
    GENERAL_BANDPASS;

in an interrupt routine that runs once every 60th of a second. I’ll also have to show them how to use the IntervalTimer in the Teensyduino software development kit to set up the interrupt routine.

Leave a Comment »

No comments yet.

RSS feed for comments on this post. TrackBack URI

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Blog at WordPress.com.

%d bloggers like this: