#!/usr/bin/env python
# coding: utf-8

# # Scipy Signal Processing: IIR Filter Design

# In[1]:


import addutils.toc ; addutils.toc.js(ipy_notebook=True)


# In[2]:


from addutils import css_notebook
import bokeh.plotting as bk
css_notebook()


# In[3]:


bk.output_notebook()


# ## 1 Filter Specification

# 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:
# * [Signal Processing (scipy.signal)](http://docs.scipy.org/doc/scipy/reference/signal.html#module-scipy.signal)

# A filter is usually defined in the frequency domain.  The following figure depicts the frequency domain specification:

# <a name="arrayobject"><img src="files/utilities/filter_design.jpg"  style="max-width: 100%"/>

# The specification figure illustrates a low-pass filter but the terminology applies to all filter types, lowpass, highpass, bandpass, and stopband.
# 
# * **Passband (Wp):** This is the frequency range which we desire to let the signal through with minimal attenuation. In the scipy functions this is in normalized frequency, 1 > Wp > 0.
# * **Stopband (Ws):** This is the frequency range which the signal should be attenuated. Again, in scipy this is in normalized frequency, 1 > Ws > 0.
# * **Passband ripple (Rp, gpass):* The max variation in the passband, in decibals.
# * **Stopband attenuation (As, gstop):** The max level in the stopband or in other words the min attenuation in the stopband, in decibals.
# 
# 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:

# <a name="arrayobject"><img src="files/utilities/IIR_filters_table.jpg"  style="max-width: 100%"/>

# ## 2 Filter Design

# In[4]:


import numpy as np
from scipy import signal
from scipy.signal import filter_design as fd

import bokeh.plotting as bk
bk.output_notebook()


# In[5]:


# 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')


# ## 3 Frequency and Phase Response

# To define plot colors, any [legal HTML name for colors](http://en.wikipedia.org/wiki/Web_colors) can be used:

# In[6]:


plotcolors = ['DarkRed', 'Crimson', 'Red',
              'Olive', 'Green', 'SeaGreen',
              'SteelBlue', 'Blue', 'Navy',
              'Fuchsia', 'DarkViolet', 'Indigo']


# Plotting:
# 
# * `a, b` are respectively numerator and denominator of the filter
# * `w` are the frequencies at which the frequency response `'h'` is computed

# In[7]:


xmn1, 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:

# In[8]:


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)


# ## 4 Step and Impulse Response

# We can have a look to the inpulse and step response for any filter (bessel in this example):

# In[9]:


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](<http://www.add-for.com/IT>) for more tutorials and updates.
# 
# This work is licensed under a <a rel="license" href="http://creativecommons.org/licenses/by-sa/4.0/">Creative Commons Attribution-ShareAlike 4.0 International License</a>.