% 第八章: 例 8.20
%            雙線性變換
%                  橢圓低通濾波器設計
%
% 數字濾波器指標:
wp = 0.2*pi;                         % 數字通帶頻率(弧度)
ws = 0.3*pi;                         % 數字阻帶頻率(弧度)
Rp = 1;                              % 通帶波動(dB)
As = 15;                             % 阻帶衰減(dB) 

% 模擬原型指標的頻率逆映射
T = 1; Fs = 1/T;                     % 置 T=1
OmegaP = (2/T)*tan(wp/2);            % 預修正原型通帶頻率
OmegaS = (2/T)*tan(ws/2);            % 預修正原型阻帶頻率
ep = sqrt(10^(Rp/10)-1);             % 通帶波動參數
Ripple = sqrt(1/(1+ep*ep));          % 通帶波動
Attn = 1/(10^(As/20));               % 阻帶衰減

% 模擬橢圓原型濾波器計算:
[cs,ds] = afd_elip(OmegaP,OmegaS,Rp,As);
%%*** 橢圓濾波器階次=  3

% 雙線性變換:
 [b,a] = bilinear(cs,ds,T);
[C,B,A] = dir2cas(b,a)
%%C = 0.1214
%%B = 1.0000   -1.4211    1.0000
%%    1.0000    1.0000         0
%%A = 1.0000   -1.4928    0.8612
%%    1.0000   -0.6183         0

% 繪圖
figure(1); subplot(1,1,1)
[db,mag,pha,grd,w] = freqz_m(b,a);
subplot(2,2,1); plot(w/pi,mag); title('幅度響應')
xlabel(''); ylabel('|H|'); axis([0,1,0,1.1])
set(gca,'XTickMode','manual','XTick',[0,0.2,0.3,1]);
set(gca,'YTickmode','manual','YTick',[0,Attn,Ripple,1]); grid
subplot(2,2,3); plot(w/pi,db); title('幅度(dB)');
xlabel('頻率（單位:pi）'); ylabel('decibels'); axis([0,1,-40,5]);
set(gca,'XTickMode','manual','XTick',[0,0.2,0.3,1]);
set(gca,'YTickmode','manual','YTick',[-50,-15,-1,0]); grid
set(gca,'YTickLabelMode','manual','YTickLabels',['50';'15';' 1';' 0'])
subplot(2,2,2); plot(w/pi,pha/pi); title('相位響應')
xlabel(''); ylabel('單位:pi'); axis([0,1,-1,1]);
set(gca,'XTickMode','manual','XTick',[0,0.2,0.3,1]);
set(gca,'YTickmode','manual','YTick',[-1,0,1]); grid
subplot(2,2,4); plot(w/pi,grd); title('群延遲')
xlabel('頻率（單位:pi）'); ylabel('樣本'); axis([0,1,0,15])
set(gca,'XTickMode','manual','XTick',[0,0.2,0.3,1]);
set(gca,'YTickmode','manual','YTick',[0:5:15]); grid