«

»

Print this Post

Signal Filtering with Python

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/

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

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-21-signal-filtering-with-python/

39 comments

  1. Anonymous

    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

  2. Gopal Koduri

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

  3. CoreDistance

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

  4. Jphn

    “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!

  5. Kamil

    Thank you, great article!

  6. Danny

    Thanks for the article, sure it’ll come in handy cleaning some of my (very noisy) data.

  7. Ken

    Be careful with this method of “filtering.” This kind of direct manipulation of the spectrum results in time-aliasing when you take the inverse FFT. The result is not truly the original signal with its high frequencies filtered out, because the time-aliasing leaves the signal littered with artifacts.

    The easiest way to do it properly would be to do a lowpass filter in the time-domain.

    You can do a similar filter in the frequency domain (sort of like you’re doing), but the time-domain data and a filter impulse response needs to be properly zero-padded before you FFT them to avoid circular convolution (which results in time-aliasing).

  8. Anonymous

    Thanks for the code snippet. I’ve got a few suggestions for you, though:
    1) When I use the line “bp=fft[:]“, modifying bp still alters fft. I’ve used “bp = fft.copy()” instead.
    2) When zeroing out FFT terms, you need to make the zeros symmetric around the mid-point of the FFT spectrum. The inverse FFT will have a non-negligible imaginary component if the frequency spectrum isn’t symmetric, even though the input data is strictly real.
    3) To zero out fft coefficients, it’s easier to write bp[10:-10] = 0
    4) To normalize the data, add the line bp *= real(fft.dot(fft))/real(bp.dot(bp)). fft.dot(fft) gives the total power in the signal, so this redistributes the power in the stop band more or less equally to the pass band.
    5) It’s perhaps better to change the variable name ‘fft’ to something else so that it doesn’t conflict with numpy.fft.fft. If you’re working in iPython, numpy.fft.fft is automagically imported into the global namespace, so trying to use the fft function after running your script could cause a few headaches.

  9. Ahmet

    Thanks for the information, but I was wondering something.
    Maybe you can help me out with a problem. I have a signal of 0.1 seconds length which consists of 10.000 datapoints. First I averaged the datapoints 40 times such that it becomes a set of 250 datapoints. Even though this is a bit rough, the signal is still clearly visible.
    When I make an FFT of this averaged signal the frequency resolution is too low; the frequency step is 10 Hz.
    To increase the frequency resolution I decreased the averaging to 10 times (1000 datapoints), but I noticed that this does not matter for the frequency resolution; it stays 10 Hz.

    So the frequency resolution is only determined by the measurement time (0.1s) and will not increase if I increase the number of datapoints during this 0.1 second…

    My question is how I can increase the freq. resolution…

  10. Doug

    Thank you for this excellent example of the power of numpy, scipy, and pylab!

    1 more note. In graph (G) you label the right hand side peak as “static”. If you zoom in on the peak with:

    pylab.axis([9840,9900,None,None])

    You’ll see that it’s the symmetric part of your sine waves.

  11. Michael

    I was thinking to make an ECG. I’m a real amateur. I don’t want to visualize my heartbeat but I want to transfer the signal of my heartbeat to a digital signal to use in a a microprocessor. I don’t wanne use a computer at all.
    Can you pleas tell me how I could do it on a simple way? I wanne use the minimum of filters.

    sorry fore my bad English.

  12. Elliot

    Thank you very much.
    This was very helpful for me processing ultrasound signals.

  13. Laboratordentar.Info

    What i do not understood is if truth be told how you are no longer really a lot more
    well-favored than you may be right now. You are so intelligent.
    You know thus considerably with regards to this topic, produced me for
    my part believe it from numerous various angles.
    Its like men and women don’t seem to be fascinated unless it’s one thing to do with Girl gaga!
    Your personal stuffs nice. At all times take care of it up!

  14. Bogdan Hlevca

    Great code, however, it has a serious weakness as it does not explain how the frequency value is selected for filtering. It can be inferred from the code, but there is a more elegant way of doing that, using the numpy functions instead of scipy functions.
    Another problem that you have is the resulting value is half in amplitude. You need to multiply by 2 to get the expected value.

    For example a filer function can be defined as:


    import numpy as np

    def fft_bandpassfilter(data, fs, lowcut, highcut):
    fft = np.fft.fft(data)
    n = len(data)
    timestep = 1.0 / fs
    freq = np.fft.fftfreq(n, d = timestep)
    bp = fft[:]
    for i in range(len(bp)):
    if freq[i] >= highcut or freq[i] < lowcut:
    bp[i] = 0

    # must multipy by 2 to get the correct amplitude
    ibp = 2 * sp.ifft(bp)
    return ibp

    where fs = sample frequency, data ids your timeseries and the other two the cutoff frequencies.

    1. Stu

      Hi @Bogdan Hlevca
      I’m trying to get fft_bandpassfilter working, is lowcut and highcut in HZ ?

      + frequencies within this range will be kept in ?

      For sample frequency, I guess 44010 should work if I’m sampling pc microphone ?

      Cheers
      S

  15. nightclub greater london

    The main reason I observe this blog is since I believe you do usually give a slightly distinct slant on things
    to numerous other internet sites so high five from me!
    ! !

  16. Boiler repair

    Thanks Scott when i started reading I didnt think it would be for me but ended up being quite interested – thanks!

  17. Niels

    I would be nice if you warned readers about the flaws in this article. By removing the correct part of the FFT results, no scaling at all is needed. That’s the great thing about spectral filtering (al least in 1D, and preferably on something with cyclic boundaries): all the information beyond the cut-off frequency is maintained (plot some spectra of the data before and after filtering and you will see the result).

    Anyhow, to warn readers: BE CAREFULL WITH THIS ARTICLE! NICE ATTEMPT, BUT WRITER IS BY FAR NO SPECIALIST, EXAMPLE CODE CONTAINS SERIOUS FLAWS!

  18. online payday loan

    Wonderful article! We will be linking to this
    great post on our site. Keep up the great writing.
    online payday loan

  19. study skills

    I’m gone to tell my little brother, that he should also pay a visit this web site on regular basis to take updated from latest reports.

  20. Http://Genwiki.Nl/

    Hey just ωаntеd to give you a quіck hеads up and let yοu κnow а few of thе pictures aren’t loading properly. I’m not sure why but I think its a linκing isѕue.
    I’ve tried it in two different web browsers and both show the same results.

  21. design expert Central london

    You really were on the money with a fantastic posting with a bit of wonderful information

  22. nail fungus

    Since the admin of this site is working, no doubt very soon it
    will be well-known, due to its feature contents.

  23. Brigette

    I was pretty pleased to find this great site. I wanted
    to thank you for ones time for this wonderful read!
    ! I definitely enjoyed every part of it and i also have you
    saved to fav to see new stuff in your web site.

  24. Archeage gold

    I am only commenting to let you know of the fine experience my cousin’s princess developed reading your site. She noticed several things, with the inclusion of what it is like to possess a great helping style to let many people clearly have an understanding of selected specialized topics. You truly did more than visitors’ expected results. Thanks for showing such informative, dependable, explanatory and even cool tips on your topic to Evelyn.

  25. Brooke

    If you desire to get much from this paragraph then you have to apply such methods to your
    won webpage.

  26. tradnonptbetback.seesaa.net

    Real nice layout and fantastic subject matter

    1. Anonymous

      Thankyou for sharing this project.

  27. creer un site internet gratuitement

    It’s going to be end of mine day, buut before finish I am reading this fantastic
    article to improve myy know-how.

    Review my blog :: creer un site internet gratuitement

  28. agence webdesign

    I will immediately grab your rss feed as I can’t to find your email subscription hyperlink
    or newsletter service. Do you’ve any? Please permit
    me know so that I may just subscribe. Thanks.

    Check out my weeb site – agence webdesign

  29. Jessika

    Highly descriptive post, I liked that a lot.
    Will there be a part 2?

    my web-site – diy home hydroponics (Jessika)

  30. controle Curriculum

    Haave you ever thought about incdluding
    a little bit more than just your articles? I mean,
    what you say is fundamental and everything. Nevertheless think about if you
    added some great graphics or video clips to give your posts
    more, “pop”! Your content is excellent but with pics and videos, tthis site could undeniably be oone of the best in its field.

    Wonderful blog!

    Also visit my homepage; controle Curriculum

  31. compliantness

    That is a really good tip especially to those fresh
    to the blogosphere. Short but very accurate info…
    Thank you for sharing this one. A must read post!

  32. Johne232

    Some genuinely great information, Glad I discovered this. Good teaching is onefourth preparation and threefourths theater. by Gail. ddgecfeadgaf

  33. hay day cheats youtube

    Today, allows examine about frontiervsecrets.com from Tony ‘Tdub’ Sanders and how you may
    be assisted by it. Bejeweled Blitz has over 10 million participants.
    Simply don’t disturb me while I am chatting, okay?

  34. Alan

    Guitar Hero is incredibly popular – specifically
    with the youngster audience. In certain methods it doesn’t seem so terrible.
    40% surpasses percent. You’ve adventure , military, and westerns.

  35. Jesus

    An individual-hung window features a fixed sash together with it, as
    the lower sash is portable, fundamentally. Eliminate any unnecessary files from your own Laptop.
    Your choices in blinds and curtains are countless.

  36. movie 2k

    Our popular life is adjusting. The normal illustration I wish to say may be the automotive
    for people. There is a great example to examine competitive in golf the same as strength training.

  37. http://optimizaresite.bravesites.com/

    În – fapt , ea a fost semnificativ mai mic decât anual,
    de când am stabilit primul meu raport on-line . Verificați din nou mai
    târziu pentru mai mult pe cel mai bun mod de a utiliza de relații publice pentru a
    opera un optimizarea motorului de căutare de vehicul .

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>