### Circuits vs. Software

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

Would I rather design circuits or software? I’m a software guy (likely due to my lack of working knowledge of circuits) so I’d rather record noisy signals and write software to eliminate the noise, rather than assembling circuits to eliminate the noise for me. In the case of my DIY ECG machine, I’d say I’ve done a great job of eliminating noise via the software route. Most DIY ECG circuits on the net use multiple op-amps and diodes to do this, and have a hardware-based band-pass filter to eliminate frequencies around 60 Hz. Instead of all that fancy stuff, I made a super-crude circuit (a single op-amp and two resisters) to record my ECG. It was INCREDIBLY noisy! So, how did I clean it up with software? I’ll tell you.

The first step in removing electrical noise is classifying it. Most of the noise in my signal were overlapping sine waves caused by my electrodes picking up signals not from my body. This was determined by simply close-up observation of the original trace. Since this sine-type interference is consistant, power-spectral analysis could be applied to determine the frequencies of the noise so I could block them out. I used the fast Fourier transform algorithm (FFT) on the values from my wave file to generate a plot of noise level with respect to frequency (noise was seen as sharp peaks). I manually blocked-out certain regions of the FFT trace that I thought were noise-related (colored bands, all the values of which were set to zero) – a process known as band-pass filtering (something which is possible to do electronically with a more complicated circuit) – then performed an inverse FFT algorithm on the trace. The result was a trace with greatly reduced noise and after a moving window smoothing algorithm was applied the signal was better than ever. Below are some figures. Note that I recorded the raw WAV file with “sound recorder” (not GoldWave) and did all of the processing (including band-pass filtering) within Python.

What an awesome chart! On top we have the power spectral analysis with band-stop filters applied at the colored regions. Below is the trace of the original WAV file (light gray), the inverse-FFT-filtered trace (dark gray), and the smoothed inverse-FFT-filtered trace (black) – the final ECG signal I intend to use.

This is a magnified view of a few heartbeats after the inverse FFT filtering (multi-band-blocking) process was applied. Not bad! Not bad at all…

And, of course, the updated code:

``` import wave, struct, numpy, pylab, scipy

fname='./success3.wav'

"""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))
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

"""a do-everything function to get usable, smoothed data from a WAV."""
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))
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.savefig('out.png',dpi=100)
pylab.show()
print "COMPLETE"
```

### DIY ECG Trial 1

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

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 lol), fed the amplified signal into my sound card, and recorded it as a WAV. The signal is INCREDIBLY 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 consecutative moving window averages (20-point window, performed on the condensed data). The result was a voltage reading that had most of the random interference oscillations removed and, behold, a BEAUTIFUL ECG signal!!! 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. Nifty, huh?

I swear, this is legitimate data that I recorded myself tonight! 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. Incredible, huh? You can clearly see the Q, R, S, and T spikes (described in the previous couple entries)! 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. Until then, here’s some real data!

This is a zoomed-in view of a representative QRST wave of mine. Amazing, huh? I’m proud of my ticker! =o)

This is the output of the python script I wrote tonight

And here’s the code I used: note that it relies on the WAV file I recorded – if you want it just ask me for it and I’ll email it to you – it’s about 12 sec long. Oh yeah, not all of the functions in this script are used to create the image (such as integration calculations), but I left ’em in because you may find them useful.

``` import wave, struct, numpy, pylab, scipy

"""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))
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

"""a do-everything function to get usable, smoothed data from a WAV."""
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

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"
```

### ECG Success!

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

I kept working on my homemade ECG machine (I had to change the values of some of the resisters) and it looks like I’m getting some valid signals! By recording the potential using my sound card (microphone hole = a nice analog to digital converter that every PC has) I was able record my ECG with sound recording software, smooth it, and this is what it looks like. Pretty cool huh?

This was based on a circuit I made using a single op-amp (A LM324 from RadioShack \$1.49). Basically the op-amp just amplifies micropotential generated by my heart and outputs it in such a way that I can connect it to a standard headphone jack to plug into my microphone hole. The signal is very noisy though. I’m thinking about making the intricate circuit (with 6 op-amps) to produce a better signal-to-noise ratio, but first I’ll try coding my way out of the noise.

### DIY ECG Attempt 1: Failure

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

So I followed-through on yesterday’s post and actually tried to build an ECG machine. I had a very small amount of time to work on it, so instead of building the fancy circuit (with 6 band-pass filtered op-amps and diodes posted in the previous entry) I built the most crude circuit that would theoretically work. I used one of the 4 available op-amps from a LM324 chip (pictured to the left) . I was working late at night, and I’m quite colorblind, so I had to take a gamble (as usual) with the resistors I used. Resistors are color coded with bands that represent their resistance. But, since I’m highly colorblind, red, orange, and black all look the same to me. So I can’t be sure I go the resistances right. I’d check it with my digital multimeter but I seem to have lost it and my wife doesn’t remember seeing it recently. Blast! Anyway, I built the sucker, hooked it up to my sound card, and made electrodes by soldering wires to pennies. After a good lick, I attached the pennies to my chest with tape and tried recording. Every time the pennies made contact with my skin, I would see noise on the trace, but I couldn’t seem to isolate a strong heartbeat signal. This is what I saw and the circuit I build to see it:

When I have more time I’ll locate my multimeter and build the full circuit described below. I’ll update everyone on my progress later. Wish me luck!

Perhaps this project will be working soon. How cool is that? Many techno-savvy people have made these DIY ECG machines, but no where on the net do these people describe how to interpret the data. Since I’m planning on building it, testing it, recording ECG data, and processing/analyzing it, I’ll have something truly unique on the internet. Perhaps it will be worthy of mention on Hack-a-day when it’s complete! How cool would that be?

### DIY ECG?

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

Last night my wife put her head on my chest while we were watching a movie. A minute or two later I felt a light sinking feeling in my upper chest, and my wife looked up at me in horror. “Your heart stopped beating!” I assured her that everything was okay (it quickly resumed), and that it happens all the time. I feel the sinking feeling often, know it’s because my heart is briefly beating irregularly, and assume it’s normal. After all, your heart isn’t a robot, it’s a living organ doing the best it can. It’s never perfectly regular, and presumably everybody has momentary irregularities, they just don’t notice them. When I got in bed I began wondering how regular irregular heartbeats are. What would the chances be that I have some kind of arrhythmia? I’ve had a checkup not too long ago by a family practice physician who used a stethoscope on my back to listen to my heartbeat, and he didn’t notice anything. Then again, how often do general practice doctors detect subtle arrhythmia?
I know that whatever problem I have is likely too small to cause any serious troubles, but at the same time I’m becoming obsessed as to determining exactly what my problem is. How many times a day does my heart skip beats? What about nighttime? If only there were some way to record heartbeat data, then I could analyze it and determine the severity of my problem. But wait, data? That would be hours of heartbeat recordings… that means… YES!!!! HOME MADE DO-IT-YOURSELF HARDWARE! WRITING SOFTWARE TO ANALYZE LARGE AMOUNT OF DATA! SUPER-COOL SCIENCY STUFF THAT I CAN MAKE MYSELF, COMBINING BASIC ELECTRONICS WITH BIOMEDICAL DATA ANALYSIS AND, HARAY, PROGRAMMING WITH PYTHON!
Naturally, my thoughts began to overwhelm my reality as soon as Python entered the scene. I wondered how I could use my PC to record my heartbeat, without spending much money on hardware, and only using software I write myself. I pondered this on the way to work this morning, and came up with two possible methods:
Method 1: acoustical recordings. This would be the easiest way to record my heartbeat. I could tape a stethoscope to my chest, insert a small microphone in the earpiece, connect the microphone to my PC, and record sound data for several hours. Theoretically it would work, but it would be highly prone to noise from breathing, and I would have to lay perfectly still to avoid noise caused by movements. The data (trace) would have to be smoothed, processed with a band-pass filter (to eliminate interference), and heartbeats could be calculated. However, this would only give me heart beat time information…
Method 1: electrical recordings. This would be a little more complicated, but generate much more information. I could record the electrical activity of my heart, and the charts would look like the cool electrocardiograms (ECGs) that you see on TV shows and stuff. I would have to build some circuity to amplify the 1mV heartbeat signal recorded from electrodes taped to my chest, but then again what’s more fun than ghetto-building circuits! (I apologize for any people reading this blog who actually live in the ghetto and enjoy building circuits.) I did a little Googling and found that similar things have been done before with a handful of diodes, resisters, and 6 op-amps. I think I’m going to follow the guide on this page and build the circuit seen below:

Schematic of a crude ECG circuit

Supposedly, the data I can obtain looks something like the image at the bottom of this blog entry. I’d attach 3 electrodes to my body (chest, arm, and leg), hook them up to my little circuit, then connect to circuit to my PCs sound card. I’d record the trace (maybe while I sleep?) and analyze it with Python, Numpy, Scipy, and Matplotlib (gosh I love Python). There are several websites which demonstrate how to build DIY ECG recording devices, but none of these actually ANALYZE the data they obtain! Hopefully I could fill this little niche on the internet. We’ll see what happens. I have my thesis to work on, and a whole bunch of other stuff on my plate right now.

A sample recording from the circuit pictured above

The terms assigned to different parts of the heartbeat

UPDATE: I found an extremely crude ECG circuit which I can make from parts I already have at my house. It has tons of noise, but maybe I can filter it out? Perhaps I’ll try this tonight? [ponders]