fs = 1e6;
fp = 20e3;
Qp = 100;

T = 2*pi*fp/fs;
r = 1 - pi*fp/fs/Qp;
cT2 = 2 * cos(T);

z = tf('z', 1/fs);

H = (z^-2 - z^-1 * cT2 + 1) / (r^2 * z^-2 - z^-1 * r * cT2 + 1);

flin = linspace(0, fs/2, 1001);
wlin = 2*pi*flin;
f = flin/fs;
Mlin = bode(H, wlin);
Mlin = (Mlin(1,:));
Mlin(Mlin < 1e-4) = 1e-4;

figure(1);
plot(flin, Mlin);

H1a = 1 - z^-1 * cT2;
H1b = 1 - z^-1 * cT2 + z^-2;
H1c = 1 - z^-1 * cT2 + z^-2 / ( 1 - z^-1 * r * cT2 );
M1a = bode(H1a, wlin);   M1a = 20*log10(M1a(1,:));
M1b = bode(H1b, wlin);   M1b = 20*log10(M1b(1,:));
M1c = bode(H1c, wlin);   M1c = 20*log10(M1c(1,:));

figure(2);
plot(f, M1a, f, M1b, f, M1c, f, 20*log10(Mlin));
legend('Partial H_{1a}', 'Partial H_{1b}', 'Partial H_{1c}', 'Filter H_1', 4);
xlabel('Frequency  [ f / f_s ]');
ylabel('Amplitude  [ dB ]');
fixfig;

H2a = 1 / ( 1 - z^-1 * r * cT2 );
H2b = 1 / ( 1 - z^-1 * r * cT2 + r^2 * z^-2);
H2c = 1 - z^-1 * cT2 / ( 1 - z^-1 * r * cT2 + r^2 * z^-2);
M2a = bode(H2a, wlin);   M2a = 20*log10(M2a(1,:));
M2b = bode(H2b, wlin);   M2b = 20*log10(M2b(1,:));
M2c = bode(H2c, wlin);   M2c = 20*log10(M2c(1,:));

figure(3);
plot(f, M2a, f, M2b, f, M2c, f, 20*log10(Mlin));
legend('Partial H_{2a}', 'Partial H_{2b}', 'Partial H_{2c}', 'Filter H_2', 4);
xlabel('Frequency  [ f / f_s ]');
ylabel('Amplitude  [ dB ]');
fixfig;