clc;
clear all;


%% without prewarp
[z, p, k] = buttap(3);
[num, den] = zp2tf(z, p, k);
wc = 15;
[num, den] = lp2lp(num, den, wc);
T = 0.1;
[numd, dend] = bilinear(num, den, 1/T)
w = -pi:2*pi/300:pi;
z = exp(j*w);
H = (numd(1)*z.^3 + numd(2)*z.^2 + numd(3)*z + numd(4)) ./ (dend(1)*z.^3 + dend(2)*z.^2 + dend(3)*z + dend(4));
subplot(2,1,1);
plot(w, abs(H));
title('without prewarp');


%% with prewarp
[z, p, k] = buttap(3);
[num, den] = zp2tf(z, p, k);
wp = 2*tan(wc*T/2)/T;
[num, den] = lp2lp(num, den, wp);
T = 0.1;
[numd, dend] = bilinear(num, den, 1/T)
w = -pi:2*pi/300:pi;
z = exp(j*w);
H = (numd(1)*z.^3 + numd(2)*z.^2 + numd(3)*z + numd(4)) ./ (dend(1)*z.^3 + dend(2)*z.^2 + dend(3)*z + dend(4));
subplot(2,1,2);
plot(w, abs(H));
title('with prewarp');