Visualising transfer functions in the Z-domain

I’m fortunate to share an office with DSP guru Richard Hayes who, over the last week or so, has given me a couple of brilliant 20-minute tutorials on the Z-transform which have really helped to clear up some points that have been confusing me for years. In fact, I’m ashamed to say that my knowledge of the Z-transform has always been a bit shaky, trailing far behind my modest level of understanding of the Fourier and Laplace transforms.

I first encountered all of these transforms as an undergraduate engineer in UCD back in the 1990s, but I really only learned them off well enough at that time to get me through the exams. Although I could sort of understand the basic principle of each transform, my intuitive understanding of the mathematics was totally inadequate. It wasn’t until I was doing my masters in UCD’s Rehabilitation Engineering lab around 1999 that I first developed the deeper intuitive grasp of these techniques which one requires in order to apply them with confidence, or to apply them creatively. My point of entry to real, meaningful understanding was the Discrete Fourier Transform (DFT) which I learned inside out by writing my own (very inefficient) C implementation. Although the DFT can be thought of as a special case of the Z-transform, I remained mostly baffled by the Z-domain even after that. In some ways, the Z-Transform equation looks innocuous enough:

X(z) = \mathscr{Z}\left\{{x[n]}\right\} = \sum\limits_{n=-\infty}^{\infty} x[n]z^{-n}

However, trying to actually visualise what happens when a signal (or system) is transformed from the time domain into the Z-domain is (for most people, including me) very challenging.

What Richard was working on last week was a MATLAB program to generate a 1-page summary of discrete-time system analysis, which included a 3D plot of a system transfer function plotted as a surface in the Z-domain. He talked me through his MATLAB code line-by-line, which really helped me to picture the meaning of the surface. Also, partly by accident, he had generated a spectacular illustration of the Z-domain region of convergence (his surface shot straight up at the boundary of this region).

Although I’ve seen lots of Z-domain surface plots of this type before, I hadn’t ever generated one myself programmatically, or read through a program that generated one. Inspired by Richard’s example, I thought I’d have a go at it myself over the weekend, programming in Python instead of MATLAB. The result is shown below:

Animated Z-domain surface of a transfer function

The code I used to generate the above animation is shown below. Please note that this code is very far from being efficient. I originally wrote it to generate a single image, but then subsequently hacked it to produce 360 images, so it recalculates the exact same Z-domain surface over and over again, once for each image.

# - Python program to visualise the Z-domain
# surface of a system transfer function. The program
# generates 360 separate images, each viewing the same
# surface from a different angle (1 degree increments).
# The output files are 000.pbm, 001.pbm, 002.pbm, etc.
# Written by Ted Burke
# Last updated 14-1-2013

import sys
import numpy
import math

for degrees in range(360):
	# convert angle to radians
	angle = math.pi*degrees/180.0
	# Array of pixels for plot
	p = numpy.zeros((400,400),int)

	# Poles and zeros of transfer function
	p1 = complex(0.8,0.5)
	p2 = complex(0.8,-0.5)
	z1 = complex(1.0, 0.0)
	z2 = complex(-1.0, 0.0)

	# Calculate a and b coefficients from poles and zeros
	a = numpy.convolve([1,-p1],[1,-p2])
	b = numpy.convolve([1,-z1],[1,-z2])

	# origin point and i,j,k unit vectors for 3D plot
	origin = numpy.array([200.0, 120.0])
	# Horizontal major axis 80, Vertical minor axis 30
	i = numpy.array([80.0*math.cos(angle), 30.0*math.sin(angle)])
	j = numpy.array([-80.0*math.sin(angle), 30.0*math.cos(angle)])
	k = numpy.array([0.0, 30.0])

	# Calculate magnitude of H over a square area on the Z-plane
	for im in numpy.arange(-1.5, 1.5, 0.05):
		for re in numpy.arange(-1.5, 1.5, 0.05):
			# Next z value
			z = complex(re,im)
			# Calculate value of H(z) for this z value
			H = (b[0]*z*z + b[1]*z + b[2])/(a[0]*z*z + a[1]*z + a[2])
			H_magn = abs(H)
			H_arg = math.atan2(H.imag, H.real)
			# Calculate corresponding coordinates on raster image
			point = origin + re*i + im*j + H_magn*k
			#point = origin + re*i + im*j + H_arg*k
			px,py = int(point[0]),int(point[1])
			# If point lies within image boundary, colour it in
			if px>=0 and px<400 and py>=0 and py<400:
				p[py,px] = 1

	# Output binary PBM image to a numbered filename
	print('Writing {0:03d}.pbm'.format(degrees))
	of = open('{0:03d}.pbm'.format(degrees), 'wb')
	# Write PBM file header
	of.write("# Z domain transfer function surface\n")
	of.write("400 400\n")
	# Write pixel data
	for y in range(399,-1,-1):
		for x in range(0,400,8):
			# Pack 8 pixels into a byte
			byte = 0
			for n in range(8):
				byte = (byte << 1) + p[y,x+n]
	# Close output file

The image files produced by the program above are called "000.pbm", "001.pbm", "002.pbm" and so on up to "359.pbm". Each file shows the same surface viewed from a different angle. The angle between successive views is 1 degree. I used ImageMagick’s convert command to combine the 360 separate PBM images into a single animated GIF, as follows:

convert -delay 2 *.pbm surf.gif

That command took a couple of minutes to run, but it did finally produce the animated GIF shown earlier in this post.

This entry was posted in Uncategorized and tagged , , , , , , , , , , . Bookmark the permalink.

2 Responses to Visualising transfer functions in the Z-domain

  1. Memo Pérez says:

    Hi Ted … Do you have the visualization of Riemann Zeta function?
    Congratulations and have a nice day !

Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s