Gas station without pumps

2016 October 5

Broken bike seat

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

Yesterday was not a good day for me.

First, I spent most of the day struggling with the homework for the control-theory class I’m sitting in on. The course is dual listed as an undergrad and grad course, with shared lectures but different homework and projects. The undergrad part of the homework was straight-forward, and I finished it Monday night, but the two additional problems for the grad students were tough.  One of them had a simple “engineering” solution that I got quickly by formal manipulation of the formulæ, but I could not justify some of the steps, since they involved a integral that was not finite.  The other problem was not difficult, but involved a rather tedious amount of algebra to linearize the system—the professor had done the linearization in lecture notes,  and we were just supposed to check it for the homework, but he’d made an error in algebra, so I had to redo the whole thing.

Late in the afternoon, I decided to take a break and replace the sump pump that had failed sometime in the past couple of weeks.  Originally I was going to disassemble the pump and see if the problem was repairable (I think that the switch for the float is not turning on reliably, possibly from corroded contacts), but I decided that I could do that later to have a spare pump, meanwhile getting a working sump pump.  (My house is built over a seep where an aquifer comes to the surface, and the water table is about 3 inches below the surface—during wet years, the water table is sometimes right at the surface.)

I put the old pump in my panniers and headed down to hardware store, when my bike seat suddenly failed.  I tried riding for a block with the failed seat and gave up and returned home.  The failure was right at edge of the block that holds the horizontal crossbar at the front of the seat:

Here is a view from the front showing the tubing displaced vertically from where it belongs.

Here is a view from the front showing the tubing displaced vertically from where it belongs.

A closer view shows a very clean break right at the surface of the block that clamps around the tube.

A closer view shows a very clean break right at the surface of the block that clamps around the tube.

I probably should have had some warning about the imminent failure, as the bike has been creaking a bit more than usual when I pedal for the past several months, but I was never able to track down the creaking. I’m not sure I could have seen the crack that was probably propagating, since it was flush with clamp block.

The seat on my Longbikes Vanguard is not a standard, off-the-shelf component, so I’m probably going to have to custom order a new seat from the manufacturer (who no longer make the Vanguard model, so probably has no spare seats built) and wait weeks or months for one to be built.

I got my old upright bike down from the garage wall, inflated the tires, adjusted one of my panniers to fit the different rack, and headed off to the hardware store, carrying the old pump in the pannier. At the hardware store, I could not find a sump pump with the same outlet size as the old one (they all had bigger outlets). I needed to match, in order to hook the sump pump up to the existing plumbing. Luckily, they did have a reducer that would adjust for the difference.

After buying the pump, I went out to my bike and realized that I couldn’t fit both pumps into one pannier—in fact the new boxed sump pump wouldn’t fit into the pannier even by itself. Normally I carry a bungee cord or two for strapping stuff onto my rear rack, but those were left on the other bike. So I had to go back into the hardware store to buy some new bungee cords—not a big deal, but an irritation.

The bike was a bit wobbly on the way home—I’d forgotten how much difference a high center of gravity makes on an upright bike—and the bike has much twitchier steering than my recumbent anyway—but I got home without incident.

On getting home, I immediately attached the pluming to the new sump pump and lowered it into the sump. Let me correct that—I tried to lower it into the sump, but it wouldn’t fit. The pump was a couple of inches wider than the old pump and though the hole at the top was more than wide enough, it narrowed significantly where the bottom of the foundation for the house spread out, and the remaining hole was simply too small for the new pump. This was particularly frustrating for me, as I was meeting my wife downtown for dinner in less than an hour, and I was going to have to walk rather than bike, so I only had about 10 minutes to come up with a fix.

I then remembered something that should have occurred to me much earlier—I had another one of the small sump pumps in a different sump in the back garden. Quickly pulling it out and attaching the plumbing got the main sump working again (though I still need to recheck the plumbing for leaks). And it turned out that the garden sump was wide enough to accept the new pump—problem solved!

I cleaned up, grabbed a backpack so I could do some shopping after dinner, and walked down to the library to meet my wife. After the stresses of the day, I felt the need for comfort food, so we went to Betty’s Noodles, a hole-in-the-wall Chinese restaurant in the bus station. This restaurant has taken over the niche that Little Shanghai used to fill of providing cheap, tasty Chinese fast food (noodles and rice bowls).  I had ma-po tofu over Chow Fun noodles, which went a long way to de-stress me.  Going to Mission Hill Creamery for a plum sorbet cone afterwards helped also.

On the walk home, a couple blocks before we got home, I realized that I had not done my shopping! I decided not to go back downtown, but to do without my chocolate soymilk for a couple of days, until I can go shopping again.

This morning I finished the homework and submitted it. I’m still a bit bothered about the inverse Laplace transform problem that  can be formally solved but that ends up with a function that doesn’t have a Laplace transform, but I’m pretty sure I did what was expected. After turning in the homework, I realized that there was a possible different interpretation of part of the linearization question than what I did, so I queried the professor about what he really meant.  (The homework isn’t due for a week, so if there is a clarification needed, he can get it to the grad students before the homework is due.)

The TA does not grade my homework, since I’m just auditing, but I’m doing the homework using Python instead of Matlab, so I’m sharing it with the TA and professor anyway, so they can see whether it would be worth switching to free tools.

Currently, the scipy.signal package and matplotlib seem as easy to use at Matlab, but there is no equivalent of SIMULINK, which the professor is relying on for students doing simulations.  I can do the simulations in Python, but setting them up is all text-based, and requires thinking explicitly about the state vector, rather than having a GUI that does all the setup for you.

I bicycled up to campus today on my old upright, after adjusting my other pannier to fit the rack.  I had forgotten how uncomfortable an upright bike is.  This evening my neck and shoulders are sore, and I have chafing on the inside of my thigh.  I really hope I can get the recumbent seat replaced quickly, so that I can go back to riding comfortably!  It might even be worth taking the seat to a local frame-builder and finding out whether they could replace the tube, even if only for a temporary fix. (Although most of the bike is chrome-moly steel, the seat appears to be all aluminum tubing.)

2016 September 5

Improved spectral response for dehumidifier noise

Filed under: Data acquisition — gasstationwithoutpumps @ 14:29
Tags: , , , , , ,

In Dehumidifier, I provided spectral analysis of the sound from the 30-pint Whynter RPD-321EW Energy Star Portable Dehumidifier, using my BitScope USB oscilloscope to record the sound and averaging the absolute value of the FFT over 100 short traces. Today I decided to improve the FFT analysis by using the scipy.signal.welch() function to split a longer trace into windows, apply a Dolph-Chebyshev window-shaping function, and average the resulting power spectra. The main improvement here is the use of the window-shaping function.

I used PteroDAQ with a Teensy LC board to record at 15kHz, which was as fast as I could go with 4× averaging communicating with my old MacBook Pro laptop. Based on the Bitscope recordings, though, that should include all the interesting signal, which appears to be limited to about 5.4kHz. I could only record for about 330,000 samples to avoid overrunning the USB buffer and dropping samples, so I also tried 12kHz, which allowed me to run longer. The analysis looked much the same, with no apparent reduction in the noise floor, so I’ll only present the 15kHz results. (Note: when I get a new laptop this fall, I should be able to run at this sampling rate for much longer, which should allow smoothing out more of the noise.)

The dehumidifier looks very much like 1/f noise, with a big spike around 1050Hz. The spikes at 60Hz, 120Hz, 240Hz, and 360Hz may be electrical noise picked up by the sound amplifier, rather than acoustic noise.

The dehumidifier looks very much like 1/f noise, with a big spike around 1050Hz. The spikes at 60Hz, 120Hz, 240Hz, and 360Hz may be electrical noise picked up by the sound amplifier, rather than acoustic noise.

Using Welch’s method and a proper window function results in a much cleaner spectrum than before. I used a 216 FFT window for the analysis, since I had plenty of data, which gives a fairly fine resolution of about 0.23Hz. The 60Hz spike is appearing at 58.36Hz, which is a little troubling, as the PteroDAQ crystal and the power company both have better frequency control than a 2.7% error. The 120Hz, 240Hz, and 360Hz spikes are at most one bin off from where they should be, so I don’t understand the shift at 60Hz.

I checked my code by testing with a recording of a triangle wave, and found that Welch’s method with an 80dB Dolph-Chebyshev window did an excellent job of removing the spectral smear that I was getting by just averaging rectangular FFT windows (my first approach). The noise floor for Welch’s method on those measurements was pretty much flat across the entire spectrum (maybe rising at higher frequencies), so the 1/f plot here is not an artifact of the analysis, but a real measurement of the acoustic noise. Using 60dB window resulted in a higher noise floor, but 100dB and larger did not lower the noise floor, so I stuck with 80dB as a good parameter value. It might be useful to use other values if the FFT window width is changed or the noise in the data is different.

I’ve included the source code for the program, so that others can use it. It does require that SciPy be installed, in order to use the scipy.signal package.

#! /usr/bin/env python3

"""Sun Sep  4 21:40:28 PDT 2016 Kevin Karplus

Reads from stdin a white-space separated file 
whose first column is a timestamp (in seconds) and 
whose subsequent columns are floating-point data values.

May contain multiple traces, separated by blank lines.

Outputs a tab-separated file with two columns:
        the frequencies
        the Fourier transform

"""

from __future__ import print_function, division

import argparse
import sys
from math import sqrt
import numpy as np
try:
    from future_builtins import zip
except ImportError:     # either before Python2.6 or one of the Python 3.*
    try:
        from itertools import izip as zip
    except ImportError: # not an early Python, so must be Python 3.*
        pass

import scipy
import scipy.signal

class Trace:
    """One contiguous trace consisting of two 1D numpy arrays of the same length
        times:  timestamps in seconds
        values: arbitrary floats
      and sampling_freq: the sampling frequency in Hz
    """
    def __init__(self,times,values,sampling_freq=None):
        self.times=np.array(times)
        self.values=np.array(values)
        assert self.times.shape==self.values.shape
        if sampling_freq is None and len(times)>1:
            sampling_freq = (len(times)-1)/(times[-1]-times[0])
        self.sampling_freq=sampling_freq
    def __len__(self):
        return len(self.values)

def read_trace(file, column=2, sampling_freq=None, echo_comments_to=None):
    """Generator that yields traces from tab-separated rows of data.
    Each row contains a timestamp in the first column and any number of subsequent
    columns, one of which is selected to be the values for the trace.
    
    New traces are started whenever the timestamp in the first column decreases.
    
    If sampling_freq is None, the sampling_freq for the
    trace is deduced from the timestamps.
    
    Comments and blank lines are echoed to a file-like object,
    echo_comments_to, if it is not None.
    """
    times=[]    # time values (in seconds) from frst column
    values=[]   # values from the specified column
    
    for line in file:
        line=line.strip()
        if line.startswith('#') or len(line)==0:
            if echo_comments_to is not None:
                print(line, file=echo_comments_to)
            continue
        fields=line.split()
        if len(fields)>=column:
            time=float(fields[0])
            if times and time<times[-1]:
                yield Trace(times=times, values=values, sampling_freq=sampling_freq)
                times=[]
                values=[]
            times.append(time)
            values.append(float(fields[column-1]))
    if times:
        yield Trace(times=times, values=values, sampling_freq=sampling_freq)


def posint(opt):
    """returns int(opt) if that is positive,
            otherwise raises an argparse.ArgumentTypeError
    """
    try:
        x = int(opt)
    except:
        x=None
    if x is None or x<=0:
        raise argparse.ArgumentTypeError("{} is not a positive integer".format(opt))
    return x

def parse_options(argv):
    """Parse the options and return what argparse does:
        a structure whose fields are the possible options
    """
    parser = argparse.ArgumentParser( description= __doc__, formatter_class = argparse.ArgumentDefaultsHelpFormatter )
    parser.add_argument("--sampling_freq","-s", type=float,
            default=None,
            help="""Sampling frequency (in Hz). 
            If none specified, then estimate sampling frequency from the timestamps.
            """)
    parser.add_argument("--window", "-w", type=posint,
            default=(1<<16),
            help="""Max window size (in samples) for FFT.
            Longer traces will be split up into windows of this size
            """)
    parser.add_argument("--column" , "-c", type=posint, default=2,
           help="""What column should be Fourier transformed?""")
    parser.add_argument("--welch" , type=float, default=None,
           help="""Use Welch's method to produce averaged spectrum from windows. 
           Value is dB difference between main lobe and sidelobes for Dolph-Chebyshev window. 
           80 seems reasonable with 16-bit data (60 produces a higher noise floor, without noticeably narrowing the spikes).
           Note, sqrt is used to get amplitude spectrum, rather than power spectrum. """)
    options=parser.parse_args(argv[1:])
    return options


def main(argv):
    """parse the command line arguments and do FFTs
    """
    options = parse_options(argv)
    sampling_freq= options.sampling_freq

    # input file may contain separate time traces.
    # Each new trace stats when the time value is less than the previous time.
    
    for trace in read_trace(sys.stdin, column=options.column, 
                        sampling_freq=options.sampling_freq, 
                        echo_comments_to=sys.stdout):
        print("# using sampling frequency {:.6f} Hz for FFT of trace with {} samples".format(trace.sampling_freq,len(trace)))


        fft_len=64      # length of the fft to do
        while fft_len<len(trace) and fft_len<options.window:
            fft_len *= 2
        
        if options.welch is not None:
            welch_f,welch_spectrum = scipy.signal.welch(trace.values, trace.sampling_freq, nperseg=fft_len, scaling="spectrum", window=("chebwin",options.welch))
            print("")       # separate the traces
            for f,w in zip(welch_f,welch_spectrum):
                print("{0:.7f}\t{1:.6f}".format(f, sqrt(w)*(fft_len//2)))
        else:
            for start in range(0,len(trace),fft_len):
                window=trace.values[start:start+fft_len]

                # remove DC component, so that zero-padding makes sense
                # WARNING: This in-place change could introduce errors if
                # overlapping windows are allowed.
                dc = sum(window)/len(window)
                window -= dc

                fft = scipy.fftpack.fft(window, n=fft_len)
                fft[0] += dc

                print("")       # separate the traces
                print("# Window width {} = {:.3f}s".format(fft_len,fft_len/trace.sampling_freq))
                for i,f in enumerate(fft[0:fft_len//2]):
                    print("{0:.7f}\t{1:.6f}".format(i*trace.sampling_freq/fft_len, f))


if __name__ == "__main__" :
    sys.exit(main(sys.argv))

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.

%d bloggers like this: