Systems & Signals practical, investigating z-domain transfer functions.

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

Prac4.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.pyplot as plt
from scipy import signal
from pylab import unwrap
from plot_zplane import zplane

def H(theta_1 = 3 * np.pi / 4.0,
      theta_2 =     np.pi / 4.0,
      r_1     = 0.95,
      r_2     = 0.95):
  return np.array([[1, -2 * np.cos(theta_1) * r_1, r_1**2],
          [1, -2 * np.cos(theta_2) * r_2, r_2**2]])

def helper1bc(num, den):
  # Calculate Magnitude response
  freqs, response = signal.freqz(num, den)

  # Calculate phase response
  # Unwrap uses small deltas to tell when the result of arctan has been wrapped
  phase = unwrap(np.arctan2(np.imag(response), np.real(response)))

  # Calculate impulse response
  impulse = np.zeros(50)
  impulse [0] = 1
  t = np.arange(50)
  iresponse = signal.lfilter(num, den, impulse)

  # 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]))
  phase = np.hstack((phase, phase[::-1]))
  
  # Plot Magnitude
  plt.subplot(311)
  plt.xlabel("Frequency (Hz)")
  plt.ylabel("Gain (DB)")
  plt.title("Frequency Response")
  plt.plot(freqs, mag)

  # Plot Phase
  plt.subplot(312)
  plt.xlabel("Frequency (Hz)")
  plt.ylabel("Phase (radians)")
  plt.title("Phase Response")
  plt.plot(freqs, phase)

  # Plot response
  plt.subplot(313)
  plt.ylabel("Amplitude")
  plt.xlabel("n (samples)")
  plt.title("Impulse Response")
  plt.stem(t, iresponse)

  #Show
  plt.subplots_adjust(hspace=0.5)
  plt.show()

# Question 1

def q1a():
  # zplane takes two arguments, the numerator and the denominator, each being lists
  # star extracts elements from a list and uses them as arguments to a function
  zplane(*H())

def q1b():
  for r_1 in [0.0, 0.5, 0.8, 1.0, 1.05]:
    # Get transfer function
    num, den = H(r_1=r_1)
    helper1bc(num, den)

def q1c():
  for r_2 in [0.0, 0.5, 0.8, 1.0, 1.05]:
    # Get transfer function
    num, den = H(r_2=r_2)
    helper1bc(num, den)

def q1d():
  for theta_1 in np.pi * np.array([0.0, 0.25 , 0.5, 0.75, 1]):
    # Get transfer function
    num, den = H(theta_1=theta_1)
    helper1bc(num, den)
 
def q1e():
  for theta_2 in np.pi * np.array([0.0, 0.25 , 0.5, 0.75, 1]):
    # Get transfer function
    num, den = H(theta_2=theta_2)
    helper1bc(num, den)

def q1f():
  num, den = H(r_2=1.0)
  t = np.arange(100)
  for i, w in enumerate(np.array([0.22, 0.25, 0.27]) * np.pi):
    inp = np.sin(t * w)
    response = signal.lfilter(num, den, inp)
    plt.subplot(3, 1, (i + 1))
    plt.ylabel("Amplitude")
    plt.xlabel("n (samples)")
    plt.title("w = " + str(w) +" (rad/sample)")
    plt.stem(t, response)

  #Show
  plt.subplots_adjust(hspace=0.5)
  plt.show()

# Question 2
def q2():
  num = [.0038, .0001, .0051, .0001, .0038]
  den = [1., -3.2821, 4.236, -2.5275, .5865]
  # Q2a
  helper1bc(num, den)
  
  # Q2b
  t = np.arange(100)
  for i, w in enumerate(np.array([0.02, 0.1, 0.25, 0.4]) * np.pi):
    inp = np.sin(t * w)
    response = signal.lfilter(num, den, inp)
    plt.subplot(4, 1, (i + 1))
    plt.ylabel("Amplitude")
    plt.xlabel("n (samples)")
    plt.title("w = " + str(w) +" (rad/sample)")
    plt.stem(t, response)
 
  #Show
  plt.subplots_adjust(hspace=0.5)
  plt.show()

  # Q2c
  zplane(np.array(num), np.array(den))

q1a()
q1b()
q1c()
q1d()
q1e()
q1f()
q2()