巴特沃斯高通数字滤波器

数字信号处理课程设计

题目巴特沃斯高通数字滤波器

老师陈忠泽老师

学院电气工程学院

班级电子信息工程0 81班

学号20084470110

姓名何依阳

二0一一年五月

目录：

一、IIR数字高通滤波器的设计

1、数字滤波器的概述

2、数字滤波器的设计步骤

3、设计方法

4、IIR巴特沃斯数字高通滤波器的实例计算

二、软件仿真工具及实现环境简介

1、计算机辅助设计方法

2、 MATLAB直接设计IIR巴特沃斯数字高通滤波器

三、滤波器结构对数字滤波器性能指标的影响分析

1、 IIR系统的基本网络结构

（1）

（2）级联型

（3）

四、有限字长运算在网络结构中对数字滤波器的影响

1 、运算量化效应对数字滤波器的影响

2 、参数的字长对数字滤波器性能指标的影响

2.1 、系数量化对数字滤波器的影响

五、运用MATLAB的辅助工具FDATOOL画出系统函数图像

六、设计心得

IIR 数字高通滤波器的设计

一、IIR 数字高通滤波器的设计 1、数字滤波器的概述

所谓数字滤波器，是指输入、输出均为数字信号，通过数值运算处理改变输入信号所含频率成分的相对比例，或者滤除某些频率成分的数字器件或程序。

2、 数字滤波器的设计步骤

设计一个IIR 数字滤波器主要包括下面5个步骤： （1） 确定滤波器要求的规范指标。

（2） 选择合适的滤波器系数的计算（如图一流程图所示）。 （3） 用一个适当的结构来表示滤波器（实现结构）。 （4） 有限字长效应对滤波器性能的影响分析。 （5） 用软件或硬件来实现滤波器。

本次设计的IIR 数字滤波器系数的计算是根据已知的模拟滤波器的特性转换到等价的数字滤波器。本次设计用双线性变换法得到数字滤波器。而且，双线性变换法得到的数字滤波器保留了模拟滤波器的幅度响应特性。

确定数字巴特沃斯

高通滤波器指标

推导出归一化模拟巴特沃斯低通滤波器指

计算出归一化模拟巴特沃斯低通滤波去归一化推导出模拟巴特沃斯高通滤波器

双线性变换推导出数字巴特沃斯高通

图一 流程图

3 、设计方法

频率变换法设计思想：

1 、从归一化模拟低通原型出发，先在模拟域内经频率变换成为所需类型的模拟滤波器；然后进行双线性变换，由 S 域变换到 Z 域，而得到所需类型的数字滤波器。

2 、先进行双线性变换，将模拟低通原型滤波器变换成数字低通滤波器；然 后在 Z 域内经数字频率变换为所需类型的数字滤波器。

4 、IIR 巴特沃斯数字高通滤波器的实例计算 （1）设计要求

设计一个巴特沃斯数字高通滤波器,要求通带截止频率0.32p ωπ

= rad ,通带衰

减不大于1dB,阻带截止频率

0.16p ωπ

=rad ，阻带衰减不小于15dB

（2）确定数字高通的技术指标：

通带截止频率

0.32p ωπ

= rad ，

通带最小衰减 1

p a = dB

阻带截止频率

0.16p ωπ

= rad ，

阻带最大衰减 15

s a = dB

（3）将巴特沃斯高通数字滤波器的技术指标转换成巴特沃斯高通模拟滤波器的设计指标：令T=2s ，预畸变校正得到的模拟边界频率：

1

tan 0.5494/,12ph p p rad s dB ωαΩ===； 1

tan

0.2566/,152sh s s rad s dB ωαΩ===。

（4）模拟低通滤波器的技术指标计算如下：对通带边界频率归一化，由于本设计的为1dB 截止频率，所以

1,1p c p dB

λλα===

把

1

p λ=和-

sh

Ω带入

p ph

s sh

λλ-Ω=

Ω求得归一化巴特沃斯低通滤波器的

阻带截止频率为：

2.141

p ph

s sh

λλ-Ω=

=Ω ，

15s α= dB

（5）设计归一化巴特沃斯模拟滤波器G （p ）。

0.10.1101

30.6228

10.8757

101

0.2589S P

sp

k αα-===-：

2.141s

sp p

λλλ=

=：

lg 3.1349

lg sp sp

k N λ=

=：

所以取N=4，根据巴特沃斯归一化低通滤波器参数表(见附录)可得归一化模拟低通原型系统函数G(p)为：

4321

() 2.613 3.4142 2.61311G p p p p p =

++++

2

2

1

(0.76541)( 1.84781)p p p p =++++ （6）利用频率变换公式

p p ph

s λΩ=

将G(p)转换成模拟高通

()

HP H s ：

p ()(p |

p ph

a s

H s G λΩ=

=）

4

4

32234

S +2.613S+3.4142S +2.6131S +S ph ph ph ph =ΩΩΩΩ

把

0.5494

ph Ω=代入此式可得：

4

234

0.09110.4333 1.0305 1.4356s s s s s =++++

用双线性变换法将模拟高通

()

a H s 转换成数字高通()H z ：

令

z

z T S 1

1112--+

-=

,

111234

1234

1210.3647 1.4587 2.1881 1.45870.3647()()|1 2.0578 1.85450.78950.1331a z s z

z z z z H z H s z z z z -----------=+-+-+==-+-+

二 、 软件仿真工具及实现环境简介 1、 计算机辅助设计方法

在优秀科技应用软件MATLAB 的信号处理工具箱中提供了一整套模拟，数字滤波器的设计命令和运算函数，方便准确，简单容行使得设计人员除了可按上述传统设计步骤快速的进行较复杂高阶选频滤波器的计算、分析外，还可通过原型变换直接进行各种典型数字滤波器设计，即应用MATLAB 设计工具从模拟原型直接变换成满足原定频域指标要求有数字滤波器。

2、 MATLAB 直接设计IIR 巴特沃斯数字高通滤波器

MATLAB 编程如下： fs=5000; wp=800*2/fs; ws=400*2/fs; rp=1; rs=15; Nn=128;

[N,wn]=buttord(wp,ws,rp,rs) [b,a]=butter(N,wn,'high') freqz(b,a,Nn,fs)

MATLAB 运算结果如下：

N = 4

wn = 0.2388

b = 0.3647 -1.4587 2.1881 -1.4587 0.3647

a = 1.0000 -2.0578 1.8545 -0.7895 0.1331

三、滤波器结构对数字滤波器性能指标的影响分析

1 、 IIR系统的基本网络结构

IIR系统的基本网络结构有三种，即直接型、级联型和并联型。（1）

N阶差分方程如下：

∑∑==-

+

-

=

M

i

N

i

i

i

i

n

y

a

i

n

x

b

n

y

01

) (

)

(

)

(

对应的系统函数为

∑∑=-=--=

N

i i

i M

i i i

z a z b z H 1

01)(

设M=N=2，按照差分方程可以直接画出网络结构如下图（a ）所示。图中第一部分系统函数用)(1z H 表示，第二部分用)(2z H 表示，那么)()()(21z H z H z H =，当然也可以写成)()()(12z H z H z H =,按照该式，相当于将下图(a)中两部分流图交换位置，如下图（b ）所示。该图中节点变量w1=w2，因此前后两部分的延时支路可以合并，形成如下图 (c)所示的网络结构流图，我们将下图 (c)所示的这类流图称为IIR 直接型网络结构。 M=N=2时的系统函数为

2

2112

21101)(------++=

z a z a z b z b b z H

对照下图 (c)的各支路的增益系数与)(z H 分母分子多项式的系数可见，可以直

接按照)(z H 画出直接型结构流图。

IIR 网络直接型结构

由bz 和az 写出数字滤波器系统函数：

1234

1234

0.3647 1.4587 2.1881 1.45870.3647()1 2.0578 1.85450.78950.1331z z z z H z z z z z ---------+-+=

-+-+

有H （z ）写出差分方程如下：

() 2.0578(1) 1.8545(2)0.7895(3)0.1331(4)y n y n y n y n y n =---+---

0.3647() 1.4587(1) 2.1881(2) 1.4587(3)0.3647(4)x n x n x n x n x n +--+---+-

直接型网络结构如下图：

直接型网络结构图

在后面的分析中我们将发现，直接型系统对滤波器的性能控制作用不明显，极点对系数的变化不灵敏，易出现不稳定或较大误差，而且运算的累积误差较大。因此，在设计时一般不选用直接型。

（2）级联型

级联型结构是将系统传递函数H(z)写成具有实系数的二阶节的乘积。将分子和分母多项式分解为各自的根，然后将一对复数共轭根（或者任意两个实数根）组合成二阶多项式。

在直接型表示的系统函数H(z)中，分子、分母均为多项式，且多项式的系

z 1

z -x(n)

y(n)

0.3642.0578

-1.8545

0.7895

-0.1331

-0.1458

2.1881

1.4587

0.3647

z 1

z -z 1

z -z 1

z

-

数一般为实数。现将分子、分母多项式分别进行因式分解，得到：

∏∏=-=---=N

r r

M

r r

z d

z C A

z H 1111)

1()1()(

上式中, A 是常数; Cr 和dr 分别表示H(z)的零点和极点。由于多项式的系数是实数，Cr 和dr 是实数或者是共轭成对的复数，将共轭成对的零点（极点）放在一起，形成一个二阶多项式，其系数仍为实数；再将分子、分母均为实系数的二阶多项式放在一起，形成一个二阶网络

)

(z H j 。

)

(z H j 如下式：

2

j 21j 12

j 211j 0j j 1)(------++=

z a z a z z z H βββ

上式中

)

(z H j 表示一个一阶或二阶的数字网络的子系统函数，每个

)

(z H j 的网络

结构均采用前面介绍的直接型网络结构，如下图所示，H(z)则由k 个子系统级联构成。

一阶和二阶直接型网络结构

级联型的MATLAB 的表示与实现：

将数字滤波器系统函数H （z ）：

1234

1234

0.3647 1.4587 2.1881 1.45870.3647()1 2.0578 1.85450.78950.1331z z z z H z z z z z

---------+-+=-+-+

用MATLAB 转换成级联型的程序如下：

bz=[0.3647 -1.4587 2.1881 -1.4587 0.3647]; az=[1.0000 -2.0578 1.8545 -0.7895 0.1331]; [S,G]=tf2sos(bz,az) freqz(bz,az)

说明：[S,G]=tf2sos(bz,az):实现直接型到级联型的变换。B 和A 分别为直接型系统函数的分子和分母多项式系数向量。返回L 级二阶级联型结构的系数矩阵S 和增益常数G 。

MATLAB 运算结果如下：

S = 1.0000 -2.1807 1.1988 1.0000 -0.8976 0.2272 1.0000 -1.8190 0.8341 1.0000 -1.1602 0.5859 G = 0.3647

由S 和G 写出数字滤波器级联型系统函数：

1212

12121 2.1807 1.19881 1.8190.8341()0.3647()()10.89760.22721 1.16020.5859z z z z H z z z z z

---------+-+=-+-+

级联型网络结构如下图：

级联型网络结构图

在级联型结构中，每一级分子的系数确定一对零点，分母的系数确定一对极点，因为子网络的零极点也即整体网络的零极点，所以整个系统的零极点都可以准确的由每一级的系数来调整和控制，这样便于调整滤波器的频率响应性能，其灵敏度特性优于直接型和正准型结构。其次，级联结构具有最少的存储器。并联

-0.2961

x(n)

y(n)

-0.8976

-0.2272

-2.1807

1.1988

z 1

z -

z 1

z

-

-1.0003

-0.12378

0.06189

z 1

z -

z 1

z -

支路的极点也是整个网络的极点，而并联支路的零点却不是整个网络的零点，因此并联网络能独立的调整系统的极点位置，但不能控制零点。并联结构的灵敏度由于直接型和正准型，运算累积误差比级联型小。

（3）

如果将级联形式的H(z)展成部分分式形式，则得到：

)

()()()(21z H z H z H z H k +++=

对应的网络结构为这k 个子系统并联。上式中，Hi(z)通常为一阶网络或二阶网

络，网络系统均为实数。二阶网络的系统函数一般为

2

2111

101)(-----+=

z z z z H i i i i i ααββ

式中，i 0β、i 1β、i 1α和i 2αi 都是实数。如果i 1β=i 2α=0，则构成一阶网络。由上

式，其输出Y(z)表示为

)

()()()()()()(21z X z H z X z H z X z H z Y k +++=

上式表明将x(n)送入每个二阶（包括一阶）网络后，将所有输出加起来得到输

出y(n)。

在并联型结构中，每一个一阶网络决定一个实数极点，每一个二阶网络决定一对共轭极点，因此调整极点位置方便，但调整零点位置不如级联型方便。另外，各个基本网络是并联的，产生的运算误差互不影响，不像直接型和级联型那样有误差积累，因此，并联形式运算误差最小。由于基本网络并联，可同时对输入信号进行运算，因此并联型结构与直接型和级联型比较，其运算速度最高。

四、有限字长运算在网络结构中对数字滤波器的影响 1 、 运算量化效应对数字滤波器的影响

在实现数字滤波器时，将遇到相乘与求和运算。在定点制运算中，每一次乘法运算之后都要作一次舍入（或截留）处理，研究定点实现相乘运算的流图如下图所示。图（a ）表示无限精度乘积y(n)；图（b ）表示有限精度乘积)(n y

∧

,]

[?表示

舍入运算。采用统计分析方法时，可以将舍入误差作为独立噪声e(n)叠加在信号上，如图(c)所示。

)

(n x )

(n y a

)(n x )(n y a

[]

?)(?n y

)(n x a

)(?n y

)

(n e

a 理想相乘

b 实际相乘的非线性流图

c 统计模型的线性流图 定点相乘运算的流图表示

显然，采用统计分析方法后，实际的输出可以表示为： )()()(?n e n y n y

+＝

对于舍入处理，e(n)的均值为零，方差为：

12

2

2q e

=σ

现在以一个一阶IIR 巴特沃斯数字高通滤波器为例来讨论分析方法。表示其输入与输出关系的差分方程为

)()1()(n x n ay n y +-=， n ≥0

式中|a|<1。它含有乘积项)

1(-?n y a ，这将引入一个舍入噪声，其统计分析流图示于下图。

一阶IIR 滤波器的幅频特性

整个系统可以当作线性系统来处理。输出噪声

)

(n e f 是由噪声源e(n)造成的输出

误差，可由量化误差通过线性系统的方法求得输出噪声。由于e(n)叠加在输入端，因此

)

()()()()(n u a n e n h n e n e n f *=*=

式中h(n)=anu(n)是一阶系统的单位脉冲响应,由式

∑∑∑∞=∞

=

∞

==-=00

2

222)

()()()(m l m e

e

f

m h

l m l h m h σ

δδσ （

2f

σ为输出噪声方差）和

z

dz z H z H j c e

f

)()(2112

2-??=πσσ

可求得输出噪声的方差

∑

∞

==0

222m m e

f a σ

σ

式中H(z)为一阶IIR 系统的传递函数，即

a z z

z H -=

)(

有此可以求得

)1(32)1(1211222222

2a a q a b

e f -=-=-=-σσ，（b

q -=2,）

由此可见，字长b 越大，数字滤波器输出端的噪声越小。

2 、参数的字长对数字滤波器性能指标的影响 2.1 、系数量化对数字滤波器的影响

系统对输入信号进行处理时需要若干参数或者称为系数，这些系数都要存储在有限位数的寄存器中，因此存在系数的量化效应。系数的量化效应误差直接影响系统函数的零、极点位置，如果发生了偏移，会使系统的频率响应偏移理论设计的频率响应，不满足实际需要。

下面分析系数量化误差对极零点位置的影响。如果极零点位置改变了，严重时不仅IIR 系统的频率响应会发生变化，还会影响系统的稳定性。因此研究极点位置的改变更加重要。为了表示系数量化对极点位置的影响，引入极点位置灵敏度的概念，所谓极点灵敏度, 是指每个极点对系数偏差的敏感程度。相应的还有

零点位置灵敏度，分析方法相同。

现分析一个N 阶直接型结构的IIR 滤波器的传递函数

)

()(1)(1

z A z B z a z b z H N

k k

k M

k k k

=

-=

∑∑=-=-

上式表示了一个N 阶直接型结构的IIR 数字滤波器的系统函数，该滤波器的极点

都在单位圆内聚集在z=1附近。系数ak 和bk 必须用有限位二进制数进行量化，存储在有限长的寄存器中，经过量化后的系数误差为△ar 和△br ，量化后的系数用r a ∧

和r b ∧

，即

r r r a a a

?+=? r r r b b b ?+=?

则实际的系统函数可表示为:

∑∑=-=--=N

k k k M

k k

k z a

z b z H

10?1?)(?

显然，系数量化后的频率响应已不同于原来设计的频率响应。用直接型结构来实现该滤波器时，系数ak 和bk 都将直接出现在信号流程图中，其中ak 影响着极点的位置。当由于系数量化误差使一个极点从单位圆内移动到单位圆上或单位圆外时，滤波器的稳定性即受到破坏。所以，只要有一个系数由于量化产生很微小的误差，就有可能使系统失去稳定。反馈支路的阶次N 越高，使滤波器失去稳定的系数量化误差的绝对值就越小，则越容易使滤波器变得不稳定。

设滤波器的传输函数)(z H ，系数ak 和bk 经舍入量化后为r r r a a a ?+=?和

r r r b b b ?+=?，这里△ar 和△br 是量化误差。

分母多项式有N 个极点，用i

p （i=1,2,… N ）表示。这样，实际的滤波器

的传输函数为:

∏∑∑∑=-=-=-=-?+-=

-=

N

i M

k k k

N

k k k M

k k

k

z pi pi z b z a

z

b z H 11

1

0])(1[??1?)(

上式中，i p ?是第i 个极点的偏差，称为极点误差，它应该和各个系数偏差都有关，它和各系数偏差的关系用下式表示:

N

i a a p p k N

k k

i

i ,,1,

1 =???=?∑=

上式中，k i a p ??的大小决定着系数k a 的偏差k a ?对极点偏差i p ?的影响程度。k i

a p ??越大，k a ?对i p ?影响也越大。称k i

a p ??为极点i p 对系数变化的灵敏度。

五 、运用MATLAB 的辅助工具FDATOOL 画出系统函数图像

系统函数系数向量经过1位二进制舍入量化前后幅频特性和相位响应

系统函数系数向量经过4位二进制舍入量化前后幅频特性和相位响应

系统函数系数向量经过1位二进制舍入量化前后零极点特性图

系统函数系数向量经过4位二进制舍入量化前后零极点特性图

由MATLAB运算1位和4位二进制舍入量化后结果可以看出，因为系数的量化，使极点位置发生变化，系数量化的相对误差（p-pq）不到10%，极点位置的相对误差ap达到了70%。问题不但是数量的变化，算一下极点的模，可以发现所有根的模都变大了，说明量化后的极点离单位圆稍远一些，如系统函数系数向量经过1和4位二进制舍入量化前后的零极点特性图如图（c）、（d）所示。这致使数字高通滤波器的幅频特性降低，运用MATLAB的辅助工具FDATOOL画出1位和4位量化墙后的幅频特性曲线分别如图（a）、（b）中的实线和虚线所示，这说明由于系数量化效应，使极点位置发生了变化，从而改变了原来设计的频率响应特性。另外，我们还注意到4位二进制舍入量化后极点都在单位圆内部，但在系统函数系数向量经过1位二进制舍入量化后，原来较小的系数相对误差变化较大，使滤波器性能偏离原设计指标要求，使本来稳定的系统变成了不稳定滤波器。

从以上分析可以看出，系数的量化效应误差直接影响系统函数的零、极点位置，如果发生了偏移，会使系统的频率响应偏移理论设计的频率响应，不满足实际需要。因此，在设计滤波器的时候应尽量选择合理的系数量化，以减小量化误差，以便设计出较为理想的实用滤波器。

六、设计心得：

在课设之前，我对 MATLAB 软件，特别是滤波器设计中的函数基本上处于一种模糊状态。但是通过与同学不断的交流，最后完成了这次课设，对滤波器的设计有了比较清楚的了解。

这次课设的完成，真的不太容易。我的理论基础不是很扎实，所以完成起来要比其他同学费劲。在课程设计的过程中，我学到了很多东西，比如设计滤波器的一些基本函数的用法，各种模拟滤波器的特性，设计滤波器的一些基本方法。但更为重要的是，我对于解决一个问题的思路更加清晰，找到了属于自己的方法。当然，在设计的过程中，不可能避免的遇到了很多问题，如刚开始思路比较混乱，没有明确的方向。也由于一些基本的概念了解的不够清楚，比如模拟频带变换，我理解为是做模拟滤波器，而老师的课程设计题上有一句话似乎要求是要做数字滤波器。所以，在这个问题上，走了很多弯路。但通过与同学交流终于弄清楚数字滤波器而不是模拟滤波器。模拟频带变换只不过是一种变化方式或者说一种变化过程。是由模拟低通变为模拟高通，再由模拟高通变为数字高通的过程。

总的来说，这次课程设计让我对 MATLAB 有了更深刻的了解，对数字滤波器的设计流程有了大致的了解，掌握了一些设计滤波器的基本方法，提高了理论用于实践的能力，掌握了更多专业相关的使用知识与技能。同时，也暴露了我很多的不足，在以后的学习中，将进一步发现并克服缺点。