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.

Last night I finished building my DIY ECG as a prototype (I finally got the circuit off the breadboard and onto a plastic sheet). This is a similar circuit to the one used to record data from the last entry (resister values are now identical to the crude circuit described several posts ago). I left-in the crude band-pass filter (made by grounding my primary electrode sensor through a 0.1µF capacitor) because it seemed to help a great deal, and wasn’t hard to implement. I picked up all of my parts (including the LM324 quad op-amp microchip) at RadioShack. Of note, the quad-op-amp is overkill because I’m only using one of the 4 op-amps. Theoretically I could add 3 more electrodes to this circuit (which would allow for multi-sensor recording) but this would require multiple microphone jacks, which isn’t very common. I guess I could use 2 microphone jacks, and differentiate right/left channels.

I made the prototype by drilling holes in a small rectangular piece of a non-conductive plastic material. I picked up a stack of these rectangular sections for a quarter at a local electrical surplus store and they’re perfect for prototyping. The two green wires coming out the left side attach to a power supply (either a plugged in AC->DC transformer, 3 or 4 AA batteries, or even a 9V should work). The blue wires on the right attach to the electrodes I put on my chest. The black wires go to a headphone-jack which I plug into the microphone hole of my PC to record the signal.

This is the back of the device which shows my crummy soldering. I’m a molecular biologist not an electrical engineer. The white/yellow wires correspond to the left/right channels of the microphone connector. I only use the left one (white), but attached the right channel (yellow) to the op-amp just in case I decide to add another sensor later – this is not required.

Here’s the full device: You can see the circuit (note its small size – easy to mount inside of a tictac box or something) with the green wires leading to a power supply, black cable being the microphone connector, and the blue wires leading to electrodes made… from… Fanta… cans…? Yes, in the spirit of rigging electronics (my specialty) I found that surprisingly nice chest electrodes can be made from aluminum soda cans! If you go this route, cut them delicately so you don’t get metal shards in your skin like I did at first. Also, note that you have to firmly scrape each side of the aluminum to get the insulating waxy-plastic stuff off or it just won’t work. I guess it’s coated with something to prevent the soda from tasting like metal. Aluminum rapidly transfers heat and it’s nearly impossible to solder leads onto these pads, so I wrapped a wire (tipped with a bead of solder) with the edge of the aluminum several times and crushed the heck out of it with pliers and it seems to stay on well and make a good connection. Also, before taping these onto your skin, it helps to put a conductive goo on it to make the connection better. I added skin moisturizer to each electrode and taped the gooey electrode directly onto my chest.

I recorded ~20 minutes of data last night with this rig and it looked pretty good. I went to analyze it with Python and it kept crashing! The python script I gave you previously loads the entire wave file into an array of numbers, but with a 20 minute wave file (at over 40,000 samples a second) it is too big for memory. I wrote an updated wave loader which loads large wave files in parts which is much more efficient. It also performs the condensation method at load time. Basically, it loads 100 data points (or whatever deg is set to), averages them, and adds this value to a list. The result is a smoothed trace with a resolution of 400 Hz instead of 40 kHz. I’d test this on the wave file I recorded last night but that’s on my laptop which is in the car and I’ve got to get back to work. Here’s that function:

 def loadWav(fname,deg=100):
     global hz
     w=wave.open(fname)
     nchannel, width, rate, length, comptype, compname = w.getparams()
     print "[%s]
 rate: %d Hz
 frames: %d
 length: %.02f sec" %
           (fname, rate, length, float(length)/rate)
     hz=rate/deg
     chunks=int(length/deg)
     data=[]
     for i in range(chunks):
         if i%7777==0:
             print "processing chunk %d of %d (%0.02f%%)" %
                   (i,chunks,100.0*i/chunks)
         vals = struct.unpack("%sh" %deg,w.readframes(deg))
         data.append(sum(vals)/float(deg))
     print "complete!"
     return data




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.

Would I rather design circuits or software? I’m a software guy (or at lesat I know more about software than circuits) so I’d rather record noisy signals and write software to eliminate the noise, rather than assembling circuits to eliminate the noise in hardware. In the case of my DIY ECG machine, I’d say I’ve done a surprisingly good job of eliminating noise using software. Most DIY ECG circuits on the net use multiple op-amps and filters to do this. Instead of all that fancy stuff, I made a crude circuit (a single op-amp and two resisters) that is capable of record my ECG and filtered it in software. The output is pretty good!

The first step in removing noise is understanding it. Most of the noise in my signal came from sine waves caused by my electrodes picking up radiated signals in the room. Since this type of interference is consistent through the entire recording, power-spectral analysis could be applied to determine the frequencies of the noise so I could selectively block them out. I used the fast Fourier transform algorithm (FFT) on the values to generate a plot of the spectral components of my signal (mostly noise) seen as sharp peaks. I manually band-stopped certain regions of the spectrum that I thought were noise-related (colored bands). This is possible to do electronically with a more complicated circuit, but is interesting to do in software. I think performed an inverse FFT on the trace. The result was a trace with greatly reduced noise. After a moving window smoothing algorithm was applied the signal was even better! Note that I recorded the WAV file with “sound recorder” (not GoldWave) and did all of the processing (including band-pass filtering) within Python.

The ECG came out better than expected! The graph above shows the power spectral analysis with band-stop filters applied at the colored regions. Below is the trace of the original signal (light gray), the inverse-FFT-filtered trace (dark gray), and the smoothed filtered trace (black) – the final ECG signal I intend to use.

This is a magnified view of a few heartbeats. It looks pretty good! Here’s the code I used to do all the calculations:

 import wave, struct, numpy, pylab, scipy

 fname='./success3.wav'

 def readwave(wavfilename):
     """load raw data directly from a WAV file."""
     global rate
     w=wave.open(wavfilename,'rb')
     (nchannel, width, rate, length, comptype, compname) = w.getparams()
     print "[%s] %d HZ (%0.2fsec)" %(wavfilename, rate, length/float(rate))
     frames = w.readframes(length)
     return numpy.array(struct.unpack("%sh" %length*nchannel,frames))

 def shrink(data,deg=100):
     """condense a linear data array by a multiple of [deg]."""
     global rate
     small=[]
     print "starting with", len(data)
     for i in range(len(data)/deg):
         small.append(numpy.average(data[i*deg:(i+1)*deg]))
     print "ending with", len(small)
     rate = rate/deg
     #return small[40000:50000]
     return small

 def normalize(data):
     """make all data fit between -.5 and +.5"""
     data=data-numpy.average(data)
     big=float(max(data))
     sml=float(min(data))
     data=data/abs(big-sml)
     data=data+float(abs(min(data)))-.47
     return data

 def smooth(data,deg=20,expand=False):
     """moving window average (deg = window size)."""
     for i in range(len(data)-deg):
         if i==0: cur,smooth=sum(data[0:deg]),[]
         smooth.append(cur/deg)
         cur=cur-data[i]+data[i+deg]
     if expand:
         for i in range(deg):
             smooth.append(smooth[-1])
     return smooth

 def smoothListGaussian(list,degree=10,expand=False):
     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)
     if expand:
         for i in range((degree*2)-1):
             smoothed.append(smoothed[-1])
     return smoothed

 def goodSmooth(data):
     #data=smooth(fix,20,True)
     data=smooth(fix,100,True)
     #data=smooth(fix,20,True)
     return data

 def makeabs(data):
     """center linear data to its average value."""
     for i in range(len(data)): data[i]=abs(data[i])
     return data

 def invert(data):
     """obviously."""
     for i in range(len(data)): data[i]=-data[i]
     return data

 def loadwav(fname):
     """a do-everything function to get usable, smoothed data from a WAV."""
     wav=readwave(fname)
     wav=shrink(wav)
     wav=invert(wav)
     wav=smooth(wav)
     wav=smooth(wav,10)
     wav=normalize(wav)
     Xs=getXs(wav)
     return Xs,wav

 def getXs(datalen):
     """calculate time positions based on WAV frequency resolution."""
     Xs=[]
     for i in range(len(datalen)):
         Xs.append(i*(1/float(rate)))
     print len(datalen), len(Xs)
     return Xs

 def integrate(data):
     """integrate the function with respect to its order."""
     inte=[]
     for i in range(len(data)-1):
         inte.append(abs(data[i]-data[i+1]))
     inte.append(inte[-1])
     return inte

 def getPoints(Xs,data,res=10):
     """return X,Y coordinates of R peaks and calculate R-R based heartrate."""
     pXs,pYs,pHRs=[],[],[]
     for i in range(res,len(data)-res):
         if data[i]>data[i-res]+.1 and data[i]>data[i+res]+.1:
             if data[i]>data[i-1] and data[i]>data[i+1]:
                 pXs.append(Xs[i])
                 pYs.append(data[i])
                 if len(pXs)>1:
                     pHRs.append((1.0/(pXs[-1]-pXs[-2]))*60.0)
     pHRs.append(pHRs[-1])
     return pXs,pYs,pHRs

 def bandStop(fft,fftx,low,high,show=True):
     lbl="%d-%d"%(low,high)
     print "band-stopping:",lbl
     if show:
         col=pylab.cm.spectral(low/1200.)
         pylab.axvspan(low,high,alpha=.4,ec='none',label=lbl,fc=col)
         #pylab.axvspan(-low,-high,fc='r',alpha=.3)
     for i in range(len(fft)):
         if abs(fftx[i])>low and abs(fftx[i])<high :
             fft[i]=0
     return fft

 def getXs(data):
     xs=numpy.array(range(len(data)))
     xs=xs*(1.0/rate)
     return xs

 def clip(x,deg=1000):
     return numpy.array(x[deg:-deg])

 pylab.figure(figsize=(12,8))
 raw = invert(shrink(readwave(fname),10))
 xs = getXs(raw)
 fftr = numpy.fft.fft(raw)
 fft = fftr[:]
 fftx= numpy.fft.fftfreq(len(raw), d=(1.0/(rate)))

 pylab.subplot(2,1,1)
 pylab.plot(fftx,abs(fftr),'k')

 fft=bandStop(fft,fftx,30,123)
 fft=bandStop(fft,fftx,160,184)
 fft=bandStop(fft,fftx,294,303)
 fft=bandStop(fft,fftx,386,423)
 fft=bandStop(fft,fftx,534,539)
 fft=bandStop(fft,fftx,585,610)
 fft=bandStop(fft,fftx,654,660)
 fft=bandStop(fft,fftx,773,778)
 fft=bandStop(fft,fftx,893,900)
 fft=bandStop(fft,fftx,1100,max(fftx))
 pylab.axis([0,1200,0,2*10**6])
 pylab.legend()
 pylab.title("Power Spectral Analysis",fontsize=28)
 pylab.ylabel("Power",fontsize=20)
 pylab.xlabel("Frequency (Hz)",fontsize=20)

 pylab.subplot(2,1,2)
 pylab.title("Original Trace",fontsize=28)
 pylab.ylabel("Potential",fontsize=20)
 pylab.xlabel("Time (sec)",fontsize=20)
 pylab.plot(clip(xs),clip(raw),color='.8',label='1: raw')

 fix = scipy.ifft(fft)
 pylab.plot(clip(xs),clip(fix)+5000,color='.6',label='2: band-stop')
 pylab.plot(clip(xs),clip(goodSmooth(fix))-5000,'k',label='3: smoothed')
 pylab.legend()
 pylab.title("Band-Stop Filtered Trace",fontsize=28)
 pylab.ylabel("Potential",fontsize=20)
 pylab.xlabel("Time (sec)",fontsize=20)

 pylab.subplots_adjust(hspace=.5)
 pylab.savefig('out.png',dpi=100)
 pylab.show()
 print "COMPLETE"




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.

I’ve succeeded in building my own electrocardiograph (ECG) to record the electrical activity of my own heart! Briefly, I built a micropotential amplifier using an op-amp and attached it to makeshift electrodes on my chest (pennies and shampoo), fed the amplified signal into my sound card, and recorded it as a WAV. The signal is very noisy though. I was able to do a great job at removing this noise using band/frequency filters in GoldWave (audio editing software designed to handle WAV files). I band-blocked 50-70 Hz (which removed the oscillations from the 60 Hz AC lines running around my apartment). I then wrote the Python code (at the bottom of this entry) to load this WAV file as a single list of numbers (voltage potentials). I performed a data condensation algorithm (converting 100 points of raw WAV data into a single, averaged point, lessening my processing load by 100x), followed by two consecutive moving window averages (20-point window, performed on the condensed data). The result was a voltage reading that had most of the noise removed and a beautiful ECG signal emerged! I also tossed in some code to determine the peak of the R wave, label it (yellow dot), and use the inverse R-R time distance to calculate heart rate.

This is my actual ECC signal as record by a circuit similar to the one in the previous entry, recorded through my sound card, and processed with the Python script below. You can start to see the Q, R, S, and T components. I can’t wait to solder-up a prototype (it’s currently breadboarded) and try to analyze hours of data rather than just a few seconds. I’ll take pictures of this device soon.

And here’s the code I used: note that it relies on the WAV file I recorded. This code has extra functions not required to produce the image above, but I left them in in case they may be useful.

 import wave, struct, numpy, pylab, scipy

 def readwave(wavfilename):
     """load raw data directly from a WAV file."""
     global rate
     w=wave.open(wavfilename,'rb')
     (nchannel, width, rate, length, comptype, compname) = w.getparams()
     print "[%s] %d HZ (%0.2fsec)" %(wavfilename, rate, length/float(rate))
     frames = w.readframes(length)
     return numpy.array(struct.unpack("%sh" %length*nchannel,frames))

 def shrink(data,deg=100):
     """condense a linear data array by a multiple of [deg]."""
     global rate
     small=[]
     print "starting with", len(data)
     for i in range(len(data)/deg):
         small.append(numpy.average(data[i*deg:(i+1)*deg]))
     print "ending with", len(small)
     rate = rate/deg
     return small

 def normalize(data):
     """make all data fit between -.5 and +.5"""
     data=data-numpy.average(data)
     big=float(max(data))
     sml=float(min(data))
     data=data/abs(big-sml)
     data=data+float(abs(min(data)))-.47
     return data

 def smooth(data,deg=20):
     """moving window average (deg = window size)."""
     for i in range(len(data)-deg):
         if i==0: cur,smooth=sum(data[0:deg]),[]
         smooth.append(cur/deg)
         cur=cur-data[i]+data[i+deg]
     return smooth

 def makeabs(data):
     """center linear data to its average value."""
     for i in range(len(data)): data[i]=abs(data[i])
     return data

 def invert(data):
     """obviously."""
     for i in range(len(data)): data[i]=-data[i]
     return data

 def loadwav(fname='./wavs/bandpassed.wav'):
     """a do-everything function to get usable, smoothed data from a WAV."""
     wav=readwave(fname)
     wav=shrink(wav)
     wav=invert(wav)
     wav=smooth(wav)
     wav=smooth(wav,10)
     wav=normalize(wav)
     Xs=getXs(wav)
     return Xs,wav

 def getXs(datalen):
     """calculate time positions based on WAV frequency resolution."""
     Xs=[]
     for i in range(len(datalen)):
         Xs.append(i*(1/float(rate)))
     print len(datalen), len(Xs)
     return Xs

 def integrate(data):
     """integrate the function with respect to its order."""
     inte=[]
     for i in range(len(data)-1):
         inte.append(abs(data[i]-data[i+1]))
     inte.append(inte[-1])
     return inte

 def getPoints(Xs,data,res=10):
     """return X,Y coordinates of R peaks and calculate R-R based heartrate."""
     pXs,pYs,pHRs=[],[],[]
     for i in range(res,len(data)-res):
         if data[i]>data[i-res]+.1 and data[i]>data[i+res]+.1:
             if data[i]>data[i-1] and data[i]>data[i+1]:
                 pXs.append(Xs[i])
                 pYs.append(data[i])
                 if len(pXs)>1:
                     pHRs.append((1.0/(pXs[-1]-pXs[-2]))*60.0)
     pHRs.append(pHRs[-1])
     return pXs,pYs,pHRs

 Xs,Ys=loadwav()
 px,py,pHR=getPoints(Xs,Ys)

 pylab.figure(figsize=(12,6))
 pylab.subplot(2,1,1)
 #pylab.axhline(color='.4',linestyle=':')
 pylab.plot(Xs,Ys,'b-')
 pylab.plot(px,py,'y.')
 pylab.axis([None,None,-.6,.6])
 pylab.title("DIY Electrocardiogram - Trial 1",fontsize=26)
 pylab.ylabel("Normalized Potential",fontsize=16)
 #pylab.xlabel("Time (sec)")
 ax=pylab.axis()
 pylab.subplot(2,1,2)
 pylab.plot(px,pHR,'k:')
 pylab.plot(px,pHR,'b.')
 pylab.axis([ax[0],ax[1],None,None])
 pylab.ylabel("Heart Rate (BPM)",fontsize=16)
 pylab.xlabel("Time (seconds)",fontsize=16)
 pylab.savefig("test.png",dpi=120)
 pylab.show()
 print "DONE"




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 discrete functions to perform data smoothing in python. I found this Gaussian smoothing function to be most useful:

def smoothListGaussian(list, degree=5):
    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)
    return smoothed 

Provide a list and it will return a smoother version of the data. The Gaussian smoothing function I wrote is leagues better than a moving window average method, for reasons that are obvious when viewing the chart below. Surprisingly, the moving triangle method appears to be very similar to the Gaussian function at low degrees of spread. However, for huge numbers of data points, the Gaussian function should perform better.

import pylab
import numpy

def smoothList(list, strippedXs=False, degree=10):
    if strippedXs == True:
        return Xs[0:-(len(list)-(len(list)-degree+1))]
    smoothed = [0]*(len(list)-degree+1)
    for i in range(len(smoothed)):
        smoothed[i] = sum(list[i:i+degree])/float(degree)
    return smoothed

def smoothListTriangle(list, strippedXs=False, degree=5):
    weight = []
    window = degree*2-1
    smoothed = [0.0]*(len(list)-window)
    for x in range(1, 2*degree):
        weight.append(degree-abs(degree-x))
    w = numpy.array(weight)
    for i in range(len(smoothed)):
        smoothed[i] = sum(numpy.array(list[i:i+window])*w)/float(sum(w))
    return smoothed

def smoothListGaussian(list, strippedXs=False, degree=5):
    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)
    return smoothed

### DUMMY DATA ###
data = [0]*30  # 30 "0"s in a row
data[15] = 1  # the middle one is "1"

### PLOT DIFFERENT SMOOTHING FUNCTIONS ###
pylab.figure(figsize=(550/80, 700/80))
pylab.suptitle('1D Data Smoothing', fontsize=16)
pylab.subplot(4, 1, 1)
p1 = pylab.plot(data, ".k")
p1 = pylab.plot(data, "-k")
a = pylab.axis()
pylab.axis([a[0], a[1], -.1, 1.1])
pylab.text(2, .8, "raw data", fontsize=14)
pylab.subplot(4, 1, 2)
p1 = pylab.plot(smoothList(data), ".k")
p1 = pylab.plot(smoothList(data), "-k")
a = pylab.axis()
pylab.axis([a[0], a[1], -.1, .4])
pylab.text(2, .3, "moving window average", fontsize=14)
pylab.subplot(4, 1, 3)
p1 = pylab.plot(smoothListTriangle(data), ".k")
p1 = pylab.plot(smoothListTriangle(data), "-k")
pylab.axis([a[0], a[1], -.1, .4])
pylab.text(2, .3, "moving triangle", fontsize=14)
pylab.subplot(4, 1, 4)
p1 = pylab.plot(smoothListGaussian(data), ".k")
p1 = pylab.plot(smoothListGaussian(data), "-k")
pylab.axis([a[0], a[1], -.1, .4])
pylab.text(2, .3, "moving gaussian", fontsize=14)
# pylab.show()
pylab.savefig("smooth.png", dpi=80)

This data needs smoothing. Below is a visual representation of the differences in the methods of smoothing.

The degree of window coverage for the moving window average, moving triangle, and Gaussian functions are 10, 5, and 5 respectively. Also note that (due to the handling of the “degree” variable between the different functions) the actual number of data points assessed in these three functions are 10, 9, and 9 respectively. The degree for the last two functions represents “spread” from each point, whereas the first one represents the total number of points to be averaged for the moving average.