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.

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

Digital filters: z-transform and finite-impulse response filters

Filed under: Uncategorized — gasstationwithoutpumps @ 00:32
Tags: , , , ,

My son decided to try implementing a loudness sensor on the KL25Z board—last night he designed a circuit to interface an electret microphone to the A-to-D on the KL25Z, wired it up, wrote a small test program, and got an LED to be controlled by the loudness on the microphone.  He still needs to adjust the gain of his amplifier (the sensitivity was too low) and improve his loudness detection algorithm.  Currently the system only responds to very loud noises, but he’d like it to be able to pick out the beat in music even if the music is being played fairly softly.

I think that the best way to do what he wants is with digital filters (and some non-linear computations, like absolute value or squaring), but he doesn’t have any filter theory.  I’m going to try writing a few blog posts to cover the basics of implementing some simple digital filters, so that he can try out his ideas.

This post is going to cover the fundamental mathematical tool underlying most digital signal processing: the z-transform (which is related to the Laplace and Fourier transforms in continuous signal processing).

If we have a semi-infinite sequence of numbers, x_{0}, x_{1}, \ldots, then the z-transform of the sequence X(z) is \sum_{k=0}^{\infty} x_{k} z^{-k}.  For example, for the constant sequence x_{k}=c, the z-transform is X(z) = c \sum_{k=0}^{\infty} z^{-k} = \frac{1}{1-z}.  For the sinusoidal function with angular frequency \omega radians per sample, amplitude A, and phase \phi, x_{k} = A e^{i(\omega k +\phi)}, the z-transform is X(z) = \frac{Ae^{i\phi}}{1- e^{-i\omega}z}.

The z-transform is a linear transform—that is, if we take any two sequences x_{k} and y_{k} with z-transforms X(z) and Y(z), then the z-transform of a x_{k} +b y_{k} is a X(z)+ b Y(z).

One thing that is cool about z-transforms is that they can easily represent what happens if we delay a signal.  Consider the two sequences x_{0}, x_{1}, x_{2}, \ldots and 0, x_{0}, x_{1}, \ldots .  If the z-transform of the first sequences is X(z), then the z-transform of the second sequence is just z^{-1} X(z).

A linear filter consists just of multiplications by constants, delay elements, and additions.  The simplest filters (mathematically) are finite-impulse response filters:

Typical finite-impulse response (FIR) filter

Typical finite-impulse response (FIR) filter

The output y_{k} is computed as y_{k}=a_{0} x_{k} + a_{1}x_{k-1} + a_{2} x_{k-2} + \cdots a_{5} x_{k-5}, and the z-transform is Y(z) = a_{0} X(z) + a_{1} X(z) z^{-1} + a_{2} X(z) z^{-2} + \cdots a_{5} X(z) z^{-5}. We can simplify this to Y(z)/X(z) = a_{0} + a_{1} z^{-1} + a_{2} z^{-2} + \cdots a_{5} z^{-5}, a simple polynomial in z^{-1}. The function Y(z)/X(z) is known as the transfer function of the filter. The transfer function is sometimes written as H(z).

The impulse response of a filter is the output it provides when the input is 1 for the first sample and 0 thereafter. The z-transform of an impulse is the constant function 1, so the z-transform of the impulse response of a filter is just the transfer function for the filter. For the FIR filter in the picture above, the impulse response is a_{0}, a_{1}, a_{2}, a_{3}, a_{4}, a_{5}, 0, 0, \ldots. Because the response is zero after a while, this is a finite impulse response.

When we are designing or analyzing filters, we often want to think about what they do to sinusoids, since we know that adding, multiplying, and delaying sinusoids of a given frequency e^{i \omega t} produces another sinusoid, but with a different amplitude and phase A e^{i (\omega t + \phi)}.

How do we convert the transfer function of the filter into the phasor A e^{i \phi}?

Note that delaying a sinusoid by 1 time unit is the same as subtracting \omega from its phase, so multiplying by the z-transform by z^{-1} must be the same as multiplying by e^{-i\omega} or setting z=e^{i\omega}.

So we can get the phasor for a given angular frequency just by plugging z=e^{i\omega} into the transfer function! This makes plotting the amplitude and phase response of a digital filter very simple.

Let’s look at a very simple filter averaging the last 4 samples, so y_{k}= 1/4 (x_{k} + x_{k-1} + x_{k-2} + x_{k-3}), and H(z)= 1/4 (1+z^{-1}+z^{-2}+z^{-3}). If we plot the gain (the absolute value of Y/X) as a function of frequency f, we can see that this is a low-pass filter:

This is a low-pass filter, with the half-power point, where the gain is 0.707, around 0.11536 cycles per sample, and gain goes to zero at frequencies of 1/4 and 1/2 (ω=π/2 or π), but there is a substantial extra peak around a frequency of 0.36614 cycles/sample. It is possible to design low-pass FIR filters with better characteristics than this one!

This is a low-pass filter, with the half-power point, where the gain is 0.707, around 0.11536 cycles per sample, and gain goes to zero at frequencies of 1/4 and 1/2 (ω=π/2 or π), but there is a substantial extra peak around a frequency of 0.36614 cycles/sample. It is possible to design low-pass FIR filters with better characteristics than this one!

The phase change is small for small frequencies (as the output is only delayed about 1.5 samples from the input), but does funny things in the higher frequencies.

The phase change is small for small frequencies (as the output is only delayed about 1.5 samples from the input), but does funny things in the higher frequencies.

Here is the gnuplot source code for the gain plot above:

set title "Averaging 4 adjacent samples"
set key bottom left
set samples 10000

set ylabel "Gain (y/x)"
set logscale y
set yrange [0.009:1.1]

set xlabel "frequency [cycles/sample]"
set logscale x
set xrange [0.001:0.5]

# transfer function
H(z)= (1+z**(-1)+z**(-2)+z**(-3))/4

j=sqrt(-1)
amplitude(omega) = abs(H(exp(j*omega)))
phase(omega) = imag(log(H(exp(j*omega)))) * 180/pi  # phase in degrees

plot amplitude(2*pi*x) notitle

I rarely use FIR filters (other than when I’m being very lazy and just averaging a few samples for downsampling), but usually use infinite-response filters, which I’ll leave for a subsequent post.

%d bloggers like this: