%《MATLAB在電子信息課程中的應用》(第三版)第七章例7.33程序q733
% 用雙線性變換法設計高通和帶通巴特沃斯數字濾波器
% 電子工業出版社出版  陳懷琛 吳大正 高西全合著 2001年10月

%用雙線性變換法設計數字高通和帶通濾波器
clear;close all
T=1;wch=pi/2;			% T為采樣間隔，wch為數字高通3dB截止頻率
wlc=0.35*pi;wuc=0.65*pi;		% wlc,wuc:數字帶通3dB截止頻率
B=1;A=[1,2.6131,3.4142,2.6131,1];
[h,w]=freqs(B,A,512);			% 求原歸一化模擬濾波器的頻率響應
subplot(3,2,1);plot(w,20*log10(abs(h)));	% 畫模擬濾波器幅頻特性
grid;axis([0,10,-90,0])
xlabel('ω/π');ylabel('模擬低通幅度(dB)')
%(1)設計高通
omegach=2*tan(wch/2)/T; 		% 預畸變求模擬高通3dB截止頻率
[Bhs,Ahs]=lp2hp(B,A,omegach);		% 模擬域低通轉換為高通系數
[Bhz,Ahz]=bilinear(Bhs,Ahs,1/T);	% 模擬轉換為數字高通系數向量
[h,w]=freqz(Bhz,Ahz,512);			% 求并畫數字濾波器幅頻特性
subplot(3,2,3);plot(w/pi,20*log10(abs(h)));
grid;axis([0,1,-150,0])
xlabel('ω/π');ylabel('數字高通幅度(dB)')
%(2)設計帶通
omegalc=2*tan(wlc/2)/T; 	% 預畸變求濾波器通帶低端截止頻率
omegauc=2*tan(wuc/2)/T;		% 預畸變求濾波器通帶高端截止頻率
wo=sqrt(omegalc*omegauc);Bw=omegauc-omegalc;
[Bbs,Abs]=lp2bp(B,A,wo,Bw);	%		% 模擬域低通轉換為帶通系數
[Bbz,Abz]=bilinear(Bbs,Abs,1/T);	% 模擬轉換為數字帶通系數向量
[h,w]=freqz(Bbz,Abz,512);			% 求并畫數字濾波器幅頻特性
subplot(3,2,4);plot(w/pi,20*log10(abs(h)));
grid;axis([0,1,-150,0])
xlabel('ω/π');ylabel('數字帶通幅度(dB)')