plotzmag.m

biomedical stuffs...

function dummy=plotzmag(nz,dz,rmin,rmax,thetamin,thetamax,nres)
%	plotzmag(nz,dz,rmin,rmax,thetamin,thetamax,nres)
%	rmin, rmax: min and max values for radius
%	thetamin, thetamax: min and max values for angle (usually -pi and pi)
%	nz = vector of numerator polynomial coefficients
%	dz = vector of denominator polynomial coefficients
%	nz and dz are in ascending powers of z^(-1)
%	nres = number of points for r (also used for theta)
%
%	figure(1) is full scale;  figure(2)  has z-axis magnified
%
%	NOTE: Relative heights of peaks at poles can be misleading
%	due to the sampling of the z-plane imposed by the grid
%	of nres X nres points.
%
%	NOTE: Use 'rotate3d' command to rotate plot after it is drawn!
%
ln=length(nz);ld=length(dz);i=sqrt(-1);
dr=(rmax-rmin)/(nres-1);dth=(thetamax-thetamin)/(nres-1);
zmag=zeros(nres,nres);
num_zn=zeros(1,ln);den_zn=zeros(1,ld);
for m=1:ln, num_zn(m)=1-m;end
for m=1:ld, den_zn(m)=1-m;end
disp('The poles of your Z-transform are:')
roots(dz)
disp('The zeros of your Z-transform are:')
roots(nz)
for j=0:nres-1
   r=rmin+j*dr;
   for k=0:nres-1
      th=thetamin+k*dth;
      zv=r*exp(th*i);
      xarg(k+1,j+1)=r*cos(th);yarg(k+1,j+1)=r*sin(th);
      zmag(k+1,j+1)=abs(sum(nz.*(zv.^num_zn))/sum(dz.*(zv.^den_zn)));
   end
end
zgain=10*(nz(1)/dz(1));
if nz(1)==0, zgain=10;end
zmagpl=min(zmag,zgain*ones(size(zmag)));
figure(1),surf(xarg,yarg,zmag)
xlabel('Re{z}');ylabel('Im{z}');zlabel('Mag(X(z))')
colormap('white')
figure(2),dummy=surf(xarg,yarg,zmagpl);
xlabel('Re{z}');ylabel('Im{z}');zlabel('Mag(X(z))')
colormap('white')