## Designing a Butterworth low-pass filter with SciPy

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

Advertisements

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

Congratulations.

Yeison Cardona09/04/2011 at 01:27

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

congratulations

Andre09/07/2011 at 20:58

[…] 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 […]

Matplotlib font manager error « Azitech01/09/2011 at 13:03

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.

dirk20/01/2012 at 12:43

You are correct!

Bogdan18/01/2013 at 21:14

Thank you, I made the change.

azimout21/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?)

timo04/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

Bogdan04/02/2013 at 19:03

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

Thanks.

mutaharchalmersMutahar24/02/2012 at 17:57

Multiplying b 1000 does some strange scaling. I found a good normalizing approach is

b *= numpy.sum(a)/numpy.sum(b)

Tom12/09/2012 at 16:00

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

Bogdan18/01/2013 at 19:56

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

azimout21/01/2013 at 10:16

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

Matus20/05/2012 at 12:21

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

azimout20/05/2012 at 14:38

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!

Juanlu00112/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?

azimout15/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!

Juanlu00115/10/2013 at 18:15