Filter design¶

Illustrate some signal filtering function from scipy.signal (http://docs.scipy.org/doc/scipy-0.14.0/reference/signal.html)

Similar to Matlab IIR filter design functions (http://mathworks.com/help/signal/ug/iir-filter-design.html)

1) Using the filter design functions¶

In [1]:
import numpy as np
import scipy.signal as sig

• iirdesign: http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.signal.iirdesign.html
• iirfilter: http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.signal.iirfilter.html
In [2]:
sig.iirdesign?


Example:

In [18]:
b,a = sig.iirdesign(wp = 0.001, ws=0.1,
gpass=3, gstop=40,
analog=False,
ftype='butter')
b, a

Out[18]:
(array([ 0.00157206,  0.00157206]), array([ 1.        , -0.99685589]))


2) Step response¶

Compute the step response, using lfilter (http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.signal.lfilter.html)

In [19]:
N = 2*10**3

x = np.zeros(N)
x[N//10:] = 1

y = sig.lfilter(b,a, x)

In [20]:
np.ones(3) * 2 +1

Out[20]:
array([ 3.,  3.,  3.])


Plot the step response¶

In [21]:
import matplotlib.pyplot as plt
%matplotlib inline

In [22]:
t = np.arange(N)
plt.plot(t, x, 'k')
plt.plot(t, y, 'b', linewidth=2)
plt.ylim(-0.1, 1.5)
plt.grid()


3) Interactive design¶

dict of available filter types

In [23]:
ftypes = {
'butter': 'Butterworth',
'cheby1': 'Chebyshev I',
'cheby2': 'Chebyshev II',
'ellip': 'Cauer/elliptic',
#'bessel': 'Bessel/Thomson' # order selection not available for Bessel
}


groupe the above code in one single function:

In [24]:
def plot_step_resp(wp_log = -3, ws_log=-1,
gpass=3, gstop=40,
ftype='butter'):
# 1) Design the filter
wp = 10**wp_log
ws = 10**ws_log
b,a = sig.iirdesign(wp, ws, gpass, gstop, False, ftype)
N_ord = max(len(b), len(a))
# 2) Compute the impulse response
y = sig.lfilter(b,a, x)
# 3) Plot
t = np.arange(N)
plt.plot(t, x, 'k')
plt.plot(t, y, 'b', lw=2)
plt.ylim(-0.1, 1.5)
plt.grid()
plt.title('{} filter of order {}, with\n wp=$10^{{{:.1f}}}$, ws=$10^{{{:.1f}}}$, Gp=-{:.1f} dB, Gs=-{:.1f} dB'.format(
ftypes[ftype], N_ord, wp_log, ws_log, gpass, gstop))


Test the function

In [26]:
plot_step_resp(wp_log = -3, ws_log= -2,
gpass=3, gstop=40,
ftype='butter')


Interactivity with IPython widgets¶

cf. widget doc: http://nbviewer.ipython.org/github/jvns/ipython/blob/master/examples/Interactive%20Widgets/Index.ipynb

In [28]:
from IPython.html.widgets import interact


For each parameter, we give an interval of allowed values

In [29]:
interact(plot_step_resp,
wp_log=[-3, -1.],
ws_log=[-2., 0],
gpass=[0.5,6], gstop=[10, 100], ftype=ftypes.keys());


4) Extra: plot the poles and zeros¶

compute the poles and zeros of the transfer function:

tf2zpk http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.signal.tf2zpk.html

In [77]:
sig.tf2zpk?

In [123]:
def plot_roots(wp_log = -3, ws_log=-1,
gpass=3, gstop=40,
ftype='butter', zoom=False):
# 1) Design the filter
wp = 10**wp_log
ws = 10**ws_log
b,a = sig.iirdesign(wp, ws, gpass, gstop, False, ftype)
N_ord = max(len(b), len(a))
# 2) Compute the zeros and poles
z, p, k = sig.tf2zpk(b, a)
# 3) Plot
theta = np.linspace(0, np.pi*2, 500)
circle = np.exp(1j*theta)
plt.plot(circle.real, circle.imag, 'k')
plt.plot(z.real, z.imag, 'c*', ms=10)
plt.plot(p.real, p.imag, 'rD')
ax = plt.gca()
ax.set(aspect='equal',
xlim=(-1.1, 1.1), ylim=(-1.1, 1.1),
)
if zoom:
ax.set(xlim=(0.9, 1.1), ylim=(-0.1, 0.1))
plt.title('{} filter of order {}, with\n wp=$10^{{{:.1f}}}$, ws=$10^{{{:.1f}}}$, Gp=-{:.1f} dB, Gs=-{:.1f} dB'.format(
ftypes[ftype], N_ord, wp_log, ws_log, gpass, gstop))

plot_roots(wp_log = -3, ws_log= -1,
gpass=3, gstop=40,
ftype='cheby1')

In [124]:
interact(plot_roots,
wp_log=[-3, -1.],
ws_log=[-2., 0],
gpass=[0.5,6], gstop=[10, 100], ftype=ftypes.keys());

In []: