MATLAB——IIR数字滤波器的设计

1.基础知识
1.1、数字滤波器设计的基本步骤
我们知道模拟滤波器的设计是数字滤波器的设计的基础。在学习数字信号处理的过程中，IIR数字滤波器的设计的步骤是
（1）确定采样间隔Ts或者采样频率fs。
（2）根据模拟频率和数字频率之间的关系，将所给出的数字滤波器的指标转化为模拟滤波器的指标。（Ω=ω/Ts）
（3）根据模拟滤波器的指标设计模拟滤波器。
（4）根据冲激响应不变法和双线性变换法，将H（s）转化为H（z）。

1.2、冲激响应不变法
根据$z=e^{sT}$,将S平面映射到Z平面，但不是一对一的映射。冲激响应不变法只适合用于有限带宽的滤波器设计。

1.3、双线性变换法
通过对S平面进行压缩，消除了频谱混叠的现象。但是在这个过程中产生了非线性畸变，所以说要进行校正。这个方法的S平面与Z平面的映射关系为
$s=\frac{2}{T}\ast \frac{1-z^{-1}}{1+z^{-1}}$
预畸变校正公式为
$\Omega =\frac{2}{T}\ast tan(\frac{\omega }{2})$
2.具体题目
2.1、题目1（脉冲响应不变法）
用脉冲响应不变法设计Butterworth低通数字滤波器，要求通带频率为 ，通带波纹小于1dB，阻带频率为 ，幅度衰减大于15dB，采样频率为Fs=2000Hz。
这里面有几个MATLAB自带的函数，需要解释一下
首先回顾一下模拟滤波器设计里面的函数
[n,wn]=buttord(fp,fs,rp,as,‘s’);
这个函数用于计算巴特沃斯滤波器的阶数和截止频率。它的输入参数包括通带截止频率（fp）、阻带截止频率（fs）、通带最大衰减（rp）和阻带最小衰减（as），以及滤波器的类型（'s’表示模拟滤波器，'z’表示数字滤波器）。这个函数的输出包括阶数（n）和截止频率（wn）。
[z,p,k]=buttap(n);
这个函数用于生成一个巴特沃斯滤波器的极点和零点。它的输入参数是阶数（n），它的输出参数是极点（p）和零点（z），以及一个系数（k），用于归一化滤波器的增益。
[b,a]=zp2tf(z,p,k);
zp2tf：这个函数用于将极点和零点转换为传递函数的分子和分母形式。它的输入参数是极点（p）、零点（z）和系数（k），它的输出是分子系数（b）和分母系数（a）。
变换到数字滤波器的函数
[b,a]=lp2lp(B,A,Wn)
一个低通原型变换为具有不同截止频率的低通滤波器，这里规定了截止频率Wn，输入分别是模拟滤波器分子、分母的系数，以及前面求得的截止频率。
[bz,az]=impinvar(b,a,fs);
impinvar 函数使用脉冲响应不变法将模拟滤波器变换为数字滤波器。具体地，impinvar 函数先对模拟滤波器进行离散化，然后使用离散化后的脉冲响应作为数字滤波器的脉冲响应。
注意
freqz和freqs
freqz 和 freqs 是 MATLAB 中用于绘制滤波器频率响应的函数。(j简而言之，一个离散，一个连续)
freqz 函数可以绘制数字滤波器的离散时间频率响应，其语法为：
freqz(b,a,n,fs)
其中，b 和 a 分别是数字滤波器的分子系数和分母系数，n 是计算离散时间频率响应的点数，fs 是采样频率。

freqs 函数则可以绘制模拟滤波器的频率响应，其语法为：
freqs(b,a,w)
其中，b 和 a 分别是模拟滤波器的分子系数和分母系数，w 是计算频率响应的频率点。需要注意的是，freqs 函数计算的是模拟滤波器的频率响应，因此输入的频率参数 w 是以弧度为单位的。
代码

clc;
%写出已知条件，设置采样点数
wp=0.25*pi;ws=0.4*pi;
rp=1;as=15;
fs=2000;Nn=128;
%由数字指标映射到模拟滤波器指标
WP=wp*fs;WS=ws*fs;
%利用巴特沃斯滤波器的设计方法，设计一个模拟滤波器
[N,Wn]=buttord(WP,WS,rp,as,'s');
[z,p,k]=buttap(N);
[B,A]=zp2tf(z,p,k);
%转换到数字滤波器
[b,a]=lp2lp(B,A,Wn);%从一个低通原型变换为具有不同截止频率的低通滤波器
[bz,az]=impinvar(b,a,fs);%根据冲激响应不变法得到数字滤波器
freqz(bz,az,Nn,fs);

运行结果

2.2题目2（双线性变换法）
用双线性变换法设计Butterworth低通数字滤波器，要求通带频率为 ，通带波纹小于1dB，阻带频率为 ，幅度衰减大于15dB，采样频率为Fs=100Hz。
有了上面的基础，我们来进行类比学习。

WS=2*fs*tan(ws/2);
WP=2*fs*tan(wp/2);

增加校正畸变的过程。

[bz,az]=bilinear(b,a,fs);

换成双线性不变法

freqz(bz,az,Nn,fs);

[H,w]=freqz(bz,az);
figure;
subplot(211);
plot(w/pi,abs(H));
subplot(212);
plot(w/pi,angle(H));

最后绘图部分，我是做了个对比。
第一个是直接freqz绘制的幅频和相频图像。
freqz函数计算数字滤波器的频率响应，返回的是频率响应的幅度和相位信息。其中w是频率，H是频率响应的复数值。因为通常我们更关心的是频率在$0$到$\pi$之间的变化情况，所以使用plot(w/pi,abs(H))和plot(w/pi,angle(H))将频率从弧度变为单位圆上的比例。就是我第二个绘图窗口的文件。

完整代码

wp=0.25*pi;ws=0.4*pi;
rp=1;as=15;
fs=100;
Ts=1/fs;
Nn=128;

WS=2*fs*tan(ws/2);
WP=2*fs*tan(wp/2);

[N,Wn]=buttord(WP,WS,rp,as,'s');
[z,p,k]  =buttap(N);
[B,A]=zp2tf(z,p,k);
[b,a]=lp2lp(B,A,Wn);
[bz,az]=bilinear(b,a,fs);

freqz(bz,az,Nn,fs);

[H,w]=freqz(bz,az);
figure;
subplot(211);
plot(w/pi,abs(H));
subplot(212);
plot(w/pi,angle(H));

运行结果