function pzd(b,a)
% pzd(b,a)
% produces a pole-zero diagram of B(s)/A(s)
	r=1;ar=[];br=[];
	if length(a)>1
		ar=roots(a);
		r=max(max(abs(ar),r));
	end	
% the following patch is necessary to correct some
% careless filter designs.  It eliminates
% very small leading numbers in b, which would produce
% large spurious zeros.
	c=b;s=(1e-13)*sum(abs(b));
	while abs(c(1)) < s
		c=c(2:length(c));
	end	 
	if length(c)>1
		br=roots(c);
		r=max(max(abs(br)),r);
	end
	ax=[-r r];bx=[0 0];
	plot(ax,bx,'k',bx,ax,'k')
	hold on
	plot(real(br),imag(br),'or',real(ar),imag(ar),'xb')
	hold off
	title('poles and zeros')
	axis([-1.05*r 1.05*r -1.05*r 1.05*r])
    axis square
    grid
    figure(gcf)