fun with gnuplot prg (poles and zeros)

[Martin Czech]

All this talking about poles and zeros in the last time made me dig out
this old theory stuff, and I thought how this could be visualized. 
(OH boy, this is sooo long ago.)

Indeed, you can do this with a pencil and a piece of paper. 

This is nothing new to some people, but nevertheless:

Draw an s-plane, mark your poles and zeros (don't forget, if they
have an imaginary part, there always must be a complex conjugate (sp?)
counterpart (mirror image) (not like me, I first forgot one mirrored zero
and got bizzare results... ;->). Now, how to get an impression of H(s)?
(This is evaluation on the jw axis, we can concentrate on the POSITIVE
jw axis, because everything is mirrored on the negative part.)

First |H(s)|:

You are somewhere on the (positive!) jw axis and want to know the
magnitude at this point (or frequency). Draw lines (more exact vectors)
from every pole to your point, and measure the length. Multiply all these
pole-to-jw distances. Do the same for the zeros. Divide the zero product
by the pole product. Voila!. We can derive from this: a pole close to
jw gives a maximum (a peak) at frequency of the pole (ie. imaginary
part). The smaller the real part of the pole, the more resonant is the
peak. Comming close to a zero means magnitude loss, the frequency of the
"dip" beeing the frequency of the zero (or imag part) , and the amount
of the notch is controlled by the real part of the zero. A zero in your
jw axis gives 0 amplitude, the signal will totally vanish in this case.
If you have more poles then zeros, then the magnitude will fall for large
w. If more zeros, it will rise (not physical possible).
If the number of poles and zeros are equal it will get to some limit value
and stay there (also not physical possible).

Now this was easy, lets proceed to the phase angle.


You have allready drawn all these pole and zero vectors to your jw point.
Now measure the pole vector angles away from the positive real axis.
Ie. if the pole vector points in the same direction as the positive real
axis, the angle is 0. If the pole vector points allmost in the positive
jw axis direction, the angle is allmost 90 deg. Add up all pole angles.
Same procedure for the zero angles. Now subtract the sum of pole angles
from the sum of zero angles. Viola! You've got the phase angle for your
jw point.

>From this geometrical construction, we can see that no sudden phase
change is possible, phase is allways steady and if you have e.g. four
real poles, the total phase angle will come close to -360deg for large w.
There is no sudden jump from -180deg to +180deg.  These artefacts are only
confusion that comes up with stupid programs or measurement equipment.
The only exeption is a zero pair on the jw axis, first you have 0
deg, then suddenly -180 deg. This is a very sharp and deep notch filter,
so the magnitude is also 0 , or very small around the zero.  I think,
phase angle has no real meaning for magnitudes -> 0.
So the discontinuity is not harmfull for group delay and such.

For all gnuplot friends (I think there are versions for all OS, it is
a good tool to visualize formulas or raw data) I've written some lines
that compute |H(s)| and angle{H(s)} for you.  There are functions for
single and double poles and zeros.  This is convenient, because you have
to type in double poles and zeros only once.

(Erata: zeros can't have real part = 0.0 because of division, use a very
little number instead : 0.0000001)


#constant for radian to degrees conversion
WI=360.0/2.0/pi

pr1=-1.0
pi1=1.0

nr1=0.00001
ni1=3.0


#function for magnitude of single (real) pole
ps(x,a)=1.0/sqrt(a**2+x**2)
#function for magnitude of double pole
pd(x,a,b)=1.0/sqrt(a**2+(x-b)**2)/sqrt(a**2+(x+b)**2)
#function for magnitude of single (real) zero
zs(x,a,b)=sqrt(a**2+x**2)
#function for magnitude of double zero
zd(x,a,b)=sqrt(a**2+(x-b)**2)*sqrt(a**2+(x+b)**2)

#function for angle of single (real) pole
vs(x,a)=atan(x/a)*WI
#function for magnitude of double pole
vd(x,a,b)=(atan((x-b)/a)+atan((x+b)/a))*WI
#function for angle of single (real) zero
ws(x,a)=-atan(x/a)*WI
#function for magnitude of double zero
wd(x,a,b)=-(atan((x-b)/a)+atan((x+b)/a))*WI


#magnitude example in dB, note: allways multiply!
h(x)=20*log10(pd(x,pr1,pi1)*zd(x,nr1,ni1))

#phase example, note: allways add
p(x)=vd(x,pr1,pi1)+wd(x,nr1,ni1)


set logscale x

plot [0.1:10]  p(x)
#plot [0.1:10] h(x)