ディジタルフィルタの振幅特性

リファレンスとなる理想アナログフィルタの伝達関数Gを定義
→理想アナログフィルタを逆フーリエ変換して理想アナログフィルタの時系列信号g(t)が求まる
→欲しいディジタルフィルタHのインパルス応答hkが求まる（hk=T•g(kT)の関係式より）
→hkとZ変換の関係式より、ディジタルフィルタHの伝達関数H(z)が求まる
→伝達関数H(z)にz=exp(jωT)=exp(jω)(∵T=1)を代入し、周波数応答関数H(e^(jω))を求める
→周波数応答関数H(e^(jω))の絶対値を取ることで（Matlabのabs()関数）、振幅特性が得られる。

clear; close all;

wc = pi/2;
T = 1;
N=61;
l = (N-1)/2;
h1 = [];
h2 = [];
h3 = [];
h4 = [];
w0 = pi/8;
ks = 1:N;
for kprime = ks
    k = kprime-1;
    hk1 = wc*T/pi*sin((k-l)*wc*T)/(k-l)/wc/T;
    h1(end+1) = hk1;
    hk2 = -wc*T/pi*sin((k-l)*wc*T)/(k-l)/wc/T;
    h2(end+1) = hk2;
    hk3 = 2*wc*T/pi*sin((k-l)*wc*T)/(k-l)/wc/T*cos(k-l)*w0*T;
    h3(end+1) = hk3;
    hk4 = -2*wc*T/pi*sin((k-l)*wc*T)/(k-l)/wc/T*cos(k-l)*w0*T;
    h4(end+1) = hk4;
end
h1(l+1) = wc*T/pi;
h2(l+1) = 1-wc*T/pi;
h4(l+1) = 1-2*wc*T/pi;

% subplot(231);
% stem(ks, h1 ,'o-','DisplayName','Low pass filter');
% legend();
% subplot(232);
% stem(ks, h2 ,'o-','DisplayName','High pass filter');
% legend();
% subplot(233);
% stem(ks, h1+h2 ,'o-','DisplayName','LPF+HPF');
% legend();
% 
% subplot(234);
% stem(ks, h3 ,'o-','DisplayName','Band pass filter');
% legend();
% subplot(235);
% stem(ks, h4 ,'o-', 'DisplayName','Band eliminate filter');
% legend();
% subplot(236);
% stem(ks, h3+h4 ,'o-', 'DisplayName','BPF+BEF');
% legend();
% big;


Z61 = 0;
syms z;
for n = 0:(length(h1) - 1)
    Z61 = Z61 + h1(n+1) * z^(-n);    
end
Z61 = vpa(simplify(Z61))
syms w;
H61 = subs(Z61, z, exp(i*w*T));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N=121;
ks121 = 1:N;
l = (N-1)/2;
h121 = [];
for kprime = ks121
    k = kprime-1;
    hk1 = wc*T/pi*sin((k-l)*wc*T)/(k-l)/wc/T;
    h121(end+1) = hk1;
end
h121(l+1) = wc*T/pi;

Z121 = 0;
syms z;
for n = 0:(length(h121) - 1)
    Z121 = Z121 + h121(n+1) * z^(-n);    
end
Z121 = vpa(simplify(Z121))
syms w;
H121 = subs(Z121, z, exp(i*w*T));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N=241;
ks241 = 1:N;
h241 = [];
l = (N-1)/2;
for kprime = ks241
    k = kprime-1;
    hk1 = wc*T/pi*sin((k-l)*wc*T)/(k-l)/wc/T;
    h241(end+1) = hk1;
end
h241(l+1) = wc*T/pi;

Z241 = 0;
syms z;
for n = 0:(length(h241) - 1)
    Z241 = Z241 + h241(n+1) * z^(-n);   
end
Z241 = vpa(simplify(Z241))
syms w;
H241 = subs(Z241, z, exp(i*w*T));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%


figure();
omegas = linspace(0,pi,100);

subplot(321);
plot(omegas, abs(subs(H61, w, omegas)),'DisplayName','N=61');
legend();
title('振幅特性');

subplot(322);
stem(ks, h1 ,'o-','DisplayName','N=61');
legend();
title('インパルス応答');

subplot(323);
plot(omegas, abs(subs(H121, w, omegas)),'DisplayName','N=121');
legend();
title('振幅特性');

subplot(324);
stem(ks121, h121 ,'o-','DisplayName','N=121');
legend();
title('インパルス応答');

subplot(325);
plot(omegas, abs(subs(H241, w, omegas)),'DisplayName','N=241');
legend();
title('振幅特性');

subplot(326);
stem(ks241, h241 ,'o-','DisplayName','N=241');
legend();
title('インパルス応答');

big;

・フィルタの次数Nを高めるとリプルは抑えられるが、遅延時間L=(N-1)/2*Tが増大する
・ギブス現象による遮断周波数ω=ωcでの不連続性は解消できない