The personal website of Scott W Harden
January 15th, 2009

# DIY ECG Trial 1

⚠️ Check out my newer ECG designs:

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

"""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]&gt;data[i-res]+.1 and data[i]&gt;data[i+res]+.1:
if data[i]&gt;data[i-1] and data[i]&gt;data[i+1]:
pXs.append(Xs[i])
pYs.append(data[i])
if len(pXs)&gt;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,ax,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"``````
```---
title: DIY ECG Trial 1
date: 2009-01-15 00:39:45
tags: diyECG, python, old
---

# DIY ECG Trial 1

> **⚠️ Check out my newer ECG designs:**
* [**Single op-amp ECG**](https://swharden.com/blog/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), 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.

<div class="text-center">

[![](diy_ecg2_thumb.jpg)](diy_ecg2.png)

</div>

__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.

<div class="text-center">

[![](diy_ecg1_thumb.jpg)](diy_ecg1.png)

</div>

__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.

```python
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]&gt;data[i-res]+.1 and data[i]&gt;data[i+res]+.1:
if data[i]&gt;data[i-1] and data[i]&gt;data[i+1]:
pXs.append(Xs[i])
pYs.append(data[i])
if len(pXs)&gt;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,ax,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"
```

```
January 15th, 2009

# Circuits vs. Software

⚠️ Check out my newer ECG designs:

Would I rather design circuits or software? I'm a software guy (or at least 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'

"""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]&gt;data[i-res]+.1 and data[i]&gt;data[i+res]+.1:
if data[i]&gt;data[i-1] and data[i]&gt;data[i+1]:
pXs.append(Xs[i])
pYs.append(data[i])
if len(pXs)&gt;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])&gt;low and abs(fftx[i])&lt;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"``````
```---
title: Circuits vs. Software
date: 2009-01-15 17:47:52
tags: diyECG, old, python
---

# Circuits vs. Software

> **⚠️ Check out my newer ECG designs:**
* [**Single op-amp ECG**](https://swharden.com/blog/2016-08-08-diy-ecg-with-1-op-amp/)

__Would I rather design circuits or software?__ I'm a software guy (or at least 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)](http://en.wikipedia.org/wiki/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.

<div class="text-center">

[![](diy_ecg4_thumb.jpg)](diy_ecg4.png)

</div>

__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.

<div class="text-center">

[![](diy_ecg3_thumb.jpg)](diy_ecg3.png)

</div>

__This is a magnified view of a few heartbeats__. It looks pretty good! Here's the code I used to do all the calculations:

```python
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]&gt;data[i-res]+.1 and data[i]&gt;data[i+res]+.1:
if data[i]&gt;data[i-1] and data[i]&gt;data[i+1]:
pXs.append(Xs[i])
pYs.append(data[i])
if len(pXs)&gt;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])&gt;low and abs(fftx[i])&lt;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"
``````
November 24th, 2008

# Compress Strings and Store to Files in Python

While writing code for my graduate research thesis I came across the need to lightly compress a huge and complex variable (a massive 3D data array) and store it in a text file for later retrieval. I decided to use the zlib compression library because it's open source and works pretty much on every platform. I ran into a snag for a while though, because whenever I loaded data from a text file it wouldn't properly decompress. I fixed this problem by adding the "rb" to the open line, forcing python to read the text file as binary data rather than ascii data. Below is my code, written in two functions to save/load compressed string data to/from files in Python.

``````import zlib

def saveIt(data,fname):
data=str(data)
data=zlib.compress(data)
f=open(fname,'wb')
f.write(data)
f.close()
return

def openIt(fname,evaluate=True):
f=open(fname,'rb')
f.close()
data=zlib.decompress(data)
if evaluate: data=eval(data)
return data  ``````

Oh yeah, don't forget the evaluate option in the openIt function. If set to True (default), the returned variable will be an evaluated object. For example, `[[1,2],[3,4]]` will be returned as an actual 2D list, not just a string. How convenient is that?

```---
title: Compress and Store Files in Python
date: 2008-11-24 17:48:16
tags: python, old
---

# Compress Strings and Store to Files in Python

__While writing code for my graduate research thesis__ I came across the need to lightly compress a huge and complex variable (a massive 3D data array) and store it in a text file for later retrieval.  I decided to use the [zlib](http://en.wikipedia.org/wiki/Zlib) compression library because it's open source and works pretty much on every platform.  I ran into a snag for a while though, because whenever I loaded data from a text file it wouldn't properly decompress.  I fixed this problem by adding the "rb" to the open line, forcing python to read the text file as binary data rather than ascii data.  Below is my code, written in two functions to save/load compressed string data to/from files in Python.

```python
import zlib

def saveIt(data,fname):
data=str(data)
data=zlib.compress(data)
f=open(fname,'wb')
f.write(data)
f.close()
return

def openIt(fname,evaluate=True):
f=open(fname,'rb')
f.close()
data=zlib.decompress(data)
if evaluate: data=eval(data)
return data
```

__Oh yeah, don't forget__ the evaluate option in the openIt function.  If set to True (default), the returned variable will be an evaluated object.  For example, `[[1,2],[3,4]]` will be returned as an actual 2D list, not just a string.  How convenient is that?```
November 17th, 2008

# Linear Data Smoothing in Python

⚠️ SEE UPDATED POST: Signal Filtering in Python

``````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 large 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 = *(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 = *30  # 30 "0"s in a row
data = 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, a, -.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, a, -.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, a, -.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, a, -.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.

```---
title: Linear Data Smoothing in Python
date: 2008-11-17 18:50:10
tags: old, python
---

# Linear Data Smoothing in Python

> **⚠️ SEE UPDATED POST:** [**Signal Filtering in Python**](https://swharden.com/blog/2020-09-23-signal-filtering-in-python/)

```python
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 large numbers of data points, the Gaussian function should perform better.

<div class="text-center">

[![](smooth_thumb.jpg)](smooth.png)

</div>

```python
import pylab
import numpy

def smoothList(list, strippedXs=False, degree=10):
if strippedXs == True:
return Xs[0:-(len(list)-(len(list)-degree+1))]
smoothed = *(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 = *30  # 30 "0"s in a row
data = 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, a, -.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, a, -.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, a, -.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, a, -.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.

<div class="text-center">

[![](smooth2_thumb.jpg)](smooth2.png)

</div>

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.```
August 12th, 2005

# Create Custom MiniDisc Labels

Great looking MiniDisc labels can be made at home! In this post I'll show you how I created the nice MiniDisc labels shown on this page. Links to the Photoshop template are at the bottom of the page so you can make your own labels too.

## MiniDisc Label Dimensions

Getting these labels to properly fit your MiniDisc can be harder than it seems. Your best bet is to make a practice sheet where the shapes are solid colors. Print it, cut it out, hold it up, and see if any of the sizes need adjusting. It might take a few tries, so don't be afraid to use up a few sheets of paper!

Measurement Width Height
Front `1.407` inches `2.077` inches
Spine `2.360` inches `0.113` inches

## MiniDisc Label Templates for Photoshop

If you have Photoshop you can download one of the blank MiniDisc label templates I created.

## Printing MiniDisc Labels

The labels I've been making lately have been printed on Avery's White Full Sheet Labels #8165. I got a pack of 25 of them at OfficeMax for \$9.99 (which equates to a little less than 40 cents a sheet). If you're interested in purchasing them, I recommend this brand and know I'll be buying more personally (that is, when I run out of 25 sheets (which is about 150 labels worth of paper)). The best printing settings for this type of paper are those that will place a little more ink on the page than normal. With my HP DeskJet 932C, I've been getting good results by printing with the settings normally used for textured greeting cards.

## Cutting MiniDisc Labels

Cutting your labels well is the crucial step for making great labels. This is your make or break moment! It's incredibly important that you relax and take your time when you cut these labels. The right tools help a lot too. While it is possible to cut your labels with scissors, I strongly recommend you use a double blade slicer. The one I'm using is made for photos. A razor blade is helpful when it comes to separating the sticky paper from the wax coated base sheet.

Begin by using scissors to clip each label set out of the paper. Leave around half an inch around the actual image so you can touch the paper to avoid smearing the fresh ink. Look closely at the image and locate the crosshairs at the corners. If you focus on slicing the image just inside the crosshairs you'll get far better results than if you try to cut on the edge of the colored portion. It's a lot easier to line up a blade with a thin line than it is with a contrasting edge. When you're ready, hold your paper firmly and slice it. Try your best to get clean slices, and don't be afraid to shave off a small portion of the colored area to ensure that you don't leave any white on the edges.

Notice how I cut a 0.5mm border around the edge of the image in preparation for my final trimming. This greatly improves the precision you get when making your final slices.

Pay extra attention when you cut the corner at the top left of the MiniDisc label. This little cut plays a large role in what the finished product will look like. Since it's the only irregularity in the shape of the label, the eye is naturally drawn to it. The angle needs to be as close to 45º as you can get it. That little slice is a lot smaller than you think, so when you cut it don't take too much off! The size of the slice in the edge on the recessed region of your minidisk is about twice the size that your slice should be, since your paper will border slightly inside the recessed area. This step is really, really important. Notice how I line up the edge of my label with the corners of the squares on my slicer to ensure I get a clean, even, and level 45º cut on the corner! Hitomi looks scared in this picture.

Once everything is cut, and before you remove the backing, hold up the labels to the MiniDisc you'll be putting them on to make sure everything is perfect. The labels should fit cleanly, with a small amount of room between the label and the edge of the recessed region. If you're dissatisfied with the color of MiniDisc you had in mind while designing this label, you still have a chance to make a last minute change! I'm second guessing my own color choice here, as the blue seems to go better with the orange now that I look at the printed label.

## Preparing the Labels

Perhaps the most frustrating part of your project is separating the sticky paper from the wax coated paper backing. Some label paper has convenient little creased pre-sliced edges on the back to help you separate the two, however more often than not the stupid little things are never in the right place so you're left with really cool images that can be ruined in a second by clumsy separation of the paper layers. I found that a sharp razor blade helps a lot with this step. If you carefully wedge the tip of your razor between the two layers, you can separate them enough that they can be split by hand.

In the photo below notice how I'm using that little corner slice to my advantage. The corner method works wonders with the spine as well. Just be careful that you don't cut it!

### Applying the Labels

Much like slicing your labels, this next step is a highly visible one that makes a big impact on what your final project will look like. I like to apply the big face label first because it's a little easier. It's really easy to stick a label on crooked if you're not careful, so put extra care into proper alignment. I found I get best results when I use one finger to line up and stick the top right corner of the label to the disc first, then slide a finger down to the bottom right corner, then brush to the left securing the rest. If you line up two corners (preferably ones far apart) you know you'll get a straight and even application.

When applying the spine, attach the leftmost edge of it to the correct place on the edge of the MiniDisc and hold it there with your thumb. Since it's such a narrow strip, a small portion of the sticky paper touching the plastic won't hold too firmly, so you can move it around until it's level to your liking. Once you're confident it's level, go ahead and press the other end down (gently) and rub your finger (lightly) from left to right to smooth out the stick. Next, to be sure this spine is stuck on firmly, put some pressure on it with something that won't smear it (like cloth, your shirt, or perhaps even another sheet of paper if you're careful). When done correctly, the spine labels can make your MiniDisc look awesome!

## Final Result

I originally wrote this article when I was a teenager, almost twenty years ago. I'm pretty embarrassed about the content of the original article now that I read it as an adult. In 2021 (at the age of 35) I rediscovered this article, edited-out the cringiest segments, and re-uploaded it to my blog using a redirect from the original URL. Articles like this are part of my history, a glimpse into my past, an opportunity to reflect on far I've come, and something that I think is worth preserving. -- Scott

```---
title: Custom MiniDisc Labels with Photoshop
date: 2005-08-12 13:11:59
tags: old
---

# Create Custom MiniDisc Labels

**Great looking MiniDisc labels can be made at home!** In this post I'll show you how I created the nice MiniDisc labels shown on this page. Links to the Photoshop template are at the bottom of the page so you can make your own labels too.

<div class="text-center img-medium img-border">

![](images/gigdisk.jpg)
![](images/smallstack.jpg)

</div>

## MiniDisc Label Dimensions

**Getting these labels to properly fit your MiniDisc can be harder than it seems.** Your best bet is to make a practice sheet where the shapes are solid colors. Print it, cut it out, hold it up, and see if any of the sizes need adjusting. It might take a few tries, so don't be afraid to use up a few sheets of paper!

<div align="center">

Measurement|Width|Height
---|---|---
Front|`1.407` inches|`2.077` inches
Spine|`2.360` inches|`0.113` inches

</div>

## MiniDisc Label Templates for Photoshop

If you have Photoshop you can download one of the blank MiniDisc label templates I created.

<div class="text-center img-small">

![](images/miniDiscLabel.png)
![](images/miniDiscLabelSheet.png)

</div>

## Using The Photoshop Templates

<div class="text-center">

![](images/photoshopping.gif)

</div>

## Printing MiniDisc Labels

The labels I've been making lately have been printed on [Avery's White Full Sheet Labels #8165](https://www.amazon.com/s?k=avery+8165). I got a pack of 25 of them at OfficeMax for \$9.99 (which equates to a little less than 40 cents a sheet). If you're interested in purchasing them, I recommend this brand and know I'll be buying more personally (that is, when I run out of 25 sheets (which is about 150 labels worth of paper)). The best printing settings for this type of paper are those that will place a little more ink on the page than normal. With my HP DeskJet 932C, I've been getting good results by printing with the settings normally used for textured greeting cards.

<div class="text-center">

![](images/paper_small.jpg)
![](images/printersettings.jpg)

</div>

## Cutting MiniDisc Labels

**Cutting your labels well is the crucial step for making great labels.** This is your make or break moment! It's incredibly important that you relax and take your time when you cut these labels. The right tools help a lot too. While it is possible to cut your labels with scissors, I strongly recommend you use a double blade slicer. The one I'm using is made for photos. A razor blade is helpful when it comes to separating the sticky paper from the wax coated base sheet.

<div class="text-center">

![](images/cutter_small.jpg)

</div>

**Begin by using scissors to clip each label set out of the paper.** Leave around half an inch around the actual image so you can touch the paper to avoid smearing the fresh ink. Look closely at the image and locate the crosshairs at the corners. If you focus on slicing the image just inside the crosshairs you'll get far better results than if you try to cut on the edge of the colored portion. It's a lot easier to line up a blade with a thin line than it is with a contrasting edge. When you're ready, hold your paper firmly and slice it. Try your best to get clean slices, and don't be afraid to shave off a small portion of the colored area to ensure that you don't leave any white on the edges.

<div class="text-center">

![](images/printed_small.jpg)

</div>

Notice how I cut a 0.5mm border around the edge of the image in preparation for my final trimming. This greatly improves the precision you get when making your final slices.

<div class="text-center img-medium img-border">

![](images/edge.jpg)

</div>

**Pay extra attention when you cut the corner at the top left of the MiniDisc label.** This little cut plays a large role in what the finished product will look like. Since it's the only irregularity in the shape of the label, the eye is naturally drawn to it. The angle needs to be as close to 45º as you can get it. That little slice is a lot smaller than you think, so when you cut it don't take too much off! The size of the slice in the edge on the recessed region of your minidisk is about twice the size that your slice should be, since your paper will border slightly inside the recessed area. This step is really, really important. Notice how I line up the edge of my label with the corners of the squares on my slicer to ensure I get a clean, even, and level 45º cut on the corner! Hitomi looks scared in this picture.

<div class="text-center img-medium img-border">

![](images/corner.jpg)

</div>

**Once everything is cut,** and before you remove the backing, hold up the labels to the MiniDisc you'll be putting them on to make sure everything is perfect. The labels should fit cleanly, with a small amount of room between the label and the edge of the recessed region. If you're dissatisfied with the color of MiniDisc you had in mind while designing this label, you still have a chance to make a last minute change! I'm second guessing my own color choice here, as the blue seems to go better with the orange now that I look at the printed label.

<div class="text-center img-medium img-border">

![](images/cut.jpg)

</div>

## Preparing the Labels

**Perhaps the most frustrating part of your project is separating the sticky paper from the wax coated paper backing.** Some label paper has convenient little creased pre-sliced edges on the back to help you separate the two, however more often than not the stupid little things are never in the right place so you're left with really cool images that can be ruined in a second by clumsy separation of the paper layers. I found that a sharp razor blade helps a lot with this step. If you carefully wedge the tip of your razor between the two layers, you can separate them enough that they can be split by hand.

<div class="text-center img-medium img-border">

![](images/snipcorner.jpg)

</div>

In the photo below notice how I'm using that little corner slice to my advantage. The corner method works wonders with the spine as well. Just be careful that you don't cut it!

<div class="text-center img-medium img-border">

![](images/snipspine.jpg)

</div>

### Applying the Labels

**Much like slicing your labels, this next step is a highly visible one that makes a big impact on what your final project will look like.** I like to apply the big face label first because it's a little easier. It's really easy to stick a label on crooked if you're not careful, so put extra care into proper alignment. I found I get best results when I use one finger to line up and stick the top right corner of the label to the disc first, then slide a finger down to the bottom right corner, then brush to the left securing the rest. If you line up two corners (preferably ones far apart) you know you'll get a straight and even application.

<div class="text-center img-medium img-border">

![](images/flat.jpg)

</div>

**When applying the spine,** attach the leftmost edge of it to the correct place on the edge of the MiniDisc and hold it there with your thumb. Since it's such a narrow strip, a small portion of the sticky paper touching the plastic won't hold too firmly, so you can move it around until it's level to your liking. Once you're confident it's level, go ahead and press the other end down (gently) and rub your finger (lightly) from left to right to smooth out the stick. Next, to be sure this spine is stuck on firmly, put some pressure on it with something that won't smear it (like cloth, your shirt, or perhaps even another sheet of paper if you're careful). When done correctly, the spine labels can make your MiniDisc look awesome!

<div class="text-center img-medium img-border">

![](images/prettyslim.jpg)

</div>

## Final Result

<div class="text-center img-small img-border">

![](images/prettynice.jpg)
![](images/finished_sideways.jpg)

</div>

<div class="text-center img-medium img-border">

![](images/inhand.jpg)

</div>