0.5dBチェビシェフフィルタのインパルス応答

clear; close all;

s = tf('s');

xi = 0.1;
e = sqrt(10.^(xi/10)-1);
tmp = sqrt(1+e^(-2))+e^(-1);
theta = 0:0.1:2*pi;
x = cos(theta);
y = sin(theta);

for n=1:9
    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);
        H = 1/D;
    end
    coef(D)
    subplot(3,3,k);
    [num,den]=tfdata(D);
    num = cell2mat(num);
    den = cell2mat(den);
    num = real(num);
    D = tf(num, den);
    H = 1/D;
    impulse(H);
    xlim([0 25]);
end
big;


function n = coef(D)
[n,d] = tfdata(D);
n = cell2mat(n);
n = n(2:end);
n = real(n);
end