import addutils.toc ; addutils.toc.js(ipy_notebook=True)
from addutils import css_notebook
import bokeh.plotting as bk
css_notebook()
The following is an introduction on how to design an infinite impulse response (IIR) filters using the Python scipy.signal
package. This notebook is not a thorough introduction to IIR filter design.
Documentation and reference:
A filter is usually defined in the frequency domain. The following figure depicts the frequency domain specification:
The specification figure illustrates a low-pass filter but the terminology applies to all filter types, lowpass, highpass, bandpass, and stopband.
NOTE: The cutoff frequency is the -3dB point. If the cutoff frequency is required the algorithm will work to meet the -3dB point at the Wc frequency. The pass frequency, Wp, is the last, first, or first and last frequencies in the passband which the Rp value occurs (the rest are less).
The above specifications are the information that will be used and passed to the IIR filter design functions. In addition to our filter specification information we need to select our "filter" type. The following table summarizes the different IIR filter design approaches supported by the scipy.signal
design functions:
import numpy as np
from scipy import signal
from scipy.signal import filter_design as fd
import bokeh.plotting as bk
bk.output_notebook()
# Specification for our filter
Wp = 0.270 # Cutoff frequency
Ws = 0.412 # Stop frequency
Rp = 0.1 # passband maximum loss (gpass)
As = 60 # stoppand min attenuation (gstop)
Filters = {'ellip' : (), 'cheby2' : (), 'butter' : (), 'cheby1' : (), 'bessel' : ()}
# The ellip and cheby2 filter design
Filters['ellip'] = fd.iirdesign(Wp, Ws, Rp, As, ftype='ellip')
Filters['cheby2'] = fd.iirdesign(Wp, Ws, Rp, As, ftype='cheby2')
# The butter and cheby1 need less constraint spec
Rpl = Rp*10; Asl = As/4.
Filters['butter'] = fd.iirdesign(Wp, Ws, Rp, As, ftype='butter')
Filters['cheby1'] = fd.iirdesign(Wp, Ws, Rp, As, ftype='cheby1')
# The bessel max order of 8 for this cutoff, can't use
# iirdesign have to use iirfilter.
Filters['bessel'] = fd.iirfilter(8, Wp, btype='lowpass', ftype='bessel')
To define plot colors, any legal HTML name for colors can be used:
plotcolors = ['DarkRed', 'Crimson', 'Red',
'Olive', 'Green', 'SeaGreen',
'SteelBlue', 'Blue', 'Navy',
'Fuchsia', 'DarkViolet', 'Indigo']
Plotting:
a, b
are respectively numerator and denominator of the filterw
are the frequencies at which the frequency response 'h'
is computedxmn1, xmx1, ymn1, ymx1 = 0.1, 0.5, -75., 5.
xmn2, xmx2, ymn2, ymx2 = 0.1, 0.5, -22., -0.
fig1 = bk.figure(plot_width=350, plot_height=400,
title='Frequency response',
y_axis_label='Magnitude (db)',
x_axis_label='Normalized Frequency',
x_range=(xmn1, xmx1), y_range=(ymn1, ymx1))
fig2 = bk.figure(plot_width=350, plot_height=400,
title='Phase response',
x_axis_label='Normalized Frequency',
x_range=(xmn2, xmx2), y_range=(ymn2, ymx2))
for i, filtername in enumerate(Filters):
b, a = Filters[filtername]
w,h = signal.freqz(b,a)
h_dB = 20*np.log10(np.abs(h))
h_Phase = np.unwrap(np.arctan2(np.imag(h),np.real(h)))
fig1.line(w/np.max(w),h_dB, color=plotcolors[i].lower(), legend=filtername)
fig2.line(w/np.max(w),h_Phase,color=plotcolors[i].lower(), legend=filtername)
fig1.line([Wp, Wp], [ymx1, ymn1], line_dash=[5])
fig1.line([Ws, Ws], [ymx1, ymn1], line_dash=[5])
fig1.line([xmn1, xmx1], [-As, -As], line_dash=[5])
fig1.legend.location = "bottom_left"
fig2.legend.location = "bottom_left"
bk.show(bk.gridplot([[fig1, fig2]]))
We can have a better look at the Passband Ripple:
fig = bk.figure(title='Frequency Response',
x_axis_label='Normalized Frequency',
y_axis_label='Magnitude (db)',
plot_width=700, plot_height=400,
x_range=(xmn1, xmx1), y_range=(ymn1, ymx1))
for i, filtername in enumerate(Filters):
b, a = Filters[filtername]
w,h = signal.freqz(b,a)
h_dB = 20*np.log10(np.abs(h))
fig.line(w/np.max(w),h_dB, color=plotcolors[i].lower(), legend=filtername)
xmn1, xmx1, ymn1, ymx1 = 0.1, 0.3, -0.11, 0.01
fig.line([Wp, Wp], [ymx1, ymn1], line_dash=[5])
fig.line([Ws, Ws], [ymx1, ymn1], line_dash=[5])
fig.line([xmn1, xmx1], [-As, -As], line_dash=[5])
fig.legend.location = "bottom_right"
bk.show(fig)
We can have a look to the inpulse and step response for any filter (bessel in this example):
from bokeh.models.tickers import AdaptiveTicker
b, a = Filters['bessel']
impulse = np.zeros(50); impulse[0] = 1
x = np.arange(0, 50)
response = signal.lfilter(b, a, impulse)
step = np.cumsum(response)
d1 = bk.figure(title="Impulse response", plot_width=350, plot_height=400)
d1.xaxis.axis_label = 'n (samples)'
d1.yaxis.axis_label = 'Amplitude'
d1.segment(range(len(response)), 0, range(len(response)), response,
line_width=2, line_color="green")
d1.circle(range(len(response)), response,
size=5, fill_color="darkred", line_color="green", line_width=1)
d2 = bk.figure(title="Step response", plot_width=350, plot_height=400)
d2.xaxis.axis_label = 'n (samples)'
d2.yaxis.axis_label = 'Amplitude'
d2.segment(range(len(step)), 0, range(len(step)), step,
line_width=2, line_color="SteelBlue")
d2.circle(range(len(step)), step,
size=5, fill_color="orange", line_color="SteelBlue", line_width=1)
bk.show(bk.gridplot([[d1, d2]]))
Visit www.add-for.com for more tutorials and updates.
This work is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.