数字信号处理实验四 IIR滤波器设计.docx
《数字信号处理实验四 IIR滤波器设计.docx》由会员分享，可在线阅读，更多相关《数字信号处理实验四 IIR滤波器设计.docx（34页珍藏版）》请在冰豆网上搜索。
数字信号处理实验四IIR滤波器设计
实验四IIR滤波器设计
一、教学目的和任务
1．熟悉用双线性变换法设计IIR数字滤波器的原理和方法；
2．了解用脉冲响应不变法设计IIR数字滤波器的原理和方法；
3．掌握双线性变换及脉冲响应不变法设计的滤波器的频域特性，了解双线性变换法及脉冲响应不变法的特点；
4．掌握数字滤波器的计算机仿真方法。
二、实验原理介绍
IIR数字滤波器的系统函数为
的有理分式：
设计IIR滤波器的系统函数，就是要确定
的阶数N及分子分母多项式的系数
和
，使其
满足指定的频率特性。
由于模拟滤波器的设计有许多简单而严谨的设计公式和大量的图表可以利用，因此IIR滤波器设计的方法之一是：
先设计一个合适的模拟滤波器，然后将模拟滤波器通过适当的变换转换成满足给定指标的数字滤波器。
1、Butterworth模拟低通滤波器
幅度平方函数：
其中，N为滤波器的阶数，
为通带截止频率。
2．Chebyshev模拟低通滤波器
3、脉冲响应不变法原理
用数字滤波器的单位脉冲响应序列h（n）逼近模拟滤波器的冲激响应
，让h（n）正好等于
的采样值，即：
其中，T为采样间隔。
如果以
和H（z）分别表示
的拉氏变换及h（n）的Z变换，则：
4、双线性变换法原理
双线性变换法是通过两次映射采用非线性频率压缩的方法，将整个频率轴上的频率范围压缩到±π/T之间，再用
转换到z平面上，从而使数字滤波器的频率响应与模拟滤波器的频率响应相似。
5、设计IIR数字滤波器的步骤
1）确定数字滤波器的通带频率、阻带频率，通带最大衰减和阻带最小衰减。
2）计算对应的模拟低通滤波器的频率。
3）确定模拟低通滤波器的阶数N和3dB截止频率
。
4）模拟低通滤波器的系统函数H（s）。
5）由H（s）经过反归一化、脉冲响应不变法和双线性变换法确定数字低通滤波器的系统函数H（z）。
6）设计其它形式的滤波器时，由模拟低通到所需类型滤波器的频率域变换直接得到。
6*、MATLAB中用于IIR数字滤波器设计的函数
1）滤波器的特性分析
1Freqz函数：
求解数字滤波器的频率响应
[h,w]=freqz（b,a,n）：
返回数字滤波器的n点复频率响应，输入参数b和a分别是滤波器系数的分子和分母向量；输出参数h是复频率响应，w是频率点。
输入参数n默认是512。
Freqz（b,a,…）：
没有输出参数，直接在当前窗口中绘制频率响应的幅频响应和相频响应。
2Freqs函数：
求解模拟滤波器的频率响应
[h,w]=freqz（b,a,n）：
返回模拟滤波器的n点复频率响应，输入参数b和a分别是滤波器系数的分子和分母向量；输出参数h是复频率响应，w是频率点。
输入参数n默认是512。
3Abs、angle函数：
分别用于从复频域响应数据中提取幅值信息和相位信息
4Zplane函数：
绘制系统的零极点图
zplane（z,p）：
以单位圆为基准；z为系统的零点向量，图中用o表示；p为系统的极点向量，图中用x表示。
zplane（b,a）：
输入参数为系统传递函数的分子向量和分母向量。
2）确定滤波器最小阶数
函数
函数功能
[n,wn]=Buttord（wp,ws,rp,rs）
估计Butterworth滤波器阶数
[n,wn]=Cheb1ord（wp,ws,rp,rs）
估计ChebyshevⅠ型滤波器阶数
[n,wn]=Cheb2ord（wp,ws,rp,rs）
估计ChebyshevⅡ型滤波器阶数
[n,wn]=Ellipord（wp,ws,rp,rs）
估计椭圆滤波器阶数
wp：
归一化的通带截止频率；ws：
归一化的阻带截止频率
rp：
通带最大衰减量；rs：
阻带最小衰减量
n：
返回符合要求的滤波器阶数；wn：
返回滤波器的截止频率
3）模拟低通滤波器的设计
函数
函数功能
[z,p,k]=Buttap（n）
返回Butterworth滤波器的零点、极点、增益
[z,p,k]=Cheb1ap（n,rp）
返回ChebyshevⅠ型滤波器的零点、极点、增益
[z,p,k]=Cheb2ord（n,rs）
返回ChebyshevⅡ型滤波器的零点、极点、增益
[z,p,k]=Ellipap（n,rp,rs）
返回椭圆滤波器的零点、极点、增益
[b,a]=Butter（n,wn,’s’）
返回Butterworth滤波器的分子分母多项式的系数
[b,a]=Cheby1（n,rp,wn,’s’）
返回ChebyshevⅠ型滤波器的分子分母多项式的系数
[b,a]=Cheby2（n,rp,wn,’s’）
返回ChebyshevⅡ型滤波器的分子分母多项式的系数
[b,a]=Ellip（n,rp,rs,wn,’s’）
返回椭圆滤波器的分子分母多项式的系数
4）模拟滤波器的离散化
1Impinvar函数：
模拟滤波器变换成数字滤波器的脉冲响应不变法
[bz,az]=impinvar（b,a,fs）：
将模拟滤波器（b,a）变换成数字滤波器（bz,az）；输入参数fs是对模拟滤波器频率响应的采样，默认为1。
2bilinear函数：
模拟滤波器转换为数字滤波器的双线性变换法
[zd,pd,kd]=bilinear（z,p,k,fs）：
将采样零极点模型表达的模拟滤波器转换为数字滤波器。
列向量zd为零点向量，列向量pd为极点向量，kd为系统增益，fs是指定的采样频率。
[numd,dend]=bilinear（num,den,fs）：
将采用传递函数模型表达的模拟滤波器转换为数字滤波器。
[ad,bd,cd,dd]=bilinear（a,b,c,d,fs）：
将采用状态空间模型表达的模拟滤波器转换为数字滤波器。
5）直接设计IIR数字滤波器
函数
函数功能
[b,a]=Butter（n,wn）
设计巴特沃斯滤波器
[b,a]=Cheby1（n,r,wn）
设计切比雪夫Ⅰ型滤波器，r指定通带内波纹大小
[b,a]=Cheby2（n,r,wn）
设计切比雪夫Ⅱ型滤波器，r指定通带内波纹大小
[b,a]=Ellip（n,rp,rs,wn）
设计椭圆滤波器，rp指定通带内波纹最大衰减，rs指定通带内波纹的最小衰减
wn：
数字滤波器的截止频率；当wn是一个二元向量[w1,w2]时，返回一个2n阶的带通滤波器，通带为w1≤ω≤w2。
返回滤波器的分子分母多项式的系数，b为分子系数向量、a为分母系数向量（降幂排列）。
三、实验相关知识准备
1．模拟滤波器的设计函数
为了利用模拟滤波器设计数字滤波器，先必须设计对应的模拟滤波器，常用的模拟滤波器有：
bessel滤波器，butterworth滤波器，chebyshev
型滤波器，chebyshev
型滤波器及椭圆函数滤波器。
1》．模拟滤波器的设计函数
1）设计bessel模拟低通滤波器
[z,p,k]=besselap（n）
2）设计butterworth模拟低通滤波器
[z,p,k]=buttap（n）
3）设计chebyshev
型模拟低通滤波器
[z,p,k]=cheb1ap（n,Rp）
Rp：
通带内的波纹系数，单位分贝
4）设计chebyshev
型模拟低通滤波器
[z,p,k]=cheb2ap（n,Rs）
Rs：
阻带内的波纹系数低于通带Rs分贝
5）设计椭圆模拟滤波器
[z,p,k]=ellipap（n,Rp.Rs）
2》滤波器阶数的选择
下列函数除了能选择模拟滤波器的阶数外，同时也能选择数字滤波器的阶数。
1）选择butterworth滤波器阶数
数字域：
[n,Wn]=buttord（Wp,Ws,Rp,Rs）
模拟域：
[n,Wn]=buttord（Wp,Ws,Rp,Rs,’s’）
2）选择chebyshev
型滤波器阶数
数字域：
[n,Wn]=cheb1ord（Wp,Ws,Rp,Rs）
模拟域：
[n,Wn]=cheb1ord（Wp,Ws,Rp,Rs,’s’）
3）选择chebyshev
型滤波器阶数
数字域：
[n,Wn]=cheb2ord（Wp,Ws,Rp,Rs）
模拟域：
[n,Wn]=cheb2ord（Wp,Ws,Rp,Rs,’s’）
4）选择椭圆滤波器阶数
数字域：
[n,Wn]=ellipord（Wp,Ws,Rp,Rs）
模拟域：
[n,Wn]=ellipord（Wp,Ws,Rp,Rs,’s’）
注意：
n：
返回符合要求性能指标的数字滤波器或模拟滤波器的最小阶数
Wn：
滤波器的截至频率（即3db频率）
Wp：
通带的截至频率，Ws：
阻带的截至频率，单位rad/s。
且均为规一化频率，即
。
1对应π弧度。
频率规一化：
信号处理工具箱中使用的频率为莱奎斯特频率，根据香农定理，它为采样频率的一半，在滤波器设计中的截止频率均使用莱奎斯特频率进行规一化。
规一化频率转换为角频率，则将规一化频率乘以π。
如果将规一化频率转换为Hz，则将规一化频率乘以采样频率的一半。
3》．模拟频率变换
1）低通到低通模拟滤波器变换
[bt,at]=lp2lp（b,a,Wo）
lp2lp函数将截至频率为1rad/s（规一化截至频率）的模拟低通滤波器原型变换成截至频率为Wo的低通滤波器。
2）低通到带通模拟滤波器变换
[bt,at]=lp2bp（b,a,Wo,Bw）
lp2bp函数可将截止频率为1rad/s的模拟低通滤波器原型变换成具有指定带宽Bw和中心频率Wo的带通滤波器。
其中Wo=sqrt（W1*W2），Bw=W2-W1，W2高端截止频率，W1低端截止频率。
3）低通到高通模拟滤波器变换
[bt,at]=lp2hp（b,a,Wo）
4）低通到带阻模拟滤波器变换
[bt,at]=lp2bs（b,a,Wo,Bw）
5）besself模拟滤波器设计
[b,a]=besself（n,Wn）
[b,a]=besself（n,Wn,’ftype’）
[z,p,k]=besself（n,Wn）
[z,p,k]=besself（n,Wn,’ftype’）
2．由模拟滤波器变换成等效的数字滤波器
方法一：
双线性变换法实现模拟滤波器到数字滤波器的变换
双线性变换将S域映射成Z域，从而将模拟滤波器变换成等效的数字滤波器
[zd,pd,kd]=bilinear（z,p,k,Fs）为零极点增益表示的bilinear函数
其中z,p,k为S域传递函数的零点、极点、和增益，Fs为取样频率，zd,pd,kd为双线性变换后Z域传递函数的零点、极点和增益。
[numd,dend]=bilinear（num,den,Fs,Fp）为传输函数表示的bilinear函数
方法二：
冲激响应不变法实现模拟到数字滤波器的变换
[bz,az]=impinvar（b,a,Fs）
[bz,az]=impinvar（b,a）采用默认值为1hz的Fs。
小结：
设计IIR数字滤波器的一般步骤：
1）把给出的数字滤波器的性能指标转换为模拟滤波器的性能指标
2）根据转换后的性能指标，通过滤波器阶数选择函数，来确定滤波器的最小阶数N和固有频率Wn
3）由最小阶数N得到低通滤波器原型
4）由固有频率Wn把模拟低通滤波器原型转换为低通、高通、带通、带阻滤波器
5）运用脉冲响应不变法或双线性变换法把模拟滤波器转换成数字滤波器
MATLAB信号处理工具箱提供了几个用于直接设计IIR数字滤波器的函数，这些函数把上面这些复杂的步骤融合为一个整体，为设计滤波器带来了极大的方便。
3．直接设计IIR数字滤波器
1》Butterworth模拟和数字滤波器设计
数字域：
[b,a]=butter（n,Wn）可设计出截止频率为Wn的n阶butterworth滤波器
[b,a]=butter（n,Wn,’ftype’）当ftype=high时，可设计出截止频率为Wn的高通滤波器；当ftype＝stop时，可设计出带阻滤波器
[z,p,k]=butter（n,Wn）
[zp,k]=buter（n,Wn,’ftype’）
[A,B,C,D]=butter（n,Wn）
[A,B,C,D]=butter（n,Wn,’ftype’）
模拟域：
[b,a]=butter（n,Wn,’s’）可设计出截止频率为Wn的n阶模拟butterworth滤波器，
其余形式类似于数字域的。
2》chebyshev
型滤波器（通带等波纹）设计
数字域：
[b,a]=cheby1（n,Rp,Wn）可设计出n阶chebyshevI滤波器,其截止频率由Wn确定，通带内的波纹由Rp确定
[b,a]=cheby1（n,Rp,Wn,’ftype’）当ftype=high时，可设计出截止频率为Wn的高通滤波器；当ftype＝stop时，可设计出带阻滤波器
[z,p,k]=cheby1（n,Rp,Wn）
[zp,k]=cheby1（n,Rp,Wn,’ftype’）
[A,B,C,D]=cheby1（n,Rp,Wn）
[A,B,C,D]=cheby1（n,Rp,Wn,’ftype’）
模拟域：
[b,a]=cheby1（n,Rp,Wn,’s’）可设计出截止频率为Wn的n阶chebyshevI型模拟滤波器，其余形式类似于数字域的。
3》chebyshev
型滤波器（阻带等波纹）设计
数字域：
[b,a]=cheby2（n,Rs,Wn）可设计出n阶chebyshevI滤波器,其截止频率由Wn确定，阻带内的波纹由Rs确定
[b,a]=cheby2（n,Rs,Wn,’ftype’）当ftype=high时，可设计出截止频率为Wn的高通滤波器；当ftype＝stop时，可设计出带阻滤波器
[z,p,k]=cheby2（n,Rs,Wn）
[zp,k]=cheby2（n,Rs,Wn,’ftype’）
[A,B,C,D]=cheby2（n,Rs,Wn）
[A,B,C,D]=cheby2（n,Rs,Wn,’ftype’）
模拟域：
[b,a]=cheby2（n,Rs,Wn,’s’）可设计出截止频率为Wn的n阶chebyshev
型模拟滤波器，其余形式类似于数字域的。
4》椭圆滤波器设计
数字域：
[b,a]=ellip（n,Rp,Rs,Wn）可设计出n阶椭圆滤波器,其截止频率由Wn确定，通带内的波纹由Rp确定,阻带内的波纹由Rs确定。
[b,a]=ellip（n,Rp,Rs,Wn,’ftype’）当ftype=high时，可设计出截止频率为Wn的高通滤波器；当ftype＝stop时，可设计出带阻滤波器
[z,p,k]=ellip（n,Rp,Rs,Wn）
[zp,k]=ellip（n,Rp,Rs,Wn,’ftype’）
[A,B,C,D]=ellip（n,Rp,Rs,Wn）
[A,B,C,D]=ellip（n,Rp,Rs,Wn,’ftype’）
模拟域：
[b,a]=ellip（n,Rp,Rs,Wn,’s’）可设计出截止频率为Wn的n阶椭圆型模拟滤波器，其余形式类似于数字域的。
5》基于幅度特性设计IIR数字滤波器
[b,a]=yulewalk（n,f,m）可得到以b,a表示的n阶IIR滤波器，其幅频特性逼近由频率矢量f和幅度矢量m所给定的特性。
矢量f和n的长度必须相同。
6》基于频率特性设计IIR数字滤波器
[b,a]=invfreqz（h,w,nb,na）可得到以b,a表示的IIR滤波器，其幅频特性逼近由频率矢量w和复频率响应矢量h所给定的特性。
矢量w和h的长度必须相同，nb,na分别为指定的分子和分母多项式的阶数。
四、实验内容和步骤
1、IIR数字滤波器经典设计法
IIR数字滤波器经典设计法的一般步骤为：
（1）根据给定的性能指标和方法不同，首先对设计性能指标中的频率指标，如数字边界频率进行变换，转换后的模拟频率指标作为模拟滤波器原型设计的性能指标。
（2）估计模拟滤波器最小阶数和截止频率，利用MATLAB工具函数buttord、cheb1ord、cheb2ord、ellipord等。
（3）设计模拟低通滤波器原型。
利用MATLAB工具函数buttap、cheb1ap、cheb2ap、ellipap等。
（4）由模拟原型低通滤波器经频率变换获得模拟滤波器（低通、高通、带通、带阻等），利用MATLAB工具函数lp2lp、lp2hp、lp2bp、lp2bs。
（5）将模拟滤波器离散化获得IIR数字滤波器，利用MATLAB工具函数bilinear或impinvar。
后面的几个步骤在前面的章节中已经论述，这里主要介绍第一步，关于设计性能指标的转换。
设计IIR滤波器时，给出的性能指标通常分数字指标和模拟指标两种。
数字性能指标给出通带截止频率
，阻带起始频率
，通带波纹Rp，阻带衰减Rs等。
数字频率
和
的取值范围为0~
，单位弧度。
而MATLAB工具函数常采用归一化频率，
和
的取值范围为0~1，对应于0～
，此时需进行转换。
模拟性能指标给出通带截止频率
，阻带起始频率
，通带波纹Rp，阻带衰减Rs等。
模拟频率
和
单位为弧度/秒（rad/s）。
MATLAB信号处理工具箱中，设计性能指标的转换应根据不同设计方法进行不同处理。
下面就举例介绍这些方法。
冲激（脉冲）响应不变法采用
来将
和
变换为
和
【例1】用脉冲响应不变法设计一个Butterworth低通数字滤波器，使其特征逼近一个低通Butterworth模拟滤波器的下列性能指标：
通带截止频率
，通带波纹Rp小于3dB，阻带边界频率为
，阻带衰减大于15dB，采样频率Fs=10000Hz。
假设一个信号
，其中f1=1000Hz,f2=4000Hz。
试将原信号与通过该滤波器的输出信号进行比较。
%Samp6_5
Wp=2000*2*pi;Ws=3000*2*pi;%滤波器截止频率
Rp=3;Rs=15;%通带波纹和阻带衰减
Fs=10000;%采样频率
Nn=128;%调用freqz所用的频率点数
[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]=impinvar（b,a,Fs）;%运用脉冲响应不变法得到数字滤波器的传递函数
figure
（1）
[H,f]=freqz（bz,az,Nn,Fs）;%求解数字滤波器的幅频特性和相频特性
subplot（2,1,1）,plot（f,20*log10（abs（H）））
xlabel（'频率/Hz'）;ylabel（'振幅/dB'）;gridon;
subplot（2,1,2）,plot（f,180/pi*unwrap（angle（H）））
xlabel（'频率/Hz'）;ylabel（'相位/^o'）;gridon;
figure
（2）
f1=1000;f2=4000;%输入信号的频率
N=100;%数据长度
dt=1/Fs;n=0:
N-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（'输入信号'）%绘制输入信号
y=filtfilt（bz,az,x）;%用函数filtfilt对输入信号进行滤波
y1=filter（bz,az,x）;%用filter函数对输入信号滤波
subplot（2,1,2）,plot（t,y,t,y1,':
'）,title（'输出信号'）,xlabel（'时间/s'）
legend（'filtfilt','filter'）%加图例
由图可知，在小于2000Hz处的衰减小于3dB，而在大于3000Hz处衰减大于15dB，满足滤波器的设计指标。
由图6-5可见滤波器对含有1000Hz和4000Hz频率成分的信号进行了滤波，滤除了4000Hz的信号。
由程序的输出还可以看出，采用filtfilt函数，输出的1000Hz信号（实线）与输入1000Hz的信号相位一致，即经过滤波后并没有改变信号波形形状。
而运用filter函数滤波后（虚线）有一些延迟，改变了信号的形状。
图6-4例6-5所设计Butterworth滤波器的频率响应。
上图：
幅频特性；下图：
相频特性
图设计滤波器的输入和输出信号
输出信号中实线表示运用filtfilt函数的输出，虚线表示运用filter函数的输出
【例2】用脉冲响应不变法设计Butterworth低通数字滤波器，要求通带频率为
，通带波纹小于1dB，阻带在
内，幅度衰减大于15dB，采样周期T=0.01s。
假设一个信号
，其中f1=5Hz,f2=30Hz。
试将原信号与经过该滤波器的输出信号进行比较。
注意该题给出的通带边界频率和阻带边界频率均为数字频率，因此设计时首先要将其转换为模拟频率。
%Samp6_6
wp=0.2*pi;ws=0.3*pi;Rp=1;Rs=15;%数字滤波器截止频率、通带波纹和阻带衰减
T=0.01;Nn=128;%采样间隔
Wp=wp/T;Ws=ws/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]=impinvar（b,a,1/T）;%脉冲响应不变法设计数字滤波器传递函数
figure
（1）
[H,f]=freqz（bz,az,Nn,1/T）;%输出幅频响应和相频响应
subplot（2,1,1）,plot（f,20*log10（abs（H）））;
xlabel（'频率/Hz'）;ylabel（'振幅/dB'）;gridon;
subplot（2,1,2）,plot（f,180/pi*unwrap（angle（H）））
xlabel（'频率/Hz'）;ylabel（'相位/^o'）;gridon;
figure
（2）
f1=5;f2=30;%输入信号含有的频率
N=100;%数据点数
n=0:
N-1;t=n*T;%时间序列
x=sin（2*pi*f1*t）+0.5*cos（2*pi*f2*t）;%输入信号
subplot（2,1,1）,plot（t,x）,title（'输入信号'）
y=filtfilt（bz,az,x）;%对信号进行滤波
subplot（2,1,2）,plot（t,y）,title（'输出信号'）,xlabel（'时间/s'）
该例所要求的通带频率为
，而该滤波器的采样间隔为0.01s，采样频率为100Hz，即2
对应于100Hz，则0.2
对应于10Hz，即该滤波器的通带范围为0~10Hz。
观看图，其通带范围最大衰减小于1dB，与该分析一致。
题意要求阻带在
内，幅度衰减
升级会员 