双一次Z変換によるチェビシェフIIRフィルタの設計

clear; close all;

syms w;
chebyshev = @(c1,c0) 2*w*c1 - c0;
c0 = 1;
c1 = w;
cs = [c0 c1];
for n=2:10
    cn_prev = cs(n-1);
    cn = cs(n);
    cn_next = chebyshev(cn, cn_prev);
    cs(end+1) = cn_next;
end

n = 8;
xis = [0.1 0.5 1.0];
syms s z;
T = 1e-4; %s
for idx=1:length(xis)
    xi = xis(idx);
    e = sqrt(10.^(xi/10)-1);

    tmp = sqrt(1+e^(-2))+e^(-1);
    D = 1;
    for k=1:n
        sinhv = 1/2*(tmp^(1/n)-tmp^(-1/n));
        coshv = 1/2*(tmp^(1/n)+tmp^(-1/n));
        sigma_k = -sin((2*k-1)/(2*n)*pi)*sinhv;
        omega_k = cos((2*k-1)/(2*n)*pi)*coshv;
        sk = sigma_k + i*omega_k;
        D = D*(s-sk);
    end

    G = vpa(1/D);

    H_gain = 1/(1+e^2*cs(n+1)^2);
    omega = linspace(0,2,100);
    figure(1);
    subplot(211);
    hold on;
    xistr = sprintf('ξ=%.3gdB',xi);
    plot(omega, subs(H_gain, w, omega),'DisplayName',xistr);
    legend();
    
    subplot(212);
    hold on;
    s_by_z = 2/T*(1-z^(-1))/(1+z^(-1));
    Hz = subs(G, s, s_by_z);
    z_exp_jwt = exp(j*(omega)*T);
    Hz_abs = abs(subs(Hz,z,z_exp_jwt));    
    plot(omega, Hz_abs/max(Hz_abs),'DisplayName',xistr);
    
    hold on;
end
nstr=sprintf('Chebyshev filter (n=%.3g)',n);
subplot(211); legend; title(nstr);
nstr=sprintf('等価IIRフィルタ (n=%.3g)',n);
subplot(212); legend; title(nstr);

big;

ゲインが変わってしまっているのが気になる