Azimout's Linux weblog

Designing a Butterworth low-pass filter with SciPy

with 17 comments

from scipy import signal
import math, numpy
from matplotlib import pyplot

# some constants
samp_rate = 20
sim_time = 60
nsamps = samp_rate*sim_time
cuttoff_freq = 0.1

fig = pyplot.figure()

# generate input signal
t = numpy.linspace(0, sim_time, nsamps)
freqs = [0.1, 0.5, 1, 4]
x = 0
for i in range(len(freqs)):
 x += numpy.cos(2*math.pi*freqs[i]*t)
time_dom = fig.add_subplot(232)
pyplot.plot(t, x)
pyplot.title('Filter Input - Time Domain')

# input signal spectrum
xfreq = numpy.fft.fft(x)
fft_freqs = numpy.fft.fftfreq(nsamps, d=1./samp_rate)
pyplot.loglog(fft_freqs[0:nsamps/2], numpy.abs(xfreq)[0:nsamps/2])
pyplot.title('Filter Input - Frequency Domain')
pyplot.text(0.03, 0.01, "freqs: "+str(freqs)+" Hz")

# design filter
norm_pass = cuttoff_freq/(samp_rate/2)
norm_stop = 1.5*norm_pass
(N, Wn) = signal.buttord(wp=norm_pass, ws=norm_stop, gpass=2, gstop=30, analog=0)
(b, a) = signal.butter(N, Wn, btype='low', analog=0, output='ba')
print("b="+str(b)+", a="+str(a))

# filter frequency response
(w, h) = signal.freqz(b, a)
pyplot.loglog(w, numpy.abs(h))
pyplot.title('Filter Frequency Response')
pyplot.text(2e-3, 1e-5, str(N)+"-th order Butterworth filter")

# filtered output
#zi = signal.lfiltic(b, a, x[0:5], x[0:5])
#(y, zi) = signal.lfilter(b, a, x, zi=zi)
y = signal.lfilter(b, a, x)
pyplot.plot(t, y)
pyplot.title('Filter output - Time Domain')

# output spectrum
yfreq = numpy.fft.fft(y)
pyplot.loglog(fft_freqs[0:nsamps/2], numpy.abs(yfreq)[0:nsamps/2])
pyplot.title('Filter Output - Frequency Domain')

butterworth filter

Written by azimout

15/03/2011 at 16:01

Posted in Tricks

17 Responses

Subscribe to comments with RSS.

  1. Very nice great code, I’ll have to study it and understand it a little.

    Yeison Cardona

    09/04/2011 at 01:27

  2. tanks for your code….it’s a very good for my project.


    09/07/2011 at 20:58

  3. […] to work on my machine stopped working at some point. It’s the script presented in the post Designing a Butterworth low-pass filter with SciPy. The last few lines of the traceback […]

  4. Nice code, insightful plot, but you are computing your filters incorrectly: norm_pass = 2*math.pi*cuttoff_freq/samp_rate should be norm_pass = cutoff_freq / (samp_rate/2). The matlab/scipy docs say that the input frequency is normalized with respect to the nyquist frequency (half the sample rate), so that 1 is identical to the nyquist frequency.


    20/01/2012 at 12:43

    • You are correct!


      18/01/2013 at 21:14

      • I found the doc saying: normalized from 0 to 1 (1 corresponds to pi radians / sample)
        So norm_pass = 2*math.pi*cuttoff_freq/samp_rate is correct (again?)


        04/02/2013 at 13:58

      • You are interpreting the doc incorrectly. First of all freq sample is actually the Nyquist frequency and the code above divides by 2 and this is correct. Second if both the freq and sample are in Hz there is no point in multiplying with PI as they are normalized just by dividing them

        The following function will give you what buttord expects.

        def convert_hertz(freq, sample):
        # convert frequency in hz to units of pi rad/sample
        return freq * 2.0 / sample


        04/02/2013 at 19:03

      • Thank you, I made the change.


        21/01/2013 at 10:16

  5. Why do you multiply the filter coefficients b by 1000? Is this a quirk of the Scipy.signal code?


    • Multiplying b 1000 does some strange scaling. I found a good normalizing approach is
      b *= numpy.sum(a)/numpy.sum(b)


      12/09/2012 at 16:00

      • It should be no scaling at all. That line needs to be removed


        18/01/2013 at 19:56

      • Thank you, I made the change. Need to redo the plotting…


        21/01/2013 at 10:16

  6. the filtered data are time-shifted to the right, does anyone know how to correct this time-shift?


    20/05/2012 at 12:21

    • Dear Matus, this is the filter’s delay and is inevitable in a filter, no way to fix it…


      20/05/2012 at 14:38

  7. Someone tried this code but there is a BadCoefficients warning and the plots are different:

    Running the code with older versions of SciPy didn’t solve the problem. Perhaps you didn’t redo the plotting after the code updates? Could you please give it a try? Thanks in advance!


    12/10/2013 at 14:17

    • Sorry, but I’ve moved away from this line of work a while ago, I fear I cannot help you.
      Best of luck anyway! Maybe another one of the readers of this blog can help you?


      15/10/2013 at 13:56

      • Too bad, we’ll keep working on it, but maybe a warning that the code does not produce the figures / post outdated or something would be appropriate. Thanks!


        15/10/2013 at 18:15

Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s


Get every new post delivered to your Inbox.

%d bloggers like this: