function plotZTP(b,a,nu,tmax)
% plotZTP(b,a,nu,tmax)
% Display Z transform pairs: plot
% (1) pole zero diagram,
% (2) frequency response
% (3) pulse response on [-tmax,tmax]
% h(k) <--> H(z)=(z^nu)*B(z)/A(z),
% B(z)=b0+b1*z^(-1)+...+bm*z^(-m), etc.
% b=[b0,b1,...,bm],a=[1,a1,a2,...,an]
% b0,bm,an should all be nonzero
% This function requires two other 
% homebrew functions, 'gpfx.m' and 'pzd.m' 
	T=max(30,floor(tmax));
	
	subplot(2,2,2) % ** pole-zero **
	pzd(b,a)
	hold on
	L=100;theta=(2*pi/L)*[0:L];z=exp(j*theta);
	plot(real(z),imag(z),'b')
	axis('off')
	axis('square')
	hold off
	
	subplot(2,2,1) % ** magnitude of H **
	L=512;[H,W]=freqz(b,a,L);H=exp(j*nu*W).*H;
	mag=abs(H);mx=max(mag);
	plot((1/L)*[0:L-1],mag,'b');
	axis([0 1 0 mx+.05])
	title('magnitude of H(z)')
	grid
	
	subplot(2,2,3) % ** angle of H **
	mu=angle(H)/pi;% give preference to pi over -pi
	plot((1/L)*[0:L-1],mu-sign(mu+.999)+1,'r');
	axis([0 1 -1.1 1.1])
	title('phase/pi of H(z)')
	grid
	
	subplot(2,2,4) % ** h(k) **
    [bp,ap,bn,an]=gpfx(b,a,'abs(z)<1'); % separate
	h=zeros(1,2*T+1);
	if max(length(bp),length(ap))>=1; %causal part
		h([T+1:2*T+1])=filter(bp,ap,[1,zeros(1,T)]);
	end	
	if max(length(bn),length(an))>1; %anticausal part
		[bn,an,nn]=z_rev(bn(2:length(bn)),an,-1);
		hn=filter(bn,an,[1,zeros(1,T)]);
		h(1-nn:T+1-nn)=h(1-nn:T+1-nn)+flipud(hn')';
	end	
	mn=min(h);[hm,tm]=max(h);d=.1*(hm-mn+.01);mn=mn-d;mx=hm+2*d;
	t0=max(1,2*nu+1);t2=2*T+min(1,2*nu+1);t=[t0:t2];t1=(t0+t2)/2;
	stem(t,h(t),'k')
	axis([t0 t2 mn mx])
	axis('off')
	hold on;plot([t1 t1],[mn mx],'w');
	title('time domain sequence \{h[k]\}')
	dt=3;plot(tm+dt*[-1 1],[hm,hm],'c',[tm tm],[mn 0],'c')
	text(tm+1+dt,hm+d,num2str(hm))
	text(tm-1,mn-d,num2str(tm-t1))
	hold off