Prac5 Examining filters

plot_zplane.py

#
# Copyright (c) 2011 Christopher Felton
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.
#

# The following is derived from the slides presented by
# Alexander Kain for CS506/606 "Special Topics: Speech Signal Processing"
# CSLU / OHSU, Spring Term 2011.

import numpy as np
import matplotlib.pyplot as plt
from  matplotlib import patches
from matplotlib.figure import Figure
from matplotlib import rcParams
    
def zplane(b,a,filename=None):
    """Plot the complex z-plane given a transfer function.
    """

    # get a figure/plot
    ax = plt.subplot(111)

    # create the unit circle
    uc = patches.Circle((0,0), radius=1, fill=False,
                        color='black', ls='dashed')
    ax.add_patch(uc)

    # The coefficients are less than 1, normalize the coeficients
    if np.max(b) > 1:
        kn = np.max(b)
        b = b/float(kn)
    else:
        kn = 1

    if np.max(a) > 1:
        kd = np.max(a)
        a = a/float(kd)
    else:
        kd = 1
        
    # Get the poles and zeros
    p = np.roots(a)
    z = np.roots(b)
    k = kn/float(kd)
    
    # Plot the zeros and set marker properties    
    t1 = plt.plot(z.real, z.imag, 'go', ms=10)
    plt.setp( t1, markersize=10.0, markeredgewidth=1.0,
              markeredgecolor='k', markerfacecolor='g')

    # Plot the poles and set marker properties
    t2 = plt.plot(p.real, p.imag, 'rx', ms=10)
    plt.setp( t2, markersize=12.0, markeredgewidth=3.0,
              markeredgecolor='r', markerfacecolor='r')

    ax.spines['left'].set_position('center')
    ax.spines['bottom'].set_position('center')
    ax.spines['right'].set_visible(False)
    ax.spines['top'].set_visible(False)

    # set the ticks
    r = 1.5; plt.axis('scaled'); plt.axis([-r, r, -r, r])
    ticks = [-1, -.5, .5, 1]; plt.xticks(ticks); plt.yticks(ticks)

    if filename is None:
        plt.show()
    else:
        plt.savefig(filename)
    

    return z, p, k

Prac5.py

# Copyright Robert Ketteringham 2013
# Licensed under WTFPL
# With help from:
#   http://www.scipy.org/Cookbook/FIRFilter
#   http://mpastell.com/2009/11/05/iir-filter-design-with-python-and-scipy/
import numpy as np
import matplotlib.pylab as plt
from scipy import signal
from plot_zplane import zplane

# Setup
x0 = lambda t: np.cos(50 * 2 * np.pi * t)
x1 = lambda t: np.cos(63 * 2 * np.pi * t)
x2 = lambda t: np.abs(x0(t))
x = lambda t: x1(t) + x2(t)
fs = 1000
p = 1. / fs
t = np.arange(0, 1, p)

# Question 1
def Q1():
    plt.subplot(321)
    plt.title("Magnitude of X1")
    plt.stem(t[:150], np.abs(x1(t[:150])))
    plt.subplot(322)
    plt.title("Spectra of X1")
    freq = np.fft.fft(x1(t))
    plt.plot(np.abs(freq))
    
    plt.subplot(323)
    plt.title("Magnitude of X2")
    plt.stem(t[:150], np.abs(x2(t[:150])))
    plt.subplot(324)
    plt.title("Spectra of X2")
    freq = np.fft.fft(x2(t))
    plt.plot(np.abs(freq))
    
    plt.subplot(325)
    plt.title("Magnitude of X")
    plt.stem(t[:150], np.abs(x(t[:150])))
    plt.subplot(326)
    plt.title("Spectra of X")
    freq = np.fft.fft(x(t))
    plt.plot(np.abs(freq))
    
    plt.show()

def Q3(num, den):  
    # Prototype Filter
    freqs, response = signal.freqz(num, den)
    # Normalize axis
    freqs = freqs / np.pi # rad -> Hz
    mag = 20 * np.log10(abs(response)) # DB
    # Append reversed results (for symmetry)
    # Normal list addition doesn't work because numpy tries to add the values of arrays
    freqs = np.hstack((freqs, freqs[-1] + freqs))
    mag = np.hstack((mag, mag[::-1]))
  
    # Plot Magnitude
    plt.xlabel("Frequency (Hz)")
    plt.ylabel("Gain (DB)")
    plt.title("Frequency Response")
    plt.plot(freqs, mag)
    plt.show()

prototype = [[1, -1], [1, -.95]]
comb = [[1, 3.331e-016, 8.327e-016, -3.331e-016, -4.996e-015, -3.997e-015, 5.551e-016, 1.443e-015, 5.551e-017, -1.11e-016, -1],
        [1, 0, -8.327e-016, 2.22e-016, 2.22e-016, -1.443e-015, 2.109e-015, 1.388e-015, -1.638e-015, -2.22e-016, -0.5987]]

def Q3a():
    Q3(*prototype)

def Q3b():
    Q3(*comb)

def Q4a():
    zplane(*prototype)

def Q4b():
    zplane(*comb)

def Q5():
    inp = x(t)
    response = signal.lfilter(comb[0], comb[1], inp)
    #print response
    plt.subplot(211)
    plt.title("Magnitude of Y")
    plt.stem(t[:150], response[:150])
    plt.subplot(212)
    plt.title("Spectra of Y")
    freq = np.fft.fft(response)
    plt.plot(np.abs(freq))
    plt.show()

Q5()