バタワースフィルタの振幅特性・位相特性

clear; close all;

s =  tf('s');
Ns = [5,8,10];
syms w;

for idx = 1:length(Ns)
    n = Ns(idx);
    D = 1;
    H_gain = 1;
    H_phi = 0;  
    for jdx=1:n
        j = jdx;
        sigma_j = -sin((2*j-1)/(2*n)*pi);
        omega_j = cos((2*j-1)/(2*n)*pi);
        sj = sigma_j + i*omega_j;
        D = D*(s-sj);
        H_gain = H_gain*1/sqrt(sigma_j^2+(omega_j-w)^2);
        H_phi = H_phi - atan((omega_j-w)/sigma_j);
    end
    H = 1/D;
    nstr = sprintf('n=%d',n);
    figure(1);
    hold on;
    omega = 0:0.05:2;
    plot(omega, subs(H_gain, w, omega),'DisplayName',nstr);
    figure(2);
    semilogy(omega, subs(H_gain, w, omega),'DisplayName',nstr);
    hold on;
    figure(3);
    hold on;
    plot(omega, subs(H_phi, w, omega),'DisplayName',nstr);
    coef(D)
end
figure(1); legend;title('振幅特性');
figure(2); legend;title('振幅特性(dB)');
figure(3); legend;title('位相特性');

Ns = [15,20,40];
syms w;

for idx = 1:length(Ns)
    n = Ns(idx);
    D = 1;
    H_gain = 1;
    H_phi = 0;  
    for jdx=1:n
        j = jdx;
        sigma_j = -sin((2*j-1)/(2*n)*pi);
        omega_j = cos((2*j-1)/(2*n)*pi);
        sj = sigma_j + i*omega_j;
        D = D*(s-sj);
        H_gain = H_gain*1/sqrt(sigma_j^2+(omega_j-w)^2);
        H_phi = H_phi - atan((omega_j-w)/sigma_j);
    end
    H = 1/D;
    nstr = sprintf('n=%d',n);
    figure(4);
    hold on;
    omega = 0:0.05:2;
    plot(omega, subs(H_gain, w, omega),'DisplayName',nstr);
    figure(5);
    omega = 0:0.05:2;
    semilogy(omega, subs(H_gain, w, omega),'DisplayName',nstr);
    hold on;
    figure(6);
    hold on;
    plot(omega, subs(H_phi, w, omega),'DisplayName',nstr);
    coef(D)
end
figure(4); legend;title('振幅特性');
figure(5); legend;title('振幅特性(dB)');
figure(6); legend;title('位相特性');

big;

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