Raspberry Pi RF Frequency Counter

I build a lot of RF circuits, and often it’s convenient to measure and log frequency with a computer. Previously I’ve built standalone frequency counters, frequency counters with a PC interface, and even hacked a classic frequency counter to add USB interface (twice, actually). My latest device uses only 2 microchips to provide a Raspberry Pi with RF frequency measurement capabilities. The RF signal clocks a 32-bit counter SN74LV8154 ($1.04 on Mouser) connected to a 16-bit IO expander MCP23017 ($1.26 on Mouser) accessable to the Raspberry Pi (via I²C) to provide real-time frequency measurements from a python script for $2.30 in components! Well, plus the cost of the Raspberry Pi. All files for this project are on my GitHub page.

img_8773

The entire circuit is only two microchips! I have a few passives to clean up the RF signal (the RF input is loaded with a 1k resistor to ground, decoupled through a series 100 nF capacitor, and balanced at VCC/2 through a voltage divider of two 47k resistors), but if the measured signal is already a strong square wave they could be omitted. The circuit requires a gate pulse which typically will be 1 pulse per second (1PPS) and can be generated by dividing-down a 32.768kHz oscillator, a spare pin on a microcontroller, a fancy 1PPS time reference, or like in my case a GPS module (Neo-6M) with 1PPS output to provide an extremely accurate gate.

schem

The connections are intuitive! The I2C address is 0x20 when A0, A1, and A2 are grounded. GPB(1-4) control the register select of the counter, and GPA(0-7) reads each bit of the selected register. The whole thing is controlled from Python, but could be trivially written in any language.

img_8777

Here’s a quick summary describing how the code works: First I send bytes to address 0 and 1 to set all pins of GPIO A as inputs, and GPIO B as outputs. Note that only 4 of 8 pins are used for the output, so technically 4 extra pins could be used for things like blinking LEDs or controlling other devices. I then set the register select pins by sending a value to 0x13 (GPIO B), and read the entire GPIO A bus (INTCAPB, 0x18). For address details, consult the datasheet. I do this 4 times (1 for each byte of the 32-bit counter), do a little math to turn it into a frequency value, and compare the current value with the last value and take the difference to display as the measured frequency.

screenshot

An advantage of this continuously running mode is that no clock cycles are lost, so a gate which accidentally fires a bit early due to jitter and cuts-off a cycle will compensate for it on a subsequent read. This is shown above, as a very stable 10MHz frequency reference is measured with this method. A “slow” 1PPS clock tick causes a reading slightly higher, compensated-for by the next reading being slightly lower. In this way, clock sources which are extremely accurate but suffer from low precision (like GPS time sources) are able to maximize the long-term measurement of frequency. Combining this frequency measurement technique with the ability to generate an analog voltage with a Raspberry Pi will allow me to perform some interesting experiments with a voltage controlled crystal oscillator.

Useful Links:

 


     

Generating Analog Voltage with Raspberry Pi

I recently had the need to generate analog voltages from the Raspberry PI, which has rich GPIO digital outputs but no analog outputs. I looked into the RPi.GPIO project which can create PWM (which I wanted to smooth using a low pass filter to create the analog voltage), but its output on the oscilloscope looked terrible! It stuttered all over the place, likely because the duty is continuously under software control. I ended up solving my problem with a MCP4921 12-bit DAC chip (about $1.50 on eBay). It’s controlled via SPI, and although I could have written a python program to bit-bang its protocol with RPi.GPIO I realized I could write directly to the Raspberry Pi SPI device using the echo command. Dividing 3.3V into 12-bits (4096) means that I can control voltage in steps of less than 1mV each, right from the bash console!

img_8696

Video: The Problem (RPi PWM jitters)

Video: My Solution (SPI DAC)

Hardware Connection

There’s very little magic in how the microchip is connected to the Pi. It’s a straight shot to its SPI bus! Here’s a quick drawing showing which pins to connect. Check your device against the Raspberry Pi GPIO pinout diagram for different devices.

img_8701-1

Controlling the DAC with a Bus Pirate

Before I used a Raspberry Pi to control the DAC chip, I tested it out with a Bus Pirate. I don’t have a lot of pictures of the project, but I have a screenshot of a serial console used to send commands to the chip. One advantage of the Bus Pirate is that I can type bytes in binary, which helps to see the individual bits. I don’t have this ability when I’m working in the bash console.

serial

I’m less familiar with the Bus Pirate, but this was a good opportunity to get to know it a little better. It look me a long time (requiring I pull out the logic analyzer) to realize that I had to manually enable/disable the chip-select line, using the “[” and “]” commands. When I set up the SPI mode (command m5) I told it to use active low, but I wasn’t sure how to reverse the active level of the chip-select commands, so I just did ]this[ instead of [this] and it worked great.

frompi

This is the signal probed when it was controlled by the Raspberry Pi, but it looked essentially identical when values were sent via the Bus Pirate. The only difference is there was an appreciable delay between the “]” commands and each of the bytes. It worked fine though.

Controlling the DAC with Console Commands

Once the hardware was configured, the software was trivial. I could control analog voltages by sending two properly-formatted bytes to the SPI hardware device. Importantly, you must use raspi-config to enable SPI.

# set analog voltage to minimum value (about 0V)
echo -ne "\x30\x00" > /dev/spidev0.0 # minimum

# set analog voltage to something a little higher
echo -ne "\x30\xAB" > /dev/spidev0.0 

# set analog voltage to maximum value (about 3.3V)
echo -ne "\x3F\xFF" > /dev/spidev0.0

Helpful Links:

 


     

Realtime Audio Visualization in Python

Python’s “batteries included” nature makes it easy to interact with just about anything… except speakers and a microphone! As of this moment, there still are not standard libraries which which allow cross-platform interfacing with audio devices. There are some pretty convenient third-party modules, but I hope in the future a standard solution will be distributed with python. I appreciate the differences of Linux architectures such as ALSA and OSS, but toss in Windows and MacOS in the mix and it gets to be a huge mess. For Linux, would I even need anything fancy? I can run “cat file.wav > /dev/dsp” from a command prompt to play audio. There are some standard libraries for operating system specific sound (i.e., winsound), but I want something more versatile. The official audio wiki page on the subject lists a small collection of third-party platform-independent libraries. After excluding those which don’t support microphone access (the ultimate goal of all my poking around in this subject), I dove a little deeper into sounddevice and PyAudio. Both of these I installed with pip (i.e., pip install pyaudio)

For a more modern, cleaner, and more complete GUI-based viewer of realtime audio data (and the FFT frequency data), check out my Python Real-time Audio Frequency Monitor project.

I really like the structure and documentation of sounddevice, but I decided to keep developing with PyAudio for now. Sounddevice seemed to take more system resources than PyAudio (in my limited test conditions: Windows 10 with very fast and modern hardware, Python 3), and would audibly “glitch” music as it was being played every time it attached or detached from the microphone stream. I tried streaming, but after about an hour I couldn’t get clean live access to the microphone without glitching audio playback. Furthermore, every few times I ran this script it crashed my python kernel! I very rarely see this happening. iPython complained: “It seems the kernel died unexpectedly. Use ‘Restart kernel’ to continue using this console” and I eventually moved back to PyAudio. For a less “realtime” application, sounddevice might be a great solution. Here’s the minimal case sounddevice script I tested with (that crashed sometimes). If you have a better one to do live high-speed audio capture, let me know!

import sounddevice #pip install sounddevice

for i in range(30): #30 updates in 1 second
    rec = sounddevice.rec(44100/30)
    sounddevice.wait()
    print(rec.shape)

Here’s a simple demo to show how I get realtime microphone audio into numpy arrays using PyAudio. This isn’t really that special. It’s a good starting point though. Note that rather than have the user define a microphone source in the python script (I had a fancy menu system handling this for a while), I allow PyAudio to just look at the operating system’s default input device. This seems like a realistic expectation, and saves time as long as you don’t expect your user to be recording from two different devices at the same time. This script gets some audio from the microphone and shows the values in the console (ten times).

import pyaudio
import numpy as np

CHUNK = 4096 # number of data points to read at a time
RATE = 44100 # time resolution of the recording device (Hz)

p=pyaudio.PyAudio() # start the PyAudio class
stream=p.open(format=pyaudio.paInt16,channels=1,rate=RATE,input=True,
              frames_per_buffer=CHUNK) #uses default input device

# create a numpy array holding a single read of audio data
for i in range(10): #to it a few times just to see
    data = np.fromstring(stream.read(CHUNK),dtype=np.int16)
    print(data)

# close the stream gracefully
stream.stop_stream()
stream.close()
p.terminate()

01

I tried to push the limit a little bit and see how much useful data I could get from this console window. It turns out that it’s pretty responsive! Here’s a slight modification of the code, made to turn the console window into an impromptu VU meter.

import pyaudio
import numpy as np

CHUNK = 2**11
RATE = 44100

p=pyaudio.PyAudio()
stream=p.open(format=pyaudio.paInt16,channels=1,rate=RATE,input=True,
              frames_per_buffer=CHUNK)

for i in range(int(10*44100/1024)): #go for a few seconds
    data = np.fromstring(stream.read(CHUNK),dtype=np.int16)
    peak=np.average(np.abs(data))*2
    bars="#"*int(50*peak/2**16)
    print("%04d %05d %s"%(i,peak,bars))

stream.stop_stream()
stream.close()
p.terminate()

The results are pretty good! The advantage here is that no libraries are required except PyAudio. For people interested in doing simple math (peak detection, frequency detection, etc.) this is a perfect starting point. Here’s a quick cellphone video:

I’ve made realtime audio visualization (realtime FFT) scripts with Python before, but 80% of that code was creating a GUI. I want to see data in real time while I’m developing this code, but I really don’t want to mess with GUI programming. I then had a crazy idea. Everyone has a web browser, which is a pretty good GUI… with a Python script to analyze audio and save graphs (a lot of them, quickly) and some JavaScript running in a browser to keep refreshing those graphs, I could get an idea of what the audio stream is doing in something kind of like real time. It was intended to be a hack, but I never expected it to work so well! Check this out…

Here’s the python script to listen to the microphone and generate graphs:

import pyaudio
import numpy as np
import pylab
import time

RATE = 44100
CHUNK = int(RATE/20) # RATE / number of updates per second

def soundplot(stream):
    t1=time.time()
    data = np.fromstring(stream.read(CHUNK),dtype=np.int16)
    pylab.plot(data)
    pylab.title(i)
    pylab.grid()
    pylab.axis([0,len(data),-2**16/2,2**16/2])
    pylab.savefig("03.png",dpi=50)
    pylab.close('all')
    print("took %.02f ms"%((time.time()-t1)*1000))

if __name__=="__main__":
    p=pyaudio.PyAudio()
    stream=p.open(format=pyaudio.paInt16,channels=1,rate=RATE,input=True,
                  frames_per_buffer=CHUNK)
    for i in range(int(20*RATE/CHUNK)): #do this for 10 seconds
        soundplot(stream)
    stream.stop_stream()
    stream.close()
    p.terminate()

Here’s the HTML file with JavaScript to keep reloading the image… 

<html>
<script language="javascript">
function RefreshImage(){
document.pic0.src="03.png?a=" + String(Math.random()*99999999);
setTimeout('RefreshImage()',50);
}
</script>
<body onload="RefreshImage()">
<img name="pic0" src="03.png">
</body>
</html>

Here’s the result! I couldn’t believe my eyes. It’s not elegant, but it’s kind of functional!

Why stop there? I went ahead and wrote a microphone listening and processing class which makes this stuff easier. My ultimate goal hasn’t been revealed yet, but I’m sure it’ll be clear in a few weeks. Let’s just say there’s a lot of use in me visualizing streams of continuous data. Anyway, this class is the truly terrible attempt at a word pun by merging the words “SWH”, “ear”, and “Hear”, into the official title “SWHear” which seems to be unique on Google. This class is minimal case, but can be easily modified to implement threaded recording (which won’t cause the rest of the functions to hang) as well as mathematical manipulation of data, such as FFT. With the same HTML file as used above, here’s the new python script and some video of the output:

import pyaudio
import time
import pylab
import numpy as np

class SWHear(object):
    """
    The SWHear class is made to provide access to continuously recorded
    (and mathematically processed) microphone data.
    """

    def __init__(self,device=None,startStreaming=True):
        """fire up the SWHear class."""
        print(" -- initializing SWHear")

        self.chunk = 4096 # number of data points to read at a time
        self.rate = 44100 # time resolution of the recording device (Hz)

        # for tape recording (continuous "tape" of recent audio)
        self.tapeLength=2 #seconds
        self.tape=np.empty(self.rate*self.tapeLength)*np.nan

        self.p=pyaudio.PyAudio() # start the PyAudio class
        if startStreaming:
            self.stream_start()

    ### LOWEST LEVEL AUDIO ACCESS
    # pure access to microphone and stream operations
    # keep math, plotting, FFT, etc out of here.

    def stream_read(self):
        """return values for a single chunk"""
        data = np.fromstring(self.stream.read(self.chunk),dtype=np.int16)
        #print(data)
        return data

    def stream_start(self):
        """connect to the audio device and start a stream"""
        print(" -- stream started")
        self.stream=self.p.open(format=pyaudio.paInt16,channels=1,
                                rate=self.rate,input=True,
                                frames_per_buffer=self.chunk)

    def stream_stop(self):
        """close the stream but keep the PyAudio instance alive."""
        if 'stream' in locals():
            self.stream.stop_stream()
            self.stream.close()
        print(" -- stream CLOSED")

    def close(self):
        """gently detach from things."""
        self.stream_stop()
        self.p.terminate()

    ### TAPE METHODS
    # tape is like a circular magnetic ribbon of tape that's continously
    # recorded and recorded over in a loop. self.tape contains this data.
    # the newest data is always at the end. Don't modify data on the type,
    # but rather do math on it (like FFT) as you read from it.

    def tape_add(self):
        """add a single chunk to the tape."""
        self.tape[:-self.chunk]=self.tape[self.chunk:]
        self.tape[-self.chunk:]=self.stream_read()

    def tape_flush(self):
        """completely fill tape with new data."""
        readsInTape=int(self.rate*self.tapeLength/self.chunk)
        print(" -- flushing %d s tape with %dx%.2f ms reads"%\
                  (self.tapeLength,readsInTape,self.chunk/self.rate))
        for i in range(readsInTape):
            self.tape_add()

    def tape_forever(self,plotSec=.25):
        t1=0
        try:
            while True:
                self.tape_add()
                if (time.time()-t1)>plotSec:
                    t1=time.time()
                    self.tape_plot()
        except:
            print(" ~~ exception (keyboard?)")
            return

    def tape_plot(self,saveAs="03.png"):
        """plot what's in the tape."""
        pylab.plot(np.arange(len(self.tape))/self.rate,self.tape)
        pylab.axis([0,self.tapeLength,-2**16/2,2**16/2])
        if saveAs:
            t1=time.time()
            pylab.savefig(saveAs,dpi=50)
            print("plotting saving took %.02f ms"%((time.time()-t1)*1000))
        else:
            pylab.show()
            print() #good for IPython
        pylab.close('all')

if __name__=="__main__":
    ear=SWHear()
    ear.tape_forever()
    ear.close()
    print("DONE")

I don’t really intend anyone to actually do this, but it’s a cool alternative to recording a small portion of audio, plotting it in a pop-up matplotlib window, and waiting for the user to close it to record a new fraction. I had a lot more text in here demonstrating real-time FFT, but I’d rather consolidate everything FFT related into a single post. For now, I’m happy pursuing microphone-related python projects with PyAudio.

UPDATE: Displaying a single frequency

Use Numpy’s FFT() and FFTFREQ() to turn the linear data into frequency. Set that target and grab the FFT value corresponding to that frequency. I haven’t tested this to be sure it’s working, but it should at least be close…

import pyaudio
import numpy as np
np.set_printoptions(suppress=True) # don't use scientific notation

CHUNK = 4096 # number of data points to read at a time
RATE = 44100 # time resolution of the recording device (Hz)
TARGET = 2100 # show only this one frequency

p=pyaudio.PyAudio() # start the PyAudio class
stream=p.open(format=pyaudio.paInt16,channels=1,rate=RATE,input=True,
              frames_per_buffer=CHUNK) #uses default input device

# create a numpy array holding a single read of audio data
for i in range(10): #to it a few times just to see
    data = np.fromstring(stream.read(CHUNK),dtype=np.int16)
    fft = abs(np.fft.fft(data).real)
    fft = fft[:int(len(fft)/2)] # keep only first half
    freq = np.fft.fftfreq(CHUNK,1/RATE)
    freq = freq[:int(len(freq)/2)] # keep only first half
    assert freq[-1]>TARGET, "ERROR: increase chunk size"
    val = fft[np.where(freq>TARGET)[0][0]]
    print(val)

# close the stream gracefully
stream.stop_stream()
stream.close()
p.terminate()

UPDATE: Display peak frequency

If your goal is to determine which frequency is producing the loudest tone, use this function. I also added a few lines to graph the output in case you want to observe how it operates. I recommend testing this script with a tone generator, or a YouTube video containing tones of a range of frequencies like this one.

import pyaudio
import numpy as np
import matplotlib.pyplot as plt

np.set_printoptions(suppress=True) # don't use scientific notation

CHUNK = 4096 # number of data points to read at a time
RATE = 44100 # time resolution of the recording device (Hz)

p=pyaudio.PyAudio() # start the PyAudio class
stream=p.open(format=pyaudio.paInt16,channels=1,rate=RATE,input=True,
              frames_per_buffer=CHUNK) #uses default input device

# create a numpy array holding a single read of audio data
for i in range(10): #to it a few times just to see
    data = np.fromstring(stream.read(CHUNK),dtype=np.int16)
    data = data * np.hanning(len(data)) # smooth the FFT by windowing data
    fft = abs(np.fft.fft(data).real)
    fft = fft[:int(len(fft)/2)] # keep only first half
    freq = np.fft.fftfreq(CHUNK,1/RATE)
    freq = freq[:int(len(freq)/2)] # keep only first half
    freqPeak = freq[np.where(fft==np.max(fft))[0][0]]+1
    print("peak frequency: %d Hz"%freqPeak)

    # uncomment this if you want to see what the freq vs FFT looks like
    #plt.plot(freq,fft)
    #plt.axis([0,4000,None,None])
    #plt.show()
    #plt.close()

# close the stream gracefully
stream.stop_stream()
stream.close()
p.terminate()

 


     

Epoch Timestamp Hashing

I was recently presented with the need to rename a folder of images based on a timestamp. This way, I can keep saving new files in that folder with overlapping filenames (i.e., 01.jpg, 02.jpg, 03.jpg, etc.), and every time I run this script all images are prepended with a timestamp. I still want the files to be sorted alphabetically, which is why an alphabetical timestamp (rather than a random hash) is preferred.

  • At first I considered a long date such as 2014-04-19-01.jpg, but that adds so much text!
    …also, it doesn’t include time of day.
  • If I include time of day, it becomes 2014-04-19-09-16-23-01.jpg
  • If I eliminate dashes to shorten it, it becomes hard to read, but might work 140419091623-01.jpg
  • If I use Unix Epoch time, it becomes 1397912944-01.jpg

The result I came up with uses base conversion and a string table of numbers and letters (in alphabetical order) to create a second-respecting timestamp hash using an arbitrary number of characters. For simplicity, I used 36 characters: 0-9, and a-z. I then wrote two functions to perform arbitrary base conversion, pulling characters from the hash. Although I could have nearly doubled my available characters by including the full ASCII table, respecting capitalization, I decided to keep it simple. The scheme goes like this:

  • Determine the date / time: 19-Apr-2014 13:08:55
  • Create an integer of Unix Epoch time (seconds past Jan 1, 1970):  1397912935
  • Do a base conversion from a character list: n4a4iv
  • My file name now becomes n4a4iv-01.jpg – I can accept this!
    and when I sort the folder alphabetically, they’re in order by the timestamp

I can now represent any modern time, down to the second, with 6 characters. Here’s some example output:

19-Apr-2014 13:08:55 <-> 1397912935 <-> n4a4iv
19-Apr-2014 13:08:56 <-> 1397912936 <-> n4a4iw
19-Apr-2014 13:08:57 <-> 1397912937 <-> n4a4ix
19-Apr-2014 13:08:58 <-> 1397912938 <-> n4a4iy
19-Apr-2014 13:08:59 <-> 1397912939 <-> n4a4iz
19-Apr-2014 13:09:00 <-> 1397912940 <-> n4a4j0
19-Apr-2014 13:09:01 <-> 1397912941 <-> n4a4j1
19-Apr-2014 13:09:02 <-> 1397912942 <-> n4a4j2
19-Apr-2014 13:09:03 <-> 1397912943 <-> n4a4j3
19-Apr-2014 13:09:04 <-> 1397912944 <-> n4a4j4

Interestingly, if I change my hash characters away from the list of 36 alphanumerics and replace it with just 0 and 1, I can encode/decode the date in binary:

19-Apr-2014 13:27:28 <-> 1397914048 <-> 1010011010100100111100111000000
19-Apr-2014 13:27:29 <-> 1397914049 <-> 1010011010100100111100111000001
19-Apr-2014 13:27:30 <-> 1397914050 <-> 1010011010100100111100111000010
19-Apr-2014 13:27:31 <-> 1397914051 <-> 1010011010100100111100111000011
19-Apr-2014 13:27:32 <-> 1397914052 <-> 1010011010100100111100111000100
19-Apr-2014 13:27:33 <-> 1397914053 <-> 1010011010100100111100111000101
19-Apr-2014 13:27:34 <-> 1397914054 <-> 1010011010100100111100111000110
19-Apr-2014 13:27:35 <-> 1397914055 <-> 1010011010100100111100111000111
19-Apr-2014 13:27:36 <-> 1397914056 <-> 1010011010100100111100111001000
19-Apr-2014 13:27:37 <-> 1397914057 <-> 1010011010100100111100111001001

Here’s the code to generate / decode Unix epoch timestamps in Python:

hashchars='0123456789abcdefghijklmnopqrstuvwxyz'
#hashchars='01' #for binary

def epochToHash(n):
  hash=''
  while n>0:
    hash = hashchars[int(n % len(hashchars))] + hash
    n = int(n / len(hashchars))
  return hash

def epochFromHash(s):
  s=s[::-1]
  epoch=0
  for pos in range(len(s)):
    epoch+=hashchars.find(s[pos])*(len(hashchars)**pos)
  return epoch

import time
t=int(time.time())
for i in range(10):
  t=t+1
  print(time.strftime("%d-%b-%Y %H:%M:%S", time.gmtime(t)),
              "<->", t,"<->",epochToHash(t))

     

Simple DIY ECG + Pulse Oximeter (version 2)

UPDATE: An improved ECG design was posted in August, 2016.
Check out: https://www.swharden.com/wp/2016-08-08-diy-ecg-with-1-op-amp/

Of the hundreds of projects I’ve shared over the years, none has attracted more attention than my DIY ECG machine on the cheap posted almost 4 years ago. This weekend I re-visited the project and made something I’m excited to share!  The original project was immensely popular, my first featured article on Hack-A-Day, and today “ECG” still represents the second most searched term by people who land on my site. My gmail account also has had 194 incoming emails from people asking details about the project. A lot of it was by frustrated students trying to recreate the project running into trouble because it was somewhat poorly documented. Clearly, it’s a project that a wide range of people are interested in, and I’m happy to revisit it bringing new knowledge and insight to the project. I will do my best to document it thoroughly so anyone can recreate it!

The goal of this project is to collect heartbeat information on a computer with minimal cost and minimal complexity.  I accomplished this with fewer than a dozen components (all of which can be purchased at RadioShack). It serves both as a light-based heartbeat monitor (similar to a pulse oximeter, though it’s not designed to quantitatively measure blood oxygen saturation), and an electrocardiogram (ECG) to visualize electrical activity generated by heart while it contracts. Let’s jump right to the good part – this is what comes out of the machine:

That’s my actual heartbeat. Cool, right? Before I go into how the circuit works, let’s touch on how we measure heartbeat with ECG vs. light (like a pulse oximeter).  To form a heartbeat, the pacemaker region of the heart (called the SA node, which is near the upper right of the heart) begins to fire and the atria (the two top chambers of the heart) contract. The SA node generates a little electrical shock which stimulated a synchronized contraction. This is exactly what defibrillators do when a heart has stopped beating. When a heart attack is occurring and a patient is undergoing ventricular fibrillation, it means that heart muscle cells are contracting randomly and not in unison, so the heart quivers instead of pumping as an organ. Defibrillators synchronize the heart beat with a sudden rush of current over the heart to reset all of the cells to begin firing at the same time (thanks Ron for requesting a more technical description).  If a current is run over the muscle, the cells (cardiomyocytes) all contract at the same time, and blood moves. The AV node (closer to the center of the heart) in combination with a slow conducting pathway (called the bundle of His) control contraction of the ventricles (the really large chambers at the bottom of the heart), which produce the really large spikes we see on an ECG.  To measure ECG, optimally we’d place electrodes on the surface of the heart. Since that would be painful, we do the best we can by measuring voltage changes (often in the mV range) on the surface of the skin. If we amplify it enough, we can visualize it. Depending on where the pads are placed, we can see different regions of the heart contract by their unique electrophysiological signature. ECG requires sticky pads on your chest and is extremely sensitive to small fluctuations in voltage. Alternatively, a pulse oximeter measures blood oxygenation and can monitor heartbeat by clipping onto a finger tip. It does this by shining light through your finger and measuring how much light is absorbed. This goes up and down as blood is pumped through your finger. If you look at the relationship between absorbency in the red vs. infrared wavelengths, you can infer the oxygenation state of the blood. I’m not doing that today because I’m mostly interested in detecting heart beats.

For operation as a pulse oximeter-type optical heartbeat detector (a photoplethysmograph which produces a photoplethysmogram), I use a bright red LED to shine light through my finger and be detected by a phototransistor (bottom left of the diagram). I talk about how this works in more detail in a previous post. Basically the phototransistor acts like a variable resistor which conducts different amounts of current depending on how much light it sees. This changes the voltage above it in a way that changes with heartbeats. If this small signal is used as the input, this device acts like a pulse oximeter.

For operation as an electrocardiograph (ECG), I attach the (in) directly to a lead on my chest. One of them is grounded (it doesn’t matter which for this circuit – if they’re switched the ECG just looks upside down), and the other is recording. In my original article, I used pennies with wires soldered to them taped to my chest as leads. Today, I’m using fancier sticky pads which are a little more conductive. In either case, one lead goes in the center of your chest, and the other goes to your left side under your arm pit. I like these sticky pads because they stick to my skin better than pennies taped on with electrical tape. I got 100 Nikomed Nikotabs EKG Electrodes 0315 on eBay for $5.51 with free shipping (score!). Just gator clip to them and you’re good to go!

In both cases, I need to build a device to amplify small signals. This is accomplished with the following circuit. The core of the circuit is an LM324 quad operational amplifier.  These chips are everywhere, and extremely cheap. It looks like Thai Shine sells 10 for $2.86 (with free shipping). That’s about a quarter each. Nice!  A lot of ECG projects use instrumentation amplifiers like the AD620 (which I have used with fantastic results), but these are expensive (about $5.00 each). The main difference is that instrumentation amplifiers amplify the difference between two points (which reduces noise and probably makes for a better ECG machine), but for today an operational amplifier will do a good enough job amplifying a small signal with respect to ground. I get around the noise issue by some simple filtering techniques. Let’s take a look at the circuit.

This project utilizes one of the op-amps as a virtual ground. One complaint of using op-amps in simple projects is that they often need + and – voltages. Yeah, this could be done with two 9V batteries to generate +9V and -9V, but I think it’s easier to use a single power source (+ and GND). A way to get around that is to use one of the op-amps as a current source and feed it half of the power supply voltage (VCC), and use the output as a virtual ground (allowing VCC to be your + and 0V GND to be your -). For a good description of how to do this intelligently, read the single supply op amps web page. The caveat is that your signals should remain around VCC/2, which can be done if it is decoupled by feeding it through a series capacitor. The project works at 12V or 5V, but was designed for (and has much better output) at 12V. The remaining 3 op-amps of the LM324 serve three unique functions:

STAGE 1: High gain amplifier. The input signals from either the ECG or pulse oximeter are fed into a chain of 3 opamp stages. The first is a preamplifier. The output is decoupled through a series capacitor to place it near VCC/2, and amplified greatly thanks to the 1.8Mohm negative feedback resistor. Changing this value changes initial gain.

STAGE 2: active low-pass filter. The 10kOhm variable resistor lets you adjust the frequency cutoff. The opamp serves as a unity gain current source / voltage follower that has high input impedance when measuring the output f the low-pass filter and reproduces its voltage with a low impedance output. There’s some more information about active filtering on this page. It’s best to look at the output of this stage and adjust the potentiometer until the 60Hz noise (caused by the AC wiring in the walls) is most reduced while the lower-frequency component of your heartbeat is retained. With the oximeter, virtually no noise gets through. Because the ECG signal is much smaller, this filter has to be less aggressive, and this noise is filtered-out by software (more on this later).

STAGE 3: final amplifier with low-pass filter. It has a gain of ~20 (determined by the ratio of the 1.8kOhm to 100Ohm resistors) and lowpass filtering components are provided by the 22uF capacitor across the negative feedback resistor. If you try to run this circuit at 5V and want more gain (more voltage swing), consider increasing the value of the 1.8kOhm resistor (wit the capacitor removed). Once you have a good gain, add different capacitor values until your signal is left but the noise reduced. For 12V, these values work fine. Let’s see it in action!

Now for the second half – getting it into the computer. The cheapest and easiest way to do this is to simply feed the output into a sound card! A sound card is an analog-to-digital converter (ADC) that everybody has and can sample up to 48 thousand samples a second! (overkill for this application) The first thing you should do is add an output potentiometer to allow you to drop the voltage down if it’s too big for the sound card (in the case of the oximeter) but but also allow full-volume in the case of sensitive measurements (like ECG). Then open-up sound editing software (I like GoldWave for Windows or Audacity for Linux, both of which are free) and record the input. You can do filtering (low-pass filter at 40Hz with a sharp cutoff) to further eliminate any noise that may have sneaked through. Re-sample at 1,000 Hz (1kHz) and save the output as a text file and you’re ready to graph it! Check it out.

Here are the results of some actual data recorded and processed with the method shown in the video. let’s look at the pulse oximeter first.

That looks pretty good, certainly enough for heartbeat detection. There’s obvious room for improvement, but as a proof of concept it’s clearly working. Let’s switch gears and look at the ECG. It’s much more challenging because it’s signal is a couple orders of magnitude smaller than the pulse oximeter, so a lot more noise gets through. Filtering it out offers dramatic improvements!

Here’s the code I used to generate the graphs from the text files that GoldWave saves. It requires Python, Matplotlib (pylab), and Numpy. In my case, I’m using 32-bit 2.6 versions of everything.

# DIY Sound Card ECG/Pulse Oximeter
# by Scott Harden (2013) http://www.SWHarden.com

import pylab
import numpy

f=open("light.txt")
raw=f.readlines()[1:]
f.close()

data = numpy.array(raw,dtype=float)
data = data-min(data) #make all points positive
data = data/max(data)*100.0 #normalize
times = numpy.array(range(len(data)))/1000.0
pylab.figure(figsize=(15,5))
pylab.plot(times,data)
pylab.xlabel("Time Elapsed (seconds)")
pylab.ylabel("Amplitude (% max)")
pylab.title("Pulse Oximeter - filtered")
pylab.subplots_adjust(left=.05,right=.98)
pylab.show()

Future directions involve several projects I hope to work on soon. First, it would be cool to miniaturize everything with surface mount technology (SMT) to bring these things down to the size of a postage stamp. Second, improved finger, toe, or ear clips (or even taped-on sensors) over long duration would provide a pretty interesting way to analyze heart rate variability or modulation in response to stress, sleep apnea, etc. Instead of feeding the signal into a computer, one could send it to a micro-controller for processing. I’ve made some darn-good progress making multi-channel cross-platform USB option for getting physiology data into a computer, but have some work still to do. Alternatively, this data could be graphed on a graphical LCD for an all-in-one little device that doesn’t require a computer. Yep, lots of possible projects can use this as a starting point.

Notes about safety: If you’re worried about electrical shock, or unsure of your ability to make a safe device, don’t attempt to build an ECG machine. For an ECG to work, you have to make good electrical contact with your skin near your heart, and some people feel this is potentially dangerous. Actually, some people like to argue about how dangerous it actually is, as seen on Hack-A-Day comments and my previous post comments. Some people have suggested the danger is negligible and pointed-out that it’s similar to inserting ear-bud headphones into your ears. Others have suggested that it’s dangerous and pointed-out that milliamps can kill a person. Others contest that pulses of current are far more dangerous than a continuous applied current. Realists speculate that virtually no current would be delivered by this circuit if it is wired properly. Rational, cautionary people worried about it reduce risk of accidental current by applying bidirectional diodes at the level of the chest leads, which short any current (above 0.7V) similar to that shown here. Electrically-savvy folks would design an optically decoupled solution. Intelligent folks who abstain from arguing on the internet would probably consult the datasheets regarding ECG input protection. In all cases, don’t attach electrical devices to your body unless you are confident in their safety. As a catch-all, I present the ECG circuit for educational purposes only, and state that it may not be safe and should not be replicated  There, will that cover me in court in case someone tapes wires to their chest and plugs them in the wall socket?

LET ME KNOW WHAT YOU THINK! If you make this, I’m especially interested to see how it came out. Take pictures of your projects and send them my way! If you make improvements, or take this project further, I’d be happy to link to it on this page. I hope this page describes the project well enough that anyone can recreate it, regardless of electronics experience. Finally, I hope that people are inspired by the cool things that can be done with surprisingly simple electronics. Get out there, be creative, and go build something cool!


     

AVR Programming in Linux

It is not difficult to program ATMEL AVR microcontrollers with linux, and I almost exclusively do this because various unofficial (inexpensive) USB AVR programmers are incompatible with modern versions of windows (namely Windows Vista and Windows 7). I am just now setting-up a new computer station in my electronics room (running Ubuntu Linux 12.04), and to make it easy for myself in the future I will document everything I do when I set-up a Linux computer to program microcontrollers.

Install necessary software

sudo apt-get install gcc-avr avr-libc uisp avrdude

Connect the AVR programmer
This should be intuitive for anyone who has programmed AVRs before. Visit the datasheet of your MCU, identify pins for VCC (+), GND (-), MOSI, MISO, SCK, and RESET, then connect them to the appropriate pins of your programmer.

Write a simple program in C
I made a file “main.c” and put the following inside. It’s the simplest-case code necessary to make every pin on PORTD (PD0, PD1, …, PD7) turn on and off repeatedly, sufficient to blink an LED.

#include <avr/io.h>
#include <util/delay.h>

int main (void)
{
 DDRD = 255; // MAKE ALL PORT D PINS OUTPUTS

 while(1) {
  PORTD = 255;_delay_ms(100); // LED ON
  PORTD = 0;  _delay_ms(100); // LED OFF
 }

 return 0;
}

Compile the code (generate a HEX file)

avr-gcc -w -Os -DF_CPU=2000000UL -mmcu=atmega8 -c -o main.o main.c
avr-gcc -w -mmcu=atmega8 main.o -o main
avr-objcopy -O ihex -R .eeprom main main.hex

note that the arguments define CPU speed and chip model – this will need to be customized for your application

Program the HEX firmware onto the AVR

sudo avrdude -F -V -c avrispmkII -p ATmega8 -P usb -U flash:w:main.hex

note that this line us customized based on my connection (-P flag, USB in my case) and programmer type (-c flag, AVR ISP mkII in my case)

When this is run, you will see something like this:


avrdude: AVR device initialized and ready to accept instructions

Reading | ################################################## | 100% 0.01s

avrdude: Device signature = 0x1e9307
avrdude: NOTE: FLASH memory has been specified, an erase cycle will be performed
         To disable this feature, specify the -D option.
avrdude: erasing chip
avrdude: reading input file "main.hex"
avrdude: input file main.hex auto detected as Intel Hex
avrdude: writing flash (94 bytes):

Writing | ################################################## | 100% 0.04s

avrdude: 94 bytes of flash written

Your Program should now be loaded! Sit back and watch your LED blink.

TIP 1: I’m using a clone AVRISP mkII from Fun4DIY.com which is only $11 (shipped!). I glued it to a perf-board with a socket so I can jumper its control pins to any DIP AVR (80 pins or less). I also glued a breadboard for my convenience, and added LEDs (with ground on one of their pins) for easy jumpering to test programs. You can also build your own USB AVR ISP (free schematics and source code) from the USBtinyISP project website.

TIP 2: Make a shell script to run your compiling / flashing commands with a single command. Name them according to architecture. i.e., “build-m8” or “build-t2313”. Make the first line delete all .hex files in the directory, so it stalls before it loads old .hex files if the compile is unsuccessful. Make similar shell scripts to program fuses. i.e., “fuse-m8-IntClock-8mhz”, “fuse-m8-IntClock-1mhz”, “fuse-m8-ExtCrystal”.

TIP 3: Use a nice text editor well suited for programming. I love Geany.


     

Multichannel USB Analog Sensor with ATMega48

Sometimes it’s tempting to re-invent the wheel to make a device function exactly the way you want. I am re-visiting the field of homemade electrophysiology equipment, and although I’ve already published a home made electocardiograph (ECG), I wish to revisit that project and make it much more elegant, while also planning for a pulse oximeter, an electroencephalograph (EEG), and an electrogastrogram (EGG). This project is divided into 3 major components: the low-noise microvoltage amplifier, a digital analog to digital converter with PC connectivity, and software to display and analyze the traces. My first challenge is to create that middle step, a device to read voltage (from 0-5V) and send this data to a computer.

This project demonstrates a simple solution for the frustrating problem of sending data from a microcontroller to a PC with a USB connection. My solution utilizes a USB FTDI serial-to-usb cable, allowing me to simply put header pins on my device which I can plug into providing the microcontroller-computer link. This avoids the need for soldering surface-mount FTDI chips (which gets expensive if you put one in every project). FTDI cables are inexpensive (about $11 shipped on eBay) and I’ve gotten a lot of mileage out of mine and know I will continue to use it for future projects. If you are interested in MCU/PC communication, consider one of these cables as a rapid development prototyping tool. I’m certainly enjoying mine!

It is important to me that my design is minimalistic, inexpensive, and functions natively on Linux and Windows without installing special driver-related software, and can be visualized in real-time using native Python libraries, such that the same code can be executed identically on all operating systems with minimal computer-side configuration. I’d say I succeeded in this effort, and while the project could use some small touches to polish it up, it’s already solid and proven in its usefulness and functionality.

This is my final device. It’s reading voltage on a single pin, sending this data to a computer through a USB connection, and custom software (written entirely in Python, designed to be a cross-platform solution) displays the signal in real time. Although it’s capable of recording and displaying 5 channels at the same time, it’s demonstrated displaying only one. Let’s check-out a video of it in action:

This 5-channel realtime USB analog sensor, coupled with custom cross-platform open-source software, will serve as the foundation for a slew of electrophysiological experiments, but can also be easily expanded to serve as an inexpensive multichannel digital oscilloscope. While more advanced solutions exist, this has the advantage of being minimally complex (consisting of a single microchip), inexpensive, and easy to build.

 To the right is my working environment during the development of this project. You can see electronics, my computer, microchips, and coffee, but an intriguingly odd array of immunological posters in the background. I spent a couple weeks camping-out in a molecular biology laboratory here at UF and got a lot of work done, part of which involved diving into electronics again. At the time this photo was taken, I hadn’t worked much at my home workstation. It’s a cool picture, so I’m holding onto it.

Below is a simplified description of the circuit schematic that is employed in this project. Note that there are 6 ADC (analog to digital converter) inputs on the ATMega48 IC, but for whatever reason I ended-up only hard-coding 5 into the software. Eventually I’ll go back and re-declare this project a 6-channel sensor, but since I don’t have six things to measure at the moment I’m fine keeping it the way it is. RST, SCK, MISO, and MOSI are used to program the microcontroller and do not need to be connected to anything for operation. The max232 was initially used as a level converter to allow the micro-controller to communicate with a PC via the serial port. However, shortly after this project was devised an upgrade was used to allow it to connect via USB. Continue reading for details…

Below you can see the circuit breadboarded. The potentiometer (small blue box) simulated an analog input signal.

The lower board is my AVR programmer, and is connected to RST, SCK, MISO, MOSI, and GND to allow me to write code on my laptop and program the board. It’s a Fun4DIY.com AVR programmer which can be yours for $11 shipped! I’m not affiliated with their company, but I love that little board. It’s a clone of the AVR ISP MK-II.

As you can see, the USB AVR programmer I’m using is supported in Linux. I did all of my development in Ubuntu Linux, writing AVR-GCC (C) code in my favorite Linux code editor Geany, then loaded the code onto the chip with AVRDude.

I found a simple way to add USB functionality in a standard, reproducible way that works without requiring the soldering of a SMT FTDI chip, and avoids custom libraries like V-USB which don’t easily have drivers that are supported by major operating systems (Windows) without special software. I understand that the simplest long-term and commercially-logical solution would be to use that SMT chip, but I didn’t feel like dealing with it. Instead, I added header pins which allow me to snap-on a pre-made FTDI USB cable. They’re a bit expensive ($12 on ebay) but all I need is 1 and I can use it in all my projects since it’s a sinch to connect and disconnect. Beside, it supplies power to the target board! It’s supported in Linux and in Windows with established drivers that are shipped with the operating system. It’s a bit of a shortcut, but I like this solution. It also eliminates the need for the max232 chip, since it can sense the voltages outputted by the microcontroller directly.

The system works by individually reading the 10-bit ADC pins on the microcontroller (providing values from 0-1024 to represent voltage from 0-5V or 0-1.1V depending on how the code is written), converting these values to text, and sending them as a string via the serial protocol. The FTDI cable reads these values and transmits them to the PC through a USB connection, which looks like “COM5” on my Windows computer. Values can be seen in any serial terminal program (i.e., hyperterminal), or accessed through Python with the PySerial module.

As you can see, I’m getting quite good at home-brewn PCBs. While it would be fantastic to design a board and have it made professionally, this is expensive and takes some time. In my case, I only have a few hours here or there to work on projects. If I have time to design a board, I want it made immediately! I can make this start to finish in about an hour. I use a classic toner transfer method with ferric chloride, and a dremel drill press to create the holes. I haven’t attacked single-layer SMT designs yet, but I can see its convenience, and look forward to giving it a shot before too long.

Here’s the final board ready for digitally reporting analog voltages. You can see 3 small headers on the far left and 2 at the top of the chip. These are for RST, SCK, MISO, MOSI, and GND for programming the chip. Once it’s programmed, it doesn’t need to be programmed again. Although I wrote the code for an ATMega48, it works fine on a pin-compatible ATMega8 which is pictured here. The connector at the top is that FTDI USB cable, and it supplies power and USB serial connectivity to the board.

If you look closely, you can see that modified code has been loaded on this board with a Linux laptop. This thing is an exciting little board, because it has so many possibilities. It could read voltages of a single channel in extremely high speed and send that data continuously, or it could read from many channels and send it at any rate, or even cooler would be to add some bidirectional serial communication capabilities to allow the computer to tell the microcontroller which channels to read and how often to report the values back. There is a lot of potential for this little design, and I’m glad I have it working.

Unfortunately I lost the schematics to this device because I formatted the computer that had the Eagle files on it. It should be simple and intuitive enough to be able to design again. The code for the microcontroller and code for the real-time visualization software will be posted below shortly. Below are some videos of this board in use in one form or another:

Here is the code that is loaded onto the microcontroller:

#define F_CPU 8000000UL
#include <avr/io.h>
#include <util/delay.h>

void readADC(char adcn){
		//ADMUX = 0b0100000+adcn; // AVCC ref on ADCn
		ADMUX = 0b1100000+adcn; // AVCC ref on ADCn
		ADCSRA |= (1<<ADSC); // reset value
        while (ADCSRA & (1<<ADSC)) {}; // wait for measurement
}

int main (void){
    DDRD=255;
	init_usart();
    ADCSRA = 0b10000111; //ADC Enable, Manual Trigger, Prescaler
    ADCSRB = 0;

    int adcs[8]={0,0,0,0,0,0,0,0};

    char i=0;
	for(;;){
		for (i=0;i<8;i++){readADC(i);adcs[i]=ADC>>6;}
		for (i=0;i<5;i++){sendNum(adcs[i]);send(44);}
		readADC(0);
		send(10);// LINE BREAK
		send(13); //return
		_delay_ms(3);_delay_ms(5);
	}
}

void sendNum(unsigned int num){
	char theIntAsString[7];
	int i;
	sprintf(theIntAsString, "%u", num);
	for (i=0; i < strlen(theIntAsString); i++){
		send(theIntAsString[i]);
	}
}


void send (unsigned char c){
	while((UCSR0A & (1<<UDRE0)) == 0) {}
	UDR0 = c;
}

void init_usart () {
	// ATMEGA48 SETTINGS
	int BAUD_PRESCALE = 12;
	UBRR0L = BAUD_PRESCALE; // Load lower 8-bits
	UBRR0H = (BAUD_PRESCALE >> 8); // Load upper 8-bits
	UCSR0A = 0;
	UCSR0B = (1<<RXEN0)|(1<<TXEN0); //rx and tx
	UCSR0C = (1<<UCSZ01) | (1<<UCSZ00); //We want 8 data bits
}

Here is the code that runs on the computer, allowing reading and real-time graphing of the serial data. It’s written in Python and has been tested in both Linux and Windows. It requires *NO* non-standard python libraries, making it very easy to distribute. Graphs are drawn (somewhat inefficiently) using lines in TK. Subsequent development went into improving the visualization, and drastic improvements have been made since this code was written, and updated code will be shared shortly. This is functional, so it’s worth sharing.

import Tkinter, random, time
import socket, sys, serial

class App:

	def white(self):
		self.lines=[]
		self.lastpos=0

		self.c.create_rectangle(0, 0, 800, 512, fill="black")
		for y in range(0,512,50):
			self.c.create_line(0, y, 800, y, fill="#333333",dash=(4, 4))
			self.c.create_text(5, y-10, fill="#999999", text=str(y*2), anchor="w")
		for x in range(100,800,100):
			self.c.create_line(x, 0, x, 512, fill="#333333",dash=(4, 4))
			self.c.create_text(x+3, 500-10, fill="#999999", text=str(x/100)+"s", anchor="w")

		self.lineRedraw=self.c.create_line(0, 800, 0, 0, fill="red")

		self.lines1text=self.c.create_text(800-3, 10, fill="#00FF00", text=str("TEST"), anchor="e")
		for x in range(800):
			self.lines.append(self.c.create_line(x, 0, x, 0, fill="#00FF00"))

	def addPoint(self,val):
		self.data[self.xpos]=val
		self.line1avg+=val
		if self.xpos%10==0:
			self.c.itemconfig(self.lines1text,text=str(self.line1avg/10.0))
			self.line1avg=0
		if self.xpos>0:self.c.coords(self.lines[self.xpos],(self.xpos-1,self.lastpos,self.xpos,val))
		if self.xpos<800:self.c.coords(self.lineRedraw,(self.xpos+1,0,self.xpos+1,800))
		self.lastpos=val
		self.xpos+=1
		if self.xpos==800:
			self.xpos=0
			self.totalPoints+=800
			print "FPS:",self.totalPoints/(time.time()-self.timeStart)
		t.update()

	def __init__(self, t):
		self.xpos=0
		self.line1avg=0
		self.data=[0]*800
		self.c = Tkinter.Canvas(t, width=800, height=512)
		self.c.pack()
		self.totalPoints=0
		self.white()
		self.timeStart=time.time()

t = Tkinter.Tk()
a = App(t)

#ser = serial.Serial('COM1', 19200, timeout=1)
ser = serial.Serial('/dev/ttyUSB0', 38400, timeout=1)
sock = socket.socket(socket.AF_INET, socket.SOCK_DGRAM)
sock.setsockopt(socket.SOL_SOCKET, socket.SO_BROADCAST, 1)

while True:
	while True: #try to get a reading
		#print "LISTENING"
		raw=str(ser.readline())
		#print raw
		raw=raw.replace("n","").replace("r","")
		raw=raw.split(",")
		#print raw
		try:
			point=(int(raw[0])-200)*2
			break
		except:
			print "FAIL"
			pass
	point=point/2
	a.addPoint(point)

If you re-create this device of a portion of it, let me know! I’d love to share it on my website. Good luck!


     

Quantifying University Network Frustrations

I’m sitting in class frustrated as could be. The Internet in this room (D3-3 in the dental tower of Shands Hospital at UF) is unbelievably annoying. For some reason, everything runs fine, then functionality drops to unusable levels. Downloading files (i.e., PDFs of lectures) occurs at about 0.5kb/s (wow), and Internet browsing is hopeless. At most, I can connect to IRC and enjoy myself in #electronics, #python, and #linux. I decided to channel my frustration into productivity, and wrote a quick Python script to let me visualize the problem.out

Notice the massive lag spikes around the time class begins. I think it’s caused by the retarded behavior of windows update and anti-virus software updates being downloaded on a gazillion computers all at the same time which are required to connect to the network on Windows machines. Class start times were 8:30am, 9:35am, and 10:40am. Let’s view it on a logarithmic scale:out2

Finally, the code. It’s two scripts. One pings a website (kernel.org) every few seconds and records the ping time to “pings.txt”, and the other graphs the data. Here are the two scripts:

import socket, time, os, sys, re

def getping():
	pingaling = os.popen("ping -q -c2 kernel.org")
	sys.stdout.flush()
	while 1:
		line = pingaling.readline()
		if not line: break
		line=line.split("n")
		for part in line:
			if "rtt" in part:
				part=part.split(" = ")[1]
				part=part.split('/')[1]
				print part+"ms"
				return part

def add2log(stuff):
	f=open("pings.txt",'a')
	f.write(stuff+",")
	f.close()

while 1:
	print "pinging...",
	stuff="[%s,%s]"%(time.time(),getping())
	print stuff
	add2log(stuff)
	time.sleep(1)
import pylab, time, datetime, numpy

def smoothTriangle(data,degree,dropVals=False):
	triangle=numpy.array(range(degree)+[degree]+range(degree)[::-1])+1
	smoothed=[]
	for i in range(degree,len(data)-degree*2):
		point=data[i:i+len(triangle)]*triangle
		smoothed.append(sum(point)/sum(triangle))
	if dropVals:
		print "smoothlen:",len(smoothed)
		return smoothed
	#smoothed=[smoothed[0]]*(degree+degree/2)+smoothed
	#while len(smoothed)<len(data):smoothed.append(smoothed[-1])
	while len(smoothed)<len(data):smoothed=[None]+smoothed+[None]
	if len(smoothed)>len(data):smoothed.pop(-1)
	return smoothed

print "reading"
f=open("pings.txt")
raw=eval("[%s]"%f.read())
f.close()

xs,ys,big=[],[],[]
for item in raw:
	t=datetime.datetime.fromtimestamp(item[0])
	maxping=20000
	if item[1]>maxping or item[1]==None:
		item[1]=maxping
		big.append(t)
	ys.append(float(item[1]))
	xs.append(t)

#print xs
#raw_input("WAIT")

print "plotting"
fig=pylab.figure(figsize=(10,7))
pylab.plot(xs,ys,'k.',alpha=.1)
pylab.plot(xs,ys,'k-',alpha=.1)
pylab.plot(xs,smoothTriangle(ys,15),'b-')
pylab.grid(alpha=.3)
pylab.axis([None,None,None,2000])
#pylab.semilogy()
#pylab.xlabel("time")
pylab.ylabel("latency (ping kernel.org, ms)")
pylab.title("D3-3 Network Responsiveness")
fig.autofmt_xdate()
#pylab.show()
pylab.savefig('out.png')
pylab.semilogy()
pylab.savefig('out2.png')
fig.autofmt_xdate()
print "done"

     

VD Labs makes it big debut

The VD Labs webpage has been published! I hope that the new VD Labs page will be a single location where I can link to descriptions and downloads of useful radio, audio analysis, and QRSS-related software. It will eventually be the home of the next (recoded-from-scratch) version of QRSS VD, but let’s not get too far ahead of ourselves!

>>> VD Labs Webpage

Since I ran out of steam from working so much on QRSS VD, I didn’t think I’d be publishing mush more “useful” software, but this one blind-sighted me. People on the Knights QRSS mailing list were talking about dividing QRSS transmissions into images which line up with the period of the transmitters repeated messages and projecting the images together in an attempt to average-out the noise, and boost the signal. It’s a simple idea, and it’s the basis behind how a lot of poor imaging devices can improve their output clarity by software (MRI anyone?). I was overwhelmed by dental school obligations the last few weeks, and it pained me so much to read what people were doing (or at least trying to do) and having to sit it out. Now that I have a free day (yay for weekends!) I sat down and wrote some code. I introduce VD Labs QRSS Stitcher and QRSS Stacker!
vd labs flyer

Converting Argo captures into continuous images:

example output:

stitched

Doing the same thing, with ultra-narrow images:

File produced:

stacked_narrow

Using QRSS Stacker to project images:

Another example output:

stacked_stitched

Screenshots:

vd labs qrss stacker
vd labs qrss stitcher

DOWNLOAD:

very soon…


     

Epic Failure 1 Year in the Making

My expression is completely flat right now. I simply cannot believe I’m about to say what I’m preparing to say. I spent nearly a year cracking large prime numbers. In short, I took-on a project I called The Flowering N’th Prime Project, where I used my SheevaPlug to generate a list of every [every millionth] prime number. The current “golden standard” is this page where one can look-up the N’th prime up to 1 trillion. My goal was to reach over 1 trillion, which I did just this morning! I was planning on being the only source on the web to allow lookups of prime numbers greater than 1 trillion. flowering_primes

However, when I went to look at the logs, I realized that the software had a small, fatal bug in it. Apparently every time the program restarted (which happened a few times over the months), although it resumed at its most recent prime number, it erased the previous entries. As a result, I have no logs below N=95 billion. In other words, although I reached my target this morning, it’s completely irrelevant since I don’t have all the previous data to prove it. I’m completely beside myself, and have no idea what I’m going to do. I can start from the beginning again, but that would take another YEAR. [sigh]

So here’s the screw-up. Apparently I coded everything correctly on paper, but due to my lack of experience I overlooked the potential for multiple appends to occur simultaneously. I can only assume that’s what screwed it up, but I cannot be confident. Honestly, I still don’t know specifically what the problem is. All in all, it looks good to me. Here is the relevant Python code.

def add2log(c,v):
	f=open(logfile,'a')
	f.write("%d,%dn"%(c,v))
	f.close()

def resumeFromLog():
	f=open('log.txt')
	raw=f.readlines()[-1]
	f.close()
	return eval("["+raw+"]")

For what it’s worth, this is what remains of the log file:

953238,28546251136703
953239,28546282140203
953240,28546313129849
...
1000772,30020181524029
1000773,30020212566353
1000774,30020243594723