数字滤波器设计指南：脉冲响应不变法和双线性变换法

%---------digital filter's guide line：数字滤波器的设计指南
clear all;  %清除工作空间中的所有变量
clc;        %清除命令窗口中的所有内容
wp=0.2*pi;              %通带截止频率
ws=0.3*pi;              %阻带截止频率
rp=1;                   %通带最大衰减
rs=15;                  %阻带最小衰减
fs=1000;                %采样频率
Ts=1/fs;                %采样周期
%------------analog filter's guide line：模拟滤波器的设计指南----------------
omegap=wp*fs;                   %脉冲响应不变法中的通带截止频率
omegas=ws*fs;                   %脉冲响应不变法中的阻带截止频率
omegap2=(2/Ts)*tan(wp/2);       %双线性变换法中的通带截止频率
omegas2=(2/Ts)*tan(ws/2);       %双线性变换法中的阻带截止频率
%-----AF Design(impluse response invariant)：脉冲响应不变法设计------------------------
[N,omegac]=buttord(omegap,omegas,rp,rs,'s');     %计算Butterworth滤波器的阶数和截止频率
[ba,aa]=butter(N,omegac,'s');                   %计算Butterworth滤波器的系数
[bd,ad]=impinvar(ba,aa,fs);                     %脉冲响应不变法
[H,W]=freqz(bd,ad,1000);                        %计算滤波器的频率响应
h=freqz(bd,ad,[wp,ws]);                         %计算滤波器在通带和阻带的增益
%------------direct DF IIR design：直接IIR设计---------  
[N1,omegac1]=buttord(wp/pi,ws/pi,rp,rs);        %计算Butterworth滤波器的阶数和截止频率
[bd1,ad1]=butter(N1,omegac1);                   %计算Butterworth滤波器的系数
[H1,W1]=freqz(bd1,ad1,1000);                    %计算滤波器的频率响应
h1=freqz(bd1,ad1,[wp,ws]);                      %计算滤波器在通带和阻带的增益

%-------------------AF Design(bilinear transform)：双线性变换法设计-----------------
[N2,omegac2]=buttord(omegap2,omegas2,rp,rs,'s'); %计算Butterworth滤波器的阶数和截止频率
[ba2,aa2]=butter(N2,omegac2,'s');                %计算Butterworth滤波器的系数
[bd2,ad2]=bilinear(ba2,aa2,fs);                  %双线性变换法
[H2,W2]=freqz(bd2,ad2,1000);                     %计算滤波器的频率响应
h2=freqz(bd2,ad2,[wp,ws]);                       %计算滤波器在通带和阻带的增益
%-------------------------draw figure：绘制图形-------------  
figure(1)   %创建第1个图形窗口
subplot(1,3,1)  %将第1个子图放在1行3列的第1个位置上
plot(W/pi,20*log(abs(H))); %绘制频率响应曲线
hold on;    %保持原有图像并添加新图像
stem([wp,ws]/pi,20*log(abs(h)),'o:'); %绘制通带和阻带的增益
grid on;    %添加网格线
axis([0,1,-100,0]); %设置坐标轴范围
title('脉冲响应不变法'); %设置子图标题
subplot(1,3,2)  %将第2个子图放在1行3列的第2个位置上
plot(W1/pi,20*log(abs(H1)));   %绘制频率响应曲线
hold on;    %保持原有图像并添加新图像
stem([wp,ws]/pi,20*log(abs(h1)),'o:'); %绘制通带和阻带的增益
grid on;    %添加网格线
axis([0,1,-100,0]); %设置坐标轴范围
title('直接设计IIR'); %设置子图标题
subplot(1,3,3)  %将第3个子图放在1行3列的第3个位置上
plot(W2/pi,20*log(abs(H2)));   %绘制频率响应曲线
hold on;    %保持原有图像并添加新图像
stem([wp,ws]/pi,20*log(abs(h2)),'o:'); %绘制通带和阻带的增益
grid on;    %添加网格线
axis([0,1,-100,0]); %设置坐标轴范围
title('双线性变换法'); %设置子图标题