Gas station without pumps

2017 August 6

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.

 

2017 July 9

Digital tone detection with Goertzel’s algorithm

Filed under: Robotics,Uncategorized — gasstationwithoutpumps @ 00:28
Tags: , ,

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

\frac{S(z)}{X(z)} = \frac{1}{1 - 2 \cos(\omega_0) z^{-1} + z^{-2}}  = \frac{1}{\left(1 - e^{+j \omega_0} z^{-1}\right)\left(1 - e^{-j \omega_0} z^{-1}\right)}~,  

where \omega_0 = 2\pi f_0/f_s is the center frequency in radians, with f_s 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):

y[n] = e^{j 2 \pi \frac{k}{N}} s[N-1] - s[N-2] ~.

I chose a sampling rate of 16kHz, so that my 2kHz center frequency is \omega_0=\pi/4, and the coefficient for z^{-1} is \sqrt{2}.

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.

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.

2015 May 14

Blood pressure lab

Despite fairly poor prelab homework turned in, the first half of the blood pressure lab went well.  After seeing how poorly students were doing on breaking down the problem into pieces (perhaps the main transferable engineering skill I’m trying to get them to develop), I ended up giving them more explicit instructions on the board at the beginning of lab:

  1. calculate sensor voltage difference for 100mmHg with 3.3v power
  2. measure sensor voltage difference for 100mmHg with 3.3v power (also 0mmHg and -100mmHg)
  3. determine upper and lower inputs of voltage for instrumentation amp INA126P from the data sheet, using worst-case rather than typical specs (“worst-case” meaning the smallest remaining voltage range)
  4. use Vout-Vref = G(V+ – V) to determine maximum gain to avoid clipping if input swing is +-180mmHg (+-24kPa)
  5. compute needed gain resistor, wire it up (and virtual ground)
  6. measure voltage at output of instrumentation amp at 0, +100mmHg, -100mmHg
  7. compute gain needed in second stage to get maximum range (without clipping) at final output
  8. wire up op amp and measure final output voltage at 0, +100mmHg, -100mmHg
  9. What is Vout as function of pressure?
  10. record with the PteroDAQ a blood pressure measurement with pressure slowly decaying from 180mmHg down to 40mmHg (not too slowly, or your hand will get swollen).  Check for clipping at high end.  Check that you are using nearly the full range. Check that pulsations are visible when plotting the data.
  11. Use bandpass-filter.py to filter the first channel of the recording (later channels will be discarded)

I may have to put some version of these instructions in the book, though this sort of hand-holding is precisely what I’m trying to cut out in the “descaffolding”.  I’m afraid we’re training a generation of technicians rather than engineers—they’re good at following very explicit instructions, but not so good at breaking problems down into smaller problems.

With these explicit instructions, most of the students managed to get breadboard versions of the pressure sensor amplifiers working. I may have to help out bench 4, as it turned out that their pressure sensor seems to have a 0.7mV offset (which is pretty big—way out of spec).  They’ll have to decide whether to change benches to get a different sensor, compensate for the sensor offset electronically, or compensate for it in the post-processing of the data.  Any of these solutions would be acceptable, but they aren’t all equally easy.

The students needed less help than in previous years in the lab, so I think that having the students struggle with the prelabs, even if they don’t get the answers right, is helping make the lab time more efficient—they only have to get past a couple of misunderstandings, rather than trying to learn all the material for the first time in lab, as so many did the last couple of years.

In lecture on Wednesday, I went over blood pressure waveforms defining pulse rate, systolic pressure, and diastolic pressure, and talking about the frequency ranges of the pulse rate. I then explained to them how the filter program was run (many students still don’t know about the “<” and “>” conventions for standard in and standard out on command lines). I also showed the gnuplot trick that allows using standard out from a program in place of a file in a plot command:
plot '< python bandpass-filter.py < pressure.data' using 1:3 with lines

I did not explain how digital filters worked, but I did say why I chose Bessel filters (to preserve as much of the time-domain structure of the signal as possible).  In response to a question I also explained the effect that choosing 5th order filters had (the rolloff as f5 or f–5, rather than f1 or f–1 as with a first-order RC filter). I also explained that the computation required more and more precise numbers as the order got higher, and that 5th-order was a good tradeoff between needed precision and fast rolloff.

One thing that I didn’t get to was explaining that “filtfilt” does the filtering twice: once with time going forward and again with time running backwards. The time reversal cancels a lot of the distortion in the time domain (so the choice of Bessel filters is not crucial), but doing two passes also doubles the order of the filter, so that the rolloff is really f10 or f–10.

I did remember to tell students that they needed to have the scipy package installed in order to run the filter program, and that if their python was installed from python.org that they could probably just run “pip install scipy”. At least one student in the class is using the Anaconda installation of python, which already has scipy installed.

At the end of the lecture I had only 10 minutes left, so I did not get into the internals of instrumentation amplifiers (needed for the EKG lab at the end of the quarter) nor transimpedance amplifiers (needed for next week’s lab). Instead I covered the voltmeter impedance measurements I made last week, explaining how I did the measurements, how I did the fitting, and what the results were.  In particular, I mentioned that swapping the sets of leads changed the behavior, so the extra capacitance (beyond the 100pF of the meter itself) appears to be coming from the leads.  I sent the data files and gnuplot script to them via e-mail, after one student requested them.

2014 August 19

Two papers submitted

I got two papers with a co-author finished and submitted today.

One of them is based on the segmenter algorithms that I wrote up in this blog (with some improved parameterization): Segmenting noisy signals from nanopores, Segmenting filtered signals, and Segmenting noisy signals revisited.

The other one is an HMM paper that uses the automatic segmentation and an HMM that models the DNA motor behavior to do the same analysis that was done by hand in an earlier paper “Error Rates for Nanopore Discrimination Among Cytosine, Methylcytosine, and Hydroxymethylcytosine Along Individual DNA Strands” in PNAS.  The automatic HMM method was more accurate than the hand curation and feature extraction followed by sophisticated machine learning methods like support vector machines and random forests—I think that we were limited before by the biases of the hand curation and feature extraction.

Unfortunately, we were not able to do side-by-side comparison on exactly the same data, as the original data for the PNAS paper was lost in a 2-drive failure on a RAID that had not been properly backed up.

The paper writing did not take much of my time this summer, as my co-author did a lot of the writing and kept me from procrastinating too much.

I’ve also been working on my book, tentatively titled Applied Electronics for Bioengineers, but without a co-author to keep me honest, I’ve been doing a lot of procrastinating.  Yesterday I had planned to show the waveform of the gate signal on a nFET being used for pulse-width modulation, to show the “flat spot” that results from feedback through the gate-drain capacitance.  I fussed with that for quite a while, but never got a recording I liked the looks of—so ended up getting nothing written on the book.  I’ll probably try again tomorrow, but with a backup plan of working on a different section if that one gets too frustrating, or of writing some of the text that will go around the figure, once I get a figure I can use.

The book is currently 215 pages, but that includes the front matter and back matter.  The core of the book is only about 188 pages, with 74 figures.  I probably need a bunch more figures (schematics, graphs, photos, …), but those are what take the longest to produce.  There’s a little over a month left until classes start, but I don’t see myself finishing by then.

Next Page »