%《數字信號處理》第六章 例6.4.2 設計程序：ep642.m
% 電子工業出版社出版  高西全  丁玉美  闊永紅合著 2006年8月
% 用雙線性變換法設計DF

T=1;Fs=1/T;
wpz=0.2;wsz=0.3;
wp=2*tan(wpz*pi/2);ws=2*tan(wsz*pi/2);rp=1;rs=1 %預畸變校正轉換指標
[N,wc]=buttord(wp,ws,rp,rs,'s');    %設計過渡模擬濾波器
[B,A]=butter(N,wc,'s');
%繪制模擬濾波器的損耗函數曲線
fk=0:1/512:1;wk=2*pi*fk;
Hk=freqs(B,A,wk);
subplot(2,2,1);
plot(fk,20*log10(abs(Hk)));grid on
xlabel('f(Hz)');ylabel('幅度(dB)')
axis([0,0.9,-100,5]);title('(a)         ');
[Bz,Az]=bilinear(B,A,Fs);      %用雙線性變換法轉換成數字濾波器
wk=0:pi/512:pi;
Hz=freqz(Bz,Az,wk);
subplot(2,2,2);
plot(wk/pi,20*log10(abs(Hz)));grid on;
xlabel('\omega/\pi');ylabel('幅度(dB)')
axis([0,1,-100,5]);title('(b)           ');
[N,wc]=buttord(wpz,wsz,rp,rs);      %調用buttord和butter直接設計數字濾波器   
[Bz,Az]=butter(N,wc);
wk=0:pi/512:pi;
Hz=freqz(Bz,Az,wk);
subplot(2,2,4);
plot(wk/pi,20*log10(abs(Hz)));grid on;
xlabel('\omega/\pi');ylabel('幅度(dB)')
axis([0,1,-100,5]);title('(b)           ');