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()
Very nice great code, I’ll have to study it and understand it a little.
Congratulations.
Yeison Cardona
09/04/2011 at 01:27
tanks for your code….it’s a very good for my project.
congratulations
Andre
09/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 « Azitech
01/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.
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
Why do you multiply the filter coefficients b by 1000? Is this a quirk of the Scipy.signal code?
Thanks.
mutaharchalmersMutahar
24/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)
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
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
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
hello, would you please tell me why you chose gpass=2, gstop=30 for the filter?!
杨薇
25/03/2017 at 15:11