第１種チェビシェフフィルタの位相特性

clear; close all;

s = tf('s');
syms w;


xis = [0.1 0.5 1.0];

n = 5;
for idx=1:length(xis)
    xi = xis(idx);
    e = sqrt(10.^(xi/10)-1);
    tmp = sqrt(1+e^(-2))+e^(-1);
    sinhv = 1/2*(tmp^(1/n)-tmp^(-1/n));
    coshv = 1/2*(tmp^(1/n)+tmp^(-1/n));

    H_gain = 1;
    H_phi = 0;
    for k=1:n
        sigma_k = -sin((2*k-1)/(2*n)*pi)*sinhv;
        omega_k = cos((2*k-1)/(2*n)*pi)*coshv;
        H_gain = H_gain*1/sqrt(sigma_k^2+(omega_k-w)^2);
        H_phi = H_phi - atan((omega_k-w)/sigma_k);
    end
    nstr = sprintf('xi=%.3gdB',xi);
    
    omega = 0:0.05:2;
    figure(1);
    plot(omega, rad2deg(subs(H_phi, w, omega)),'DisplayName',nstr);
    hold on;
end
nstr = sprintf('位相特性 (n=%.3g)',n);
figure(1); legend;title(nstr);

n = 8;
for idx=1:length(xis)
    xi = xis(idx);
    e = sqrt(10.^(xi/10)-1);
    tmp = sqrt(1+e^(-2))+e^(-1);
    sinhv = 1/2*(tmp^(1/n)-tmp^(-1/n));
    coshv = 1/2*(tmp^(1/n)+tmp^(-1/n));

    H_gain = 1;
    H_phi = 0;
    for k=1:n
        sigma_k = -sin((2*k-1)/(2*n)*pi)*sinhv;
        omega_k = cos((2*k-1)/(2*n)*pi)*coshv;
        H_gain = H_gain*1/sqrt(sigma_k^2+(omega_k-w)^2);
        H_phi = H_phi - atan((omega_k-w)/sigma_k);
    end
    nstr = sprintf('xi=%.3gdB',xi);
    
    omega = 0:0.05:2;
    figure(2);
    plot(omega, rad2deg(subs(H_phi, w, omega)),'DisplayName',nstr);
    hold on;
end

nstr = sprintf('位相特性 (n=%.3g)',n);
figure(2); legend;title(nstr);


n = 10;
for idx=1:length(xis)
    xi = xis(idx);
    e = sqrt(10.^(xi/10)-1);
    tmp = sqrt(1+e^(-2))+e^(-1);
    sinhv = 1/2*(tmp^(1/n)-tmp^(-1/n));
    coshv = 1/2*(tmp^(1/n)+tmp^(-1/n));

    H_gain = 1;
    H_phi = 0;
    for k=1:n
        sigma_k = -sin((2*k-1)/(2*n)*pi)*sinhv;
        omega_k = cos((2*k-1)/(2*n)*pi)*coshv;
        H_gain = H_gain*1/sqrt(sigma_k^2+(omega_k-w)^2);
        H_phi = H_phi - atan((omega_k-w)/sigma_k);
    end
    nstr = sprintf('xi=%.3gdB',xi);
    
    omega = 0:0.05:2;
    figure(3);
    plot(omega, rad2deg(subs(H_phi, w, omega)),'DisplayName',nstr);
    hold on;
end

nstr = sprintf('位相特性 (n=%.3g)',n);
figure(3); legend;title(nstr);

big;