Azitech

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')
pyplot.grid(True)

# input signal spectrum
xfreq = numpy.fft.fft(x)
fft_freqs = numpy.fft.fftfreq(nsamps, d=1./samp_rate)
fig.add_subplot(233)
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")
pyplot.grid(True)

# 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)
fig.add_subplot(131)
pyplot.loglog(w, numpy.abs(h))
pyplot.title('Filter Frequency Response')
pyplot.text(2e-3, 1e-5, str(N)+"-th order Butterworth filter")
pyplot.grid(True)

# 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)
fig.add_subplot(235)
pyplot.plot(t, y)
pyplot.title('Filter output - Time Domain')
pyplot.grid(True)

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

pyplot.show()

butterworth filter

About these ads

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

    Yeison Cardona

    09/04/2011 at 01:27

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

    Andre

    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.

    dirk

    20/01/2012 at 12:43

    • You are correct!

      Bogdan

      18/01/2013 at 21:14

      • Thank you, I made the change.

        azimout

        21/01/2013 at 10:16

      • 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?)

        timo

        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

        Bogdan

        04/02/2013 at 19:03

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

    Thanks.

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

      Tom

      12/09/2012 at 16:00

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

        Bogdan

        18/01/2013 at 19:56

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

        azimout

        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?

    Matus

    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…

      azimout

      20/05/2012 at 14:38

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

    https://github.com/scipy/scipy/issues/2980

    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!

    Juanlu001

    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?

      azimout

      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!

        Juanlu001

        15/10/2013 at 18:15


Leave a Reply

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

WordPress.com Logo

You are commenting using your WordPress.com 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

Follow

Get every new post delivered to your Inbox.

%d bloggers like this: