双线性变换法设计数字切比雪夫带通IIR滤波器.docx

双线性变换法设计数字切比雪夫带通IIR滤波器.docx

《双线性变换法设计数字切比雪夫带通IIR滤波器.docx》由会员分享，可在线阅读，更多相关《双线性变换法设计数字切比雪夫带通IIR滤波器.docx（17页珍藏版）》请在冰点文库上搜索。

双线性变换法设计数字切比雪夫带通IIR滤波器

摘　要

在进行DSP系统设计时,往往先采用MATLAB软件等对算法进行仿真,确定最佳算法和参数。

利用MATLAB的信号处理工具箱可以直接设计数字滤波器,也可以建立模拟原型,离散化设计数字滤波器。

本文介绍了IIR数字带通滤波器的设计原理、步骤以及在MATLAB中的实现方法，并能应用到实际的信号处理中。

关键词:

IIR数字滤波器，数字信号处理（DSP），MATLAB，仿真

Abstract

WhentheDSPsystemisdesigned,weoftensimulatethealgorithmanddecidethebestalgorithmandparametersontheMATLABsoftware.WecandesignadigitalfiltereitherusingtheMATLABToolboxdirectlyorcreatingaanalogfilterprototypeandscateringthedigitalfilter.ThispaperintroducesthedesignprincipleofIIRdigitalbandpassfilter,thestepsaswellasthemethodsinMATLAB,anditalsocanbeappliedtotheactualsignalprocessing.

Keywords:

IIRdigitalfilters,Digitalsignalprocessing（DSP）,MATLAB,Simulation

目　录

Matlab课程设计

——双线性变换法设计数字切比雪夫带通IIR滤波器

1绪论

在现代通信系统中，由于信号中经常混有各种复杂成分，所以很多信号的处理和分析都是基于滤波器而进行的。

数字滤波器可以理解为是一个计算程序或算法，将代表输入信号的数字时间序列转化为代表输出信号的数字时间序列，并在转化过程中，使信号按预定的形式变化。

数字滤波器有多种分类，根据数字滤波器冲激响应的时域特征，可将数字滤波器分为两种，即无限长冲激响应（IIR）滤波器和有限长冲激响应（FIR）滤波器。

但是，传统的数字滤波器的设计使用繁琐的公式计算，改变参数后需要重新计算，从而在设计滤波器尤其是高阶滤波器时工作量很大。

利用MATLAB信号处理箱（SignalProcessingToolbox）可以快速有效地实现数字滤波器的设计与仿真。

   本文设计一个IIR数字带通滤波器。

IIR数字滤波器具有无限宽的冲激响应，与模拟滤波器相匹配，所以IIR滤波器的设计可以采取在模拟滤波器设计的基础上进一步变换的方法。

其设计方法主要有经典设计法、直接设计法和最大平滑滤波器设计法。

在对滤波器实际设计时，整个过程的运算量是很大的。

   设计完成后对已设计的滤波器的频率响应要进行校核，要得到幅频、相频响应特性，运算量也是很大的。

平时所要设计的数字滤波器，阶数和类型并不一定是完全给定的，很多时候要根据设计要求和滤波效果不断地调整，以达到设计的最优化。

在这种情况下，滤波器设计就要进行大量复杂的运算，单纯的靠公式计算和编制简单的程序很难在短时间内完成。

利用MATLAB强大的计算功能进行计算机辅助设计，可以快速有效地设计数字滤波器，大大地简化了计算量。

2IIR数字滤波器设计的原理与方法

2.1IIR数字滤波器设计的原理

IIR数字滤波器具有无限持续时间冲激响应,需要用递归模型来实现,其系统函数为:

（2.1）

设计IIR滤波器的任务就是寻求一个因果、物理上可实现的系统函数H（z），使其频率响应

满足所希望得到的频域指标，即符合给定的通带截止频率、阻带截止、通带衰减和阻带衰减.。

利用冲激响应不变法设计数字滤波器时可能会导致频域混叠现象，为了克服这一问题，需要找到由s平面到z平面的另外的映射关系，这种关系应保证：

1）s平面的整个jΩ轴仅映射为z平面单位圆上的一周；

2）若G（s）是稳定的，由G（s）映射得到的H（z）也应该是稳定的；

3）这种映射是可逆的，既能由G（s）得到H（z），也能由H（z）得到G（s）；

4）如果G（j0）=1，那么

＝1。

双线性Ｚ变换满足以上4个条件的映射关系，其变换公式为

（2.2）

双线性Ｚ变换的基本思路：

首先将整个Ｓ平面压缩到一条从－π/Ts变换到2π/Ts的横带里，然后通过标准的变换关系

将横带变换到整个Ｚ平面上去，这样就得到了Ｓ平面与Ｚ平面间的一一对应的单值关系。

图2.1双线性变换法S平面到Z平面的映射关系

2.2IIR数字滤波器设计的基本方法

IIR数字滤波器的设计方法有两类,一类是借助于模拟滤波器的设计方法设计出模拟滤波器,利用冲激响应不变法或双线性变换法转换成数字滤波器,再用硬件或软件实现；另一类是直接在频域或时域中进行设计,设计时需要计算机作辅助工具。

随着MATLAB软件尤其是MATLAB的信号处理工作箱的不断完善,不仅数字滤波器的计算机辅助设计有了可能,而且还可以使设计达到最优化。

IIR数字滤波器设计的基本步骤如下:

（1）根据任务,确定性能指标。

在设计一个滤波器之前,首先根据工程实际的需要确定滤波器的技术指标如：

边界频率：

ωp，ωs，ωc；阻带最小衰减As和通带最大衰减Rp；

（2）将数字滤波器的技术指标转换成模拟滤波器指标。

利用冲激响应不变法与双线性变换法进行频率间的转换，主要是边界频率Wp与Ws的转换。

（3）用模拟滤波器设计方法得到模拟滤波器的传输函数Ha（s）；可借助巴特沃斯（Butterworth）滤波器、切比雪夫（Chebyshev）滤波器、椭圆（Cauer）滤波器、贝塞尔（Bessel）滤波器等，这些滤波器都有严格的设计公式、现成的曲线和图表供设计人员使用。

（4）映射实现。

利用双线性变换法将模拟滤波器Ha（s）转换成数字滤波器H（z）。

（5）用有限精度算法实现这个系统函数H（z）（包括选择运算结构、选择合适的字长、有效数字处理方法）。

（6）用适当的软、硬件技术实现。

包括采用通用计算机软件、数字滤波器硬件或者软硬件结合，确定DF采用的结构将会影响其精度、稳定性、经济性及运算速度等很多重要性质。

3IIR带通滤波器的MATLAB设计

IIR带通滤波器的设计框图如下：

图3.1IIR带通滤波器的设计框图

3.1IIR带通滤波器的设计流程

图3.2IIR带通滤波器的设计流程

本文设计的IIR带通滤波器是从低通变换过来的，利用的是双线性变换以及切比雪夫II滤波器的原型，其具体的设计流程为上图所示。

首先根据题目要求确定带通滤波器的技术指标，先要进行频率的预畸变，并且归一化频率，再设计出切比雪夫II模拟低通滤波器，并求出其阶数等相关参数。

其次利用双线性变换法设计数字带通滤波器，，再调用函数进行双线性变换，并求出分子、分母的系数向量。

最后通过画图求出其幅频响应、相频响应、幅度特性曲线与零极点，并画出波形图。

最后进行验证，看所设计的滤波器能否达到要求的指标，若能达到，则说明该滤波器设计符合要求。

3.2IIR带通滤波器的设计步骤

（1）根据设计流程，首先确定所要设计的数字带通滤波器的相关指标：

①　通带截止频率wp1=0.4π，wp2=0.6π，通带最大衰减Rp=2dB；

②　阻带截止频率ws1=0.2π和ws2=0.8π，阻带最小衰减Rs=30dB；

③　取样间隔T=0.1s。

其实现程序如下（程序中pi代表π）：

Ts=0.1;Fs=1/Ts;％取样周期或频率

Rp=2;％通带最大衰减

Rs=30;％阻带最小衰减

wp1=0.4*pi;％通带、阻带上、下限截止频率

wp2=0.6*pi;

ws1=0.2*pi;

ws2=0.8*pi;

（2）频率的预畸变。

双线性变换中无法避免的一个问题即是频率的非线性偏移，因为数字频率的最大值为π，而模拟频率可以向无穷延伸，两者之间又要保持一一对应的映射关系。

双线性变换中的模拟角频率

与数字角频率

之间的关系为：

（3.1）

表明S平面与Z平面是单值的一一对应关系，即频率轴是单值变换关系。

虽然避免了脉冲响应不变法的频率响应的混叠现象，但是经过变换后，得到的幅频响应特性各分段边缘频率不能保持原来的比例关系，必须通过预修正加以校正。

做法是将数字频率

按

＝2/T*tan（w/2）的关系，变成模拟频率

，利用这组做过修正的模拟频率来设计模拟带通滤波器作为模拟原型。

Wp1=（2/T）*tan（wp1/2）;

Wp2=（2/T）*tan（wp2/2）;

Wp=[Wp1,Wp2];％模拟滤波器的通带截止频率

Ws1=（2/T）*tan（ws1/2）;

Ws2=（2/T）*tan（ws2/2）;

Ws=[Ws1,Ws2];％模拟滤波器的阻带截止频率BW=Ws2-Ws1;％模拟滤波器的带宽

Omegaw0=sqrt（Ws1*Ws2）;　％模拟滤波器的中心频率

（2）设计切比雪夫模拟低通滤波器。

％求模拟低通滤波器的阶数与边缘频率

[N,OmegaC]=cheb2ord（Wp,Ws,Rp,Rs,'s'）;

％求切比雪夫II型模拟低通滤波器的零、极点与增益

[z0,p0,k0]=cheb2ap（N,Rs）;

利用函数[N,OmegaC]=cheb2ord（Wp,Ws,Rp,Rs,’s’），通过给定滤波器的技术指标Wp、Ws、Rp、Rs，求得滤波器的阶数N与边缘频率OmegaC。

Wp、Ws、与OmegaC均在［0,1］区间归一化，以π弧度为单位。

利用函数[z,p,k]=cheb2ap（N,Rs），来设计一个阶数为N，阻带波动为Rs的归一化切比雪夫II型原型滤波器，得到左半平面零极点。

数组Z中返回零点，数组P中返回极点，并且返回增益K。

（3）设计归一化的模拟原型带通滤波器：

%求原型滤波器的分子系数

AnalogB=k0*real（poly（z0））;

%求原型滤波器的分母系数

AnalogA=real（poly（p0））;

%模拟低通到模拟带通的分子、分母系数的变换

[BandB,BandA]=lp2bp（AnalogB,AnalogA,Omegaw0,BW）;

%双线性变换：

模拟带通与数字带通的分子分母系数的变换

[DigitalB,DigitalA]=bilinear（BandB,BandA,Fs）;

%变为二阶节级联结构

[sos,G]=tf2sos（DigitalB,DigitalA）;

利用函数p=poly（A）来计算模拟滤波器的分子、分母系数向量，因其为实数，因此用real（）函数取其实部，即可得模拟滤波器的分子、分母系数向量。

这两个函数实现的功能可以用函数[Bs,As]=zp2tf（z,p,k）直接求得传递函数的分子、分母系数向量。

利用函数[BandB,BandA]=lp2bp（AnalogB,AnalogA,Omegaw0,BW），将模拟域的低通变为带通，并且得到模拟带通滤波器的分子、分母系数向量，Omegaw0取为中心频率，BW为带宽。

利用函数[DigitalB,DigitalA]=bilinear（BandB,BandA,Fs），双线性变换为数字带通滤波器的指标，如分子、分母的系数向量。

函数[sos,G]=tf2sos（DigitalB,DigitalA），即把z变换传递函数的直接形式转换成级联形式。

需要注意的是，这个函数是针对以z的负幂排列的多项式开发的，虽然可以推广到s域，但连续系统传递函数是按s的正幂排列的，要使两者一致，关键是使分子、分母系数向量同长，两序列中各元素的幂次排列一致

（4）求数字带通滤波器的幅频、相频特性、及其群延迟

％求数字带通滤波器的幅频特性

[Hz,Wz]=freqz（DigitalB,DigitalA,1024,'whole'）;

％将数字带通滤波器的幅频特性转化为分贝表示

dbHz=20*log10（（abs（Hz）+eps）/max（abs（Hz）））;

％求数字带通滤波器的相频特性

φ＝angle（Hz）

％求数字带通滤波器的群延迟特性

grd=grpdelay（DigitalB,DigitalA,Wz）;

函数[Hz,Wz]=freqz（DigitalB,DigitalA,1024,'whole'）可以求数字带通滤波器的幅频特性，而其幅度（即模值）的最大值可以归一化为1，则其模值（单位为dB）即可以用公式表示为dbHz=20*log10（（abs（Hz）+eps）/max（abs（Hz）））。

函数φ＝angle（Hz）可求得其相频特性，而对于一个滤波器来说，要满足其线性相位，则其群延迟

要为一常数，则其相位特性必须为一直线，即满足

。

因此，利用群延迟函数grd=grpdelay（DigitalB,DigitalA,Wz）可以判断所设计的滤波器是否是线性相位，如果不符合，可以更改参数加以较正或者用其他方法重新设计，从而方便了设计。

（5）作图

程序见附录。

以上是根据切比雪夫II型滤波器设计的带通滤波器，如果以切比雪夫I型滤波器设计，则只需改变几个参数则可，其程序如下所示（只附改变的参数与程序）。

％模拟滤波器的带宽

BW=Wp2-Wp1;

％模拟滤波器的中心频率

Omegaw0=sqrt（Wp1*Wp2）;

％求切比雪夫Ｉ型模拟低通滤波器的阶数与边缘频率

[N,OmegaC]=cheb1ord（Wp,Ws,Rp,Rs,'s'）;

％求切比雪夫I型模拟低通滤波器的零、极点与增益

[z0,p0,k0]=cheb1ap（N,Rp）;

4IIR带通滤波器的仿真结果及波形

根据上述对IIR切比雪夫I、II带通滤波器的设计过程，在Matlab软件中得到其幅频响应、相频响应、零极点图以及群延迟特性曲线如下图所示。

图4.1ChebyshevI型IIR带通滤波器特性

图4.2ChebyshevII型IIR数字带通滤波器特性

5IIR带通滤波器的仿真结果分析

根据上述的仿真波形，可以看出：

（1）由ChebyshevI型、ChebyshevII型低通原型变换成带通模型的幅频特性可以用分贝形式表示。

对于ChebyshevI型的幅频响应在通带内为等纹波衰减、阻带内为单调减小的，且通带（归一化）截止频率在波形上显示为：

0.4与0.6，阻带（归一化）截止频率在波形上显示为：

0.2与0.8，与设计要求基本一致。

对于ChebyshevII型的幅频响应来说，在通带内单调减小、阻带内为等纹波衰减的，而通带、阻带（归一化）截止频率为：

0.4与0.6，0.2与0.8，均满足设计要求。

（2）此次设计的IIR带通滤波器的阶数可由函数[N,OmegaC]=cheb1ord（Wp,Ws,Rp,Rs,'s'）

或[N,OmegaC]=cheb2ord（Wp,Ws,Rp,Rs,'s'）求得，N=3。

则从ChebyshevI型带通滤波器的幅频响应曲线上看，当N＝3时，在

（归一化后为0--1）范围内通带波动只有两个波谷，阻带则单调减小，这与理论内容一致。

对于ChebyshevII型带通滤波器来说，其幅频响应曲线的通带波动在

（归一化后为0--1）范围内阻带波动只有两个波峰，通带则单调减小，这与理论内容一致。

因此，对于滤波器阶数这一指标来说，满足要求。

（3）对于分贝化后的幅度特性，ChebyshevI型带通滤波器的通带最大衰减大约为2dB，阻带衰减最小为30dB，而ChebyshevII型带通滤波器的阻带最小衰减大约为30dB，通带最小衰减大约为2dB，这与带通滤波器的技术指标相同。

因此，在衰减性能上，所设计的滤波器达到了要求。

（4）对于相频响应来说，ChebyshevI型带通滤波器在通带内的相频响应曲线接近为一条平滑曲线，在阻带内则存在衰减、畸变。

ChebyshevII型带通滤波器在通带截止频率以内的相频响应曲线几乎为一条直线，但在通带截止频率处有很大角度的转折，使相频特性在通带内产生了的畸变，而在阻带内存在很大的畸变。

因此，对于双线性变换法，其相位特性得不到满足，必须用其它方法加以校正或采用其它方法重新设计滤波器，使之满足线性相位。

（5）群延迟是衡量一个滤波器或是整个系统性能指标的一个重要参数。

对于利用这两种低通滤波器原型模型设计的这两种带通滤波器，带通滤波器的延迟是同等带宽的低通滤波器延迟的两倍。

这个结果是由低通特性向带通特性变换引起的，n阶低通滤波器传递函数总是变换为2n阶带通滤波器传递函数，其滤波器的群延迟特性的中心频率与理论上的很接近。

（6）零、极点图则反映Z域的传递函数，因为从图上得知，ChebyshevI型带通滤波器在Ｚ域单位圆内的传递函数存在3对极点，其每一对则关于实轴对称，而且存在1对零点，且关于虚轴对称。

ChebyshevII型带通滤波器在Ｚ域单位圆内的传递函数存在3对极点，其每一对则关于实轴对称，而且存在3对零点，且关于虚轴对称。

从理论上分析，可以求出其传递函数。

综上所述，由于软件的精度要求以及所调用函数的特性，在误差允许范围内，上述的滤波器特性的仿真结果基本符合题目的要求，并且在一定的程度上设计的指标远远好于题目的要求。

因此，所设计的滤波器达到了题目的要求。

6总结

此次Matlab课程设计是在一定的理论基础之上进行的，在先修课程《信号与系统》与《数字信号处理》中，大量有Matlab设计方面的知识，所以做起来还比较容易，而且经过了自己的亲身实践，学到了许多实践方面的知识。

首先，在信号滤波系统中，有时因为模拟滤波器阶数太高，硬件占用空间太大为某些仪器的实现设置了障碍，而对于一些窄带情况下的低通滤波器用模拟手段往往很难实现。

在这些情况下，数字滤波器将会是一个很好的解决办法。

MATLAB信号处理工具箱提供了丰富而简便的设计、实现FIR和IIR的方法，使原来繁琐的程序设计简化成函数的调用，特别是滤波器的表达方式和滤波器之间的相互转换显得十分简便。

其次，IIR数字滤波器的设计和模拟滤波器的设计有着紧密的关系。

通常要先设计出适当的模拟滤波器，再通过一定的频带变换把它转换成为所需的数字IIR滤波器。

此外，任何数字信号处理系统中也还不可避免地用到模拟滤波器，比如A/D变换器前的抗混叠滤波器和D/A变换器后的平滑滤波器，因此模拟滤波器设计也是很重要的。

最后，在比较设计滤波器的方法上应该明确其技术指标以及某些参数的实际意义。

比如本文用双线性变换法设计数字带通滤波器时，必须先将频率归一化，并且进行频率预畸变，然后设计模拟滤波器，再利用频率变换法将模拟低通变为模拟带通，最后经过双线性变换法将模拟带通变换为数字带通。

如果不进行频率预畸变，那么设计出来的带通滤波器的幅频特性与相频特性将会产生很严重的畸变（如图5、图6所示，为设计时没有进行频率预畸变的特性），使设计的结果不满足给定的要求，在实际中会造成很严重的危害。

在做本次课程设计的过程中，我深深地感受到了自己所学到知识的有限，明白了只学好课本上的知识是不够的，要通过图书馆和互联网等各种渠道来扩充自己的知识。

在实验过程中我们曾经遇到过问题。

但是从中我们学习到了如何对待遇到的困难，进一步培养了我们一丝不苟的科学态度和不厌其烦的耐心。

所有的这些心得会对我以后的学习和工作有帮助作用，忠心感谢学校给我们提供这次实验机会。

图6.1　ChebyshevI型IIR数字带通滤波器特性（无频率预畸变）

图6.2ChebyshevII型IIR数字带通滤波器特性（无频率预畸变）

参考文献

［1］数字信号处理教程——MATLAB释义与实现（第二版）.陈怀琛编著.北京.电子工业出版社,2008.10

［2］MATLAB及在理工课程中的应用指南.陈怀琛编著.西安.电子科技大学出版社,2000

［3］数字信号处理原理与实现.刘泉,阙大顺编著.北京.电子工业出版社,2005.8

［4］数字信号处理教程（第二版）.北京.清华大学出版社,1997

［5］数字信号处理及其MATLAB实现.VinayK.Ingle编著.电子工业出版社,1998

［6］DigitalSignalProcessingLaboratoryUsingMatlab.SanjitK.Miltra编著.McGraw-Hill出版社,2000

附录：

原程序

％所设计的数字滤波器的指标

Ts=0.1，Fs=1/Ts，Rp=2，Rs=30;

wp1=0.4*pi，wp2=0.6*pi;

ws1=0.2*pi，ws2=0.8*pi;

％频率的预畸变

Wp1=（2/T）*tan（wp1/2）;

Wp2=（2/T）*tan（wp2/2）;

Wp=[Wp1,Wp2];％模拟滤波器的通带截止频率

Ws1=（2/T）*tan（ws1/2）;

Ws2=（2/T）*tan（ws2/2）;

Ws=[Ws1,Ws2];％模拟滤波器的阻带截止频率BW=Ws2-Ws1;％模拟滤波器的带宽

%BW=Wp2-Wp1;

Omegaw0=sqrt（Ws1*Ws2）;　％模拟滤波器的中心频率

%Omegaw0=sqrt（Wp1*Wp2）;

％求模拟低通滤波器的阶数与边缘频率

[N,OmegaC]=cheb2ord（Wp,Ws,Rp,Rs,'s'）;

%[N,OmegaC]=cheb1ord（Wp,Ws,Rp,Rs,'s'）

％求切比雪夫模拟低通滤波器的零、极点与增益

[z0,p0,k0]=cheb2ap（N,Rs）;

%[z0,p0,k0]=cheb1ap（N,Rp）;

%设计归一化的模拟原型带通滤波器

%求原型滤波器的分子系数

AnalogB=k0*real（poly（z0））;

%求原型滤波器的分母系数

AnalogA=real（poly（p0））;

%模拟低通到模拟带通的分子、分母系数的变换

[BandB,BandA]=lp2bp（AnalogB,AnalogA,Omegaw0,BW）;

%双线性变换：

模拟带通与数字带通的分子分母系数的变换

[DigitalB,DigitalA]=bilinear（BandB,BandA,Fs）;

%变为二阶节级联结构

[sos,G]=tf2sos（DigitalB,DigitalA）;

%求数字带通滤波器的幅频、相频特性、及其群延迟

％求数字带通滤波器的幅频特性

[Hz,Wz]=freqz（DigitalB,DigitalA,1024,'whole'）;

％将数字带通滤波器的幅频特性转化为分贝表示

dbHz=20*log10（（abs（Hz）+eps）/max（abs（Hz）））;

％求数字带通滤波器的相频特性

φ＝angle（Hz）

％求数字带通滤波器的群延迟特性

grd=grpdelay（DigitalB,DigitalA,Wz）;

％作图

subplot（2,3,1）;plot（Wz/pi,abs（Hz））;title（'幅频响应'）;

xlabel（''）,ylabel（'幅度:

|Hz|'）;axis（[0,1,0,1.1]）;

set（gca,'XTickMode','manual','XTick',[0