%PROGRAM p205
%bessel lowpass prototype
n=0:0.01:2;
N=10;   
[z,p,k]=besselap(N);
[b,a]=zp2tf(z,p,k);%Zero-pole to transfer function conversion.
[H,w]=freqs(b,a,n);%returns the complex frequency response vector H of the filter b/a:
magH2=(abs(H)).^2;
phaH=unwrap(angle(H));
phaH=phaH*180/pi;
subplot(211);
plot(w,magH);
axis([0 2 0 1]);%AXIS([XMIN XMAX YMIN YMAX]) sets scaling for the x- and y-axes on the current plot    
xlabel('w/wc');
ylabel('|H(jw)|^2');
title(['N=',num2str(N)]);
grid
%display phrase_freq figure
subplot(212)    
plot(w,phaH);    
xlabel('w/wc');
ylabel('phase');
grid