Beacon detector board

I’m planning to sit in on CMPE 118/L (Mechatronics) this fall, and so I started looking over some of the material for the course from previous years.  One exercise involved designing a “beacon detector” that signals when it sees an infrared signal modulated with a 2kHz square wave. The exercise calls for an all-analog solution (phototransistor, transimpedance amplifier, active filter, rectifier, peak detector, …), which I plan to do when the time comes, but I got intrigued by the idea of doing an almost purely digital design.  That was the motivation for the Goertzel filter blog post.

I decided to take the digital design further and make a beacon detector that not only detects the 2kHz IR beacon, but also indicates what direction it is in.  To do this, I wanted 8 phototransistors around a circle, with one every 45°.  One can get a crude estimate of the angle from just which detector gets the strongest signal, but with wide-angle sensors one should be able to get finer estimates by looking at the ratio of the signals from two adjacent channels.

Because I was deliberately going for minimal analog design as an exercise, I used just a phototransistor and a pullup resistor for each channel, and a Teensy LC board for all the processing.  I chose SFH325FA side-look phototransistors, because they provide a nice, cubical package that I thought would make them easier to align.  They are surface-mount components, but with a 2.54mm pitch for the two terminals, they aren’t much harder to solder than through-hole parts.

For testing, I soldered one of the SFH325FA phototransistors to a pair of male header pins. This allowed me to experiment with different pullup resistor sizes in different lighting conditions and at different distances from an IR emitter.

My experimentation with different pullup resistors indicated that I could not operate in full sunlight, no matter what size pullup resistor I used—when the pullup was small enough that the phototransistor had sufficient voltage across it in full sunlight, the signal from the IR beacon was too small to be useful.  With AC-coupled amplifiers, I could have made it work, but I was committed to nearly pure digital solution.  If I had sunlight nearby, but the phototransistor itself was shaded, then I could go up to 22kΩ for the pullup without problems. The limiting factor was then that a very close beacon could saturate the detector, making it hard to determine power ratios.  I ended up choosing 22kΩ as my pullup, as I wanted to detect beacons from up to 2m away.

With a 22kΩ pullup and a strong signal, I can get a very clean “shark’s fin” waveform, because of the capacitance of the phototransistor acting with the resistor as a low-pass filter.  This low-pass filtering helps remove aliasing artifacts from the sampling by the analog-to-digital filter.

Behind each phototransistor I placed an LED, to indicate which channel had the maximum signal.  I charlieplexed the LEDs, so that I needed only 4 I/O pins for the 8 LEDs (I could encode up to 12 LEDs with 4 charlieplex wires). I used WP710A10SGC green LEDs for the indicators, because I had a lot of them around, but it would have been better to use an amber LED with a wide-angle, diffuse case, as the green indicators are only clearly visible over a fairly narrow angle.  The current-limiting resistors I used would keep the current low enough even with a low-forward-voltage red LED.

I also added an RGB LED to indicate the overall state—I used a common-anode one, because my son had a lot of them with diffuse cases, which are best for an indicator that needs to be seen from any angle.  The RGB LEDs were part of an unidentified lot from China, so I don’t have a part number or specs for them.  I hooked them up to PWM pins, so that I could adjust the brightness and color as needed.

I also designed in a connector for SPI connection, so that the board could be used a slave peripheral of another microcontroller.  I used Eagle on my old laptop to design the PCB (I know that in Need to find new PC board software and PCB CAD tools I said that I would be trying to switch to KiCAD or Diptrace, but I was lazy and stuck with what I already knew how to use).

I sent the design to Seeedstudio to make the PCB—like many of the Chinese firms that cater to hobbyists, they had a sale in July of boards up to 10cm×10cm for only $5.  Of course, it cost me another $21 for shipping with DHL, because I was too impatient to wait for Chinese air mail (and not trusting enough of the reliability). I chose Seeedstudio because they had the best user interface to their design-rule check, and at least as good pricing as anyone else.

The boards arrived quite quickly—only 8 days from order to delivery. Here is what they look like:

The back of the board has some documentation to remind me what is connected where.

The front of the board has relatively little documentation.

The boards are 8cm×8cm, fitting comfortably within both the cheap manufacturing limit and the Eagle free-version size limit.  I made them this big so that the Teensy LC board would fit without interfering with mounting holes (for M2 and M3 screws) and a central hole for attaching the board to a servo arm.

Getting a good right angle for the surface-mount phototransistors took some practice:

My first solder joint for the SFH325FA phototransistor tilted the transistor substantially, because the solder underneath formed a wedge. I reworked it a couple of times and finally got something sort of acceptable.

By the fifth phototransistor, I had worked out a technique that seemed to hold the phototransistor cleanly at a right angle. There is probably not much solder underneath the phototransistor, but the large wedge behind the phototransistor should provide sufficient mechanical strength.

Here are pictures of the populated board:

The back of the board does not look very different after being populated. I put one M3 screw in the bottom M3 screw hole, in order to see how compatible the screw hole design is with the screw head.

The top of the board, when populated shows how close channels 0 and 7 come to the Teensy board.

Here the board is on (powered through the USB cable), detecting an IR beacon about a meter away in the direction of channel 4. The green LED by the channel indicates the direction, and the green RGB LED indicates a strong signal.

I put the header pins on the pullup resistors, so that I could record the signal seen by the analog-to-digital converters.  Unfortunately, the signal is much noisier than I expected:

The large spikes are at 15kHz, and probably correspond to noise injected by the ADC sampling.

I turned off the beacon and looked just at the noise spikes, using the 10× probe on the oscilloscope to get better bandwidth.  I sampled at 100MHz and averaged 1000 traces to get a clean view of the signal (the triggering was set far enough from the background level that only fairly large spikes were captured, but close enough that the 1000 frames were gathered as quickly as possible).

The pulses are very short (2–3 µs), so almost certainly correspond to the charging of the sampling capacitor. After looking at this waveform, I changed my sample time on the ADC to 2µs (it had been 1µs) to reduce the noise in what the ADC reports.

This post has gotten more than long enough, though I still have a couple of things that should be added: a schematic diagram of the board and a plot of the reported angle vs actual angle.

The schematics from Eagle are really ugly, so I’ve been thinking about redrawing them with SchemeIt—perhaps that can be a subsequent post.

The reported angle vs. actual angle plot is going to require building a jig that allow me to set the angle precisely.  I tried jury-rigging something out of Lego today, but the result had too much slop—the board could not be reliably set level with the IR emitter at a fixed angle.  The OED-EL-1L2 5mm IR emitter does fit in the middle hole of a Lego Technic brick, so I using Lego to align the IR emitter was attractive. I may have to make something out of wood and MDF (medium-density fiberboard) for a more solid test jig.  The M3 screw holes will allow me to attach the board firmly to MDF with standoffs, and I can drill holes in Lego bricks to make a stand for the IR emitter.

 

Digital tone detection with Goertzel’s algorithm

In Multiple-feedback bandpass active filter I built a simple active filter to be part of a tone-detection circuit (looking for 2kHz).  Although the filter worked, it did not have as much gain as I expected, its out-of-band rejection was not impressive, and centering it at 20kHz accurately would require adding a potentiometer to tune it. To make a full tone detector, I’d have to have a rectifier and low-pass filter (or some other envelope detector).

I thought that it would be simpler to do the whole thing digitally, as long as my input signal was large enough to be reasonably measured by an analog-to-digital converter.  I started designing a digital filter, with the intent of implementing the same block diagram: bandpass filter, rectifier, and low-pass filter.  Although I came up with a fairly simple filter that looked like it would work well, my son suggested a different approach.  If I don’t need the tone detection to be super-fast, I can process the data a block at a time, and use Goertzel’s algorithm to compute the power at the frequency of interest.  By playing with the block length and windowing functions, I can shape the detection to have different bandwidths and shapes, but with a far steeper rolloff than with the active filter.

Goertzel’s algorithm is a way to compute one coefficient of a discrete Fourier transform (DFT) efficiently with only real-number computation.  In fact, with judicious scaling, I was able to do the computation with just 32-bit integer arithmetic, which can be done very quickly on microcontrollers (like the K20 processor used in the Teensy3.2 boards).

The per-sample computation is just a digital filter with transfer function

where is the center frequency in radians, with being the sampling frequency.  This filter is not quite stable (poles right on the unit circle), so we can’t run it continuously without danger of overflow.  However, if we run it for just N samples, it essentially computes one term of the DFT (with a little post-processing magic):

I chose a sampling rate of 16kHz, so that my 2kHz center frequency is , and the coefficient for is .

The algorithm I used is

sprev = 0
sprev2 = 0
    
for i in range(N):
    x = int(round(2**15 * (1+cos(freq/fs*2*pi*i))))  # full-scale unsigned 16-bit sine wave for testing
    x = x* window[i]//2**14
    s = x + sprev*coeff//scale - sprev2
    sprev2 = sprev
    sprev = s

scaled_sp2 = sprev2//N
scaled_sp = sprev//N
scaled_power =  scaled_sp2 * scaled_sp2 + scaled_sp * scaled_sp - coeff * scaled_sp * scaled_sp2//scale

The square-root of two for the coefficient in the recurrence was represented by coeff=181 and scale=128. The block length (N) was chosen to be 160, so that a block would be processed every 10ms. A 100Hz update rate is fairly typical for systems that control mechanical devices (50Hz is often good enough), and is quick enough not to be a perceptible delay.

I also kept track of the maximum absolute value for s, to make sure that I never exceeded 32 bits anywhere in the computation.
Here is a plot of the scaled_power output, as a function of the sine-wave frequency:

With a rectangular window, the detector has a bandwidth of about 89 Hz, for a Q of 22.5, but the sidelobes are pretty big.

To suppress the sidelobes, I prescaled the signals by a Slepian window, which maximizes the energy in the mainlobe. The Slepian window was computed by the SciPy signal-processing package, scaled to be full-range unsigned 16-bit integers.  The width of the mainlobe is somewhat adjustable—I chose a setting that gave me a somewhat lower Q (wider bandwidth), so that the filter could tolerate some mistuning of the center frequency in the source.

The sidelobes are so small that they get rounded to 0 or 1, while 2kHz is about 384e6, a ratio of about 85dB. The half-power bandwidth is about 226Hz, for a Q of about 8.85.

I could reduce the Slepian width to 0.1, getting a bandwidth of 192Hz (Q of 10.4) and still have the sidelobes fully suppressed. Trying to get higher Q with a block length of 160 starts getting sidelobes, but increasing N allows smaller bandwidths.

I suppose that my next step is to implement this algorithm in C on a Teensy 3.2, connected to an ADC, and see how well it performs in the real world. It looks very promising.

Digital filter lecture

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 and introduced the z-transform .  I explained that this was a linear transform, so that if we multiplied by a constant, we would multiply 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 , which is delayed by one tick, has the z-transform .

I claimed (without proof) that linear filters consisted of delay elements, multiplication by constants, and addition, so that the z-transforms of the input () and output () were related by a transfer function: , and that 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/)

I showed them how this filter had , 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/)

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 . Somewhere in the lecture I mentioned that the poles of (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 , giving a brief explanation of the advantages of using exponentials rather than trig functions and reminding them of the popular abbreviation using angular frequency instead of frequency .

I claimed, without proof, that the response of a filter to a sinusoid with angular frequency was just . 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.

I showed them why any bandpass biquad has essentially the same numerator: we want to have a zero at DC (frequency 0, so at ) and at the Nyquist frequency , so the numerator is always a multiple of .

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.

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.

Optical pulse monitor software

I have previously written several posts about doing an optical pulse monitor as a lab exercise in the applied electronics course or as a project in the freshman design seminar. Most of what I’ve written up as has been the electronics design, with recording using PteroDAQ and post-processing using filters written in Python for display using matplotlib or gnuplot.

I decided today to try writing the software to do the pulse monitor on the Teensy and report the pulse rate out the serial port.  This meant sampling at a fixed rate, doing digital bandpass filtering to remove drift and noise that is not eliminated by the electronics, detecting rising edges in the filtered output (with hysteresis to increase noise immunity), and converting the measured period to beats per minute.  I used the LED on the board to report the pulses detected (flashing the light synchronized with the pulse).

I also added some tests to see if the pulse has been lost, and only report the pulse rate if there have been enough pulses properly detected since the last loss so that a smoothed pulse rate can be reported.

I output the time, the unsmoothed pulse rate, and the smoothed pulse rate to the serial monitor, and I can save that output to a file for plotting:

The movement artifacts around 20s and 40s are greatly reduced by either averaging or median filtering, but median filtering provides a more stable output.
At 70s, I removed my finger from the sensor for about 2s, so that the pulse was lost and had to be reacquired.

 

 

 

Before anyone asks, the pulse rate of 44 bpm is normal for me when resting—this is not likely to be a timing bug.

I’ll also want to try the 240×320 RGB ILI9341 TFT display that I bought, so that I can display the pulse (and maybe the pulse waveform) without needing a USB connection for anything but power.

Reducing 60Hz interference

My optical pulse monitor (with a phototransistor, a transimpedance amplifier, a high-pass filter, and a second stage to increase gain) was having problems with capacitively coupled 60Hz interference. Because the combined gain of the two stages is over 700MΩ at the gain setting I usually use, even very tiny currents at the input cause large output signals. On the design I had in the lab to demo the method, I was getting 60Hz interference that was almost a large as the pulse signal (though grounding myself reduced it by a factor of about 4).

I have been relying on a low-pass filter in the transimpedance amplifier (a capacitor in parallel with the feedback resistor that sets the gain) to reduce the 60Hz interference, and I realized last week that I could set the corner frequency much lower, from around 15Hz down to about 2.2 Hz. The signals I’m interested in are blood pulses, so the interesting signals are about 0.3–3Hz. I don’t need to keep the detailed shape of the pulse waveform, so a little attenuation at the highest and lowest parts of the passband is a good tradeoff for reducing the 60Hz interference over 6-fold.

The transimpedance amplifier I wired up with a rather low corner frequency for the low-pass filter. The filtering is built into the transimpedance stage, by putting a capacitor in parallel with the resistor, so that the gain of the amplifier is (R||Zc).

Here is a plot of the gain of the minimum gain of the two-stage amplifier (with pot set to be 10kΩ):

The narrow bandpass of the new design peaks around 0.95Hz (59.4 bpm), but has adequately high gin over the entire range of likely pulse rates. The old design gave more high-frequency detail to the pulse shape, but at the expense of too much 60Hz pickup.

Although the minimum gain for the amplifier is about 55MΩ, I usually need to set the gain higher (up to around 700MΩ) to get a clear pulse signal to record with the PteroDAQ data acquisition system.