11:37:13 pm on 9/6/10
Menu
» Home
» About Scott
» VD Labs
» QRSS VD
» Old Stuff
» Archive
» Contact

Categories
» C/C++
» Circuitry
» DIY ECG
» General
» Linux
» Microcontrollers
» Molecular Biology
» My Website
» PHP
» Prime Numbers
» Python
» Radio
» UCF Lab
» Everything
Writings
» MD Labels
» Streamrip
» AIM Thoughts
» WindowsXP?
» Partitioning
» CD/DVD Repair
» Monitor Info
» CRT Deflection
» Venomcrack
» Flash Thing
» Heart/Brain
» Diabetes
» Triops
» Biomed

Friends
» Fred
» Kyle W
» Nick
» Louis
» Tom
» Kyle H




Archives
» September 2010
» August 2010
» July 2010
» June 2010
» May 2010
» April 2010
» March 2010
» February 2010
» January 2010
» December 2009
» September 2009
» August 2009
» July 2009
» June 2009
» May 2009
» April 2009
» March 2009
» February 2009
» January 2009
» December 2008
» November 2008
» October 2008
» September 2008
» September 2007
» December 2006
» August 2006
» January 2006
» August 2005
» July 2005
» June 2005
» May 2005
» April 2005
» March 2005
» February 2005
» January 2005
» December 2004
» November 2004
» October 2004
» September 2004
» August 2004
» July 2004
» June 2004
» May 2004
» April 2004
» March 2004
» February 2004
» January 2004
» December 2003
» November 2003
» October 2003
» September 2003
» August 2003
» July 2003
» June 2003
» May 2003
» April 2003
» March 2003
» February 2003
» January 2003
» December 2002
» November 2002
» October 2002
» September 2002
» June 2001
« DIY ECG Detected an Irregular Heartbeat
Using PHP to Create Apache-Style Access.log »


Signal Filtering with Python
1,437 words | Posted by Scott on January 21st, 2009
Scott was 23.33 years old when he wrote this!
Filed under: DIY ECG, General, Python

It’s time for a lecture. I’ve been spending a lot of time creating a DIY dlectrocardiogram and it produces fairly noisy signals. I’ve spent some time and effort researching the best ways to clean-up these signals, and the results are incredibly useful! Therefore, I’ve decided to lightly document these results in a blog entry.

Here’s an example of my magic! I take a noisy recording and turn it into a beautiful trace. See the example figure with the blue traces. How is this possible? Well I’ll explain it for you. Mostly, it boils down to eliminating excess high-frequency sine waves which are in the original recording due to electromagnetic noise. A major source of noise can be from the alternating current passing through wires traveling through the walls of your house or building. My original ECG circuit was highly susceptible to this kind of interference, but my improved ECG circuit eliminates most of this noise. However, noise is still in the trace (see the figure to the left), and it needed to be removed.

The key is the FFT (Fast Fourier Transformation) algorithm which can get pretty intimidating to research at first! I’ll simplify this process. Let’s say you have a trace with repeating sine-wave-shaped noise. The output of the FFT transformation of the signal is the breakdown of the signal by frequency. Check out this FFT trace of a noisy signal from a few posts ago (top graph). High peaks represent frequencies which are common. See the enormous peak around 60 Hz? (Hz means “per second” by the way, so 60 Hz is a sine wave that repeats 60 times a second) That’s from my AC noise. Other peaks (shown in colored bands) are other electromagnetic noise sources, such as wireless networks, TVs, telephones, and your computer processor. The heart produces changes in electricity that are very slow (the heartbeat is about 1 Hz, or 1 beat per second), so if we can eliminate all of the sine waves with frequencies higher than what we want to isolate we can get a pretty clear trace. This is called a band-stop filter (we block-out certain bands of frequencies) A band-pass filter is the opposite, where we only allow frequencies which are below (low-pass) or above (high-pass) a given frequency. By eliminating each of the peaks in the colored regions (setting each value to 0), then performing an inverse fast Fourier transformation (going backwards from frequency back to time), the result is the signal trace (seen as light gray on the bottom graph) with those high-frequency sine waves removed! (the gray trace on the bottom graph). A little touch-up smoothing makes a great trace (black trace on the bottom graph).

Here’s some Python code to get you started in cleaning-up your noisy signals! The image below is the output of the Python code at the bottom of this entry. This python file requires that test.wav (~700kb) (an actual ECG recording of my heartbeat) be saved in the same folder. Brief descriptions of each portion of the graph will follow.

(A) The original signal we want to isolate. (IE: our actual heart signal)

(B) Some electrical noise. (3 sine waves of different amplitudes and periods)

(C) Electrical noise (what happens when you add those 3 sine waves together)

(D) Static (random noise generated by a random number generator)

(E) Signal (A) plus static (D)

(F) Signal (A) plus static (D) plus electrical noise (C)

(G) Total FFT trace of (F). Note the low frequency peak due to the signal and electrical noise (near 0) and the high frequency peak due to static (near 10,000)

(H) This is a zoomed-in region of (F) showing 4 peaks (one for the original signal and 3 for high frequency noise). By blocking-out (set it to 0) everything above 10Hz (red), we isolate the peak we want (signal). This is a low-pass filter.

(I) Performing an inverse FFT (iFFT) on the low-pass iFFT, we get a nice trace which is our original signal!

(J) Comparison of our iFFT with our original signal shows that the amplitude is kinda messed up. If we normalize each of these (set minimum to 0, maximum to 1) they line up. Awesome!

(K) How close were we? Graphing the difference of iFFT and the original signal shows that usually we’re not far off. The ends are a problem though, but if our data analysis trims off these ends then our center looks great.

Here’s the code I used to make the image:

 import numpy, scipy, pylab, random  

 # This script demonstrates how to use band-pass (low-pass)  
 # filtering to eliminate electrical noise and static  
 # from signal data!  

 ##################  
 ### PROCESSING ###  
 ##################  

 xs=numpy.arange(1,100,.01) #generate Xs (0.00,0.01,0.02,0.03,...,100.0)  
 signal = sin1=numpy.sin(xs*.3) #(A)  
 sin1=numpy.sin(xs) # (B) sin1  
 sin2=numpy.sin(xs*2.33)*.333 # (B) sin2  
 sin3=numpy.sin(xs*2.77)*.777 # (B) sin3  
 noise=sin1+sin2+sin3 # (C)  
 static = (numpy.random.random_sample((len(xs)))-.5)*.2 # (D)  
 sigstat=static+signal # (E)  
 rawsignal=sigstat+noise # (F)  
 fft=scipy.fft(rawsignal) # (G) and (H)  
 bp=fft[:]  
 for i in range(len(bp)): # (H-red)  
     if i>=10:bp[i]=0  
 ibp=scipy.ifft(bp) # (I), (J), (K) and (L)  

 ################  
 ### GRAPHING ###  
 ################  

 h,w=6,2  
 pylab.figure(figsize=(12,9))  
 pylab.subplots_adjust(hspace=.7)  

 pylab.subplot(h,w,1);pylab.title("(A) Original Signal")  
 pylab.plot(xs,signal)  

 pylab.subplot(h,w,3);pylab.title("(B) Electrical Noise Sources (3 Sine Waves)")  
 pylab.plot(xs,sin1,label="sin1")  
 pylab.plot(xs,sin2,label="sin2")  
 pylab.plot(xs,sin3,label="sin3")  
 pylab.legend()  

 pylab.subplot(h,w,5);pylab.title("(C) Electrical Noise (3 sine waves added together)")  
 pylab.plot(xs,noise)  

 pylab.subplot(h,w,7);pylab.title("(D) Static (random noise)")  
 pylab.plot(xs,static)  
 pylab.axis([None,None,-1,1])  

 pylab.subplot(h,w,9);pylab.title("(E) Signal + Static")  
 pylab.plot(xs,sigstat)  

 pylab.subplot(h,w,11);pylab.title("(F) Recording (Signal + Static + Electrical Noise)")  
 pylab.plot(xs,rawsignal)  

 pylab.subplot(h,w,2);pylab.title("(G) FFT of Recording")  
 fft=scipy.fft(rawsignal)  
 pylab.plot(abs(fft))  
 pylab.text(200,3000,"signals",verticalalignment='top')  
 pylab.text(9500,3000,"static",verticalalignment='top',  
            horizontalalignment='right')  

 pylab.subplot(h,w,4);pylab.title("(H) Low-Pass FFT")  
 pylab.plot(abs(fft))  
 pylab.text(17,3000,"sin1",verticalalignment='top',horizontalalignment='left')  
 pylab.text(37,2000,"sin2",verticalalignment='top',horizontalalignment='center')  
 pylab.text(45,3000,"sin3",verticalalignment='top',horizontalalignment='left')  
 pylab.text(6,3000,"signal",verticalalignment='top',horizontalalignment='left')  
 pylab.axvspan(10,10000,fc='r',alpha='.5')  
 pylab.axis([0,60,None,None])  

 pylab.subplot(h,w,6);pylab.title("(I) Inverse FFT")  
 pylab.plot(ibp)  

 pylab.subplot(h,w,8);pylab.title("(J) Signal vs. iFFT")  
 pylab.plot(signal,'k',label="signal",alpha=.5)  
 pylab.plot(ibp,'b',label="ifft",alpha=.5)  

 pylab.subplot(h,w,10);pylab.title("(K) Normalized Signal vs. iFFT")  
 pylab.plot(signal/max(signal),'k',label="signal",alpha=.5)  
 pylab.plot(ibp/max(ibp),'b',label="ifft",alpha=.5)  

 pylab.subplot(h,w,12);pylab.title("(L) Difference / Error")  
 pylab.plot(signal/max(signal)-ibp/max(ibp),'k')  

 pylab.savefig("SIG.png",dpi=200)  
 pylab.show()  
 




This entry was posted on Wednesday, January 21st, 2009 at 4:25 pmand is filed under DIY ECG, General, Python. You can follow any responses to this entry through the RSS 2.0 feed. You can skip to the end and leave a response. Pinging is currently not allowed.



5 Responses to “Signal Filtering with Python”

Anonymous wrote the following at 03:38:42 AM on August 10th, 2009

I hope you can help or teach me how to smoothing my signal. data[0] is the realtime data and data[1] and data [2] is calculated from data [0]. because the initial data (data[0]) have a noise so that it disturb the other datas. Here is part of my program,

def update(): # Called periodically by the Tk toolkit
global data, running, timegap, ph
#s = time.time()
if running == False:
return
res = ph.get_voltage()
x = res[0] – starts[0]
v = res[1]
m = v2d(v)
data[0].append((x,m))
if len(data[0]) > 1: # plot position & calculate instantaneous velocity
dispobjects[0].delete_lines()
dispobjects[0].line(data[0], ‘black’)
vel = (data[0][-1][1] – data[0][-2][1]) / (data[0][-1][0] – data[0][-2][0])
data[1].append( (data[0][-1][0], vel) )

if len(data[1]) > 1: # plot velocity & calculate acceleration
dispobjects[1].delete_lines()
dispobjects[1].line(data[1], ‘black’)
acn = (data[1][-1][1] – data[1][-2][1]) / (data[1][0][-1] – data[1][0][-2])
data[2].append( (data[0][-1][0], acn) )

if len(data[2]) > 1: # plot acceleration
dispobjects[2].delete_lines()
dispobjects[2].line(data[2], ‘black’)
root.after(timegap, update)
if x > maxtime: running = False

nice to hear from yoyu soon

Gopal Koduri wrote the following at 08:14:11 AM on November 11th, 2009

thanks babai (a kind of *dude* in my language ;) )! I was looking for a similar thing(filtering) in python. This helped me get started :)

CoreDistance wrote the following at 06:16:30 PM on January 30th, 2010

I think the wav file is not used in the code.

Jphn wrote the following at 06:39:32 AM on May 17th, 2010

“This python file requires that test.wav (~700kb) (an actual ECG recording of my heartbeat) be saved in the same folder.”
=> false

Thanks for the code!

Kamil wrote the following at 05:13:43 PM on August 19th, 2010

Thank you, great article!

Leave a Reply




copyright © 2006 swharden@gmail.com