实验四 IIR数字滤波器的设计(1) (2)

实验四 IIR 数字滤波器的设计及网络结构

一、实验目的

1．了解IIR 数字滤波器的网络结构。

2.掌握模拟滤波器、IIR 数字滤波器的设计原理和步骤。 3．学习编写数字滤波器的设计程序的方法。

二、实验内容

数字滤波器：是数字信号处理技术的重要内容。它的主要功能是对数字信号进行处理，保留数字信号中的有用成分，去除信号中的无用成分。

1．数字滤波器的分类

滤波器的种类很多，分类方法也不同。

（1）按处理的信号划分：模拟滤波器、数字滤波器 （2）按频域特性划分；低通、高通、带通、带阻。 （3）按时域特性划分：FIR 、IIR

2．IIR 数字滤波器的传递函数及特点

数字滤波器是具有一定传输特性的数字信号处理装置。它的输入和输出均为离散的数字信号，借助数字器件或一定的数值计算方法，对输入信号进行处理，改变输入信号的波形或频谱，达到保留信号中有用成分去除无用成分的目的。如果加上A/D 、D/A 转换，则可以用于处理模拟信号。

设IIR 滤波器的输入序列为x(n)，则IIR 滤波器的输入序列x(n)与输出序列y(n)之间的关系可以用下面的方程式表示：

1

()()()

M N

i j i j y n b x n i a y n j ===-+-∑∑

（5-1）

其中，j a 和i b 是滤波器的系数，其中j a 中至少有一个非零。与之相对应的差分方程为：

10111....()()()1....M

M N

N b b z b z Y z H Z X z a z a z ----++==

++ （5-2）

由传递函数可以发现无限长单位冲激响应滤波器有如下特点： （1） 单位冲激响应h(n)是无限长的。

（2） 系统传递函数H(z)在有限z 平面上有极点存在。

（3） 结构上存在着输出到输入的反馈，也就是结构上是递归型的。

3．IIR 滤波器的结构

IIR 滤波器包括直接型、级联型和并联型三种结构：

① 直接型：优点是简单、直观。但由于系数bm 、a k 与零、极点对应关系不明显，一

个bm 或a k 的改变会影响H(z)所有零点或极点的分布，所以一方面，bm 、a k 对滤波器性能的控制关系不直接，调整困难；另一方面，零、极点分布对系数变化的灵敏度高，对有限字长效应敏感，易引起不稳定现象和较大误差。 Matlab 实现：

filter( )函数实现IIR 数字滤波器直接形式。格式为： y=filter(b,a,x)

b ，a 为差分方程输入、输出系数向量（或系统函数的分子、分母多项式，降幂），x 为输入序列，y 为输出序列。 其中，传递函数（tf ）形式

N

N M

M z

a z a z

b z b b z H ----++++++= 111101)(若

则a=[1 a 1 a 2… a N ] b=[b 0 b 1 b 2… b M ]

② 级联型：基于因式分解，将系统函数H(z)分解为因子乘积的形式。

(5-3）

级联型结构：

Matlab 实现： tf2zp( )函数用于求系统函数的零、极点和增益常数，zp2sos ( )函数则根据tf2zp( )函数结果求出各基本节系数。格式为：

[z,p,K]=tf2zp(b,a)； sos=zp2sos(z,p,K)；

b ，a 为差分方程输入、输出系数向量（系统函数的分子、分母多项式，降幂）。 其中，零极点增益形式（zp ）：

∏∏-=--=---=10

1

10

1

)

1()

1()(N k k

M i i

z

z z

z k

z H 若

则 零点向量 z=[z 1 z 2 … z M-1];

极点向量 p=[z 1,z 2,…,z N-1]

12

1212

111211()()11M

m

N N m

m k k k N

k k k k k k k b

z z z H z K K H z z z a z

ββαα---=---===++===+++∑∏∏∑

k 为系统增益。

二阶分式形式（sos ）为： 把H(z)划成二阶因式

∏∏==----++++==N

k N

k k k k k k k k z z

z z z H z H 11221

102

2110)()(αααβββ 则其二阶因式为：

??????

?

?????

??=N N N N N b b b b b a a b b b a a b b b sos 2121022

1222

120221

11211101111

③ 并联型：基于部分分式展开，将系统函数H(z)分解为部分分式和的形式。

（5-4）

并联型结构：

Matlab 实现：residue( )函数可以实现并联型结构，有两种格式：

[K,r,p]=residue(b,a); [b,a]=residue(b,a);

其中，部分分式形式：

)

(1)1()1()(1)()1(1)1()(N M n

z N M k k z

n p n r z p r z H ----+-+++-++-=

若 则极点向量 p=[p(1) p(2) … p(n)]

其对应系数向量 r=[r(1) r(2) … r(n)]

余数多项式系数向量 k=[k(1) k(2) … k(M-N+1)] 【实例5-1】已知三阶IIR 数字滤波器的系统函数

001001001211

121

()()11M m N N m m k k

k N k k k k k k k b z z H z K K H z z z a z γγαα--=---===+==+=++++∑

∑∑

∑

3212

16

36132353)(------++++=

z z z z z z H 求：①直接形式的单位采样响应h(n)；

②级联型结构的各基本节系数； ③并联型结构的部分分式系数。 解：MATLAB 源程序为 ①

b=[3,5/3,2/3];

a=[1,1/6,1/3,-1/6]; x=[1,zeros(1,50)]; y=filter(b,a,x); n=0:50;plot(n,y); ②

b=[3,5/3,2/3,0];

a=[1,1/6,1/3,-1/6]; [z,p,K]=tf2zp(b,a); sos=zp2sos(z,p,K);

③

b=[3,5/3,2/3];

a=[1,1/6,1/3,-1/6]; [K,r,p]=residue(b,a); KK1=[K(1),K(2)]; zz1=[z(1),z(2)];

[b2,a2]=residue(KK1,zz1,0);

5．IIR 数字滤波器的具体设计

（1）巴氏模拟原型滤波器的设计

巴氏模拟低通滤波器幅度平方函数为

N

C

j H 211)(???

? ??ΩΩ+=

Ω

MATLAB 工具箱函数buttap,buttord 和butter 是巴氏滤波器设计函数。 1）[Z,P,K]=buttap(N)

该格式用于计算N 阶巴氏归一化（3dB 截止频率Ωc=1）模拟低通原型滤波器系统函数的零、极点和增益。得到的系统函数为：

)

())(()

())(()(2121N N a p p p p p p z p z p z p K

p G ------=

如果要从计算得到的零、极点得到系统函数的分子和分母向量B 和A ，可以调用结

构转换函数[B,A]=zp2tf(Z,P,K)。

2）[N,wc]=buttord(wp,ws,Rp,Rs,’s ’)

该格式用于计算巴氏滤波器的阶数N 和3dB 截止频率wc 。其中： N ——滤波器的阶数。

Wc ——3dB 截止频率的归一化值。

wp 、ws ——通带、阻带边界频率的归一化值。要求：1ws 0,1wp 0≤≤≤≤，1表示数字频率π。

Rp 、Rs ——通带最大、阻带最小衰减。

s ——可选项，直接设计模拟滤波器，此时wp 和ws 均为实际模拟角频率。 3）[B,A]=butter(N,wc,’ftype ’，’s ’)

计算N 阶巴氏数字滤波器系统函数分子和分母多项式的系数向量B 和A 。

【实例5-2】已知通带截止频率kHz f p 5=，通带最大衰减dB a p 5=，阻带截止频率

kHz f s 12=，阻带最小衰减dB a s 30=，设计巴特沃斯型模拟低通滤波器。

解：MATLAB 源程序为

Wp=2*pi*5000;Ws=2*pi*12000;Rp=5;As=30; [N,Wc]=buttord(Wp,Ws,Rp,Rs,'s'); [B,A]=butter(N,Wc,'s'); freqs(B,A);

【实例5-3】设计阶数为3，5，10，15的巴氏模拟原型滤波器。并画出幅频响应曲线。 解：MATLAB 源程序为 for i=1:4;

switch i

case 1 N=3; case 2 N=5; case 3 N=10; case 4; N=15; end;

[z,p,k]=buttap(N); [b,a]=zp2tf(z,p,k); [h,w]=freqs(b,a,n); Ah=abs(h);

subplot(2,2,i),plot(w,Ah);axis([0 2 0 1]); xlabel('w/wc');ylabel('|H(jw)|.^2');

title(['filer N=',num2str(N)]);grid;end;

【实例5-4】设通带、阻带截止频率fp=0.5kHz 、fs=1.2kHz ，通带、阻带最大衰减Rp=1dB,Rs=30dB ，要求设计巴氏低通滤波器。

解：MATLAB源程序为:

>> OmegaP=2*pi*500;OmegaS=2*pi*1200;

>> Rp=1;Rs=30;

>> [N,OmegaC]=buttord(OmegaP,OmegaS,Rp,Rs,'s'); %确定阶数N

>> [z0,p0,k0]=buttap(N);%确定传递函数

>> z=z0*OmegaC;%去归一化

>> k=k0*OmegaC^N;

>> p=p0*OmegaC;

>> bs=k*real(poly(z));

>> as=real(poly(p));

>> freqs(bs,as);

（2）IIR数字滤波器的设计

模拟滤波器Ha(s) 转换成数字滤波器H(z)应满足要求： (1) 因果稳定的模拟滤波器转换成数字滤波器，仍是因果稳定的；(2)数字滤波器的频率响应模仿模拟滤波器的频响，s 平面的虚轴映射z平面的单位圆，相应的频率之间成线性关系。脉冲响应不变法和双线性变换法都满足如上要求。

①脉冲响应不变法

用数字滤波器的单位脉冲响应序列h(n)模仿模拟滤波器的冲激响应h a(t),让h(n)正好等于h a(t)的采样值，即h(n)=h a(nT)，其中T为采样间隔。

②双线性变换法

s平面与z平面之间满足以下映射关系：

1

111--+-=z z s

s 平面的虚轴单值地映射于z 平面的单位圆上，s 平面的左半平面完全映射到z 平面的单位圆内。双线性变换不存在混叠问题。

双线性变换时一种非线性变换)2/(ωtg =Ω，这种非线性引起的幅频特性畸变可通过预畸而得到校正。

以低通数字滤波器为例，将设计步骤归纳如下：

(1)确定数字滤波器的性能指标：通带临界频率f p 、阻带临界频率f s ；通带内的最大衰减R p ；阻带内的最小衰减R s ；

(2)确定相应的数字角频率，ωp=2πfp ；ωs=2πfs ； (3)计算经过预畸的相应模拟低通原型的频率，)2/(ωtg =Ω；

(4)根据Ωp 和Ωs 计算模拟低通原型滤波器的阶数N ，并求得低通原型的传递函数

H a (s)； (5)用上面的双线性变换公式代入H a (s)，求出所设计的传递函数H(z)； (6)分析滤波器特性，检查其指标是否满足要求。

1)冲激响应不变法

MATLAB 实现：采用impinvar 函数 [bz,az]=impinvar(b,a,fs,tol)

bz,az ――转换的分子分母向量 fs ――采样频率，默认值为1Hz

tol ――误差容限，表示转换后的离散系统函数是否有重复的极点。 2）双线性变换法

MATLAB 实现采用bilinear 函数： [bz, az]=bilinear(bs,as,fs)

――把模拟滤波器的传递函数转换为数字的 其中：

bs,as ——模拟滤波器的分子、分母向量

bz,az ——数字滤波器的分子、分母向量

【实例5-5】基于Butterworth 模拟滤波器原型，使用双线性转换设计数字滤波器，其中参数指标为：

通带截止频率：0.2p ωπ= 通带波动值：dB R p 1= 阻带截止频率：0.3s ωπ= 阻带波动值：dB A s 15=

解：首先确定滤波器的阶数N ，同时根据s p ΩΩ和确定c Ω＝0.5。

接着使用bilinear 进行双线性转换，最后绘制在频域上的各种图像，其源代码如

下：

％数字滤波器指标 wp=0.2*pi; ws=0.3*pi;

Rp=1;

As=15;

%将数字滤波器指标反转变化为模拟滤波器的参数

T=1;fs=1/T;

omegap=(2/T)*tan(wp/T);

omegas=(2/T)*tan(ws/T);

%butterworth原型模拟滤波器的设计

[N,Wc]=buttord(omegap,omegas,Rp,As,'s');

[B,A]=butter(N,Wc,'s');

%双线性变换

[b,a]= bilinear(B,A,fs);

%频域图像的绘制

freqz(b,a);

程序运行后，产生4阶的butterworth数字滤波器，频率响应如图5-5所示的波形。

图5-5 Butterworth数字滤波器

三、思考题

1、用脉冲响应不变法和双线性变换法设计一个数字低通滤波器，满足：wp=0.2π,Rp=1dB,ws=0.3π,Rs=15dB,采样频率为1000Hz.

四、实验报告要求

1．实验目的和实验原理

2．总结实验中的主要结论，写出收获和体会。