%impinvar_bilinear.m
%chebyshev1模擬原型濾波器 
Omegap=0.2*pi;Rp=1;     %設計指標 
Omegas=0.3*pi;As=16; 
[n1,Wn1]=cheb1ord(Omegap,Omegas,Rp,As,'s');  %獲得chebyshev1的階次和截止頻率 
[B1,A1]=cheby1(n1,0.5,Wn1,'s');       %獲得chebyshev1模擬原型濾波器 
w1=[0:500]*0.5*pi/500;              %頻率取樣點從0到0.5pi 
h1=freqs(B1,A1,w1);                %計算在各取樣點的傅立葉變換值 
subplot(3,1,1);plot(w1,abs(h1));        %繪制幅頻特性 
axis([0,0.5*pi,0,1]); 
grid; 
ylabel('模擬原型濾波器'); 

%利用脈沖響應不變法設計數(shù)字濾波器
T=0.2; %取樣周期 
[Bz,Az]=impinvar(B1,A1,1/T); 
w2=[0:500]*0.5*pi*T/500;           %頻率取樣點從0到0.5pi 
[H,W]=freqz(Bz,Az,w2);            %計算在各取樣點的傅立葉變換值 
w=W/T;                          %將數(shù)字頻率轉(zhuǎn)換為模擬頻率 
subplot(3,1,2); 
plot(w,abs(H));                    %繪制幅頻特性 
axis([0,0.5*pi,0,1]); 
grid; 
ylabel('脈沖響應不變法'); 

%利用雙線性變換法設計數(shù)字濾波器
wp=Omegap*T;                    %將模擬頻率轉(zhuǎn)換為數(shù)字頻率 
ws=Omegas*T; 
Omegap1=(2/T)*tan(wp/2);           %對模擬頻率進行預畸變 
Omegas1=(2/T)*tan(ws/2); 
[n,Wn]=cheb1ord(Omegap1,Omegas1,Rp,As,'s'); %獲得chebyshev1的階次和截止頻率 
[B,A]=cheby1(n,0.5,Wn,'s');          %獲得chebyshev1模擬原型濾波器 
[Bz,Az]=bilinear(B,A,1/T);           %利用雙線性變換法設計數(shù)字濾波器 
w2=[0:500]*0.5*pi*T/500;           %頻率取樣點從0到0.5pi 
[H,W]=freqz(Bz,Az,w2);             %計算在各取樣點的傅立葉變換值 
w=W/T;                           %將數(shù)字頻率轉(zhuǎn)換為模擬頻率 
subplot(3,1,3); 
plot(w,abs(H));                     %繪制幅頻特性 
axis([0,0.5*pi,0,1]); 
grid; 
ylabel('雙線性變換法');