«

»

Print this Post

Circuits vs. Software

Notice

I have simplified and improved my ECG machine design! Check out the new post:

http://www.swharden.com/blog/2013-04-14-simple-diy-ecg-pulse-oximeter-version-2/

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

About the author

Scott W Harden

Scott Harden has had a lifelong passion for computer programming and electrical engineering, and recently has become interested in its relationship with biomolecular sciences. He has run a personal website since he was 15, which has changed names from HardenTechnologies.com, to KnightHacker.com, to ScottIsHot.com, to its current SWHarden.com. Scott has been in college for 10 years, with 3 more years to go. He has an AA in Biology (Valencia College), BS in Cell Biology (Union University), MS in Molecular Biology and Microbiology (University of Central Florida), and is currently in a combined DMD (doctor of dental medicine) / PhD (neuroscience) program through the collaboration of the College of Dentistry and College of Medicine (Interdisciplinary Program in Biomedical Science, IDP) at the University of Florida in Gainesville, Florida. In his spare time Scott builds small electrical devices (with an emphasis on radio frequency) and enjoys writing cross-platform open-source software.

Permanent link to this article: http://www.SWHarden.com/blog/2009-01-15-circuits-vs-software/

1 comment

  1. Hiram

    I am a student and research about audio. Thanks for you blog. It is very helpful for me.

Leave a Reply

Your email address will not be published.

You may use these HTML tags and attributes: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <strike> <strong>