% ////////////巴特沃茲帶阻濾波器/////////////////////////////
% wp1=0.3*pi;wp2=0.75*pi;ws1=0.45*pi;ws2=0.65*pi;Rp=1;Rs=40
% =============雙線型變換法==============================
wp1=0.3*pi;     wp2=0.75*pi; 
ws1=0.45*pi;    ws2=0.65*pi;
Rp=1;           Rs=40;
Wp1=tan(wp1/2);   Wp2=tan(wp2/2);
Ws1=tan(ws1/2);   Ws2=tan(ws2/2);
BW=Ws2-Ws1;       W0=Ws1*Ws2;W00=sqrt(W0);
WS=1;     WP=WS*(Wp1*BW)/(W0-Wp1^2);

[N,Wn]=buttord(WP,WS,Rp,Rs,'s');
[B,A]=butter(N,Wn,'s');
[BT,AT]=lp2bs(B,A,W00,BW);
[num,den]=bilinear(BT,AT,0.5);
[h,omega]=freqz(num,den,64);

subplot(2,1,1);stem(omega/pi,abs(h));grid;axis([0 1 0 1]);
xlabel('雙線性法\omega/\pi');ylabel('|H(z)|');

subplot(2,1,2);stem(omega/pi,20*log10(abs(h)));grid;axis([0 1 -200 0]);
xlabel('雙線性法\omega/\pi');ylabel('增益.dB');

% ===============直接型==============================
figure
wp1=0.3*pi;     wp2=0.75*pi; 
ws1=0.45*pi;    ws2=0.65*pi;
Rp=1;          Rs=40;
Wp=[wp1/pi,wp2/pi];  Ws=[ws1/pi,ws2/pi];
[N,Wn]=buttord(Wp,Ws,Rp,Rs);
[B,A]=butter(N,Wn,'stop');
[h,omega]=freqz(B,A,64);

subplot(2,1,1);stem(omega/pi,abs(h));grid;
xlabel('直接法\omega/\pi');ylabel('|H(z)|');

subplot(2,1,2);stem(omega/pi,20*log10(abs(h)));grid;
xlabel('直接法\omega/\pi');ylabel('增益.dB');