IIR数字带通滤波器设计.docx
《IIR数字带通滤波器设计.docx》由会员分享，可在线阅读，更多相关《IIR数字带通滤波器设计.docx（18页珍藏版）》请在冰豆网上搜索。
IIR数字带通滤波器设计
前言
随着信息时代和数字世界的到来，数字信号处理已成为当今一门极其重要的学科和技术领域。
目前数字信号处理在通信、语音、图像、自动控制、雷达、军事、航空航天、医疗和家用电器等众多领域得到了广泛的应用。
在数字信号处理中起着重要的作用并已获得广泛应用的是数字滤波器（DF，DigitalFilter）。
数字滤波器是一种用来过滤时间离散信号的数字系统，通过对抽样数据进行数学处理来达到频域滤波的目的。
MATLAB是英文MATrixLABoratory（矩阵实验室）的缩写。
它是美国的MathWorks公司推出的一套用于科学计算和图形处理可视化、高性能语言与软件环境。
它的信号处理工具箱包含了各种经典的和现代的数字信号处理技术，是一个非常优秀的算法研究与辅助设计的工具。
在设计数字滤波器时，通常采用MATLAB来进行辅助设计和仿真。
本次基课程设计将完成一个数字切比雪夫带通IIR滤波器的设计，利用双线性变换和无限冲激响应IIR原理完成设计，并利用MATLAB进行仿真。
工程概括
1.1IIR数字滤波器工作原理
数字滤波器是一个离散时间系统，输入x（n）是一个时间序列，输出y（n）也是一个时间序列。
如数字滤波器的系统函数为H（z），其脉冲响应为h（n），则在时间域内存在下列的关系。
在z域内，输入和输出存在下列关系：
式中，X（z）、Y（z）分别为输入x（n）和输出y（n）的z变换。
同样在频率域内，输入和输出存在下列关系：
式中，
为数字滤波器的频率特性;
和
分别为x（n）和y（n）的频谱。
为数字角频率，单位rad。
通常设计
在某些频段的响应值为1，在某些频段的响应为0。
和
的乘积在频率响应为1的那些频段的值仍为
，即在这些频段的振动可以无阻碍地通过滤波器，这些频带为通带。
和
的乘积在频率响应为0的那些频段的值不管
大小如何均为零，即在这些频段里的振动不能通过滤波器，这些频带称为阻带。
为数字角频率，单位为弧度（rad），
表示模拟角频率，单位弧度/秒（rad/s）。
数字角频率
在0～
范围内。
一个合适的数字滤波器系统函数H（z）可以根据需要改变输入x（n）的频率特性。
经数字滤波器处理后的信号y（n）保留信号x（n）中的有用频率成分，去除无用频率成分。
正文
2.1数字滤波器介绍
数字滤波器是具有一定传输选择特性的数字信号处理装置,其输入、输出均为数字信号,实质上是一个由有限精度算法实现的线性时不变离散系统。
它的基本工作原理是利用离散系统特性对系统输入信号进行加工和变换,改变输入序列的频谱或信号波形,让有用频率的信号分量通过,抑制无用的信号分量输出。
数字滤波器和模拟滤波器有着相同的滤波概念,根据其频率响应特性可分为低通、高通、带通、带阻等类型,与模拟滤波器相比,数字滤波器除了具有数字信号处理的固有优点外,还有滤波精度高（与系统字长有关）、稳定性好（仅运行在0与l两个电平状态）、灵活性强等优点。
时域离散系统的频域特性:
其中
、
分别是数字滤波器的输出序列和输入序列的频域特性（或称为频谱特性）,
是数字滤波器的单位取样响应的频谱，又称为数字滤波器的频域响应。
输入序列的频谱
经过滤波后
因此，只要按照输入信号频谱的特点和处理信号的目的，适当选择
，使得滤波后的
满足设计的要求，这就是数字滤波器的滤波原理。
数字滤波器根据其冲激响应函数的时域特性，可分为两种，即无限长冲激响应（IIR）数字滤波器和有限长冲激响应（FIR）数字滤波器。
IIR数字滤波器的特征是，具有无限持续时间冲激响应，需要用递归模型
来实现，其差分方程为：
系统函数为：
设计IIR滤波器的任务就是寻求一个物理上可实现的系统函数H（z），使其频率响应H（z）满足所希望得到的频域指标，即符合给定的通带截止频率、阻带截止频率、通带衰减系数和阻带衰减系数。
2.2数字滤波器的分类
按时间域特性，数字滤波器可以分为无限冲激（脉冲）响应数字滤波器（Infiniteimpulseresponsedigitalfilter,简称IIR滤波器）和有限冲激（脉冲）响应数字滤波器（Finiteimpulseresponsedigitalfilter,简称FIR滤波器）两类。
IIR滤波器的传递函数为：
h（n）为滤波器的脉冲响应，n=0~
均有值。
M和N为分解的分子和分母多项式的系数个数。
FIR滤波器的传递函数为：
该滤波器的脉冲响应h（n）在n=0,1,…,N-1的有限个点（N个点）上有值。
式中分母
全为零时，H（z）具有全零点形式，IIR滤波器退化为FIR滤波器。
按频率特性来讲，数字滤波器和模拟滤波器一样可分为低通、高通、带通和带阻等。
数字滤波器是一个离散时间系统，在频率特性中具有周期性，因此我们讨论的频率范围仅在
的范围内，相应的归一化频率在0~1之间，
和1对应于Nyquist频率。
和模拟滤波器一样，理想数字滤波器的频率特性
在通带内必须满足：
式中，K,
均为常数。
和模拟滤波器一样，数字滤波器的设计目的是使滤波器的频率特性达到所给定的性能指标。
其性能指标也包括通带波纹Rp（dB）、阻带衰减Rs（dB）、通带边界频率（Hz）、阻带边界频率（Hz）等。
2.3脉冲响应不变法
所谓脉冲响应不变法就是使数字滤波器的脉冲响应序列h（n）等于模拟滤波器的脉冲响应ha（t）的采样值，即
式中，T为采样周期。
因此数字滤波器的系统函数H（z）可由下式求得：
Z[-]表示对[-]的内容进行Z变换，Z变换的内容请参考相应的数字信号处理教材。
如果已经获得了满足性能指标的模拟滤波器的传递函数
，求与之对应的数字滤波器的传递函数H（z）的方法是：
（1）求模拟滤波器的单位脉冲响应
。
式中，
表示对
的Laplace逆变换。
Laplace变换内容请参考高等数学的积
分变换或信号处理教材。
（2）求模拟滤波器单位冲激响应
的采样值，即数字滤波器冲激响应序列h（n）。
（3）对数字滤波器的冲激响应h（n）进行z变换，得到传递函数H（z）。
由上述方法推论出更直接地由模拟滤波器系统函数
求出数字滤波器系统函数H（z）的步骤是：
（1）利用部分分式展开将模拟滤波器的传递函数H（s）展开成：
在MATLAB中这步可通过residue函数实现。
若调用residue函数的形式为[R,P,K]=residue（a,b）形式，则将下式（传递函数形式）：
变换为：
这种形式为极点留数商向量形式，对于本节所讲的特定情况，K为空矩阵。
若为[b,a]=residue（R,P,K）则为上面调用形式的反过程。
（2）将模拟极点pk变换为数字极点
即得到数字系统的传递函数
其中T为采样间隔。
（3）将
转换为传递函数形式，在该步骤中，可采用[R,P,K]=residue（b,a）.
，MATLAB中已经提供了冲激响应不变法设计数字滤波器的函数，调用格式为：
[bz,az]=impinvar（b,a[[,Fs[，Fp]）
式中，b,a为模拟滤波器分子和分母多项式系数向量；Fs为采样频率（所滤波数据），单位Hz，缺省时为1Hz。
Fp为预畸变频率（Prewarpedfrequency）,是一个“匹配”频率,在该频率上，频率响应在变换前后和模拟频率可精确匹配。
一般设计中可以不考虑。
bz,az分别为数字滤波器分子和分母多项式系数向量。
前面已提到过，函数输入变量中的[]表示可添加也可略去的内容。
下面我们用例子说明如何使用这个函数。
（1）脉冲响应不变法将模拟滤波器
变换为数字滤波器H（z），采样周期为T=0.1s。
%Samp6_1
b=[32];a=[231];T=0.1;%模拟滤波器分子和分母多项式系数及采样间隔
[bz1,az1]=impinvar（b,a,1/T）
程序输出为：
bz1=0.3000-0.2807
az1=2.0000-3.71211.7214
在应用冲激响应不变法设计数字滤波器时要注意它的特点。
脉冲响应不变法由
这一基本关系得到数字角频率
和模拟角频率
满足
线性变换关系，T为采样间隔。
这使得
轴上每隔
便映射到z域中的单位圆一周。
如果模拟滤波器频率响应是有限带宽的话，通过变换得到的数字滤波器的频率响应非常接近于模拟滤波器的频率响应。
由于数字滤波器的频率响应是模拟滤波器频率响应的周期延拓，因此对于高通和带阻滤波器存在混叠效应，会造成频率响应失真，因此这种方法原则上只适用于有限带宽滤波器。
对于高通、带阻等滤波器，由于它们高频成分不衰减，势必产生严重的混迭失真。
双线性变换法可以弥补这方面的不足。
2.4双线性变换法
双线性变换法将s平面的整个频率轴映射到z域的一个频率周期中。
因此s平面到z平面的映射是非线性的，其单值双线性映射关系为：
式中，T为采样周期。
因此若已知模拟滤波器的传递函数
，将（6-14）式的第一式代入
即可得到数字滤波器的传递函数：
在双线性变换中，模拟角频率和数字角频率存在下面关系：
可见，模拟角频率
和数字角频率
之间的关系是非线性的。
在MATLAB中，函数bilinear采用双线性变换法实现模拟s域至数字z域的映射，直接用于模拟滤波器变换为数字滤波器。
其调用方式为：
[zd,pd,kd]=bilinear（z,p,k,Fs）
[numd,dend]=bilinear（num,den,Fs）
式中，z,p分别为模拟滤波器零点、极点列向量；k为模拟滤波器的增益；Fs为采样频率，单位Hz。
zd,pd,kd为数字滤波器的零极点和增益。
num,den分别为模拟滤波器传递函数分子和分母多项式系数向量，模拟滤波器传递函数具有下面的形式：
numd和dend分别为数字滤波器传递函数分子和分母多项式系数向量。
（2）用双线性变换法将模拟滤波器
变换为数字滤波器H（z），采样周期（间隔）T=0.1s。
%Samp6_2
b=[32];a=[231];T=0.1;%模拟滤波器分子和分母多项式的系数，采样间隔
[bz1,az1]=bilinear（b,a,1/T）%将模拟滤波器传递函数转换为数字滤波器传递函数
程序输出为：
bz1=0.07200.0046-0.0674
az1=1.0000-1.85600.8606
双线性变换法克服了脉冲响应不变法的频谱混迭问题，其幅值逼近程度好，可适用于高通、带阻等各种类型滤波器的设计。
s域和z域对应关系也简单。
缺点是频率变换的非线性导致数字滤波器与模拟滤波器在幅度和频率的对应关系上发生畸变。
但一般滤波器的幅频响应具有分段常数的特点，即滤波器允许某一频段信号通过，而不允许另外频段的信号通过的特点，故变换后这一特点仍保留，影响不大。
由数字边界频率计算模拟边界频率时，不是按线性关系进行的，这就是所谓的“预畸变”。
但如果给定预畸变频率为边界频率，经预畸变频率校正则可以保证所要设计的模拟边界频率精确映射在所要求的数字边界频率上。
2.5滤波器的特性及使用函数
（1）freqz
对于模拟滤波器，可以用freqs求解滤波器的频率响应。
与之对应的函数freqz用于求数字滤波器的频率响应，其调用格式为：
[[h,w]=]freqz（b,a,n[,’whole’]）;或[h,f]=freqz（b,a,n[,’whole’],Fs）;
式中，b,a为数字滤波器分子和分母多项式的系数，n为复数频率的响应点数，为整数，最好为2的幂，缺省时为512；Fs为采样频率，单位Hz。
如果给定该值，则f位置输出为频率Hz，若没有给定，则按角频率（Angularfrequency）给定f的频率矢量；’whole’表示返回的频率f或w值包含z平面整个单位圆频率矢量，即0~2
；缺省时，频率f或w值包含z平面上半单位圆（0~
）之间等间距n个点频率矢量。
h为复频率响应；w为n点频率向量（单位rad）；f为n点频率向量（Hz）。
函数返回值缺省时，绘制幅频响应和相频响应图。
该函数适用于下面形式的数字滤波器：
函数freqz输出的频率向量在0~
。
为了获得一个滤波器真正的相频特性图，要对相位角进行解缠绕。
为此MATLAB提供了一个函数unwrap来解决这个问题，P=unwrap（angle（H））。
（2）impz
impz用于产生数字滤波器的脉冲响应。
调用格式为：
[[h,t]=]impz（b,a[,n,Fs]）
式中，b,a分别为滤波器分子和分母多项式系数向量；n为采样点数；Fs为采样频率，缺省值为1；h为滤波器单位脉冲响应向量；t为和h对应的时间向量。
当函数输出缺省时，绘制滤波器脉冲响应图；当n缺省时，函数自动选择n值。
（3）零极点图
滤波器的零极点位置决定了滤波器稳定性和性能，因此考察滤波器的零极点的位置是分析滤波器特性的重要方面之一。
MTALAB信号处理工具箱提供绘制数字滤波器零极点位置图的工具zplane,调用格式为：
zplane（z,p）或zplane（b,a）
式中，z,p为零极点向量（为复数），b,a为滤波器分子和分母多项式的系数（为实数）。
函数在z平面绘出零点和极点。
极点用’×’表示，零点用’o’表示。
（4）群延迟
信号传输的不失真条件之一为：
滤波器相频特性是一条经过原点的直线，即
，
为常数。
但一般滤波器不满足这个条件，衡量实际滤波器相位平均延迟的物理量是群延迟。
群延迟定义为信号通过滤波器的延迟随频率变化的函数，即滤波器相频特性图上切线的负斜率：
MATLAB信号处理工具箱提供计算群延迟函数grpdelay，调用格式为：
[gd,w]=grpdelay（b,a,n[,’whole’]）
[gd,f]=grpdelay（b,a,n[,’whole’],Fs）
gd=grpdelay（b,a,w）
gd=grpdelay（b,a,f,Fs）
grpdelay
其中，gd为群延迟；其他各项意义同函数freqz,函数输出项缺省时，绘制群延迟图。
’whole’参数表示绘制包括大于Nyquist频率的一个周期的群延迟。
（5）filter函数用来实现数字滤波器对数据的滤波，函数的调用格式为：
y=filter（b,a,x）
其中，b,a分别为滤波器传递函数H（z）的分子和分母多项式系数。
x为滤波器的输入。
y为滤波器的输出。
y为与x具有相同大小的向量。
（6）filtfilt函数实现零相位前向与后向结合的滤波。
调用格式为：
y=filtfilt（b,a,x）
式中，b,a分别为滤波器传递函数H（z）的分子和分母多项式系数。
x为滤波器的输入，为值向量。
y为滤波器的输出。
该函数对序列x进行正常的正向滤波后，将滤波后的输出翻转重新用该滤波器进行滤波，第二次滤波后的输出序列的翻转即得到零相位的滤波输出。
这样就可以把延迟后的相位校正至零。
但该函数只能用于数字滤波器，FIR滤波器或IIR滤波器均能使用。
3.1设计步骤
根据以上IIR数字滤波器设计方法，下面运用双线性变换法基于MATLAB设计一个IIR带通滤波器，其中带通的中心频率为wp0=0.55π，;通带截止频率wp1=0.45π,wp2=0.65π;通带最大衰减Ap=1dB;阻带最小衰减As=40dB;阻带截止频率ws2=0.75π
（1）确定性能指标
在设计带通滤波器之前,首先根据工程实际的需要确定滤波器的技术指标:
通带截止频率wp1=0.45π,wp2=0.65π;阻带截止频率ws1=0.3π，ws2=0.75π；阻带最小衰减As=40dB和通带最大衰减Ap=1dB;中心频率wp0=0.55π。
（2）频率预畸变
用Ω=2/T*tan（w/2）对带通数字滤波器H（z）的数字边界频率预畸变,得到带通模拟滤波器H（s）的边界频率主要是通带截止频率Wp1,Wp2;阻带截止频率Ws1，Ws2的转换。
双线性变换法一般T=2s。
通带截止频率Wp1=（2/T）*tan（wp1/2）
Wp2=（2/T）*tan（wp2/2）
阻带截止频率Ws1=（2/T）*tan（ws1/2）
Ws2=（2/T）*tan（ws2/2）
（3）模拟带通性能指标转换成模拟低通性能指标
BW=Wp2-Wp1；%带通滤波器的通带宽度
W0=Wp1*Wp2；
WP=1；%归一化处理
WS=WP*（W0^2-Ws1^2）/（Ws1*BW）；
（4）模拟低通滤波器的构造
借助切比雪夫（Chebyshev）滤波器得到模拟低通滤波器的传输函数Ha（s）。
（5）模拟低通滤波器转换成模拟带通滤波器
调用lp2bp函数将模拟低通滤波器转化为模拟带通滤波器。
（6）模拟带通滤波器转换成数字带通滤波器
利用双线性变换法将模拟带通滤波器Ha（s）转换成数字带通滤波器H（z）。
（7）输入信号检验滤波器性能
输入不同频率的正弦波，观察输出波形，检验滤波器性能。
3.2程序流程图
开始
↓
读入数字滤波器技术指标
↓
将指标转换成归一化模拟低通滤波器的指标
↓
设计归一化的模拟低通滤波器阶数N和3db截止频率
↓
模拟域频率变换，将G（P）变换成模拟带通滤波器H（s）
↓
用双线性变换法将H（s）转换成数字带通滤波器H（z）
↓
输入信号后显示相关结果
↓
结束
程序流程图
3.3MATLAB程序
（1）数字滤波器部分
clear
%数字滤波器的技术指标
wp1=0.45*pi;wp2=0.65*pi;wp0=0.55*pi;
ws1=0.3*pi;ws2=0.75*pi;
Ap=1;As=40;T=2;
%带通到低通的频率变换
Wp1=（2/T）*tan（wp1/2）;Wp2=（2/T）*tan（wp2/2）;
Ws1=（2/T）*tan（ws1/2）;Ws2=（2/T）*tan（ws2/2）;Wp=（2/T）*tan（wp0/2）;
BW=Wp2-Wp1;%带通滤波器的通带宽度
W0=Wp1*Wp2;
WP=1;%归一化处理
WS=WP*（W0^2-Ws1^2）/（Ws1*BW）;
%切比雪夫模拟低通原型滤波器设计
[N,Wn]=cheb1ord（WP,WS,Ap,As,'s'）;
[B1,A1]=cheby1（N,Ap,Wn,'s'）;
%模拟低通原型滤波器幅频特性曲线（dB）
[h1,w1]=freqs（B1,A1）;
subplot（3,2,1）;plot（w1/pi,20*log（abs（h1）））;gridon;
xlabel（'w（rad）'）;
ylabel（'|H（jw）|.dB'）;
title（'模拟低通滤波器幅频特性曲线'）;
%由模拟低通原型滤波器变换为模拟带通滤波器
[B2,A2]=lp2bp（B1,A1,Wp,BW）;
[h2,w2]=freqs（B2,A2）;
%模拟带通滤波器幅频特性曲线（dB）
subplot（3,2,2）;plot（w2,20*log（abs（h2）））;axis（[0,3,-400,50]）;gridon;
xlabel（'w（rad）'）;
ylabel（'|H（jw）|.dB'）;
title（'模拟带通滤波器幅频特性曲线'）;
%双线性变换：
由模拟滤波器向数字滤波器的变换
[B3,A3]=bilinear（B2,A2,0.5）;
[h,w]=freqz（B3,A3,64）;
phz=unwrap（angle（h））;
%数字带通滤波器幅频响应曲线
subplot（3,2,3）;plot（w/pi,abs（h））;
xlabel（'w（rad）'）;ylabel（'|H（z）|'）;
title（'数字带通滤波器幅频特性曲线'）;
%数字带通滤波器幅频响应曲线（dB）
subplot（3,2,4）;plot（w/pi,20*log10（abs（h）））;axis（[-1,2,-250,50]）;
xlabel（'w（rad）'）;ylabel（'|H（z）|.dB'）;
title（'数字带通滤波器幅频特性曲线'）;
%数字带通滤波器相频特性曲线（dB）
subplot（3,2,5）;plot（w/pi,phz）;gridon;
xlabel（'w（rad）'）;
ylabel（'H（z）'）;
title（'数字带通滤波器相频特性曲线'）;
（2）输入正弦波检验性能部分
%输入正弦波波验证滤波器特性
n=0:
600;t=n/11000;
x1=2*sin（2*pi*2750*t）;%正弦波信号
figure;
subplot（121）;plot（x1）;gridon;%500Hz正弦波波形
axis（[0,10*pi,-5,5]）;
xlabel（'t（s）'）;
ylabel（'x1'）;
title（'正弦波信号'）;
y1=filter（B3,A3,x1）;%数字滤波函数输出
subplot（122）;plot（y1）;gridon;%数字滤波器输出波形
axis（[0,10*pi,-3,3]）;
xlabel（'f（hz）'）;
ylabel（'y'）;
title（'数字滤波器输出波形'）;
注:
应输入一系列不同频率的正弦波进行验证，只需将x1=2*sin（2*pi*2750*t）中的2750频率值改变即可，这里取500HZ、2750HZ和5000HZ进行验证。
3.4仿真结果
3.4.1滤波器性能仿真
源程序设计了模拟低通滤波器、模拟带通滤波器与数字带通滤波器等滤波器，对各部分滤波器的性能仿真如下，下面五个图分别为模拟低通原型滤波器幅频特性曲线、模拟带通滤波器幅频特性曲线、数字带通滤波器幅频响应曲线（有两个）、数字带通滤波器相频特性曲线，可以看到各部分滤波器波形基本满足设计要求。
3.1滤波器性能仿真
3.4.2滤波器性能验证
为了验证滤波器性能，可以输入一些列不同频率的正弦波，观察通阻状态，分别输入500HZ、2750HZ和5000HZ进行验证仿真后的结果如下：
（1）输入为500HZ的正弦波信号
3.2输入500ＨＺ正弦波波形
仿真的结果显示当频率为500Hz时经过滤波器后输出波形接近零说明频率是500Hz未达到滤波器的频率所以滤波器将其滤除。
（2）输入为2750HZ的正弦波信号
3.3输入2750ＨＺ正弦波波形
仿真的结果显示当频率为2750Hz时经过滤波器后输出波形是符合滤波器的设计说明频率是2750Hz是符合滤波器的频率所以滤波器将其完整输出。
（3）输入为5000HZ的正弦波信号
3.4输入5000ＨＺ正弦波波形
仿真的结果显示当频率为500Hz时经过滤波器后输出波形接近零说明频率是500Hz超过到滤波器的频率所以滤波器将其滤除。
3.5总结
这次课设下来，对设计带通数字滤波器的整个过程有了很好的掌握，懂得了设计滤波器的基本方法，对双线性变换法，切比雪夫滤波器有了一定了解，同时呢也熟悉了MATLAB的环境，巩固了相关知识。
最大的收获是初步了解了数字滤波器的原理及设计方法，加深了对滤波器的认识，一切从零开始，虽然没有以绝对完美结束，但在这么短的时间内能够设计成功已经出乎意料之外了，总之，收获还是很大的。
致谢
通过一个半星期的努力，首先感谢老师给我这次动手
升级会员 