无限长单位脉冲响应IIR滤波器的设计方法MATLAB

数字信号处理实验指导
实验四、 无限长单位脉冲响应(IIR)滤波器的设计方法
（一） 实验目的
加深对无限冲激响应（ IIR ）数字滤波器的常用指标和设计过程的理解。

（二） 实验内容
常用函数介绍：
1、Matlab 信号处理工具箱中提供了设计巴特沃思模拟滤波器的函数buttord 、buttap 和butter ，格式如下：
（1）[,](,,,,C P S P S N W buttord W W R R s ='')
用于计算巴特沃思模拟低通滤波器的阶N 和3dB 截止频率Wc (即本书中的符号c Ω)。

其中，Wp 和Ws 分别是滤波器的通带截止频率p Ω和阻止截止频率s Ω，单位为rad/s ；Rp 和Rs 分别是通带最大衰减系数p α和阻带最小衰减系数s α，单位为dB 。

（2）[,,]()z p G buttap N =
用于计算N 阶巴特沃思归一化（c Ω=1）模拟低通滤波器系统函数的零、极点和增益因子，返回长度为N 的向量z 和p 分别给出N 个零点和极点，G 是滤波器增益。

得到的滤波器系统函数形式如下：
1212()()()()()()()()()
a N a a N Q s s z s z s z H s G P s s p s p s p ---==--- 其中，k z 和k p 分别是向量z 和p 的第k 个元素。

如果要从零、极点得到系统函数的分子和分母多项式系数向量B 和A ，可以调用结构转换函数
（3）[,]2(,,)B A zp tf z p G =，结构转换后系统函数的形式为
111111()()()M M M a N N N
b s b s b B s H s A s a s a s a ----+++==+++ 其中，M 是向量B 的长度，N 是向量A 的长度，k k b a 和分别是向量B 和A 的第k 个元素。

（3）[,](,,,)C B A butter N W ftype s =''''
用于计算巴特沃思模拟滤波器系统函数中分子和分母多项式系数向量B 和A ，其中N 和C W 分别是滤波器的阶和3dB 截止频率c Ω，返回向量B 和A 中的元素k a 和k b 分别是上面的()a H s 表示式中的分母和分子系数。

ftype 缺省时，设计低通滤波器;ftype=high 时，设计高通滤波器；ftype=stop 时，设计带阻滤波器，此时C W 为向量[1,]c cu W W ，且ftype 缺省时，设计带通滤波器，带通的频率区间为1c cu W W 〈Ω〈。

S 缺省时，设计数字滤波器。

例如：设计一个满足下列指标要求的巴特沃思模拟低通滤波器。

指标
s rad dB s rad d s s p p /4,30,/2,1=Ω≥=Ω≤ααB
Matlab 程序：
2;4;1;30[,](,,,,[,](,,)00.0016[](p s p s C P S P S C W W R R N W buttord W W R R s B A butter N W s W H W freqs ====='')=''==； %滤波器参数
； %计算滤波器的阶和3dB 截止频率
； %计算滤波器系统函数分子、分母多项式系数：：； %在0~6rad/s 频率范围内取1001个频率点
，,,);20*log10(());(,),(//);(/)
B A W W H abs H B
plot W H grid xlabel rad s ylabel dB ='''' %计算频率向量上的滤波器频率响应 %纵坐标的单位取d on %绘制幅频响应曲线
频率（幅度
运行结果如下：
N=6
Wc=2.2496
B= 0 0 0 0 0 0 129.5917 A=1.0000 8.6916; 37.7720; 104.0667;191.1447 ;222.5973 ; 129.5917 即：
65432129.5917()8.691637.7720104.0667191.1447222.5793129.5917
a H s s s s s s s =++++++幅频响应曲线如图6.3.4所示。

2、Matlab 中设计数字滤波器的函数都是采用双线性变换法，将模拟滤波器转换为数字滤波器。

这些函数及其凋用格式如下（巴特沃思数字滤波器）：
（1）[,](,,,)C P S P S N W buttord W W R R =
该格式用于计算巴特沃思数字滤波器的阶N 和3dB 截止频率的归一化值C
W
（关于π归一化）。

调用参数Wp 和Ws 分别是数字滤波器的通带截止频率和阻带截止频率的归一化值（关于π归一化），要求01p W ≤≤和01s W ≤≤，其中1表示数字频率π（对应模拟频率/2s f ，s f 为采样频率）。

p R 和s R 分别是通带最大衰减和阻带最小衰减，单位dB 。

当s p W W ≤时，设计高通滤波器；当p W 和s W 是二元向量时，设计带通（s p W W ≤）或带阻（s p W W ≥）滤波器，这时返回参数C W 也是二元向量。

（2）[,](,,)C B A butter N W ftype =''；
该格式用于计算巴特沃思模拟滤波器系统函数中分子和分母多项式系数向量B 和A 。

调用参数N 和W C 分别是巴特沃思数字滤波器的阶和3dB 截止频率归一化
值（关于π归一化）。

当ftype 缺省时，设计低通滤波器；当ftype=high 时，设计高通滤波器；ftype=stop 时，设计带阻滤波器，此时C W 为二元向量[1,]c cu W W ，1c W 和cu W 分别是带阻滤波器的通带3dB 下，上截止频率的归一化值（关于π归一化）；C W 为向量[1,]c cu W W ，且ftype 缺省时，设计带通滤波器，带通的频率区间为1c cu W W ω〈〈。

注意设计出的带通和带阻数字滤波器是2N 阶的，这是因为带通滤波器可表示为一个N 阶低通滤波器与一个N 阶高通滤波器的级联。

由函数的返回向量B 和A 可写出数字滤波器的系统函数为：
1112111121()N
N N N N N b b z b z b z H z a a z a z a z
---+---+++++=++++ 其中，k k b a 和分别是向量B 和A 的第k 个元素。

例如：用双线性变换法设计巴特沃思数字低通滤波器，指标要求
通带截止频率0.2p w rad π=，通带最大衰减1p dB α=；
阻带截止频率0.35s w rad π=，阻带最小衰减10s dB α=
解 设计步骤如下：
（1）给出数字滤波器的指标。

0.2;0.35;10;P S S W W R ===
（2）计算巴特沃思数字滤波器的阶N 和3dB 截止频率。

[,](,,,)C P S P S N W buttord W W R R =；
（3）用双线性变换法设计巴特沃思数字低通滤波器。

[,](,)Z Z C B A butter N W =；
运行结果如下：
Z B =0.0335 0.1006 0.1006 0.0335
Z A =1.000 -1.4245 0.8827 -1.1900
设计出的数字滤波器系统函数为
123
123
0.03350.10060.10060.0335() 1.000 1.42450.88270.1900z z z H z z z z ------+++=-+- 程序如下：
0.2;0.35;10;10;[,](,,,);[,](,);P S P S C P S P S Z Z C W W R R N W buttord W W R R B A butter N W ====== %滤波器参数
%计算数字滤波器的阶和3dB 截止频率
%计算数字滤波器系统函数分子和分母多项 式的系数向量B 和A
%省去绘图部分程序
例：利用脉冲响应不变法，把系统函数为6
51)(2+++=
s s s s H a 的模拟滤波器变换成等价的数字滤波器，采样间隔S T 1.0= 手算：2132651)(2+++=+++=s s s s s s H a 2
11
12136065.05595.1108966.01.0112)(-------+--=---==z z z z e T z e T z H T T Matlab:
B=[1,1];
A=[1,5,6];
T=0.1;
Fs=1/T;
[Bz,Az]=impinvar(B,A,Fs);%用脉冲响应不变法将模拟滤波器变换成数字滤波器
运行结果：
Bz=0.1000 -0.0897
Az=1.000 -1.5595 0.6065
1、设采样周期T=250μs（采样频率f s =4kHz），用脉冲响应不变法和双线性变换法设计一个三阶巴特沃兹滤波器，其3dB边界频率为f c =1kHz。

[B,A]=butter(3,2*pi*1000,'s');
[num1,den1]=impinvar(B,A,4000);
[h1,w]=freqz(num1,den1);
[B,A]=butter(3,2/0.00025,'s');
[num2,den2]=bilinear(B,A,4000);
[h2,w]=freqz(num2,den2);
f=w/pi*2000;
plot(f,abs(h1),'-.',f,abs(h2),'-');
grid;
xlabel('频率/Hz ')
ylabel('幅值/dB')
程序中第一个butter的边界频率2π×1000，为脉冲响应不变法原型低通滤波器的边界频率；第二个butter的边界频率2/T=2/0.00025，为双线性变换法原型低通滤波器的边界频率.图3.1给出了这两种设计方法所得到的频响，
虚线为脉冲响应不变法的结果；实线为双线性变换法的结果。

脉冲响应不变法由于混叠效应，使得过渡带和阻带的衰减特性变差，并且不存在传输零点。

同时，也看到双线性变换法，在z=-1即ω=π或f=2000Hz处有一个三阶传输零点，这个三阶零点正是模拟滤波器在Ω=∞处的三阶传输零点通过映射形成的。

2、设计一数字高通滤波器，它的通带为400～500Hz，通带内容许有0.5dB 的波动，阻带内衰减在小于317Hz的频带内至少为19dB，采样频率为1,000Hz。

wc=2*1000*tan(2*pi*400/(2*1000));
wt=2*1000*tan(2*pi*317/(2*1000));
[N,wn]=cheb1ord(wc,wt,0.5,19,'s');
[B,A]=cheby1(N,0.5,wn,'high','s');
[num,den]=bilinear(B,A,1000);
[h,w]=freqz(num,den);
f=w/pi*500;
plot(f,20*log10(abs(h)));
axis([0,500,-80,10]);
grid;
xlabel('')
ylabel('幅度/dB')
图3.2给出了MATLAB计算的结果，可以看到模拟滤波器在Ω=∞处的三阶零点通过高通变换后出现在ω=0（z=1）处，这正是高通滤波器所希望得到的。

3、设计一巴特沃兹带通滤波器，其３dB边界频率分别为f2=110kHz和
f1=90kHz，在阻带f3 = 120kHz处的最小衰减大于１０dB，采样频率f s=400kHz。

w1=2*400*tan(2*pi*90/(2*400));
w2=2*400*tan(2*pi*110/(2*400));
wr=2*400*tan(2*pi*120/(2*400));
[N,wn]=buttord([w1 w2],[0 wr],3,10,'s');
[B,A]=butter(N,wn,'s');
[num,den]=bilinear(B,A,400);
[h,w]=freqz(num,den);
f=w/pi*200;
plot(f,20*log10(abs(h)));
axis([40,160,-30,10]);
grid;
xlabel('频率/kHz')
ylabel('幅度/dB')
图3.3给出了MATLAB计算的结果，可以看出数字滤波器将无穷远点的二阶零点映射为z=±1的二阶零点，数字带通滤波器的极点数是模拟低通滤波器的极点数的两倍。

4、一数字滤波器采样频率f s= 1kHz，要求滤除100Hz的干扰，其３dB的边界频率为95Hz和105Hz，原型归一化低通滤波器为
w1=95/500;
w2=105/500;
[B,A]=butter(1,[w1, w2],'stop');
[h,w]=freqz(B,A);
f=w/pi*500;
plot(f,20*log10(abs(h)));
axis([50,150,-30,10]);
grid;
xlabel('频率/Hz')
ylabel('幅度/dB')
图3.4为MATLAB的计算结果
（三）实验报告要求
1、简述实验目的和实验原理，用笔算求出设计的滤波器，并用计算机计算结果验证，总结实验中的主要结论，写出收获和体会。

2、编写教材P193(1t)和P194(11t)的MATLAB计算程序。

有源模拟带通滤波器的设计
时间：2009-08-21 10:51:10 来源：电子科技作者：张亚黄克平
滤波器是一种具有频率选择功能的电路，它能使有用的频率信号通过。

而同时抑制(或衰减)不需要传送频率范围内的信号。

实际工程上常用它来进行信号处理、数据传送和抑制干扰等，目前在通讯、声纳、测控、仪器仪表等领域中有着广泛的应用。

1 滤波器的结构及分类
以往这种滤波电路主要采用无源元件R、L和C组成，60年代以来，集成运放获得迅速发展，由它和R、C组成的有源滤波电路，具有不用电感、体积小、重量轻等优点。

此外，由于集成运放的开环电压增益和输入阻抗都很高，输出阻抗比较低，构成有源滤波电路后还具有一定的电压放大和缓冲作用。

通常用频率响应来描述滤波器的特性。

对于滤波器的幅频响应，常把能够通过信号的频率范围定义为通带，而把受阻或衰减信号的频率范围称为阻带，通带和阻带的界限频率叫做截止频率。

滤波器在通带内应具有零衰减的幅频响应和线性的相位响应，而在阻带内应具有无限大的幅度衰减。

按照通带和阻带的位置分布，滤波器通常分为低通滤波器、高通滤波器、带通
滤波器和带阻滤波器。

文中结合实例，介绍了设计一个工作在低频段的二阶有源模拟带通滤波器应该注意的一些问题。

2 二阶有源模拟带通滤波器的设计
2．1 基本参数的设定
二阶有源模拟带通滤波器电路，如图1所示。

图中R1、C2组成低通网络，R3、C1组成高通网络，A、Ra、Rb组成了同相比例放大电路，三者共同组成了具有放大作用的二阶有源模拟带通滤波器，以下均简称为二阶带通滤波器。

根据图l可导出带通滤波器的传递函数为
令s=jω，代入式(4)，可得带通滤波器的频率响应特性为
波器的通频带宽度为BW0.7=ω0／(2πQ)=f0／Q，显然Q值越高，则通频带越窄。

通频带越窄，说明其对频率的选择性就越好，抑制能力也就越强。

理想的幅频特性应该是宽度为BW0.7的矩形曲线，如图3(a)所示。

在通频带内A(f)是平坦的，而通带外的各种干扰信号却具有无限抑制能力。

各种带通滤波器总是力求趋近理想矩形特性。

然而实际设计出来的带通滤波器的幅频特性曲线，如图3(b)所示。

在工程上，定义增益自A(f0)下降3 dB(即0．707倍)时的上、下限频率之差值为通频带，用BW0.7表示。

要求其值大于有用信号的频谱宽度，保证信号的不失真传输。

综上分析可知：当有源带通滤波器的同相放大倍数变化时，既影响通带增益
A0，又影响Q值(进而影响通频带BW0.7)，而中心角频率ω0与通带增益A0无关。

2．2 实际电路设计效果分析
为了能更好的了解二阶带通滤波器在实际电路中应用的效果，设计了如图4的电路进行实验验证。

图中U1A部分为放大电路，UlB部分为二阶带通滤波器电路。

根据式(2)~式(4)，设计出了中心频率在30 kHz附近，品质因素Q为1．55，频带宽度约为19．35 kHz的二阶带通滤波器，并分别对它进行了一级到四级级联所产的电压及频率数据的记录，将记录结果绘制成电压／V~频率／kHz图，如图5所示。

从图5(a)中可以看出，随着级联次数的增加，A(f0)在逐渐变大，BW0.7也在逐渐变窄，说明其对频率的选择性越来越好，对干扰信号的抑制能力也越来越强。

除了级联能增强带通滤波器对频率的选择能力以外，另外，改变品质因素Q值的大小也能达到此效果。

众所周知，品质因素Q如果小于0，电路就会自激振荡，无法正常工作。

从图2可以看出，Q值越高，则通频带越窄，也就是说滤波器对频率的选择性就越好，对干扰信号的抑制能力也就越强，但并不是Q值越大，电路就越好越稳定。

为此，也做了如下实验，即根据式(2)～式(4)，设计出了品质因素Q分别为1．55、2．99、7．87这3种中心频率(理论值)一样的二阶带通滤波器，并分别绘制出了它们的电压/V～频率/kHz图，如图5(b)所示。

从图5(b)中可以发现，品质因素Q值越大，其A(f0)在逐渐变大，BW0.7也在逐渐变窄，但是随着Q值的增加，其中心频率也在向低频端倾斜，并且低频端上升的坡度较陡，相对于低频端，高频端下降的幅度较缓。

根据前面的分析也不难看出，Q值如果无限的大，会造成电路的自激振荡，无法正常工作。

为了确定这点，也分别测试了Q值为2．99和7．87两种带通滤波器在无信号输入情况下输出端的情况，如图6(a)，图6(b)所示。

从两个示波器的图可以看出，Q值越大，其自激的程度也就越大，当Q值达到一定数值时，自激程度
与输入信号的强度相当或者比输入信号还要强，就会影响整个电路的正常工作。

2．3 数值的选取
值得注意的是，在设计电路时，首先要根据式(3)确定带通滤波器的中心频率，因为二阶带通滤波器中的元器件比较多，相互干系也比较烦琐。

首先确定中心频率对以后的数值计算会有很大的简化。

为了方便，也可以取R1=R3=R，C1=C2=C，Ra=Rb=R’，如果想设计一
个带放大的带通滤波器，可以根据式(2)或者根据有源带通滤波器的同相放大倍数
在确定了其它数值后适当改变Ra和Rb的值得到你想要的放大倍数。

这里建议不要随意大幅度改变Ra和Rb的值，因为根据式(4)可以看出在确定了其他数值后改变Ra和Rb会影响Q值，而Q值的大小直接影响到电路的工作状态是否稳定。

此外，Q值对元器件数值的大小比较敏感，所以在选择元器件时尽量选取精度较高的器件。

3 结束语
虽然由集成运放和R、C组成的有源滤波电路，具有不用电感、体积小、重量轻，集成运放的开环电压增益和输入阻抗均很高，输出阻抗又低，构成有源滤波电路后还具有一定的电压放大和缓冲作用等优点。

但是因其品质因素Q值无法做的很大，也就导致其通频带宽度无法做的很窄，造成了该滤波器对频率的选择性不是很好，对干扰信号的抑制能力也不是很强，所以在选择设计滤波器方案的同时，要注意结合实际情况，在满足实际要求的状态下
合理选用滤波器的设计方案。

MATLAB设计模拟带通滤波器
悬赏分：200 - 解决时间：2006-12-26 13:53
设计一个模拟带通滤波器，其Ap=1dB,As=34dB,通带从200HZ到300HZ，阻带在
100HZ-400HZ，设计满足上述条件的巴特沃斯滤波器切比雪夫1 切比雪夫2 椭圆滤波器。

提问者：tryo521 - 三级最佳答案
这是我考试时候的程序，参数自己改一下就可以了
cheb1
% wp1=0.45*pi;wp2=0.65*pi;ws1=0.3*pi;ws2=0.75*pi;Rp=1;Rs=40
% =============双线型变换法========================================= wp1=0.45*pi; wp2=0.65*pi;
ws1=0.3*pi; ws2=0.75*pi;
Rp=1; Rs=40;
Wp1=tan(wp1/2); Wp2=tan(wp2/2);
Ws1=tan(ws1/2); Ws2=tan(ws2/2);
BW=Wp2-Wp1; W0=Wp1*Wp2; W00=sqrt(W0);
WP=1; WS=WP*(W0^2-Ws1^2)/(Ws1*BW);
[N,Wn]=cheb1ord(WP,WS,Rp,Rs,'s');
[B,A]=cheby1(N,Rp,Wn,'s');
[BT,AT]=lp2bp(B,A,W00,BW);
[num,den]=bilinear(BT,AT,0.5);
[h,omega]=freqz(num,den,64);
subplot(2,2,1);stem(omega/pi,abs(h));
xlabel('\omega/\pi');ylabel('|H(z)|');
subplot(2,2,2);stem(omega/pi,20*log10(abs(h)));
xlabel('\omega/\pi');ylabel('增益.dB');
% =============直接法=================================
wp1=0.45*pi; wp2=0.65*pi;
ws1=0.3*pi; ws2=0.75*pi;
Rp=1; Rs=40;
Wp=[wp1/pi,wp2/pi]; Ws=[ws1/pi,ws2/pi];
[N,Wn]=cheb1ord(Wp,Ws,Rp,Rs);
[B,A]=cheby1(N,Rp,Wn);
[h,omega]=freqz(B,A,64);
subplot(2,2,3);stem(omega/pi,abs(h));
xlabel('\omega/\pi');ylabel('|H(z)|');
subplot(2,2,4);stem(omega/pi,20*log10(abs(h)));
xlabel('\omega/\pi');ylabel('增益.dB');
%cheby2%
% wp1=0.45*pi;wp2=0.65*pi;ws1=0.3*pi;ws2=0.75*pi;Rp=1;Rs=40
% =============双线型变换法========================================= wp1=0.45*pi; wp2=0.65*pi;
ws1=0.3*pi; ws2=0.75*pi;
Rp=1; Rs=40;
Wp1=tan(wp1/2); Wp2=tan(wp2/2);
Ws1=tan(ws1/2); Ws2=tan(ws2/2);
BW=Wp2-Wp1; W0=Wp1*Wp2; W00=sqrt(W0);
WP=1; WS=WP*(W0^2-Ws1^2)/(Ws1*BW);
[N,Wn]=cheb2ord(WP,WS,Rp,Rs,'s');
[B,A]=cheby2(N,Rs,Wn,'s');
[BT,AT]=lp2bp(B,A,W00,BW);
[num,den]=bilinear(BT,AT,0.5);
[h,omega]=freqz(num,den,64);
subplot(2,2,1);stem(omega/pi,abs(h));
xlabel('\omega/\pi');ylabel('|H(z)|');
subplot(2,2,2);stem(omega/pi,20*log10(abs(h)));
axis([0 1 -100 0]);xlabel('\omega/\pi');ylabel('增益.dB');
% =============直接法=================================
wp1=0.45*pi; wp2=0.65*pi;
ws1=0.3*pi; ws2=0.75*pi;
Rp=1; Rs=40;
Wp=[wp1/pi,wp2/pi]; Ws=[ws1/pi,ws2/pi];
[N,Wn]=cheb2ord(Wp,Ws,Rp,Rs);
[B,A]=cheby2(N,Rs,Wn);
[h,omega]=freqz(B,A,64);
subplot(2,2,3);stem(omega/pi,abs(h));
xlabel('\omega/\pi');ylabel('|H(z)|');
subplot(2,2,4);stem(omega/pi,20*log10(abs(h)));
axis([0 1 -100 0]);xlabel('\omega/\pi');ylabel('增益.dB');
%椭圆%
% wp1=0.45*pi;wp2=0.65*pi;ws1=0.3*pi;ws2=0.75*pi;Rp=1;Rs=40
% =============双线型变换法========================================= wp1=0.45*pi; wp2=0.65*pi;
ws1=0.3*pi; ws2=0.75*pi;
Rp=1; Rs=40;
Wp1=tan(wp1/2); Wp2=tan(wp2/2);
Ws1=tan(ws1/2); Ws2=tan(ws2/2);
BW=Wp2-Wp1; W0=Wp1*Wp2; W00=sqrt(W0);
WP=1; WS=WP*(W0^2-Ws1^2)/(Ws1*BW);
[N,Wn]=ellipord(WP,WS,Rp,Rs,'s');
[B,A]=ellip(N,Rp,Rs,Wn,'s');
[BT,AT]=lp2bp(B,A,W00,BW);
[num,den]=bilinear(BT,AT,0.5);
[h,omega]=freqz(num,den,64);
subplot(2,2,1);stem(omega/pi,abs(h));grid;
xlabel('\omega/\pi');ylabel('|H(z)|');
subplot(2,2,2);stem(omega/pi,20*log10(abs(h)));grid;
xlabel('\omega/\pi');ylabel('增益.dB');
% =============直接法================================= wp1=0.45*pi; wp2=0.65*pi;
ws1=0.3*pi; ws2=0.75*pi;
Rp=1; Rs=40;
Wp=[wp1/pi,wp2/pi]; Ws=[ws1/pi,ws2/pi];
[N,Wn]=ellipord(Wp,Ws,Rp,Rs);
[B,A]=ellip(N,Rp,Rs,Wn);
[h,omega]=freqz(B,A,64);
subplot(2,2,3);stem(omega/pi,abs(h));grid;
xlabel('\omega/\pi');ylabel('|H(z)|');
subplot(2,2,4);stem(omega/pi,20*log10(abs(h)));grid;
xlabel('\omega/\pi');ylabel('增益.dB');。