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:

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:

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.

# # z_surf.py - 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("P4\n") 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] of.write(chr(byte)) # Close output file of.close()

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.

Hi Ted … Do you have the visualization of Riemann Zeta function?

Congratulations and have a nice day !

Hi Memo,

No, I’m afraid I’ve never programmed a visualisation for the Riemann Zeta function.

Ted