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