% Example of IIR filter design by bilinear Z transform.
% Example 7.6 on page 391 of Ifeachor and Jervis text.
% Rich Kozick, Spring 1997

Fs = 1.28e3;	% Sampling frequency
T = 1/Fs;
wdc = 2*pi*150/Fs;	% Desired normalized digital cutoff freq.
f = linspace(0, Fs/2, 201);	% Frequencies for plots, in Hz
w = 2*pi*f;

% Pre-warp the desired digital frequency
wac = (2/T) * tan(wdc/2);

% Get analog filter frequency response
% NO FREQUENCY WARPING, so analog and digital filter 
% frequency response plots have same cutoff.
wac_true = wdc*Fs;
[Ha, Wa] = freqs(wac_true^2, [1, sqrt(2)*wac_true, wac_true^2], w);
HadB = 20*log10(abs(Ha));

% Get digital filter frequency response
Hd = freqz([1, 2, 1]*0.0878, [1, -1.0048, 0.3561], f, Fs);
HddB = 20*log10(abs(Hd));

% Plot analog and digital responses on same axis for comparison
figure(1)
plot(f, HadB, '-', f, HddB, '--')
legend('Analog Filter', 'Digital Filter')
xlabel('Frequency (Hz)')
ylabel('dB')
title('Magnitude response, Bilinear Z Transform method')
axis([0 Fs/2 -40 0])
