数字信号处理(DSP)课程设计—利用Matlab设计IIR数字滤波器

数字信号处理

课

程

设

计

报

告

姓名：一.蒲钇霖 二.舒琳

学号：一.************ 二.************

学院：微固

利用Matlab设计IIR数字滤波器

一、前言

数字化、智能化和网络化是当代信息技术发展的大趋势，而数字化是智能化和网络化的基础，这也就促进了数字信号处理这门技术的飞速发展，使其不仅自成一门学科，更是以不同形式影响和渗透到其它学科。如今，数字信号处理技术与国民经济息息相关，与国家建设紧密相连，影响着并改变着我们的生产、生活方式，因此它受到了社会的普遍关注。

而在数字信号处理当中，数字滤波器又占有相当重要的位置。因为任何需要我们处理的信号都包含有噪声，而滤波是去除噪声的基本手段。数字滤波器简单的说就是一种用来过滤时间离散信号的线性时不变系统，通过对采样数据进行数学处理来达到频域滤波的目的。根据其单位冲激响应函数的时域特性可以分为两类：无限冲激响应(IIR)滤波器和有限冲激响应（FIR）滤波器。数字滤波本质上是一种运算过程，输入数字信号（数字序列）通过特定的运算转变为输出的数字序列，因此，数字滤波器的本质上是一个完成特定运算的数字计算过程。描述离散系统输入和输出关系的卷积和差分方程只是给数字信号滤波器提供运算规则，使其按照这个规则完成对输入数据的处理。

二、数字滤波器设计原理

在数字滤波过程中，我们主要讨论离散时间序列。设输入序列为x(n)，滤波器对单位采样序列δ(n)的响应为h(n)。因δ(n)在时域离散信号和系统中所起的作用相当于单位冲激函数在时域连续信号和系统中所起的作用。所以数字滤波器的序列y(n)就是以上两个序列的卷积，即：

                  (1)

同样，两个序列卷积的z变换等于各自z变换的乘积，即：

                     (2)

用代入上式，则可得到时域离散系统的频谱特性：

                   (3)

式中和分别是数字滤波器的输出序列和输入序列的频谱特性，而则是单位采样序列h(n)的频特性谱，又称为数字滤波器的频域响应。输入序列经过滤波后变为。因此，只需要按照输入信号频谱的特点和我们处理信号的目的，选择适当的频域响应，使得滤波后的输出信号频谱满足我们设计的要求。

三、课程设计目的与内容

通过完成本次课程设计，加深我们对于在数字信号处理这门学科的印象，巩固我们在课堂上学习到的理论知识和实验方法，增强我们将理论知识化为实践技巧的能力，并且让我们进一步熟悉数字信号处理中常用软件Matlab的使用方法。在这个过程中同时培养我们发现问题、分析问题以及解决问题的能力。

在这次课程设计当中，我们需要分析和研究采用不同设计方法时，滤波器最低阶数与设计指标之间的关系，通过理论分析方式推导建立起相应的简单数学模型。在此基础上，自选设计指标，分别采用双线性变换法和冲激响应不变法进行设计Butterworth低通滤波器，然后通过仿真检验所建模型的正确性和局限性。

四、具体设计步骤

1.推导滤波器最小阶数与设计指标的关系

滤波器阶数就是其系统函数的极点个数。为了避免滤波器的复杂程度与我们的要求不匹配，造成不必要的成本浪费，我们在滤波器设计前先需要确定其合适的阶数，即满足设计要求的最小阶数。下面就分别以Butterworth滤波器和Chebyshev1型滤波器为例，推导其阶数的数学模型。

①对于Butterworth滤波器，若给出滤波器的设计指标为：通带截止频率，阻带截止频率，通带最大纹波，阻带最小纹波。因为滤波器幅频特性为：

                           （4）

其中为3dB截止频率，N为滤波器阶数。所以当以及时，可得到：

                (5)

然后由与的关系式可得到：

                 (6)

由上面的(5)式和(6)式可以联立求得：

             (7)

通过上面的结果我们就可以求得滤波器阶数N为：

                        (8)

所以滤波器的最小阶数就是大于上式所求得的值的最小整数。

②对于Chebyshev1型滤波器，推导方式和上面类似，同样的给出滤波器的设计指标为：通带截止频率，阻带截止频率，通带最大纹波，阻带最小纹波。因为滤波器幅频特性为：

                  (9)

又因为通带内起伏的等纹波参数，系数，所以可以得出滤波器阶数为：

                  (10)

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

2.双线性变换法设计Butterworth低通滤波器

双线性变换法是在完成模拟系统设计后，再进行反变换，从模拟系统函数得到数字系统函数。它没有采样过程，不存在频率混叠问题，适应于各类滤波器的变换，但是变换在高频区域会出现强烈的非线性。下面就是基于Butterworth模拟原型滤波器，使用双线性变换法设计数字低通滤波器的过程。设计指标为：通带截止频率=0.2π,阻带截止频率=0.3π，通带最大纹波Rp=1dB，阻带最大纹波Rs=15dB，采样频率Fs=20000Hz。同时，为了检测滤波器的性能，我们假定一个输入信号，其中f1=1000Hz，f2=4000Hz，将该信号与通过滤波器之后产生的输出信号进行比较以测试滤波器的性能。

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);

[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;

得到的Butterworth低通滤波器幅频特性和相频特性曲线为：

由语句[N,Wn]=buttord(wp,ws,Rp,Rs,'s')我们可以用Matlab得出滤波器的最小阶数N=6。再利用上面我们推导出的式（8），将滤波器各项指标代入可以算出N=4.7433，取大于它的最小整数得到滤波器的最小阶数为5。二者并不一致，但差距只有1，说明我们推导出的计算Butterworth滤波器最小阶数的数学模型存在误差，但误差不大。又因为Matlab计算滤波器的最小阶数也是估算出来的，所以这个这个误差在接受范围之类，我们推导出的数学模型可以认为是正确的。

输入信号和经过滤波后的输出信号图像为：

3.冲激响应不变法设计Butterworth低通滤波器

冲激响应不变法是对模拟滤波器的单位冲激响应h(t)等间距采样获得数字滤波器的单位冲激响应，由此得到数字滤波器的系统函数。这次我们设定与上面不同的设计指标：通带截止频率=2000π，阻带截止频率=3000π，通带最大纹波Rp=3dB，阻带最小纹波Rs=18dB,采样频率Fs=10000Hz。为了方便比较，我们设定一样的输入信号。 

Matlab程序如下：

wp=2000*pi;

ws=3000*pi;

fs=10000;

Rp=3;

Rs=18;

[N,wn]=buttord(wp,ws,Rp,Rs,'s');

[z,p,k]=buttap(N);

[Bap,Aap]=zp2tf(z,p,k);

[b,a]=lp2lp(Bap,Aap,wn);

[bz,az]=impinvar(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;

得到的Butterworth低通滤波器幅频特性和相频特性曲线为：

通过比较两组幅频特性和相频特性曲线我们可以发现，其实双线性变换法和冲激响应不变法在设计较为简单的Butterworth低通滤波器时差别并不大，并且都拥有不错的设计效果。同样的，我们通过语句[N,wn]=buttord(wp,ws,Rp,Rs,'s')可以得到Matlab估算出的滤波器最小阶数为N=6，和上一个滤波器一样。用式(8)出来的N=5.1925，取大于它的最小整数得到滤波器的最小阶数为6，两者一致，再次验证了该数学模型的可靠性。

输入信号和经过滤波后的输出信号图像为：

4.双线性变换法设计Chebyshev1型低通滤波器

为了和Butterworth滤波器进行对比，我们采用和双线性变换法设计Butterworth低通滤波器时一样的设计指标：通带截止频率=0.2π,阻带截止频率=0.3π，通带最大纹波Rp=1dB，阻带最大纹波Rs=15dB，采样频率Fs=20000Hz。输入信号同样为：。

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]=cheb1ord(wp,ws,Rp,Rs,'s');

[z,p,k]=cheb1ap(N,Rp);

[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('Chebyshev1 LPF幅频特性');

xlabel('频率(Hz) ');ylabel('幅值(dB) ');grid on;

subplot(2,1,2),plot(w,angle(h));

title('Chebyshev1 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;

得到的Chebyshev1型低通滤波器幅频特性和相频特性曲线为：

输入信号和经过滤波后的输出信号图像为：

5.冲激响应不变法设计Chebyshev1型低通滤波器

同样的为了和冲激响应不变法设计的Butterworth低通滤波器进行对比，我们还是采用和它一样的设计指标：通带截止频率=2000π，阻带截止频率=3000π，通带最大纹波Rp=3dB，阻带最小纹波Rs=18dB,采样频率Fs=10000Hz。依然使用一样的输入信号：。

Matlab程序如下：

wp=2000*pi;

ws=3000*pi;

fs=10000;

Rp=3;

Rs=18;

[N,wn]=cheb1ord(wp,ws,Rp,Rs,'s');

[z,p,k]=cheb1ap(N,Rp);

[Bap,Aap]=zp2tf(z,p,k);

[b,a]=lp2lp(Bap,Aap,wn);

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

figure(1);

[h,w]=freqz(bz,az,N,fs);

subplot(2,1,1),plot(w,abs(h));

title('Chebyshev1 LPF幅频特性');

xlabel('频率(Hz) ');ylabel('幅值(dB) ');grid on;

subplot(2,1,2),plot(w,angle(h));

title('Chebyshev1 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;

得到的Chebyshev1型低通滤波器幅频特性和相频特性曲线为:

输入信号和经过滤波后的输出信号图像为：

五、总结

这次课程设计同时采用双线性变换法和冲激响应不变法两种方法分别设计了两种滤波器——Butterworth低通滤波器和Chebyshev1型低通滤波器。通过比较四组幅频特性和相频特性曲线，我们可以看出用双线性变换法和冲激响应不变法设计出的简单低通滤波器的频谱特性差别并不明显。由于阻带纹波较小，所以用冲激响应不变法时没有发生明显的混叠现象；而且由于采样频率有限，所以用双线性变换法时也并没有表现出很强烈的非线性。再比较四组输入信号和经过滤波后的输出信号图像我们可以发现，同一种方法设计出的IIR数字滤波器的滤波性能比较相似。另外，我们利用Butterworth滤波器和Chebyshev1型滤波器的频域响应函数与通带最大纹波以及阻带最小纹波之间的关系推导出了其滤波器最小阶数N的数学模型，然后通过Matlab仿真对该模型进行了检验，证明了它具有一定的可靠性。

六、心得体会

这是本学期数字信号处理课程的第二次课程设计，由于有了第一次的基础，所以这次的进度明显加快了很多，设计过程也进行得比较顺利。这次设计完成后数字信号处理的课程也已经结束了，所以此次课程设计可以说是对这学期学习这门课程的一个总结与回顾，也再一次巩固了我在这门课上学习到的知识，加深了我对数字信号处理的总体印象。不过由于对Matlab的掌握有限，所以和上次一样，这次设计中所使用的程序语句大部分还是参考了课件和网上查阅的资料，这点希望以后能逐步改进。