Wow, I can’t believe I took-on such a massive challenge this week! Ironically I’ve worked harder and learned more in the last 4 days than I have over the last 9 months of dental school. This is an incredible feeling of accomplishment. My program, coded from scratch, polls the sound card continuously and makes extremely large spectrographs. It’s not finished, but it’s working. I’m impressed!

gotit

The image above is the actual-size display of the image below. The image below was actually cropped to less than 3,000 pixels high, whereas the original is over 8,000 pixels high!

qrss_big

Note from Future Scott (ten years later, August 2019)

I remember this day! It was the first time I completed a big software project, and it was a great feeling. It motivated me to continue working on software for years to come. I also got a kick out of this quote from early-dental-student Scott:

“I’ve worked harder and learned more in the last 4 days than I have over the last 9 months of dental school.”





Warning: This post is several years old and the author has marked it as poor quality (compared to more recent posts). It has been left intact for historical reasons, but but its content (and code) may be inaccurate or poorly written.

Update: More modern Python data visualization code and examples can be found on GitHub: https://github.com/swharden/Python-GUI-examples

spectrogram scrollbars

My project is coming along nicely. This isn’t an incredibly robust spectrograph program, but it sure gets the job done quickly and easily. The code below will produce a real time scrolling spectrograph entirely with Python! It polls the microphone (or default recording device), should work on any OS, and can be adjusted for vertical resolution / FFT frequency discretion resolution. It has some simple functions for filtering (check out the de-trend filter!) and might serve as a good start to a spectrograph / frequency analysis project. It took my a long time to reach this point! I’ve worked with Python before, and dabbled with the Python Imaging Library (PIL), but this is my first experience with real time linear data analysis and high-demand multi-threading. I hope it helps you. Below are screenshots of the program (two running at the same time) listening to the same radio signals (mostly Morse code) with standard output and with the “de-trending filter” activated.

nofilter
filter
import pyaudio
import scipy
import struct
import scipy.fftpack

from Tkinter import *
import threading
import time, datetime
import wckgraph
import math

import Image, ImageTk
from PIL import ImageOps
from PIL import ImageChops
import time
import random
import threading
import scipy


#ADJUST RESOLUTION OF VERTICAL FFT
bufferSize=2**11
#bufferSize=2**8

#ADJUSTS AVERAGING SPEED NOT VERTICAL RESOLUTION
#REDUCE HERE IF YOUR PC CANT KEEP UP
sampleRate=24000
#sampleRate=64000

p = pyaudio.PyAudio()
chunks=[]
ffts=[]
def stream():
        global chunks, inStream, bufferSize
        while True:
                chunks.append(inStream.read(bufferSize))

def record():
        global w, inStream, p, bufferSize
        inStream = p.open(format=pyaudio.paInt16,channels=1,
                rate=sampleRate,input=True,frames_per_buffer=bufferSize)
        threading.Thread(target=stream).start()
        #stream()

def downSample(fftx,ffty,degree=10):
        x,y=[],[]
        for i in range(len(ffty)/degree-1):
                x.append(fftx[i*degree+degree/2])
                y.append(sum(ffty[i*degree:(i+1)*degree])/degree)
        return [x,y]

def smoothWindow(fftx,ffty,degree=10):
        lx,ly=fftx[degree:-degree],[]
        for i in range(degree,len(ffty)-degree):
                ly.append(sum(ffty[i-degree:i+degree]))
        return [lx,ly]

def smoothMemory(ffty,degree=3):
        global ffts
        ffts = ffts+[ffty]
        if len(ffts)< =degree: return ffty ffts=ffts[1:] return scipy.average(scipy.array(ffts),0) def detrend(fftx,ffty,degree=10): lx,ly=fftx[degree:-degree],[] for i in range(degree,len(ffty)-degree): ly.append((ffty[i]-sum(ffty[i-degree:i+degree])/(degree*2)) *2+128) #ly.append(fft[i]-(ffty[i-degree]+ffty[i+degree])/2) return [lx,ly] def graph(): global chunks, bufferSize, fftx,ffty, w if len(chunks)>0:
                data = chunks.pop(0)
                data=scipy.array(struct.unpack("%dB"%(bufferSize*2),data))
                ffty=scipy.fftpack.fft(data)
                fftx=scipy.fftpack.rfftfreq(bufferSize*2, 1.0/sampleRate)
                fftx=fftx[0:len(fftx)/4]
                ffty=abs(ffty[0:len(ffty)/2])/1000
                ffty1=ffty[:len(ffty)/2]
                ffty2=ffty[len(ffty)/2::]+2
                ffty2=ffty2[::-1]
                ffty=ffty1+ffty2
                ffty=(scipy.log(ffty)-1)*120
                fftx,ffty=downSample(fftx,ffty,2)
                updatePic(fftx,ffty)
                reloadPic()

        if len(chunks)>20:
                print "falling behind...",len(chunks)

def go(x=None):
        global w,fftx,ffty
        print "STARTING!"
        threading.Thread(target=record).start()
        while True:
                #record()
                graph()


def updatePic(datax,data):
     global im, iwidth, iheight
     strip=Image.new("L",(1,iheight))
     if len(data)>iheight:
             data=data[:iheight-1]
     #print "MAX FREQ:",datax[-1]
     strip.putdata(data)
     #print "%03d, %03d" % (max(data[-100:]), min(data[-100:]))
     im.paste(strip,(iwidth-1,0))
     im=im.offset(-1,0)
     root.update()

def reloadPic():
     global im, lab
     lab.image = ImageTk.PhotoImage(im)
     lab.config(image=lab.image)


root = Tk()
im=Image.open('./ramp.tif')
im=im.convert("L")
iwidth,iheight=im.size
im=im.crop((0,0,500,480))
#im=Image.new("L",(100,1024))
iwidth,iheight=im.size
root.geometry('%dx%d' % (iwidth,iheight))
lab=Label(root)
lab.place(x=0,y=0,width=iwidth,height=iheight)
go()

UPDATE! I’m not going to post the code for this yet (it’s very messy) but I got this thing to display a spectrograph on a canvas. What’s the advantage of that? Huge, massive spectrographs (thousands of pixels in all directions) can now be browsed in real time using scrollbars, and when you scroll it doesn’t stop recording, and you don’t lose any data! Super cool.

spectrogram scrollbars




Warning: This post is several years old and the author has marked it as poor quality (compared to more recent posts). It has been left intact for historical reasons, but but its content (and code) may be inaccurate or poorly written.

Update: More modern Python data visualization code and examples can be found on GitHub: https://github.com/swharden/Python-GUI-examples

I’m stretching the limits of what these software platforms were designed to to, but I’m impressed such a hacked-together code produces fast, functional results. The code below is the simplest case code I could create which graphs the audio spectrum of the microphone input. It seems to run fine with about 30+ FPS on my modest machine. It should work on Windows and Linux. I chose not to go with matplotlib because I didn’t think it was fast enough for my needs in this one case. Here’s what the code below looks like running:

python real time tk wav fft

NOTE that this program was designed with the intent of recording the FFTs, therefore if the program “falls behind” the real time input, it will buffer the sound on its own and try to catch up (accomplished by two layers of threading). In this way, all audio gets interpreted. If you’re just trying to create a spectrograph for simple purposes, have it only sample the audio when it needs to, rather than having it sample audio continuously.

import pyaudio
import scipy
import struct
import scipy.fftpack

from Tkinter import *
import threading
import time, datetime
import wckgraph
import math

#ADJUST THIS TO CHANGE SPEED/SIZE OF FFT
bufferSize=2**11
#bufferSize=2**8

# ADJUST THIS TO CHANGE SPEED/SIZE OF FFT
sampleRate=48100
#sampleRate=64000

p = pyaudio.PyAudio()
chunks=[]
ffts=[]
def stream():
        global chunks, inStream, bufferSize
        while True:
                chunks.append(inStream.read(bufferSize))

def record():
        global w, inStream, p, bufferSize
        inStream = p.open(format=pyaudio.paInt16,channels=1,
                rate=sampleRate,input=True,frames_per_buffer=bufferSize)
        threading.Thread(target=stream).start()

def downSample(fftx,ffty,degree=10):
        x,y=[],[]
        for i in range(len(ffty)/degree-1):
                x.append(fftx[i*degree+degree/2])
                y.append(sum(ffty[i*degree:(i+1)*degree])/degree)
        return [x,y]

def smoothWindow(fftx,ffty,degree=10):
        lx,ly=fftx[degree:-degree],[]
        for i in range(degree,len(ffty)-degree):
                ly.append(sum(ffty[i-degree:i+degree]))
        return [lx,ly]

def smoothMemory(ffty,degree=3):
        global ffts
        ffts = ffts+[ffty]
        if len(ffts)< =degree: return ffty
        ffts=ffts[1:]
        return scipy.average(scipy.array(ffts),0)

def detrend(fftx,ffty,degree=10):
        lx,ly=fftx[degree:-degree],[]
        for i in range(degree,len(ffty)-degree):
                ly.append(ffty[i]-sum(ffty[i-degree:i+degree])/(degree*2))
                #ly.append(fft[i]-(ffty[i-degree]+ffty[i+degree])/2)
        return [lx,ly]

def graph():
        global chunks, bufferSize, fftx,ffty, w
        if len(chunks)>0:
                data = chunks.pop(0)
                data=scipy.array(struct.unpack("%dB"%(bufferSize*2),data))
                #print "RECORDED",len(data)/float(sampleRate),"SEC"
                ffty=scipy.fftpack.fft(data)
                fftx=scipy.fftpack.rfftfreq(bufferSize*2, 1.0/sampleRate)
                fftx=fftx[0:len(fftx)/4]
                ffty=abs(ffty[0:len(ffty)/2])/1000
                ffty1=ffty[:len(ffty)/2]
                ffty2=ffty[len(ffty)/2::]+2
                ffty2=ffty2[::-1]
                ffty=ffty1+ffty2
                ffty=scipy.log(ffty)-2
                #fftx,ffty=downSample(fftx,ffty,5)
                #fftx,ffty=detrend(fftx,ffty,30)
                #fftx,ffty=smoothWindow(fftx,ffty,10)
                ffty=smoothMemory(ffty,3)
                #fftx,ffty=detrend(fftx,ffty,10)
                w.clear()
                #w.add(wckgraph.Axes(extent=(0, -1, fftx[-1], 3)))
                w.add(wckgraph.Axes(extent=(0, -1, 6000, 3)))
                w.add(wckgraph.LineGraph([fftx,ffty]))
                w.update()
        if len(chunks)>20:
                print "falling behind...",len(chunks)

def go(x=None):
        global w,fftx,ffty
        print "STARTING!"
        threading.Thread(target=record).start()
        while True:
                graph()

root = Tk()
root.title("SPECTRUM ANALYZER")
root.geometry('500x200')
w = wckgraph.GraphWidget(root)
w.pack(fill=BOTH, expand=1)
go()
mainloop()




Warning: This post is several years old and the author has marked it as poor quality (compared to more recent posts). It has been left intact for historical reasons, but but its content (and code) may be inaccurate or poorly written.

tk scrolling

The goal is simple: have a very large image (larger than the window) automatically scroll across a Python-generated GUI window. I already have the code created to generate spectrograph images in real time, now I just need a way to have them displayed in real time. At first I tried moving the coordinates of my images and even generating new images with create_image(), but everything I did resulted in a tacky “flickering” effect (not to mention it was slow). Thankfully I found that self.canv.move(self.imgtag,-1,0) can move a specific item (self.imgtag) by a specified amount and it does it smoothly (without flickering). Here’s some sample code. Make sure “snip.bmp” is a big image in the same folder as this script

from Tkinter import *
import Image, ImageTk

class scrollingImage(Frame):

     def go(self):
		self.canv.move(self.imgtag,-1,0)
		self.canv.update()
		self.after(100,self.go)

     def __init__(self, parent=None):
		Frame.__init__(self, parent)
		self.master.title("Spectrogram Viewer")
		self.pack(expand=YES, fill=BOTH)
		self.canv = Canvas(self, relief=SUNKEN)
		self.canv.config(width=200, height=200)
		self.canv.config(highlightthickness=0)

		sbarV = Scrollbar(self, orient=VERTICAL)
		sbarH = Scrollbar(self, orient=HORIZONTAL)

		sbarV.config(command=self.canv.yview)
		sbarH.config(command=self.canv.xview)

		self.canv.config(yscrollcommand=sbarV.set)
		self.canv.config(xscrollcommand=sbarH.set)

		sbarV.pack(side=RIGHT, fill=Y)
		sbarH.pack(side=BOTTOM, fill=X)

		self.canv.pack(side=LEFT, expand=YES, fill=BOTH)
		self.im=Image.open("./snip.bmp")
		width,height=self.im.size
		#self.canv.config(scrollregion=(0,0,width,height))
		self.canv.config(scrollregion=(0,0,300,300))
		self.im2=ImageTk.PhotoImage(self.im)
		x,y=0,0
		self.imgtag=self.canv.create_image(x,y,
		     anchor="nw",image=self.im2)
		self.go()

scrollingImage().mainloop()

Alternatively, I found a way to accomplish a similar thing with PyGame. I’ve decided not to use PyGame for my software package however, because it’s too specific and can’t be run well alongside Tk windows, and it would be insanely hard to add scrollbars to the window. However it’s extremely effective at scrolling images smoothly. Anyhow, here’s the code:

import pygame
from PIL import Image

im=Image.open("1hr_original.jpg")
graphic = pygame.image.fromstring(im.tostring(),im.size,im.mode)
screen = pygame.display.set_mode((400, 300))
clock = pygame.time.Clock()
running = 1
x,y=0,0
while running:
   clock.tick(30)
   for event in pygame.event.get():    #get user input
      if event.type == pygame.QUIT:    #if user clicks the close X
           running = 0                 #make running 0 to break out of loop
   screen.blit(graphic, (x, y))
   pygame.display.flip()   #Update screen
   x-=1




Warning: This post is several years old and the author has marked it as poor quality (compared to more recent posts). It has been left intact for historical reasons, but but its content (and code) may be inaccurate or poorly written.

I wrote a program to display extremely large images in Python using TK. It’s interesting how simple this program is, yet frustrating how long it took me to figure out.

specview

This little Python program will load an image (pretty much any format) using the Python Imaging Library (PIL, which must be installed) and allows you to see it on a scrollable canvas (in two directions) with Tkinter and ImageTk. The above screenshot is of the program viewing the image below:

1hr_original

What is that image? I won’t get ahead of myself, but it’s about 5kHz of audio from 10.140mHz which includes a popular QRSS calling frequency. The image displays an hour of data. My ultimate goal is to have it scroll in the TK window, with slide-adjustable brightness/contrast/etc.

from Tkinter import *
import Image, ImageTk

class ScrolledCanvas(Frame):
     def __init__(self, parent=None):
          Frame.__init__(self, parent)
          self.master.title("Spectrogram Viewer")
          self.pack(expand=YES, fill=BOTH)
          canv = Canvas(self, relief=SUNKEN)
          canv.config(width=400, height=200)
          canv.config(highlightthickness=0)

          sbarV = Scrollbar(self, orient=VERTICAL)
          sbarH = Scrollbar(self, orient=HORIZONTAL)

          sbarV.config(command=canv.yview)
          sbarH.config(command=canv.xview)

          canv.config(yscrollcommand=sbarV.set)
          canv.config(xscrollcommand=sbarH.set)

          sbarV.pack(side=RIGHT, fill=Y)
          sbarH.pack(side=BOTTOM, fill=X)

          canv.pack(side=LEFT, expand=YES, fill=BOTH)
          self.im=Image.open("./1hr_original.jpg")
          width,height=self.im.size
          canv.config(scrollregion=(0,0,width,height))
          self.im2=ImageTk.PhotoImage(self.im)
          self.imgtag=canv.create_image(0,0,anchor="nw",image=self.im2)

ScrolledCanvas().mainloop()
 




I wanted a way to have a bunch of Morse code mp3s on my mp3 player (with a WPM/speed that I decide and I found an easy way to do it with Linux. Rather than downloading existing mp3s of boring text, I wanted to be able to turn ANY text into Morse code, so I could copy something interesting (perhaps the news? hackaday? bash.org?). It’s a little devious, but my plan is to practice copying Morse code during class when lectures become monotonous. [The guy who teaches about infectious diseases is the most boring person I ever met, I learn nothing from class, and on top of that he doesn’t allow laptops to be out!] So, here’s what I did in case it helps anyone else out there…

Step 0: GET THE REQUIRED PROGRAMS: Make sure you have installed Python, cwtext, and lame. Now you’re ready to roll!

Step 1: PREPARE SOME TEXT: I went to Wikipedia and copy/pasted an ENTIRE article into a text file called in.txt. Don’t worry about special characters (such as ” and * and #), we’ll fix them with the following python script.

import time
f=open("out.txt")
raw=f.read()
f.close()

cmd  = """echo "TEST" | cwpcm -w 7 | """
cmd += """lame -r -m m -b 8 --resample 8 -q9 - - > text.mp3"""

import os
i=0
for chunk in raw.split("n")[5:]:
        if chunk.count(" ")>50:
                i+=1
                print "nnfile",i, chunk.count(" "), "wordsn"
		do = cmd.replace("TEST",chunk).replace("text","%02d"%i)
		print "running:",do,
		time.sleep(1)
		print "nnSTART ...",
                os.system(do)
		print "DONE"

Step 2: MAKE MP3s OF THE TEXT!:There should be a new file, out.txt, which is cleaned-up nicely. Run the following script to turn every paragraph of text with more than 50 words into an mp3 file…

f=open("out.txt")
raw=f.read()
f.close()
cmd = """echo "TEST" | cwpcm -w 13 | sox -r 44k -u -b 8 -t raw - text.wav"""
cmd+="""; lame --preset phone text.wav text.mp3; rm text.wav"""
import os
i=0
for chunk in raw.split("n")[5:]:
	if chunk.count(" ")>50:
		i+=1
		print i, chunk.count(" "), "words"
		os.system(cmd.replace("TEST",chunk).replace("text","%02d"%i))

Now you should have a directory filled with mp3 files which you can skip through (or shuffle!) using your handy dandy mp3 player. Note that “-w 13” means 13 WPM (words per minute). Simply change that number to change the speed.

Good luck with your CW practice!





Update: The DIY ECG project has had several iterations. The latest one can be viewed here: https://www.swharden.com/wp/2019-03-15-sound-card-ecg-with-ad8232/

Warning: This post is several years old and the author has marked it as poor quality (compared to more recent posts). It has been left intact for historical reasons, but but its content (and code) may be inaccurate or poorly written.

Background

ecgman

You’ve probably seen somebody in a hospital hooked up to a bunch of wires used to analyze their heartbeat. The machine that visualizes heartbeat is called an electrocardiograph, or ECG. It amplifies, measures, and records the natural electrical potential created by the heart. Note that cardiac electrical signals are different than heart sounds, which are listened to with a stethoscope. The intrinsic cardiac pacemaker system is responsible for generating these electrical signals which serve to command and coordinate contraction of the four chambers at the heart at the appropriate intervals: atria (upper chambers) first, then the ventricles (lower chambers) a fraction of a second later. Analysis of these signals reveals a wealth of information about cardiac regulation, as well insights into pathological conditions.

ecg_principle_slow

Each heartbeat produces a little squiggle with unique properties. The squiggle is called the PQRST wave. The smooth curve in the ECG (P) is caused by the stimulation of the atria via the Sinoatrial (SA) node in the right atrium. There is a brief pause, as the electrical impulse is slowed by the Atrioventricular (AV) node and Purkinje fibers in the bundle of His. The prominent spike in the ECG (the QRS complex) is caused by this step, where the electrical impulse travels through the inter-ventricular septum and up through the outer walls of the ventricles. The sharp peak is the R component, and exact heart rate can be calculated as the inverse of the R-to-R interval (RRi).

WARNING: This page documents how I made a simple ECG machine with a minimum of parts to view the electrical activity of my own heart. Feel free to repeat my experiment, but do so at your own risk.

Project Goal

The goal of this project is to generate an extremely cheap, functional ECG machine made from common parts, most of which can be found around your house. This do-it-yourself (DIY) ECG project is different than many others on the internet in that it greatly simplifies the circuitry by eliminating noise reduction components, accomplishing this via software-based data post-processing. Additionally, this writeup is intended for those without any computer, electrical, or biomedical experience, and should be far less convoluted than the suspiciously-cryptic write-ups currently available online. In short, I want to give everybody the power to visualize and analyze their own heartbeat!

The ECG of my own heart:

ecg31

Video Overview

I know a lot of Internet readers aren’t big fans of reading. Therefore, I provided an outline of the process in video form. Check out the videos, and if you like what you see read more!

Video 1/3: Introducing my ECG machine

Video 2/3: Recording my ECG

Video 3/3: Analyzing my ECG

Electrical Theory

Measurement: The electrical signals which command cardiac musculature can be detected on the surface of the skin. In theory one could grab the two leads of a standard volt meter, one with each hand, and see the voltage change as their heart beats, but the fluctuations are rapid and by the time these signals reach the skin they are extremely weak (a few millionths of a volt) and difficult to detect with simple devices. Therefore, amplification is needed.

Amplification: A simple way to amplify the electrical difference between two points is to use a operational amplifier, otherwise known as an op-amp. The gain (multiplication factor) of an op-amp is controlled by varying the resistors attached to it, and an op-amp with a gain of 1000 will take a 1 millivolt signal and amplify it to 1 volt. There are many different types of microchip op-amps, and they’re often packaged with multiple op-amps in one chip (such as the quad-op-amp lm324, or the dual-op-amp lm358n). Any op-amp designed for low voltage will do for our purposes, and we only need one.

Noise: Unfortunately, the heart is not the only source of voltage on the skin. Radiation from a variety of things (computers, cell phones, lights, and especially the wiring in your walls) is absorbed by your skin and is measured with your ECG, in many cases masking your ECG in a sea of electrical noise. The traditional method of eliminating this noise is to use complicated analog circuitry, but since this noise has a characteristic, repeating, high-frequency wave pattern, it can be separated from the ECG (which is much slower in comparison) using digital signal processing computer software!

Digitization: Once amplified, the ECG signal along with a bunch of noise is in analog form. You could display the output with an oscilloscope, but to load it into your PC you need an analog-to-digital converter. Don’t worry! If you’ve got a sound card with a microphone input, you’ve already got one! It’s just that easy. We’ll simply wire the output of our ECG circuit to the input of our sound card, record the output of the op-amp using standard sound recording software, remove the noise from the ECG digitally, and output gorgeous ECG traces ready for visualization and analysis!

Parts/Cost

I’ll be upfront and say that I spent $0.00 making my ECG machine, because I was able to salvage all the parts I needed from a pile of old circuit boards. If you need specific components, check your local RadioShack. If that’s a no-go, hit-up Digikey (it’s probably cheaper too). Also, resistor values are flexible. Use mine as a good starter set, and vary them to suit your needs. If you buy everything from Digikey, the total cost of this project would be about $1. For now, here’s a list of all the parts you need:

  • 1x low voltage op-amp LM358N $0.40
  • 1x 100kOhm resistor (brn,blk,yel) virtually free
  • 1x 1kOhm resistor (brn,blk,red) virtually free
  • 1x 0.1uF capacitor (104Z) virtually free
  • Microphone cable to get from the op-amp to your PC
  • Electrodes 3 pennies should do. ($0.03)

Making the Device

Keep in mind that I’m not an electrical engineer (I have a masters in molecular biology but I’m currently a dental student if you must know) and I’m only reporting what worked well for me. I don’t claim this is perfect, and I’m certainly open for (and welcome) suggestions for improvement. With that in mind, here’s what I did!

This is pretty much it. First off is a power source. If you want to be safe, use three AAA batteries in series. If you’re a daredevil and enjoy showing off your ghettorigging skills, do what I did and grab 5v from a free USB plug! Mua ha ha ha. The power goes into the circuit and so do the leads/electrodes connected to the body. You can get pretty good results with only two leads, but if you want to experiment try hooking up an extra ground lead and slap it on your foot. More on the electrodes later. The signal from the leads is amplified by the circuit and put out the headphone cable, ready to enter your PC’s sound card through the microphone jack!

Note how I left room in the center of the circuit board. That was intentional! I wanted to expand this project by adding a microcontroller to do some on-board, real-time analysis. Specifically, an ATMega8! I never got around to it though. Its purpose would be to analyze the output of the op-amp and graph the ECG on a LCD screen, or at least measure the time between beats and display HR on a screen. (More ideas are at the bottom of this document.) Anyway, too much work for now, maybe I’ll do it one day in the future.

ECG circuit diagram:

This is the circuit diagram. This is a classical high-gain analog differential amplifier. It just outputs the multiplied difference of the inputs. The 0.1uF capacitor helps stabilize the signal and reduce high frequency noise (such as the audio produced by a nearby AM radio station). Use Google if you’re interested in learning exactly how it works.

ECG schematic:

This is how I used my LM358N to create the circuit above. Note that there is a small difference in my board from the photos and this diagram. This diagram is correct, but the circuit in some of the pictures is not. Briefly, when I built it I accidentally connected the (-) lead directly to ground, rather than to the appropriate pin on the microchip. This required me to place a 220kOhm between the leads to stabilize the signal. I imagine if you wire it CORRECTLY (as shown in these circuit diagrams) it will work fine, but if you find it too finicky (jumping quickly from too loud to too quiet), try tossing in a high-impedance resistor between the leads like I did. Overall, this circuit is extremely flexible and I encourage you to build it on a breadboard and try different things. Use this diagram as a starting point and experiment yourself!

The Electrodes:

You can make electrodes out of anything conductive. The most recent graphs were created from wires with gator clips on them clamping onto pennies (pictured). Yeah, I know I could solder directly to the pennies (they’re copper) but gator clips are fast, easy, and can be clipped to different materials (such as aluminum foil) for testing. A dot of moisturizing lotion applied to the pennies can be used to improve conduction between the pennies and the skin, but I didn’t find this to be very helpful. If pressed firmly on the body, conduction seems to be fine. Oh! I just remembered. USE ELECTRICAL TAPE TO ATTACH LEADS TO YOUR BODY! I tried a million different things, from rubber bands to packaging tape. The bottom line is that electrical tape is stretchy enough to be flexible, sticky enough not to fall off (even when moistened by the natural oils/sweat on your skin), and doesn’t hurt that bad to peel off.

Some of the best electrodes I used were made from aluminum cans! Rinse-out a soda can, cut it into “pads”, and use the sharp edge of a razor blade or pair of scissors to scrape off the wax coating on all contact surfaces. Although a little unconformable and prone to cut skin due to their sharp edges, these little guys work great!

Hooking it Up

This part is the most difficult part of the project! This circuit is extremely finicky. The best way to get it right is to open your sound editor (In Windows I use GoldWave because it’s simple, powerful, and free, but similar tools exist for Linux and other Unix-based OSes) and view the low-frequency bars in live mode while you set up. When neither electrode is touched, it should be relatively quiet. When only the + electrode is touched, it should go crazy with noise. When you touch both (one with each hand) the noise should start to go away, possibly varying by how much you squeeze (how good of a connection you have). The whole setup process is a game between too much and too little conduction. You’ll find that somewhere in the middle, you’ll see (and maybe hear) a low-frequency burst of noise once a second corresponding to your heartbeat. [note: Did you know that’s how the second was invented? I believe it was ] Once you get that good heartbeat, tape up your electrodes and start recording. If you can’t get it no matter what you do, start by putting the ground electrode in your mouth (yeah, I said it) and pressing the + electrode firmly and steadily on your chest. If that works (it almost always does), you know what to look for, so keep trying on your skin. For short recordings (maybe just a few beats) the mouth/chest method works beautifully, and requires far less noise reduction (if any), but is simply impractical for long-term recordings. I inside vs. outside potential is less susceptible to noise-causing electrical radiation. Perhaps other orifices would function similarly? I’ll leave it at that. I’ve also found that adding a third electrode (another ground) somewhere else on my body helps a little, but not significantly. Don’t give up at this step if you don’t get it right away! If you hear noise when + is touched, your circuit is working. Keep trying and you’ll get it eventually.

Recording the ECG

This is the easy part. Keep an eye on your “bars” display in the audio program to make sure something you’re doing (typing, clicking, etc) isn’t messing up the recording. If you want, try surfing the net or playing computer games to see how your heart varies. Make sure that as you tap the keyboard and click the mouse, you’re not getting noise back into your system. If this is a problem, try powering your device by batteries (a good idea for safety’s sake anyway) rather than another power source (such as USB power). Record as long as you want! Save the file as a standard, mono, wave file.

Digitally Eliminating Noise

Now it’s time to clean-up the trace. Using GoldWave, first apply a lowpass filter at 30 Hz. This kills most of your electrical noise (> 30hz), while leaving the ECG intact (< 15Hz). However, it dramatically decreases the volume (potential) of the audio file. Increase the volume as necessary to maximize the window with the ECG signal. You should see clear heartbeats at this point. You may want to apply an auto-gain filter to normalize the heartbeats potentials. Save the file as a raw sound file (.snd) at 1000 Hz (1 kHz) resolution.

Presentation and Analysis

Now you’re ready to analyze! Plop your .snd file in the same folder as my [ecg.py script], edit the end of the script to reflect your .snd filename, and run the script by double-clicking it. (Keep in mind that my script was written for python 2.5.4 and requires numpy 1.3.0rc2 for python 2.5, and matplotlib 0.99 for python 2.5 – make sure you get the versions right!) Here’s what you’ll see!

This is a small region of the ECG trace. The “R” peak is most obvious, but the details of the other peaks are not as visible. If you want more definition in the trace (such as the blue one at the top of the page), consider applying a small collection of customized band-stop filters to the audio file rather than a single, sweeping lowpass filter. Refer to earlier posts in the DIY ECG category for details. Specifically, code on Circuits vs. Software for noise reduction entry can help. For our purposes, calculating heart rate from R-to-R intervals (RRIs) can be done accurately with traces such as this.

Your heart rate fluctuates a lot over time! By plotting the inverse of your RRIs, you can see your heart rate as a function of time. Investigate what makes it go up, go down, and how much. You’d be surprised by what you find. I found that checking my email raises my heart rate more than first-person-shooter video games. I get incredibly anxious when I check my mail these days, because I fear bad news from my new university (who knows why, I just get nervous about it). I wonder if accurate RRIs could be used to assess nervousness for the purposes of lie detection?

This is the RRI plot where the value of each RRI (in milliseconds) is represented for each beat. It’s basically the inverse of heart rate. Miscalculated heartbeats would show up as extremely high or extremely low dots on this graph. However, excluding points above or below certain bounds means that if your heart did double-beat, or skip a beat, you wouldn’t see it. Note that I just realized my axis label is wrong (it should be sec, not ms).

A Poincare Plot is a commonly-used method to visually assess heart rate variability as a function of RRIs. In this plot, each RRI is plotted against the RRI of the next subsequent beat. In a heart which beats at the same speed continuously, only a single dot would be visible in the center. In a heart which beats mostly-continuously, and only changes its rate very slowly, a linear line of dots would be visible in a 1:1 ratio. However, in real life the heart varies RRIs greatly from beat to beat, producing a small cloud of dots. The size of the cloud corresponds to the speed at which the autonomic nervous system can modulate heart rate in the time frame of a single beat.

The frequency of occurrence of various RRIs can be expressed by a histogram. The center peak corresponds to the standard heart rate. Peaks to the right and left of the center peak correspond to increased and decreased RRIs, respectively. A gross oversimplification of the interpretation of such data would be to state that the upper peak represents the cardio-inhibitory parasympathetic autonomic nervous system component, and the lower peak represents the cardio-stimulatory sympathetic autonomic nervous system component.

Taking the Fast Fourier Transformation of the data produces a unique trace whose significance is extremely difficult to interpret. Near 0Hz (infinite time) the trace heads toward ∞ (infinite power). To simplify the graph and eliminate the near-infinite, low-frequency peak we will normalize the trace by multiplying each data point by its frequency, and dividing the vertical axis units by Hz to compensate. This will produce the following graph…

This is the power spectrum density (PSD) plot of the ECG data we recorded. Its physiological interpretation is extraordinarily difficult to understand and confirm, and is the subject of great debate in the field of autonomic neurological cardiac regulation. An oversimplified explanation of the significance of this graph is that the parasympathetic (cardio-inhibitory) branch of the autonomic nervous system works faster than the sympathetic (cardio-stimulatory) branch. Therefore, the lower peak corresponds to the sympathetic component (combined with persistent parasympathetic input, it’s complicated), while the higher-frequency peak corresponds to the parasympathetic component, and the sympathetic/parasympathetic relationship can be assessed by the ratio of the integrated areas of these peaks after a complicated curve fitting processes which completely separates overlapping peaks. To learn more about power spectral analysis of heart rate over time in the frequency domain, I recommend skimming this introduction to heart rate variability website and the article on Heart Rate Variability following Myocardial Infarction (heart attack). Also, National Institute of Health (NIH) funded studies on HRV should be available from pubmed.org. If you want your head to explode, read Frequency-Domain Characteristics and Filtering of Blood Flow Following the Onset of Exercise: Implications for Kinetics Analysis for a lot of good frequency-domain-analysis-related discussion and rationalization.

Encouraging Words:

Please, if you try this don’t die. The last thing I want is to have some kid calling me up and yelling at me that he nearly electrocuted himself when he tried to plug my device directly into a wall socket and now has to spend the rest of his life with two Abraham Lincolns tattooed onto his chest resembling a second set of nipples. Please, if you try this use common sense, and of course you’re responsible for your own actions. I provide this information as a description of what I did and what worked for me. If you make something similar that works, I’ve love to see it! Send in your pictures of your circuit, charts of your traces, improved code, or whatever you want and I’ll feature it on the site. GOOD LUCK!

Fancier Circuit:

If you want to try this, go for it! Briefly, this circuit uses 6 op-amps to help eliminate effects of noise. It’s also safer, because of the diodes interconnecting the electrodes. It’s the same circuit as on [this page].

Last minute thoughts:

  • More homemade ECG information can be found on my earlier posts in the DIY ECG category, however this page is the primary location of my most recent thoughts and ideas.
  • You can use moisturizing lotion between the electrodes and your skin to increase conduction. However, keep in mind that better conduction is not always what you want. You’ll have to experiment for yourself.
  • Variation in location of electrodes will vary the shape of the ECG. I usually place electrodes on each side of my chest near my arms. If your ECG appears upside-down, reverse the leads!
  • Adding extra leads can improve grounding. Try grounding one of your feet with a third lead to improve your signal. Also, if you’re powering your device via USB power consider trying battery power – it should be less noisy.
  • While recording, be aware of what you do! I found that if I’m not well-grounded, my ECG is fine as long as I don’t touch my keyboard. If I start typing, every keypress shows up as a giant spike, bigger than my heartbeat!
  • If you get reliable results, I wonder if you could make the device portable? Try using a portable tape recorder, voice recorder, or maybe even minidisc recorder to record the output of the ECG machine for an entire day. I haven’t tried it, but why wouldn’t it work? If you want to get fancy, have a microcontroller handle the signal processing and determine RRIs (should be easy) and save this data to a SD card or fancy flash logger.
  • The microcontroller could output heart rate via the serial port.
  • If you have a microcontroller on board, why not display heart rate on a character LCD?
  • While you have a LCD on there, display the ECG graphically!
  • Perhaps a wireless implementation would be useful.
  • Like, I said, there are other, more complicated analog circuits which reduce noise of the outputted signal. I actually built Jason Nguyen’s fancy circuit which used 6 op-amps but the result wasn’t much better than the simple, 1 op-amp circuit I describe here once digital filtering was applied.
  • Arrhythmic heartbeats (where your heart screws-up and misfires, skips a beat, double-beats, or beats awkwardly) are physiological (normal) and surprisingly common. Although shocking to hear about, sparse, single arrhythmic heartbeats are normal and are a completely different ball game than chronic, potentially deadly heart arrhythmias in which every beat is messed-up. If you’re in tune with your body, you might actually feel these occurrences happening. About three times a week I feel my heart screw up a beat (often when it’s quiet), and it feels like a sinking feeling in my chest. I was told by a doctor that it’s totally normal and happens many times every day without me noticing, and that most people never notice these single arrhythmic beats. I thought it was my heart skipping a beat, but I wasn’t sure. That was my motivation behind building this device – I wanted to see what my arrhythmic beats looked like. It turns out that it’s more of a double-beat than a skipped beat, as observed when I captured a single arrhythmic heartbeat with my ECG machine, as described in this entry.
  • You can improve the safety of this device by attaching diodes between leads, similar to the more complicated circuit. Theory is that if a huge surge of energy does for whatever reason get into the ECG circuit, it’ll short itself out at the circuit level (conducting through the diodes) rather than at your body (across your chest / through your heart).
  • Alternatively, use an AC opto-isolator between the PC sound card and the ECG circuit to eliminate the possibility of significant current coming back from the PC.
  • On the Hackaday post, Flemming Frandsen noted that an improperly grounded PC could be dangerous because the stored charge would be manifest in the ground of the microphone jack. If you were to ground yourself to true ground (using a bench power supply or sticking your finger in the ground socket of an AC wall plug) this energy could travel through you! So be careful to only ground yourself with respect to the circuit using only battery power to minimize this risk.
  • Do not attempt anything on this page. Ever. Don’t even read it. You read it already! You’re sill reading it aren’t you? Yeah. You don’t follow directions well do you?

SAMPLE FILTERED RECORDING:

I think this is the same one I used in the 3rd video from my single op-amp circuit. [scottecg.snd] It’s about an hour long, and in raw sound format (1000 Hz). It’s already been filtered (low-pass filtered at 30Hz). You can use it with my code below!

CODE

print "importing libraries..."
import numpy, pylab
print "DONE"

class ECG:

    def trim(self, data,degree=100):
        print 'trimming'
        i,data2=0,[]
        while i<len(data):
            data2.append(sum(data[i:i+degree])/degree)
            i+=degree
        return data2

    def smooth(self,list,degree=15):
        mults=[1]
        s=[]
        for i in range(degree): mults.append(mults[-1]+1)
        for i in range(degree): mults.append(mults[-1]-1)
        for i in range(len(list)-len(mults)):
            small=list[i:i+len(mults)]
            for j in range(len(small)):
                small[j]=small[j]*mults[j]
            val=sum(small)/sum(mults)
            s.append(val)
        return s

    def smoothWindow(self,list,degree=10):
        list2=[]
        for i in range(len(list)):
            list2.append(sum(list[i:i+degree])/float(degree))
        return list2

    def invertYs(self):
        print 'inverting'
        self.ys=self.ys*-1

    def takeDeriv(self,dist=5):
        print 'taking derivative'
        self.dys=[]
        for i in range(dist,len(self.ys)):
            self.dys.append(self.ys[i]-self.ys[i-dist])
        self.dxs=self.xs[0:len(self.dys)]

    def genXs(self, length, hz):
        print 'generating Xs'
        step = 1.0/(hz)
        xs=[]
        for i in range(length): xs.append(step*i)
        return xs

    def loadFile(self, fname, startAt=None, length=None, hz=1000):
        print 'loading',fname
        self.ys = numpy.memmap(fname, dtype='h', mode='r')*-1
        print 'read %d points.'%len(self.ys)
        self.xs = self.genXs(len(self.ys),hz)
        if startAt and length:
            self.ys=self.ys[startAt:startAt+length]
            self.xs=self.xs[startAt:startAt+length]

    def findBeats(self):
        print 'finding beats'
        self.bx,self.by=[],[]
        for i in range(100,len(self.ys)-100):
          if self.ys[i]<15000: continue # SET THIS VISUALLY
          if self.ys[i]<self.ys[i+1] or self.ys[i]<self.ys[i-1]: continue
          if self.ys[i]-self.ys[i-100]>5000 and self.ys[i]-self.ys[i+100]>5000:
              self.bx.append(self.xs[i])
              self.by.append(self.ys[i])
        print "found %d beats"%(len(self.bx))

    def genRRIs(self,fromText=False):
        print 'generating RRIs'
        self.rris=[]
        if fromText: mult=1
        else: 1000.0
        for i in range(1,len(self.bx)):
            rri=(self.bx[i]-self.bx[i-1])*mult
            #if fromText==False and len(self.rris)>1:
                #if abs(rri-self.rris[-1])>rri/2.0: continue
            #print i, "%.03ft%.03ft%.2f"%(bx[i],rri,60.0/rri)
            self.rris.append(rri)

    def removeOutliers(self):
        beatT=[]
        beatRRI=[]
        beatBPM=[]
        for i in range(1,len(self.rris)):
            #CHANGE THIS AS NEEDED
            if self.rris[i]<0.5 or self.rris[i]>1.1: continue
            if abs(self.rris[i]-self.rris[i-1])>self.rris[i-1]/5: continue
            beatT.append(self.bx[i])
            beatRRI.append(self.rris[i])
        self.bx=beatT
        self.rris=beatRRI

    def graphTrace(self):
        pylab.plot(self.xs,self.ys)
        #pylab.plot(self.xs[100000:100000+4000],self.ys[100000:100000+4000])
        pylab.title("Electrocardiograph")
        pylab.xlabel("Time (seconds)")
        pylab.ylabel("Potential (au)")

    def graphDeriv(self):
        pylab.plot(self.dxs,self.dys)
        pylab.xlabel("Time (seconds)")
        pylab.ylabel("d/dt Potential (au)")

    def graphBeats(self):
        pylab.plot(self.bx,self.by,'.')

    def graphRRIs(self):
        pylab.plot(self.bx,self.rris,'.')
        pylab.title("Beat Intervals")
        pylab.xlabel("Beat Number")
        pylab.ylabel("RRI (ms)")

    def graphHRs(self):
        #HR TREND
        hrs=(60.0/numpy.array(self.rris)).tolist()
        bxs=(numpy.array(self.bx[0:len(hrs)])/60.0).tolist()
        pylab.plot(bxs,hrs,'g',alpha=.2)
        hrs=self.smooth(hrs,10)
        bxs=bxs[10:len(hrs)+10]
        pylab.plot(bxs,hrs,'b')
        pylab.title("Heart Rate")
        pylab.xlabel("Time (minutes)")
        pylab.ylabel("HR (bpm)")

    def graphPoincare(self):
        #POINCARE PLOT
        pylab.plot(self.rris[1:],self.rris[:-1],"b.",alpha=.5)
        pylab.title("Poincare Plot")
        pylab.ylabel("RRI[i] (sec)")
        pylab.xlabel("RRI[i+1] (sec)")

    def graphFFT(self):
        #PSD ANALYSIS
        fft=numpy.fft.fft(numpy.array(self.rris)*1000.0)
        fftx=numpy.fft.fftfreq(len(self.rris),d=1)
        fftx,fft=fftx[1:len(fftx)/2],abs(fft[1:len(fft)/2])
        fft=self.smoothWindow(fft,15)
        pylab.plot(fftx[2:],fft[2:])
        pylab.title("Raw Power Sprectrum")
        pylab.ylabel("Power (ms^2)")
        pylab.xlabel("Frequency (Hz)")

    def graphFFT2(self):
        #PSD ANALYSIS
        fft=numpy.fft.fft(numpy.array(self.rris)*1000.0)
        fftx=numpy.fft.fftfreq(len(self.rris),d=1)
        fftx,fft=fftx[1:len(fftx)/2],abs(fft[1:len(fft)/2])
        fft=self.smoothWindow(fft,15)
        for i in range(len(fft)):
            fft[i]=fft[i]*fftx[i]
        pylab.plot(fftx[2:],fft[2:])
        pylab.title("Power Sprectrum Density")
        pylab.ylabel("Power (ms^2)/Hz")
        pylab.xlabel("Frequency (Hz)")

    def graphHisto(self):
        pylab.hist(self.rris,bins=20,ec='none')
        pylab.title("RRI Deviation Histogram")
        pylab.ylabel("Frequency (count)")
        pylab.xlabel("RRI (ms)")
        #pdf, bins, patches = pylab.hist(self.rris,bins=100,alpha=0)
        #pylab.plot(bins[1:],pdf,'g.')
        #y=self.smooth(list(pdf[1:]),10)
        #x=bins[10:len(y)+10]
        #pylab.plot(x,y)

    def saveBeats(self,fname):
        print "writing to",fname
        numpy.save(fname,[numpy.array(self.bx)])
        print "COMPLETE"

    def loadBeats(self,fname):
        print "loading data from",fname
        self.bx=numpy.load(fname)[0]
        print "loadded",len(self.bx),"beats"
        self.genRRIs(True)

def snd2txt(fname):
    ## SND TO TXT ##
    a=ECG()
    a.loadFile(fname)#,100000,4000)
    a.invertYs()
    pylab.figure(figsize=(7,4),dpi=100);pylab.grid(alpha=.2)
    a.graphTrace()
    a.findBeats()
    a.graphBeats()
    a.saveBeats(fname)
    pylab.show()

def txt2graphs(fname):
    ## GRAPH TXT ##
    a=ECG()
    a.loadBeats(fname+'.npy')
    a.removeOutliers()
    pylab.figure(figsize=(7,4),dpi=100);pylab.grid(alpha=.2)
    a.graphHRs();pylab.subplots_adjust(left=.1,bottom=.12,right=.96)
    pylab.savefig("DIY_ECG_Heart_Rate_Over_Time.png");
    pylab.figure(figsize=(7,4),dpi=100);pylab.grid(alpha=.2)
    a.graphFFT();pylab.subplots_adjust(left=.13,bottom=.12,right=.96)
    pylab.savefig("DIY_ECG_Power_Spectrum_Raw.png");
    pylab.figure(figsize=(7,4),dpi=100);pylab.grid(alpha=.2)
    a.graphFFT2();pylab.subplots_adjust(left=.13,bottom=.12,right=.96)
    pylab.savefig("DIY_ECG_Power_Spectrum_Weighted.png");
    pylab.figure(figsize=(7,4),dpi=100);pylab.grid(alpha=.2)
    a.graphPoincare();pylab.subplots_adjust(left=.1,bottom=.12,right=.96)
    pylab.savefig("DIY_ECG_Poincare_Plot.png");
    pylab.figure(figsize=(7,4),dpi=100);pylab.grid(alpha=.2)
    a.graphRRIs();pylab.subplots_adjust(left=.1,bottom=.12,right=.96)
    pylab.savefig("DIY_ECG_RR_Beat_Interval.png");
    pylab.figure(figsize=(7,4),dpi=100);pylab.grid(alpha=.2)
    a.graphHisto();pylab.subplots_adjust(left=.1,bottom=.12,right=.96)
    pylab.savefig("DIY_ECG_RR_Deviation_Histogram.png");
    pylab.show();

fname='publish_05_10min.snd' #CHANGE THIS AS NEEDED
#raw_input("npress ENTER to analyze %s..."%(fname))
snd2txt(fname)
#raw_input("npress ENTER to graph %s.npy..."%(fname))
txt2graphs(fname)




Warning: This post is several years old and the author has marked it as poor quality (compared to more recent posts). It has been left intact for historical reasons, but but its content (and code) may be inaccurate or poorly written.

Update: this project is now on GitHub https://github.com/FredEckert/pySquelch

I’ve been working on the pySquelch project which is basically a method to graph frequency usage with respect to time. The code I’m sharing below listens to the microphone jack on the sound card (hooked up to a radio) and determines when transmissions begin and end. I ran the code below for 24 hours and this is the result:

1png

This graph represents frequency activity with respect to time. The semi-transparent gray line represents the raw frequency usage in fractional minutes the frequency was tied-up by transmissions. The solid blue line represents the same data but smoothed by 10 minutes (in both directions) by the Gaussian smoothing method modified slightly from my linear data smoothing with Python page.

2png

I used the code below to generate the log, and the code further below to create the graph from the log file. Assuming your microphone is enabled and everything else is working, this software will require you to determine your own threshold for talking vs. no talking. Read the code and you’ll figure out how test your sound card settings.

If you want to try this yourself you need a Linux system (a Windows system version could be created simply by replacing getVolEach() with a Windows-based audio level detection system) with Python and the alsaaudio, numpy, and matplotlib libraries. Try running the code on your own, and if it doesn’t recognize a library “aptitude search” for it. Everything you need can be installed from packages in the common repository.

#pySquelchLogger.py
import time, random, alsaaudio, audioop
inp = alsaaudio.PCM(alsaaudio.PCM_CAPTURE,alsaaudio.PCM_NONBLOCK)
inp.setchannels(2)
inp.setrate(1000)
inp.setformat(alsaaudio.PCM_FORMAT_S8)
inp.setperiodsize(100)
addToLog=""
lastLogTime=0

testLevel=False ### SET THIS TO 'True' TO TEST YOUR SOUNDCARD

def getVolEach():
        # this is a quick way to detect activity.
        # modify this function to use alternate methods of detection.
	while True:
		l,data = inp.read() # poll the audio device
		if l>0: break
	vol = audioop.max(data,1) # get the maximum amplitude
	if testLevel: print vol
	if vol>10: return True ### SET THIS NUMBER TO SUIT YOUR NEEDS ###
	return False

def getVol():
        # reliably detect activity by getting 3 consistant readings.
	a,b,c=True,False,False
	while True:
		a=getVolEach()
		b=getVolEach()
		c=getVolEach()
		if a==b==c:
			if testLevel: print "RESULT:",a
			break
	if a==True: time.sleep(1)
	return a

def updateLog():
        # open the log file, append the new data, and save it again.
	global addToLog,lastLogTime
	#print "UPDATING LOG"
	if len(addToLog)>0:
        	f=open('log.txt','a')
        	f.write(addToLog)
        	f.close()
        	addToLog=""
	lastLogTime=time.mktime(time.localtime())

def findSquelch():
        # this will record a single transmission and store its data.
	global addToLog
	while True: # loop until we hear talking
		time.sleep(.5)
		if getVol()==True:
			start=time.mktime(time.localtime())
			print start,
			break
	while True: # loop until talking stops
		time.sleep(.1)
		if getVol()==False:
			length=time.mktime(time.localtime())-start
			print length
			break
	newLine="%d,%d "%(start,length)
	addToLog+=newLine
	if start-lastLogTime>30: updateLog() # update the log

while True:
	findSquelch()

The logging code (above) produces a log file like this (below). The values represent the start time of each transmission (in seconds since epoch) followed by the duration of the transmission.

#log.txt
1245300044,5 1245300057,4 1245300063,16 1245300094,13 1245300113,4 1245300120,14 1245300195,4 1245300295,4 1245300348,4 1245300697,7 1245300924,3 1245301157,4 1245301207,12 1245301563,4 1245302104,6 1245302114,6 1245302192,3 1245302349,4 1245302820,4 1245304812,13 1245308364,10 1245308413,14 1245312008,14 1245313953,11 1245314008,6 1245314584,4 1245314641,3 1245315212,5 1245315504,6 1245315604,13 1245315852,3 1245316255,6 1245316480,5 1245316803,3 1245316839,6 1245316848,11 1245316867,5 1245316875,12 1245316893,13 1245316912,59 1245316974,12 1245316988,21 1245317011,17 1245317044,10 1245317060,6 1245317071,7 1245317098,33 1245317140,96 1245317241,15 1245317259,14 1245317277,8 1245317298,18 1245317322,103 1245317435,40 1245317488,18 1245317508,34 1245317560,92 1245317658,29 1245317697,55 1245317755,33 1245317812,5 1245317818,7 1245317841,9 1245317865,25 1245317892,79 1245317972,30 1245318007,8 1245318021,60 1245318083,28 1245318114,23 1245318140,25 1245318167,341 1245318512,154 1245318670,160 1245318834,22 1245318859,9 1245318870,162 1245319042,57 1245319102,19 1245319123,30 1245319154,18 1245319206,5 1245319214,13 1245319229,6 1245319238,6 1245319331,9 1245319341,50 1245319397,71 1245319470,25 1245319497,40 1245319540,8 1245319551,77 1245319629,4 1245319638,36 1245319677,158 1245319837,25 1245319865,40 1245319907,33 1245319948,92 1245320043,26 1245320100,9 1245320111,34 1245320146,8 1245320159,6 1245320167,8 1245320181,12 1245320195,15 1245320212,14 1245320238,18 1245320263,46 1245320310,9 1245320326,22 1245320352,27 1245320381,15 1245320398,24 1245320425,57 1245320483,16 1245320501,40 1245320543,43 1245320589,65 1245320657,63 1245320722,129 1245320853,33 1245320889,50 1245320940,1485 1245322801,7 1245322809,103 1245322923,5 1245322929,66 1245323553,4 1245324203,15 1245324383,5 1245324570,7 1245324835,4 1245325200,8 1245325463,5 1245326414,12 1245327340,12 1245327836,4 1245327973,4 1245330006,12 1245331244,11 1245331938,11 1245332180,5 1245332187,81 1245332573,5 1245333609,12 1245334447,10 1245334924,9 1245334945,4 1245334971,4 1245335031,9 1245335076,11 1245335948,16 1245335965,27 1245335993,113 1245336107,79 1245336187,64 1245336253,37 1245336431,4 1245336588,5 1245336759,7 1245337048,3 1245337206,13 1245337228,4 1245337309,4 1245337486,6 1245337536,8 1245337565,38 1245337608,100 1245337713,25 1245337755,169 1245337930,8 1245337941,20 1245337967,6 1245337978,7 1245337996,20 1245338019,38 1245338060,127 1245338192,30 1245338227,22 1245338250,15 1245338272,15 1245338310,3 1245338508,4 1245338990,5 1245339136,5 1245339489,8 1245339765,4 1245340220,5 1245340233,6 1245340266,10 1245340278,22 1245340307,7 1245340315,28 1245340359,32 1245340395,4 1245340403,41 1245340446,46 1245340494,58 1245340554,17 1245340573,21 1245340599,3 1245340604,5 1245340611,46 1245340661,26 1245340747,4 1245340814,14 1245341043,4 1245341104,4 1245341672,4 1245341896,5 1245341906,3 1245342301,3 1245342649,6 1245342884,5 1245342929,4 1245343314,6 1245343324,10 1245343335,16 1245343353,39 1245343394,43 1245343439,62 1245343561,3 1245343790,4 1245344115,3 1245344189,5 1245344233,4 1245344241,6 1245344408,12 1245344829,3 1245345090,5 1245345457,5 1245345689,4 1245346086,3 1245347112,12 1245348006,14 1245348261,10 1245348873,4 1245348892,3 1245350303,11 1245350355,4 1245350766,5 1245350931,3 1245351605,14 1245351673,55 1245351729,23 1245351754,5 1245352123,37 1245352163,21 1245352186,18 1245352209,40 1245352251,49 1245352305,8 1245352315,5 1245352321,6 1245352329,22 1245352353,48 1245352404,77 1245352483,58 1245352543,17 1245352570,19 1245352635,5 1245352879,3 1245352899,5 1245352954,4 1245352962,6 1245352970,58 1245353031,21 1245353055,14 1245353071,52 1245353131,37 1245353170,201 1245353373,56 1245353431,18 1245353454,47 1245353502,13 1245353519,106 1245353627,10 1245353647,12 1245353660,30 1245353699,42 1245353746,28 1245353776,29 1245353806,9 1245353818,21 1245353841,10 1245353853,6 1245353862,224 1245354226,4 1245354964,63 1245355029,4 1245355036,142 1245355180,148 1245355330,7 1245355338,23 1245355363,9 1245355374,60 1245355437,142 1245355581,27 1245355609,5 1245355615,2 1245355630,64 1245355700,7 1245355709,73 1245355785,45 1245355834,85 1245355925,9 1245356234,5 1245356620,6 1245356629,12 1245356643,29 1245356676,120 1245356798,126 1245356937,62 1245357001,195 1245357210,17 1245357237,15 1245357258,24 1245357284,53 1245357339,2 1245357345,27 1245357374,76 1245357452,28 1245357482,42 1245357529,14 1245357545,35 1245357582,74 1245357661,30 1245357693,19 1245357714,38 1245357758,11 1245357777,37 1245357817,49 1245357868,19 1245357891,31 1245357931,48 1245357990,49 1245358043,24 1245358082,22 1245358108,17 1245358148,18 1245358168,7 1245358179,6 1245358186,19 1245358209,17 1245358229,5 1245358240,9 1245358252,10 1245358263,6 1245358272,9 1245358296,26 1245358328,49 1245358381,6 1245358389,38 1245358453,19 1245358476,24 1245358504,21 1245358533,76 1245358628,24 1245358653,10 1245358669,105 1245358781,20 1245358808,14 1245358836,6 1245358871,61 1245358933,0 1245358936,44 1245358982,11 1245358996,25 1245359023,15 1245359040,32 1245359076,19 1245359099,13 1245359117,16 1245359138,12 1245359161,33 1245359215,32 1245359249,14 1245359272,7 1245359314,10 1245359333,36 1245359371,21 1245359424,10 1245359447,61 1245359514,32 1245359560,42 1245359604,87 1245359700,60 1245359762,23 1245359786,4 1245359791,8 1245359803,6 1245359813,107 1245359922,29 1245359953,22 1245359978,86 1245360069,75 1245360147,22 1245360170,0 1245360184,41 1245360239,15 1245360256,34 1245360301,37 1245360339,1 1245360342,28 1245360372,20 1245360394,32 1245360440,24 1245360526,3 1245360728,3 1245361011,4 1245361026,35 1245361064,137 1245361359,5 1245362172,11 1245362225,21 1245362248,51 1245362302,20 1245362334,42 1245362418,12 1245362468,7 1245362557,9 1245362817,3 1245363175,4 1245363271,4 1245363446,3 1245363539,4 1245363573,4 1245363635,1 1245363637,3 1245363740,5 1245363875,3 1245364075,4 1245364354,14 1245364370,19 1245364391,49 1245364442,34 1245364478,23 1245364502,80 1245364633,15 1245364650,8 1245364673,16 1245364691,47 1245364739,53 1245364795,39 1245364836,25 1245365353,4 1245365640,11 1245365665,5 1245365726,8 1245365778,7 1245365982,4 1245366017,13 1245366042,6 1245366487,4 1245366493,4 1245366500,4 1245366507,3 1245366622,5 1245366690,5 1245366946,4 1245366953,16 1245366975,8 1245366996,7 1245367005,7 1245367031,6 1245367040,9 1245367051,7 1245367059,23 1245367084,76 1245367166,158 1245367740,4 1245367804,3 1245367847,4 1245367887,9 1245369300,10 1245369611,12 1245370038,10 1245370374,8 1245370668,5 1245370883,5 1245370927,7 1245370945,9 1245370961,16 1245370978,414 1245371398,135 1245371535,252 1245371791,238 1245372034,199 1245372621,4 1245372890,5 1245373043,7 1245373060,9 1245373073,6 1245373081,68 1245373151,10 1245373162,49 1245373212,79 1245373300,12 1245373313,38 1245373353,20 1245373374,59 1245373435,28 1245373465,94 1245373560,11 1245373574,53 1245373629,22 1245373654,6 1245373662,334 1245373998,169 1245374176,41 1245374219,26 1245374246,51 1245374299,31 1245374332,57 1245374391,55 1245374535,4 1245374759,7 1245374769,200 1245374971,215 1245375188,181 1245375371,81 1245375455,59 1245375516,33 1245375552,19 1245375572,56 1245375629,220 1245375850,32 1245375884,26 1245375948,7 1245375964,114 1245376473,4 1245376810,13 1245378296,10 1245378950,12 1245379004,3 1245379569,4 1245379582,4 1245379615,6 1245380030,3 1245380211,4 1245380412,14 1245380727,4 1245380850,4 

This log file is only 7.3 KB. At this rate, a years’ worth of log data can be stored in less than 3MB of plain text files. The data presented here can be graphed (producing the image at the top of the page) using the following code:

#pySquelchGrapher.py
print "loading libraries...",
import pylab, datetime, numpy
print "complete"

def loadData(fname="log.txt"):
	print "loading data...",
	# load signal/duration from log file
	f=open(fname)
	raw=f.read()
	f.close()
	raw=raw.replace('n',' ')
	raw=raw.split(" ")
	signals=[]
	for line in raw:
		if len(line)<3: continue
		line=line.split(',')
		sec=datetime.datetime.fromtimestamp(int(line[0]))
		dur=int(line[1])
		signals.append([sec,dur])
	print "complete"
	return signals

def findDays(signals):
	# determine which days are in the log file
	print "finding days...",
	days=[]
	for signal in signals:
		day = signal[0].date()
		if not day in days:
			days.append(day)
	print "complete"
	return days

def genMins(day):
	# generate an array for every minute in a certain day
	print "generating bins...",
	mins=[]
	startTime=datetime.datetime(day.year,day.month,day.day)
	minute=datetime.timedelta(minutes=1)
	for i in xrange(60*60):
		mins.append(startTime+minute*i)
	print "complete"
	return mins

def fillMins(mins,signals):
	print "filling bins...",
	vals=[0]*len(mins)
	dayToDo=signals[0][0].date()
	for signal in signals:
		if not signal[0].date() == dayToDo: continue
		sec=signal[0]
		dur=signal[1]
		prebuf = sec.second
		minOfDay=sec.hour*60+sec.minute
		if dur+prebuf<60: # simple case, no rollover seconds
			vals[minOfDay]=dur
		else: # if duration exceeds the minute the signal started in
			vals[minOfDay]=60-prebuf
			dur=dur+prebuf
			while (dur>0): # add rollover seconds to subsequent minutes
				minOfDay+=1
				dur=dur-60
				if dur< =0: break
				if dur>=60: vals[minOfDay]=60
				else: vals[minOfDay]=dur
	print "complete"
	return vals

def normalize(vals):
	print "normalizing data...",
	divBy=float(max(vals))
	for i in xrange(len(vals)):
		vals[i]=vals[i]/divBy
	print "complete"
	return vals

def smoothListGaussian(list,degree=10):
	print "smoothing...",
	window=degree*2-1
	weight=numpy.array([1.0]*window)
	weightGauss=[]
	for i in range(window):
		i=i-degree+1
		frac=i/float(window)
		gauss=1/(numpy.exp((4*(frac))**2))
		weightGauss.append(gauss)
	weight=numpy.array(weightGauss)*weight
	smoothed=[0.0]*(len(list)-window)
	for i in range(len(smoothed)):
	  smoothed[i]=sum(numpy.array(list[i:i+window])*weight)/sum(weight)
	while len(list)>len(smoothed)+int(window/2):
		smoothed.insert(0,smoothed[0])
	while len(list)>len(smoothed):
		smoothed.append(smoothed[0])
	print "complete"
	return smoothed

signals=loadData()
days=findDays(signals)
for day in days:
	mins=genMins(day)
	vals=normalize(fillMins(mins,signals))
	fig=pylab.figure()
	pylab.grid(alpha=.2)
	pylab.plot(mins,vals,'k',alpha=.1)
	pylab.plot(mins,smoothListGaussian(vals),'b',lw=1)
	pylab.axis([day,day+datetime.timedelta(days=1),None,None])
	fig.autofmt_xdate()
	pylab.title("147.120 MHz Usage for "+str(day))
	pylab.xlabel("time of day")
	pylab.ylabel("fractional usage")
	pylab.show()




Warning: This post is several years old and the author has marked it as poor quality (compared to more recent posts). It has been left intact for historical reasons, but but its content (and code) may be inaccurate or poorly written.

I’m often drawn toward projects involving data analysis with Python. When I found out a fellow ham in Orlando was using his computer to stream a popular local repeater frequency over the internet I got excited because of the potential for generating data from the setup. Since this guy already has his radio connected to his PC’s microphone jack, I figured I could write a Python app to check the microphone input to determine if anyone is using the frequency. By recording when people start and stop talking, I can create a log of frequency activity. Later I can write software to visualize this data. I’ll talk about that in a later post. For now, here’s how I used Python and a Linux box (Ubuntu, with the python-alsaaudio package installed) to generate such logs.

We can visualize this data using some more simple Python code. Long term it would be useful to visualize frequency activity similarly to how I graphed computer usage at work over the last year but for now since I don’t have any large amount of data to work with. I’ll just write cote to visualize a QSO (conversation) with respect to time. It should be self-explanatory. This data came from data points displayed in the video (provided at the end of this post too).

qsographpng

And, of course, the code I used to generate the log files (seen running in video above): Briefly, this program checks the microphone many times every second to determine if its state has changed (talking/no talking) and records this data in a text file (which it updates every 10 seconds). Matplotlib can EASILY be used to graph data from such a text file.

import alsaaudio, time, audioop, datetime
inp = alsaaudio.PCM(alsaaudio.PCM_CAPTURE,alsaaudio.PCM_NONBLOCK)
inp.setchannels(1)
inp.setrate(4000)
inp.setformat(alsaaudio.PCM_FORMAT_S16_LE)
inp.setperiodsize(1)

squelch = False
lastLog = 0
dataToLog = ""

def logIt(nowSquelch):
 global dataToLog, lastLog
 timeNow = datetime.datetime.now()
 epoch = time.mktime(timeNow.timetuple())
 if nowSquelch==True: nowSquelch=1
 else: nowSquelch=0
 logLine="%s %dn"%(timeNow, nowSquelch)
 print timeNow, nowSquelch
 dataToLog+=logLine
 if epoch-lastLog>10:
 #print "LOGGING..."
 f=open('squelch.txt','a')
 f.write(dataToLog)
 f.close()
 lastLog = epoch
 dataToLog=""

while True:
 l,data = inp.read()
 if l:
 vol = audioop.max(data,2)
 #print vol #USED FOR CALIBRATION
 if vol>800: nowSquelch = True
 else: nowSquelch = False
 if not nowSquelch == squelch:
 logIt(nowSquelch)
 squelch = nowSquelch
 time.sleep(.01)

To use this code make sure that you’ve properly calibrated it. See the “vol>800” line? That means that if the volume in the microphone is at least 800, it’s counted as talking, and less than it’s silence. Hopefully you can find a value that counts as silence when the squelch is active, but as talking when the squelch is broken (even if there’s silence). This is probably best achieved with the radio outputting at maximum volume. You’ll have to run the program live with that line un-commented to view the data values live. Find which values occur for squelch on/off, and pick your threshold accordingly.

After that you can visualize the data with the following code. Note that this is SEVERELY LIMITED and is only useful when graphing a few minutes of data. I don’t have hours/days of data to work with right now, so I won’t bother writing code to graph it. This code produced the graph seen earlier in this page. Make sure matplotlib is installed on your box.

import pylab

def loadData():
 #returns Xs
 import time, datetime, pylab
 f=open('good.txt')
 raw=f.readlines()
 f.close()
 onTimes=[]
 timeStart=None
 lastOn=False
 for line in raw:
 if len(line)<10: continue
 line = line.strip('n').split(" ")
 t=line[0]+" "+line[1]
 t=t.split('.')
 thisDay=time.strptime(t[0], "%Y-%m-%d %H:%M:%S")
 e=time.mktime(thisDay)+float("."+t[1])
 if timeStart==None: timeStart=e
 if line[-1]==1: stat=True
 else: stat=False
 if not lastOn and line[-1]=="1":
 lastOn=e
 else:
 onTimes.append([(lastOn-timeStart)/60.0,
 (e-timeStart)/60.0])
 lastOn=False
 return onTimes

times = loadData()
pylab.figure(figsize=(8,3))
for t in times:
 pylab.fill([t[0],t[0],t[1],t[1]],[0,1,1,0],'k',lw=0,alpha=.5)
pylab.axis([None,None,-3,4])
pylab.title("A little QSO")
pylab.xlabel("Time (minutes)")
pylab.show()




Warning: This post is several years old and the author has marked it as poor quality (compared to more recent posts). It has been left intact for historical reasons, but but its content (and code) may be inaccurate or poorly written.

I enjoy writing Python scripts to analyze and display linear data. One of my favorite blog entries is Linear Data Smoothing with Python, developed for my homemade electrocardiogram project. I installed a program called TimeTrack.exe on my work computer. It basically logs whenever you open or close a program. The data output looks like this:

"Firefox","Prototype of a Digital Biopsy Device - Mozilla Firefox","05/19/2009  9:45a","05/19/2009  9:45a","766ms","0.0"
"Firefox","Dual-Channel Mobile Surface Electromyograph - Mozilla Firefox","05/19/2009  9:46a","05/19/2009  9:46a","797ms","0.0"
"Windows Explorer","","03/24/2008  9:30a","05/19/2009  9:48a","49d 6h 9m","20.7"
"Windows Explorer","09_04_07_RA_SA_AV","05/19/2009  8:48a","05/19/2009  8:48a","1.0s","0.0"
"Windows Explorer","Image003.jpg - Windows Picture and Fax Viewer","05/18/2009  4:03p","05/18/2009  4:03p","1.2s","0.0"

I have a 13 MB file containing lines like this which I parse, condense, analyze, and display with Python. The script finds the first and last entry time and creates a dictionary where keys are the hours between the 1st and last log lines, parses the log, determines which time block each entry belongs to, and increments the integer (value of the dictionary) for its respective key. Something similar is repeated, but with respect to days rather than hours. The result is:

compusage_white

The code I used to generate this graph is:

# This script analyzes data exported from "TimeTrack" (a free computer usage
# monitoring program for windows) and graphs the data visually.
import time, pylab, datetime, numpy

# This is my computer usage data.  Generate yours however you want.
allHours = ['2008_10_29 0', '2009_03_11 5', '2009_04_09 5', '2008_07_04 10',
'2008_12_18 9', '2009_01_30 12', '2008_09_04 7', '2008_05_17 1',
'2008_05_11 5', '2008_11_03 3', '2008_05_21 3', '2009_02_19 11',
'2008_08_15 13', '2008_04_02 4', '2008_07_16 5', '2008_09_16 8',
'2008_04_10 5', '2009_05_10 1', '2008_12_30 4', '2008_06_07 2',
'2008_11_23 0', '2008_08_03 0', '2008_04_30 4', '2008_07_28 9',
'2008_05_19 0', '2009_03_30 7', '2008_06_19 3', '2009_01_24 3',
'2008_08_23 6', '2008_12_01 0', '2009_02_23 6', '2008_11_27 0',
'2008_05_02 5', '2008_10_20 13', '2008_03_27 5', '2009_04_02 9',
'2009_02_21 0', '2008_09_13 1', '2008_12_13 0', '2009_04_14 11',
'2009_01_31 7', '2008_11_04 10', '2008_07_09 6', '2008_10_24 10',
'2009_02_22 0', '2008_09_25 12', '2008_12_25 0', '2008_05_26 4',
'2009_05_01 10', '2009_04_26 11', '2008_08_10 8', '2008_11_08 6',
'2008_07_21 12', '2009_04_21 3', '2009_05_13 8', '2009_02_02 8',
'2008_10_07 2', '2008_06_10 6', '2008_09_21 0', '2009_03_17 9',
'2008_08_30 7', '2008_11_28 4', '2009_02_14 0', '2009_01_22 6',
'2008_10_11 0', '2008_06_22 8', '2008_12_04 0', '2008_03_28 0',
'2009_04_07 2', '2008_09_10 0', '2008_05_15 5', '2008_08_18 12',
'2008_10_31 5', '2009_03_09 7', '2009_02_25 8', '2008_07_02 4',
'2008_12_16 7', '2008_09_06 2', '2009_01_26 5', '2009_04_19 0',
'2008_07_14 13', '2008_11_01 5', '2009_01_18 0', '2009_05_04 0',
'2008_08_13 10', '2009_02_27 3', '2009_01_16 12', '2008_09_18 8',
'2009_02_03 7', '2008_06_01 0', '2008_12_28 0', '2008_07_26 0',
'2008_11_21 1', '2008_08_01 8', '2008_04_28 3', '2009_05_16 0',
'2008_06_13 5', '2008_10_02 11', '2009_03_28 6', '2008_08_21 7',
'2009_01_13 6', '2008_11_25 4', '2008_06_25 1', '2008_10_22 11',
'2008_03_25 6', '2009_02_07 6', '2008_12_11 4', '2009_01_01 4',
'2008_09_15 2', '2009_02_05 12', '2008_07_07 9', '2009_04_12 0',
'2008_04_11 5', '2008_10_26 4', '2008_05_28 3', '2008_09_27 14',
'2009_05_03 0', '2008_12_23 5', '2009_05_12 10', '2008_11_14 3',
'2008_07_19 0', '2009_04_24 8', '2008_04_07 1', '2008_08_08 11',
'2008_06_04 0', '2009_05_15 12', '2009_03_23 13', '2009_02_01 10',
'2008_09_23 11', '2009_02_08 3', '2008_08_28 4', '2008_11_18 9',
'2008_07_31 7', '2008_10_13 0', '2008_06_16 9', '2009_03_27 6',
'2008_12_02 0', '2008_05_01 7', '2009_04_05 1', '2008_08_16 9',
'2009_03_15 0', '2008_04_16 6', '2008_10_17 4', '2008_06_28 5',
'2009_01_28 10', '2008_04_18 0', '2008_12_14 0', '2008_11_07 6',
'2009_04_17 7', '2008_04_14 7', '2008_07_12 0', '2009_01_15 7',
'2009_05_06 8', '2008_12_26 0', '2008_06_03 7', '2008_09_28 0',
'2008_05_25 4', '2008_08_07 8', '2008_04_26 7', '2008_07_24 1',
'2008_04_20 0', '2008_11_11 4', '2009_04_29 0', '2008_10_04 0',
'2009_05_18 9', '2009_03_18 4', '2008_06_15 8', '2009_02_13 6',
'2008_05_04 5', '2009_03_04 2', '2009_03_06 3', '2008_05_06 0',
'2008_08_27 11', '2008_04_22 0', '2009_03_26 6', '2008_03_31 9',
'2008_06_27 5', '2008_10_08 4', '2008_09_09 4', '2008_12_09 3',
'2008_05_10 0', '2008_05_14 5', '2009_04_10 0', '2009_01_11 0',
'2008_07_05 8', '2009_01_05 7', '2008_10_28 0', '2009_02_18 11',
'2009_03_10 7', '2008_05_30 3', '2008_09_05 7', '2008_12_21 6',
'2009_03_02 6', '2008_08_14 5', '2008_11_12 5', '2008_07_17 8',
'2008_04_05 6', '2009_04_22 11', '2009_05_09 0', '2008_06_06 0',
'2009_01_03 0', '2008_09_17 6', '2009_03_21 3', '2009_02_10 7',
'2008_05_08 4', '2008_08_02 0', '2008_11_16 0', '2008_07_29 12',
'2008_10_15 5', '2008_06_18 5', '2009_03_25 2', '2009_01_10 0',
'2009_04_03 5', '2008_08_22 7', '2009_03_13 11', '2008_10_19 0',
'2008_06_30 8', '2008_09_02 9', '2008_05_23 4', '2008_12_12 7',
'2008_07_10 11', '2008_11_05 8', '2008_04_12 4', '2009_04_15 7',
'2008_12_24 1', '2008_09_30 0', '2008_05_27 2', '2008_08_05 10',
'2008_04_24 6', '2009_04_27 6', '2008_07_22 3', '2008_11_09 1',
'2008_06_09 6', '2008_10_06 14', '2009_03_16 7', '2008_05_22 5',
'2009_01_29 12', '2008_11_29 4', '2008_04_09 7', '2008_08_25 12',
'2009_02_15 0', '2008_03_29 7', '2008_06_21 7', '2008_10_10 9',
'2008_05_12 6', '2009_02_16 10', '2008_09_11 11', '2008_12_07 0',
'2008_07_03 6', '2009_04_08 3', '2009_01_23 7', '2009_01_27 5',
'2008_10_30 0', '2009_03_08 0', '2009_01_21 8', '2008_12_19 0',
'2008_05_16 2', '2009_01_25 1', '2009_02_26 5', '2008_09_07 2',
'2008_04_03 1', '2008_08_12 6', '2008_04_13 10', '2008_11_02 0',
'2008_07_15 0', '2009_04_20 3', '2009_02_24 10', '2009_05_11 8',
'2008_12_31 8', '2008_04_15 7', '2008_09_19 10', '2009_01_19 0',
'2008_11_22 3', '2008_07_27 2', '2009_02_04 7', '2009_03_31 1',
'2008_05_24 3', '2008_10_01 8', '2008_06_12 6', '2009_01_12 11',
'2008_11_26 8', '2009_04_01 10', '2009_02_28 0', '2008_08_20 6',
'2008_10_21 10', '2008_06_24 4', '2008_03_26 4', '2008_12_10 0',
'2008_09_12 0', '2008_05_09 7', '2009_02_17 7', '2008_07_08 6',
'2008_10_25 5', '2009_04_13 9', '2009_05_02 0', '2008_12_22 8',
'2008_09_24 9', '2009_01_20 5', '2008_11_15 6', '2009_04_25 10',
'2008_08_11 9', '2008_04_06 8', '2008_07_20 1', '2009_03_22 3',
'2008_06_11 6', '2008_09_20 3', '2009_05_14 10', '2008_11_19 0',
'2008_08_31 2', '2009_02_09 8', '2008_10_12 0', '2008_04_25 5',
'2008_06_23 4', '2009_01_07 8', '2008_08_19 0', '2008_12_05 2',
'2008_07_01 8', '2008_10_16 6', '2009_04_06 3', '2009_03_14 5',
'2008_09_01 2', '2008_12_17 14', '2008_05_18 7', '2008_04_01 2',
'2009_04_18 0', '2008_04_17 0', '2008_07_13 0', '2008_06_02 10',
'2008_09_29 6', '2008_12_29 0', '2009_05_05 8', '2008_04_19 0',
'2009_04_30 8', '2008_08_06 4', '2008_11_20 0', '2008_07_25 6',
'2009_02_06 6', '2009_03_29 3', '2009_05_17 0', '2009_03_19 7',
'2008_10_03 1', '2008_06_14 3', '2008_05_07 5', '2008_08_26 3',
'2008_11_24 9', '2008_04_21 8', '2008_04_23 4', '2008_10_23 11',
'2008_06_26 4', '2008_03_24 8', '2008_12_08 5', '2008_09_14 2',
'2009_01_02 6', '2008_04_08 0', '2008_10_27 6', '2009_04_11 0',
'2008_07_06 0', '2008_12_20 3', '2009_04_23 6', '2008_09_26 9',
'2008_05_31 0', '2008_07_18 4', '2008_11_13 6', '2008_08_09 2',
'2008_04_04 0', '2009_03_20 5', '2008_09_22 7', '2009_05_08 9',
'2008_06_05 7', '2008_07_30 7', '2008_11_17 10', '2008_05_03 0',
'2008_08_29 3', '2009_02_11 12', '2009_01_08 8', '2008_06_17 0',
'2008_10_14 7', '2009_03_24 11', '2008_08_17 6', '2008_12_03 0',
'2009_01_09 4', '2008_05_29 5', '2008_06_29 9', '2008_10_18 5',
'2009_04_04 0', '2008_12_15 10', '2009_03_12 0', '2009_03_05 7',
'2008_05_20 4', '2008_09_03 7', '2009_03_07 8', '2009_01_14 6',
'2008_05_05 5', '2008_11_06 7', '2008_07_11 6', '2009_04_16 9',
'2009_02_20 0', '2008_12_27 0', '2009_01_17 0', '2009_05_07 7',
'2008_11_10 5', '2008_07_23 11', '2009_04_28 0', '2008_04_27 2',
'2008_08_04 0', '2009_03_01 11', '2008_10_05 0', '2008_06_08 8',
'2009_05_19 5', '2008_04_29 4', '2008_11_30 0', '2009_01_06 8',
'2009_02_12 3', '2008_08_24 2', '2009_03_03 10', '2008_10_09 6',
'2008_06_20 2', '2008_05_13 10', '2008_12_06 0', '2008_03_30 7']

def genTimes():
    ## opens  exported timetrack data (CSV) and re-saves a compressed version.
    print "ANALYZING..."
    f=open('timetrack.txt')
    raw=f.readlines()
    f.close()
    times=["05/15/2009 12:00am"] #start time
    for line in raw[1:]:
        if not line.count('","') == 5: continue
        test = line.strip("n")[1:-1].split('","')[-3].replace("  "," ")+"m"
        test = test.replace(" 0:"," 12:")
        times.append(test) #end time
        test = line.strip("n")[1:-1].split('","')[-4].replace("  "," ")+"m"
        test = test.replace(" 0:"," 12:")
        times.append(test) #start time

    times.sort()
    print "WRITING..."
    f=open('times.txt','w')
    f.write(str(times))
    f.close()

def loadTimes():
    ## loads the times from the compressed file.
    f=open("times.txt")
    times = eval(f.read())
    newtimes=[]
    f.close()
    for i in range(len(times)):
        if "s" in times[i]: print times[i]
        newtimes.append(datetime.datetime(*time.strptime(times[i],
                                        "%m/%d/%Y %I:%M%p")[0:5]))
        #if i>1000: break #for debugging
    newtimes.sort()
    return newtimes

def linearize(times):
    ## does all the big math to calculate hours per day.
    for i in range(len(times)):
        times[i]=times[i]-datetime.timedelta(minutes=times[i].minute,
                                             seconds=times[i].second)
    hr = datetime.timedelta(hours=1)
    pos = times[0]-hr
    counts = {}
    days = {}
    lasthr=pos
    lastday=None
    while pos1:counts[pos]=1 #flatten
        if not daypos in days: days[daypos]=0
        if not lasthr == pos:
            if counts[pos]>0:
                days[daypos]=days[daypos]+1
                lasthr=pos
        pos+=hr
    return days #[counts,days]

def genHours(days):
    ## outputs the hours per day as a file.
    out=""
    for day in days:
        print day
        out+="%s %in"%(day.strftime("%Y_%m_%d"),days[day])
    f=open('hours.txt','w')
    f.write(out)
    f.close()
    return

def smoothListGaussian(list,degree=7):
    ## (from an article I wrote) - Google "linear data smoothing with python".
    firstlen=len(list)
    window=degree*2-1
    weight=numpy.array([1.0]*window)
    weightGauss=[]
    for i in range(window):
     i=i-degree+1
     frac=i/float(window)
     gauss=1/(numpy.exp((4*(frac))**2))
     weightGauss.append(gauss)
    weight=numpy.array(weightGauss)*weight
    smoothed=[0.0]*(len(list)-window)
    for i in range(len(smoothed)):
     smoothed[i]=sum(numpy.array(list[i:i+window])*weight)/sum(weight)
    pad_before = [smoothed[0]]*((firstlen-len(smoothed))/2)
    pad_after  = [smoothed[-1]]*((firstlen-len(smoothed))/2+1)
    return pad_before+smoothed+pad_after

### IF YOU USE MY DATA, YOU ONLY USE THE FOLLOWING CODE ###

def graphIt():
    ## Graph the data!
    #f=open('hours.txt')
    #data=f.readlines()
    data=allHours
    data.sort()
    f.close()
    days,hours=[],[]
    for i in range(len(data)):
        day = data[i].split(" ")
        if int(day[1])<4: continue
        days.append(datetime.datetime.strptime(day[0], "%Y_%m_%d"))
        hours.append(int(day[1]))
    fig=pylab.figure(figsize=(14,5))
    pylab.plot(days,smoothListGaussian(hours,1),'.',color='.5',label="single day")
    pylab.plot(days,smoothListGaussian(hours,1),'-',color='.8')
    pylab.plot(days,smoothListGaussian(hours,7),color='b',label="7-day gausian average")
    pylab.axhline(8,color='k',ls=":")
    pylab.title("Computer Usage at Work")
    pylab.ylabel("hours (rounded)")
    pylab.legend()
    pylab.show()
    return

#times = genTimes()
#genHours(linearize(loadTimes()))
graphIt()