ch16ex7.m

% Chapter 16 Exercise 7.
P=8; wn=20; Zi = 0.65;
num      = wn^2;
den      = conv([1/P 1],[1 2*Zi*wn wn^2]);

oms      = logspace(-1,2);
s        = sqrt(-1)*oms;

%--- compute magnitude and phase of the analog signal
H        = polyval(num,s)./polyval(den,s);
magHdb   = 10*log10(H.*conj(H));
phaHdeg  = 180/pi * unwrap(angle(H));

answer   = input('Please enter sampling frequency fs, or 0 to exit --> ');
while (answer ~= 0)
   fs = answer;

   %--- compute magnitude and phase of the digital equivalent
   [numd, dend]   = bilinear(num,den,fs);
   z          = exp(s/fs);
   Hd         = polyval(numd,z)./polyval(dend,z);
   magHddb    = 10*log10(Hd.*conj(Hd));
   phaHddeg   = 180/pi * unwrap(angle(Hd));

   %--- plot analog and digital-equivalent magnitude
   subplot(2,1,1);
   semilogx(oms, magHdb, ...
            oms, magHddb,'b:');   grid;
   str        = ['Sampling frequency: ',num2str(fs), ' Herz'];
   title(str);

   %--- plot analog and digital-equivalents phase
   subplot(2,1,2);
   semilogx(oms, phaHdeg,       ...
            oms, phaHddeg,'b:');  grid

   answer     = input('Please enter sampling frequency fs, or 0 to exit --> ');
end; % of while loop

disp('The variables fs, numd, dend contain the sampling frequency');
disp('and the coefficients of the equivalent digital filter.     ');