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