% quantdemo.m
% Simple script to demonstrate effects of coefficient quantization
% on direct-form and cascade-form IIR filters
% Last revised 10/18/99 - rms

% First section:  Compare DTFTs of direct form realization with and without
% quantized coefficients



% Generating coefficients and frequency response for direct-form
% realization with unquantized coefficients
poles = [.95*exp(j*pi/4) 0.95*exp(-j*pi/4) .9*exp(j*pi/4) 0.9*exp(-j*pi/4) ];
zeros = [.9*exp(j*.95*pi/4) 0.9*exp(-j*.95*pi/4) .9*exp(1.05*j*pi/4) 0.9*exp(-1.05*j*pi/4) ];

%Uncomment the next line to check pole-zero patterns
subplot(1,1,1)
zplane0(zeros',poles') 
pause
  
for N = 12:-4:8,
  
b=real(poly(zeros));
a = real(poly(poles));
[h w] = freqz(b,a);
iN = num2str(N);
subplot(2,1,1),plot(w/pi,20*log10(abs(h)))
title(['Ideal and Direct-Form response, N = ',iN],'FontSize',14)
ylabel('|H(\omega)|','fontsize',14)
xlabel('\omega/\pi','fontsize',14)
hold on
subplot(2,1,2),plot(w/pi,angle(h))
ylabel('Angle(H(\omega))','fontsize',14)
xlabel('\omega/\pi','fontsize',14)
hold on

% Generating coefficients and frequency response for direct-form
% realization with quantized coefficients
maxa= max(abs(a)); maxb = max(abs(b));
a2 = maxa*fxquant((a/maxa),N,'round','sat');
b2 = maxb*fxquant((b/maxb),N,'round','sat');
[a' a2' b' b2']
[h2 w] = freqz(b2,a2);
subplot(2,1,1),plot(w/pi,20*log10(abs(h)))
subplot(2,1,1),plot(w/pi,20*log10(abs(h2)),'r')
subplot(2,1,2),plot(w/pi,angle(h2),'r')
hold off
pause

subplot(2,1,1), plot(0,0)
subplot(2,1,2), plot(0,0)
clear h2
subplot(1,1,1)
zplane0(zeros',poles') 

% Generating coefficients for cascade-form
% realization with unquantized coefficients
[b0,B,A] = dir2cas(b,a);

% Generating coefficients for cascade-form
% realization with quantized coefficients
maxA= max(max(abs(A))); maxB = max(max(abs(B)));
A2 = maxA*fxquant((A./maxA),N,'round','sat');
B2 = maxB*fxquant((B./maxB),N,'round','sat');
[A A2];
[B B2];


% Display results of casecade-form realizations
A2_ALL = conv(A2(1,:),A2(2,:));
B2_ALL = conv(B2(1,:),B2(2,:));

subplot(2,1,1),plot(w/pi,20*log10(abs(h)))
title(['Ideal and Cascade-Form response, N = ',iN],'FontSize',14)
ylabel('|H(\omega)|','fontsize',14)
xlabel('\omega/\pi','fontsize',14)
hold on
subplot(2,1,2),plot(w/pi,angle(h))
ylabel('Angle(H(\omega))','fontsize',14)
xlabel('\omega/\pi','fontsize',14)
hold on
[h3,w] = freqz(B2_ALL,A2_ALL);
subplot(2,1,1),plot(w/pi,20*log10(abs(h3)),'r')
subplot(2,1,2),plot(w/pi,angle(h3),'r')


pause

end