利用Matlab设计IIR数字滤波final

利用Matlab 设计IIR 数字滤波 在数字信号处理当中，数字滤波器又占有相当重要的位置。

因为任何需要我们处理的信号都包含有噪声，而滤波是去除噪声的基本手段。

数字滤波器简单的说就是一种用来过滤时间离散信号的线性时不变系统，通过对采样数据进行数学处理来达到频域滤波的目的。

根据其单位冲激响应函数的时域特性可以分为两类：无限冲激响应(IIR)滤波器和有限冲激响应（FIR ）滤波器。

数字滤波本质上是一种运算过程，输入数字信号（数字序列）通过特定的运算转变为输出的数字序列，因此，数字滤波器的本质上是一个完成特定运算的数字计算过程。

描述离散系统输入和输出关系的卷积和差分方程只是给数字信号滤波器提供运算规则，使其按照这个规则完成对输入数据的处理。

一、数字滤波器设计原理
在数字滤波过程中，我们主要讨论离散时间序列。

设输入序列为x(n)，滤波器对单位采样序列δ(n)的响应为h(n)。

因δ(n)在时域离散信号和系统中所起的作用相当于单位冲激函数在时域连续信号和系统中所起的作用。

所以数字滤波器的序列y(n)就是以上两个序列的卷积，即：
()()()k y n h k x n k +∞
=-∞=-∑ (1) 同样，两个序列卷积的z 变换等于各自z 变换的乘积，即：
()()()Y z H z X z = (2)
用j z e ω=代入上式，则可得到时域离散系统的频谱特性：
()()()j j j Y e H e X e ωωω= (3)
式中()j Y e ω和()j X e ω分别是数字滤波器的输出序列和输入序列的频谱特性，而()j H e ω则是单位采样序列h(n)的频特性谱，又称为数字滤波器的频域响应。

输入序列()j X e ω经过滤波后变为()j X e ω()j H e ω。

因此，只需要按照输入信号频谱()j X e ω的特点和我们处理信号的目的，选择适当的频域响应()j H e ω，使得滤波后的输出信号频谱()j Y e ω满足我们设计的要求。

二、具体设计步骤
1. 推导滤波器最小阶数与设计指标的关系
滤波器阶数就是其系统函数的极点个数。

为了避免滤波器的复杂程度与我们的要求不匹配，造成不必要的成本浪费，我们在滤波器设计前先需要确定其合适的阶数，即满足设计要求的最小阶数。

下面就分别以Butterworth 滤波器和Chebyshev1型滤波器为例，推导其阶数的数学模型。

①对于Butterworth 滤波器，若给出滤波器的设计指标为：通带截止频率p ω，阻带截止频率s ω，通带最大纹波()p R dB ，阻带最小纹波()s R dB 。

因为滤波器幅频特性为： 21
(1()N
c H j ωωω=+ （4）
其中c ω为3dB 截止频率，N 为滤波器阶数。

所以当p ωω=以及s ωω=时，可得到： 221
1((1()1()p s N p s N c c H j H j ωωωωωω==++、 (5)
然后由()p R dB 与()s R dB 的关系式可得到： 221022101110lg ()(101110lg ()(10p s p p R p s s R s R H j H j R H j H j ωωωω⎧⎪⎪⎪⎨⎪⎪⎪⎩
=⇒==⇒= (6)
由上面的(5)式和(6)式可以联立求得： 21010102101()10101()1011()10p p s s R p N R p c N R R s N s c ωωωωωω⎧⎪⎪⎨⎪⎪⎩
+=-⇒=-+= (7) 通过上面的结果我们就可以求得滤波器阶数N 为： 10
10101lg 101lg()p s R R p s
N ωω--= (8) 所以滤波器的最小阶数就是大于上式所求得的值的最小整数。

②对于Chebyshev1型滤波器，推导方式和上面类似，同样的给出滤波器的设计指标为：通带截止频率p ω，阻带截止频率s ω，通带最大纹波()p R dB ，阻带最小纹波()s R dB 。

因为滤波器幅频特性为： 221
()1()c N H j V ωωωε=+ (9) 又因为通带内起伏的等纹波参数/1010
1p R ε=-，系数2010p R A =，所以可以得出滤波器阶数为： 22lg(1)
lg(()1)s s p p g g N ωωωω++=+- (10)
式中22(1)/g A ε=-。

滤波器的最小阶数就是取大于上式的值的 最小整数。

2. 双线性变换法设计Butterworth 低通滤波器
双线性变换法是在完成模拟系统设计后，再进行反变换，从模拟系统函数得到数字系统函数。

它没有采样过程，不存在频率混叠问题，适应于各类滤波器的变换，但是变换在高频区域会出现强烈的非线性。

下面就是基于Butterworth 模拟原型滤波器，使用双线性变换法设计数字低通滤波器的过程。

设计指标为：通带截止频率p ω=0.2π,阻带截止频率s ω=0.3π，通带最大纹波
Rp=1dB ，阻带最大纹波Rs=15dB ，采样频率Fs=20000Hz 。

同时，为了检测滤波器的性能，我们假定一个输入信号21)()sin(2)0.5cos(2t x t f t f ππ=+，其中f1=1000Hz ，f2=4000Hz ，将该信号与通过滤波器之后产生的输出信号进行比较以测试滤波器的性能。

得到的Butterworth 低通滤波器幅频特性和相频特性曲线为：
由语句[N,Wn]=buttord(wp,ws,Rp,Rs,'s')我们可以用Matlab得出滤波器的最小阶数N=6。

再利用上面我们推导出的式（8），将滤波器各项指标代入可以算出N=4.7433，取大于它的最小整数得到滤波器的最小阶数为5。

二者并不一致，但差距只有1，说明我们推导出的计算Butterworth滤波器最小阶数的数学模型存在误差，但误差不大。

又因为Matlab计算滤波器的最小阶数也是估算出来的，所以这个这个误差在接受范围之类，我们推导出的数学模型可以认为是正确的。

输入信号和经过滤波后的输出信号图像为：
总结
本文中采用双线性变换法设计了一种滤波器——Butterworth低通滤波器。

通过比较幅频特性和相频特性曲线，我们可以看出用双线性变换法设计出的简单低通滤波器的频谱特性。

由于采样频率有限，所以用双线性变换法时也并没有表现出很强烈的非线性。

再比较输入信号和经过滤波后的输出信号图像，我们可以发现，双线性变换法设计出的IIR数字滤波器的滤
波性能还挺好的，输入信号有很多的小尖峰，但是经过滤波后，这些干扰的小尖峰都被滤除了，获得了效果很好的信号，对信号的处理是一个很大的提高。

另外，我们利用Butterworth滤波器的频域响应函数与通带最大纹波以及阻带最小纹波之间的关系推导出了其滤波器最小阶数N的数学模型，然后通过Matlab仿真对该模型进行了检验，证明了它具有一定的可靠性。

程序代码:
双线性变换法设计Butterworth低通滤波器
Matlab程序如下：
Wp=0.2*pi; %通带截止频率
Ws=0.3*pi; %阻带截止频率
fs=20000; %采样频率
T=1/fs;
Rp=1; %通带最大衰减
Rs=15; %阻带最小衰减
wp=2*tan(Wp/2)/T; %双线性映射
ws=2*tan(Ws/2)/T;
[N,Wn]=buttord(wp,ws,Rp,Rs,'s');
[z,p,k]=buttap(N); %Butterworth滤波器函数[Bap,Aap]=zp2tf(z,p,k); %滤波器类型转换
[b,a]=lp2lp(Bap,Aap,Wn); %原型低通到低通模拟滤波器的变换函数
[bz,az]=bilinear(b,a,fs); %双线性变换法函数
figure(1);
[h,w]=freqz(bz,az,N,fs); %数字滤波器的频率响应
subplot(2,1,1),plot(w,abs(h));
title('Butterworth LPF幅频特性');
xlabel('频率(Hz) ');ylabel('幅值(dB) ');grid on;
subplot(2,1,2),plot(w,angle(h)); %相频特性
title('Butterworth LPF相频特性');
xlabel('频率(Hz) ');ylabel('相位(degree) ');grid on;
figure(2);
f1=1000;f2=4000;
N1=100;
dt=1/fs;n=0:N1-1;t=n*dt;
x=sin(2*pi*f1*t)+0.5*cos(2*pi*f2*t); %含噪声的信号
subplot(2,1,1),plot(t,x);
title('输入信号'),grid on;
y=filter(bz,az,x); %使用设计的滤波器进行滤波
subplot(2,1,2),plot(t,y, 'r-');
title('输出信号') ,grid on;。