毕业论文 基于Matlab的数字滤波器的设计

基于MATLAB的数字滤波器的设计

陶雯（安庆师范学院物理与电气工程学院安徽安庆 246011）

指导老师:张杰

摘要：数字滤波器是数字信号处理的重要环节，数字滤波器可分为IIR和FIR两大类。本文介绍了IIR和FIR的基本设计原理以及在MATLAB环境下如何利用直接程序设计法、SPTOOL设计法和FDATOOL设计法给出IIR和FIR数字滤波器的设计方法和操作步骤，并给出设计设计实例及运行结果，同时利用MATLAB 环境下的仿真软件SIMULINK对所设计的滤波器进行模拟仿真，仿真结果表示设计参数设置合理。

关键词：数字滤波器，仿真，MATLAB，FIR，IIR。

1 引言

MATLAB是“矩阵实验室”（MATrix LABoratoy）的缩写，它是由美国Mathworks公司于1984年正式推出的，是一种以矩阵运算为基础的交互式程序语言，专门针对科学、工程计算及绘图的需求。MATLAB 是功能强大的科学及工程计算软件，它不但具有以矩阵计算为基础的强大数学计算和分析功能，而且还具有丰富的可视化图形表现功能和方便的程序计算能力。MATLAB的应用领域极为广泛，除数学计算和分析外，还被广泛地应用于自动控制、系统仿真、数字信号处理、图形图像分析、数理统计、人工智能、虚拟现实技术、通信工程、金融系统等领域，因此，MATLAB是面向21世纪的计算机程序设计及科学计算语言。

数字滤波器是指完成信号滤波处理功能的，用有限精度算法实现的离散线性非时变系统，其输入是一组（由模拟信号取样和量化的)数字量，其输出是经过变换或说处理的另一组数字量。因此，它本身既可以是用数字硬件装配成德一台完成给定运算的专用数字计算机，也可将所需的运算编成程序，让通用计算机来执行。数字滤波器具有稳定性高、精度高、灵活性大等突出优点。随着数字技术的发展，用数字技术实现滤波器的功能愈来愈受到人们的注意和广泛使用。这里所说的数字滤波器是指理想带通，低通等的频率选择数字滤波器。

数字滤波器设计的一个重要步骤是确定一个可实现的传输函数H(z)，这个确定传输函数H(z)的过程称为数字滤波器设计。数字滤波器的一般设计过程为：（1）按照实际需要，确定滤波器的性能要求（通常在频域内给定数字滤波的性能要求）。（2）寻找一满足预定性能要求的离散时间线性系统。（3）用有限精度的运算实现所设计的系统。（4）通过模拟，验证所设计的系统是否符合给定性能要求。

2 数字滤波器的设计

滤波器分为两种，分别为模拟滤波器和数字滤波器。数字滤波器可以理解为是一个计算程序或算法，将代表输入信号的数字时间序列转化为代表输出信号的数字时间序列，并在转化的过程中，使信号按预定的形式变化。数字滤波器有多种分类，从数字滤波器功能上分可分为低通、高通、带阻、带通滤波器，根据数字滤波器冲激响应的时域特征，可将数字滤波器分为两种，即无限长冲激响应滤波器（IIR）和有限长冲激响应滤波器（FIR）。

数字滤波器指标：一般来说，滤波器的幅频特性是分段常数的，以低通为例，在通带内逼近于1，阻带内逼近与0，实际设计的滤波器并非是锐截止的通带和阻带两个范围，两者之间总有一个过渡带。

在设计滤波器时事先给定幅频特性允许误差，在通带范围内幅度响应以误差1σ逼近于1，在阻带内幅度响应以误差2σ逼近于0。

π

σσ≤≤≤≤≤≤-w w e H w w e H r jw

c jw ,2|)(|,1|)(|11 （1）

式中wc 和wr 分别为通带边界频率和阻带边界频率，wr-wc 为过渡带。在具体的技术指标中往往用通带波动σ来表示1σ，用最小阻带衰减At 来表示2σ，其具体的对应公式这里就不详述了。 2.1 IIR 数字滤波器设计

IIR DF 的冲激响应h(n)是无限长的，其输入输出的关系为：

)()()(i n x i h n y i -=

∑+∞

-∞

= (2)

系统函数为

∑+∞

=-∞

=-=

n n n

z

n h z H )()(=

∑∑=-=--n

k k

k m

r r

r

z a z

b 1

1 (3)

设计无限长单位脉冲响应（IIR ）数字滤波器一般可有三种方法。

第一种方法，先设计一个合适的模拟滤波器，然后将其数字话，即将S 平面映射到Z 平面得到所需的数字滤波器。模拟滤波器的设计技巧非常成熟，不仅得到的是闭合形式的公式，而且设计系数已经表格化了。因此，由模拟滤波器设计数字滤波器的方法准确，简便，得到普遍采用。对于这种方法，工程上有两种常见得变换法——脉冲响应不变法及双线性变换法。

第二种方法，在Z 平面直接设计IIR 数字滤波器，给出闭合形式的公式，或者以所希望的滤波器响应作为依据，直接在Z 平面上通过多次选定极点和零点的位置，以逼近该响应。

第三种方法，利用最优化技术设计参数，选定极点和零点在Z 平面上的合适位置，在某种最优化准则意义上逼近所希望的响应。但一般不能得到滤波器的系数（即零，极点的位置）作为给定响应的闭合形式函数表达式。优化设计需要完成大量的迭代运算，这种设计法实际上也是IIR 滤波器的直接设计。

本文着重介绍由模拟滤波器设计相应的IIR 数字滤波器的方法。 （1）脉冲响应不变法

脉冲响应不变法是使数字滤波器的单位脉冲响应序列h(n)逼近模拟滤波器的冲激响应()a h t ，让h(n)正好等于()a h t 的采样值。设已有满足要求的模拟滤波器, 则可

()a H s →()a h t →()()a h n h nT =→()H z

{因为:()()a h t h n ≈的图形的图形}, 公式导出:

具体转换如下:设(以一阶极点为例)

1()N

k

a k k

A H s s s ==-∑

(4) 作拉氏反变换,得

1

1

()[()]()k N

s t a a k k h t FT H s A e u t -===∑

采样得

1

()()()k N

s nT a k k h n h nT A e u nT ===∑

作Z 变换,得

101

1()()1k k N

N

s nT

n

k

k s T n k k A H z A e

u nT z

e

z ∞--=====-∑∑∑

(5)

()H z 与()a

H s 极点关系为: k s T

k z e = (6)

一般对应关系

sT

j T j T s j j z re

z e re e e ωσΩσΩ

ω=+==→=?

,T r e T σωΩ?== (7)

图1 S 平面到Z 平面变换示意图

所以, 模拟系统稳定因果→数字系统稳定因果。

按照脉冲响应不变法，从S 平面到Z 平面的映射不是单值关系，而是先将()a H s 在S 平面沿虚轴作周期严拓，再按照映射关系将()a H s 映射到Z 平面，得到()H z ，因此，脉冲响应不变法只适用于带限的滤波器（如低通、带通）。

在Matlab 中利用M 文件impinvar 可以对模拟传输函数实行脉冲响应不变法。 （2）双线性变换法

脉冲响应不变法不适带阻和高通滤波器的设计，因为高频带为通带，前述方法易引起混频。故希望:s 平面虚轴?z 平面单位圆一周, 且应有

:0:0Ωωπ→+∞?→, :0:0Ωωπ→-∞?→-,

因为tan(/2)π±=±∞, 所以选变换

tan

2

K ω

Ω= (8)

其中K 可取任意正常数, 后面将导出2/K T =

.

j Ω

σ

O

1

ω

1

r =O s

平面

0σ<0

σ=z 平面

1

r

r e

σ=k

s ?

k

z ?

设计思路: tan 2

,,,,,,p s p s p s p s K ΩωωωααΩΩαα=========>.→设计出模拟滤波器→转化成数字滤波器.

图2 数字域频率与模拟域频率的对应关系 转化公式推导如下:

sin(/2)

tan 2cos(/2)j j jK K

ω

ωΩω==/2/2/2/211j j j j j j e e e K K e e e ωωωωωω

------==++

因只关心频率转换, 故可设s j Ω=,j z

e ω=, 则有

11

11z s K z ---=+, (称为双线性变换) (9)

所以模拟滤波器转换成数字滤波器的公式为

1

1

11()()z a s K z

H z H s ---=+= (10)

由双线性变换公式, 可得11Ks

z Ks

+=

-, 视为两复平面变换, 再由

,

j s j z re ωσΩ=+=.

可得

22

22

(1)()(1)()K K r K K σΩσΩ++=-+ (11)

从上式可得:

0σ=时,1r =, s 上虚轴?z 上单位圆周。 0σ<时,1r <, s 上左半平面? z 上单位圆内。 0σ>时,1r >, s 上右半平面?z 上单位圆外。

故若模拟滤波器稳定，则双线性变换后数字滤波器也稳定。

由于双线性变换法是一种单值映射，因此消除了频率混叠的现象。双线性变换法的缺点是模拟频率与数字频率间的非线性，这种非线性关系要求被变换的连续系统的幅度响应是分段常数型的（某一频率范围内幅度响应近似于某一常数），不然所映射出的数字频率响应相对于原来的模拟频率响应会产生变形。为解决双线性变换中的频率非线性关系，我们采用预畸的方法，即

tan

2

K ω

Ω=，其中K=2/T 。

在Matlab 中利用M 文件bilinear 可以对模拟传输函数实行双线性变换法。

-10

-5

05

10

-2

02w

2 atan(w)

Ω

ωp

ωp

Ωππ

-O

MATLAB 中IIR 数字滤波器的设计过程包括两步。第一步，根据给定指标，确定滤波器的阶数N 和频率缩放因子Wn 。第二步，利用这些参数和给定的波纹参数，确定传输函数的关系。阶数估计：利用双线性变换法设计数字滤波器时，首先要对IIR 数字滤波器的阶数进行估计，相应的M 文件为：buttord 用于巴特沃斯滤波器，cheb1ord 用于切比雪夫1型滤波器，cheb2ord 用于切比雪夫2型滤波器，ellipord 用于椭圆滤波器。滤波器的设计：对于基于双线性变换法的IIR 滤波器设计，对应于四种逼近技术（即巴特沃斯、切比雪夫1型和2型及椭圆逼近），MATLAB 工具箱中有相应的函数。特别地可以用到下面的M 文件：butter 用于巴特沃斯滤波器的设计，cheby1用于切比雪夫1型滤波器的设计，cheby2用于切比雪夫2型滤波器的设计，ellip 用于椭圆滤波器的设计。这些函数的输出可以是滤波器传输函数分子和分母的系数向量，也可以是滤波器的零极点向量和标量增益因子。同时，利用zp2tf 可以由滤波器的零极点向量和标量增益因子得到传输函数分子和分母的系数向量。相应地，利用函数zp2sos 可以得到传输函数分子和分母系数向量的二次项因子。在计算出传输函数的系数之后，可以利用M 文件freqz 来计算频率响应。 2.2 FIR 数字滤波器设计

FIR DF 的冲激响应h(n)是有限长的，M 阶FIR DF 可以表示为：

∑-=-=1

)()()(M i i n x i h n y (12)

其系统函数为:

∑-=-=1

)()(M n n z n h z H (13)

与IIR 数字滤波器的设计不同，FIR 滤波器的设计与模拟滤波器的设计没有任何联系。因此，FIR 滤波器的设计基于对指定幅度响应的直接逼近，并通常要求其具有线性相位响应。为了保证滤波器具有线性相位特性，滤波器系数必须满足条件：h(n)=±h(M-1-n)。

目前关于FIR 滤波器的设计方法主要有三种，即窗函数法，频率取样法和切比雪夫等波纹逼近的最优化设计方法。一般应用较多的是第一种和第三种方法。这是因为窗函数法比较简单，可应用现成的窗函数公式，在技术指标要求不严格的情况下市比较灵活的。最优化设计法必须借助计算机计算，但是它能得到最佳的等波纹的线性相位FIR 滤波器。目前切比雪夫等波纹的线性相位FIR 滤波器的计算机机助设计程序已经比较完善，由于采用了REMEZ 迭代算法，所以设计效率也很高，在应用中越来越占优势。 （1）窗函数法

一般设计过程总是先给定一理想的滤波器频率响应)(jw d e H ，然后设计一个FIR 滤波器，用它的频率响应jwn

M n jw

e

n h e H --=∑=

10

)()(来逼近理想的)(jw d e H 。这种逼近中最直接的方法，是在时域中用FIR 滤

波器的单位脉冲响应h(n)去逼近理想的单位脉冲响应)(n h d 。因而，先由)(jw d e H 的IDTFT 导出)(n h d

dw e e H

n h jwn jw d

d )(21)(?-

=

π

ππ

(14)

由于)(jw d e H 是矩形频率特性，故)(n h d 一定是无限长的序列，且是非因果的。然而FIR 滤波器是有限长的，所以用有限长的h(n)来逼近无限长的)(n h d ，最简单的方法是截取)(n h d 中最重要的一段，将无限长的)(n h d 截取成长度为M 的有限长序列，等效于再)(n h d 上施加了一个长度为M 的矩形窗口，更为一般的，可以用一个长度为M 的窗口函数w(n)来截取)(n h d ，即

)()()(n h n w n h d = (15)

这一方法通常称为窗函数法，窗口函数的形状及长度M 的选择是窗函数法的关键。 下面我们一低通为例，了解一下窗函数法的运用: ①提出希望频率响应函数（低通）

图3 理想低通滤波器的频响 线性相位, 具有片断特点, 即

||()0

||-?≤?

=?

<≤??j j c

d c

e H e ωτω

ωωωωπ

②算出

-1

()()d 2=

?j j n d h n H e e π

ωωπ

ωπ

1d 2--

=?c

c

j j n e e ωωτωω

ωπ

sin(())

()

-=

-c n n ωτπτ(无限长)

图4 理想低通的单位脉冲响应（无限长的一部分） ③加窗()w n ,长N , 得

()()()=d h n h n w n (*)

要线性相位, 就要()h n 关于(1)/2-N 偶对称,而()d h n 关于τ偶对称, 故要求

(1)/2=-N τ

所以要求()w n 关于(1)/2=-N τ偶对称.

O

0.25-π

|()|

j d H e ω0.25π

ω

1

π-π

10

20

30

0.5

1

102030

-0.1

0.1

0.20.3

10

20

30

00.5

1

10

20

30

-0.1

0.1

0.20.3

图5 窗函数 图6 加窗后的单位脉冲响应 再回过来检验()j H e ω是否满足精度要求.

图7 图4的脉冲响应的频响 图8 理想频响与实际频响的对比 若基本满足, 则依截取的()h n , 制硬件, 编软件.

为便于选择使用, 将5种常见的窗函数基本参数如表1所示。

表1 5种常见的窗函数基本参数

类型 窗函数的 旁瓣峰n α 过渡带宽度B ? 加窗后滤波器的 阻带最小衰减s α

rectwin -13

4π/N -21

bartlet 三角 -25 8π/N -25 hanning -31 8π/N -44 hamming -41 8π/N -53 blackman

-57

12π/N

-74

（2）频率取样法

窗口设计法事从时域出发，把理想的)(n h d 用一定形状的窗口函数截取成有限长的h(n)，以此h(n)来近似理想的)(n h d ，从而频率响应)(jw e H 也近似于理想的频率响应)(jw d e H 。我们知道一个有限长序列可以通过其频谱的相同长度的等间隔采样值准确地恢复原有的序列。频率采样法便是从频域出发，对理想的频率响应)(jw d e H 加以等间隔采样 )()|

(2k H e

H d k M

w jw

d ==

π (16)

然后以此)(k H d 作为实际FIR 滤波器的频率特性的离散样本H(k)，即

1,...2,1,0,|

)()()(2-====

M k e H k H k H k M

w jw d d π (17)

由H(k)通过IDFT 可求出有限长序列h(n)为

1

2

3

0.5

1O

0.25-π|()|j d H e ω0.25π

ω

1

π-π

()

j H e ω0

1

2

3

0.51?

1,...1,0,)(1

)(10

2-==

∑-=M n e

k H M

n h M k M

nk

j π (18)

利用M 个频率的离散样本H(k)同样可求出FIR 滤波器的系统函数H(z)及频率响应)(jw e H 。

M

j M k k M M n n

e W z W k H M

z z

n h z H π21

1

1

1)(1)()(--=----=-=--=

=∑∑，其中 (19) 令jw

e z =可得到滤波器的频率响应)(jw e H 。如果设计的是线性相位的FIR 数字滤波器，其采样值H(k)的相位的幅度一定要满足特定的约束条件，这个设计时一定要注意。 （3）最优化设计法

最优化设计法事以最佳一致逼近（最大误差最小化）理论为基础，利用雷米兹算法设计的具有等波纹特性的设计方法。具体设计步骤如下： ①对设计指标进行归一化处理。

②确定remezord 函数所需要的参数。包括归一化边界频率、各频带的幅度要求和波纹要求等。归一化边界频率总是从0开始到1结束，故只需递增列出中间的边界频率；频带幅度要求不含过渡区，个数是边界频率个数的一半加1；波纹要求是频带内幅度允许的波动要求，与分贝间的关系是：

)(log 20)1(log 20),11(

log 202101

2101110δδδ

δδ-≈+-=+--=s p R R (20) ③利用remezord 函数确定remez 所需参数。 ④调用remez 函数进行设计。

⑤利用freqz 函数验算技术指标是否满足要求。 2.3 数字滤波器类型的选择

IIR 和FIR 各有优缺点，在实际运用中如何选择它们，这里做一个简单的比较。

表2 IIR 与FIR 的比较

IIR

FIR

设计方法

利用AF 的设计图表，可简单，有效的完成设计

一般无解析的设计公式，要借助计算机程序完成

设计结果

只能得到幅频特性，相频特性

未知（缺点），如需要线性相位，需用全通网络校准，但增加滤波器的阶数和复杂性

可得到幅频特性（可以多带）和线性相位（优

点）

稳定性 有稳定性问题 极点全部在原点（永远稳定），无稳定性问题 因果性 总是满足，任何一个非因果的有限长序列，总可以通过一定的延时，转变为因果序列 结构 递归系统

非递归

运算误差 有反馈，由于运算中的四舍五入会产生极限环 一般无反馈，运算误差小 快速算法

无快速运算方法

可用FFT 减少运算量

从以上简单的比较可以得到，IIR与FIR滤波器各有所长，所以应根据实际应用要求，从多方面考虑加以选择。

3 数字滤波器的MATLAB设计

3.1 直接程序设计法

（1）IIR的直接程序设计法

例如欲设计一数字（IIR）带阻滤波器，其数字域指标为：数字阻带边缘频率分别为0.4π和0.7π，数字通带边缘频率为0.25π和0.8π，通带波动为1db最小阻带衰减为40db。

此题的MATLAB程序为：

ws=[0.4*pi 0.7*pi]; %数字阻带边缘频率

wp=[0.25*pi 0.8*pi]; %数字通带边缘频率

rp=1 ; %通带波动（db）

as=40; %阻带衰减

[n,wn]=cheb2ord(wp/pi,ws/pi,rp,as);

%根据给定指标，确定滤波器的阶数N和频率缩放因子Wn

[b,a]=cheby2(n,as,ws/pi,'stop');%返回的b,a分别为H（z)的分子、分母。

[h,w]=freqz(b,a,512);%返回的h,w分别为滤波器的频率响应及其频率

plot(w/pi,abs(h));%画出频率响应（以w/pi为横轴）

grid;

xlabel('w/pi');

ylabel('幅值');

title('频率响应');

程序运行结果为：

图9 所设计的带阻滤波器的频率响应

在设计中如果该滤波器的特性不满足要求，原有的参数必须做相应的调整，在程序中只需对参数做新的设定就可以得到所需的滤波器。接下来我们来看看此题所设计的滤波器的滤波效果：S为含有3个频率成分的信号（归一化频率（w/2π）分别为0.1、0.3、0.45），用所设计的滤波器滤除归一化频率为0.3的成分。

n=0:100;

s1=sin(pi*0.2*n);

s2=sin(pi*0.6*n);

s3=sin(pi*0.9*n);

s4=s1+s3;

s=s1+s2+s3;

sf=filter(b,a,s);

subplot(311)

stem(n,s);title('滤波前的信号');

subplot(312);

stem(n,sf);title('滤波后的信号');

subplot(313);

stem(n,s4);title('想要保留的信号');

程序运行的结果为：

图10 采用filter函数进行数字滤波前后信号比较示意图

由图可以看出，滤波后的信号与想要保留的信号基本一致（相位有些许偏差，但基本一致），所以我们可以说该滤波器基本满足了以上所提出的滤波要求。

（2）FIR的直接程序设计法

例如欲设计一个线性相位数字（FIR）带通滤波器，其数字域指标为：数字通带边界频率为0.35π和

0.65π，数字阻带边界频率为0.2π和0.8π，通带波动为1db，最小阻带衰减为60db。

① FIR数字滤波器的窗函数法

此题的MATLAB程序为：

ws1=0.2*pi;

wp1=0.35*pi;

wp2=0.65*pi;

ws2=0.8*pi;

as=60;

tr=min((wp1-ws1),(ws2-wp2));

M=ceil(11*pi/tr)+1;

%滤波器的阶数，程序运行后M=75

n=[0:1:M-1];

r=(M-1)/2;

wc1=(ws1+wp1)/2;

wc2=(wp2+ws2)/2;

hd=sin(wc2*((n-r)+eps))./(pi*((n-r)+eps))-sin(wc1*((n-r)+eps))./(pi*((n-r)+eps));

%hd为理想滤波器的脉冲响应

w_bla=(blackman(M))';

%长度为M的blackman窗

h=hd.*w_bla;

%h为滤波器的实际脉冲响应

stem(n,h);

title('滤波器的实际单位脉冲响应');

freqz(h,1,512);

title('幅度响应和相位响应');

图11 所设计的滤波器的实际单位脉冲响应

由上图可知滤波器的实际脉冲响应h是偶对称的，即h(n)=h(M-1-n),故该滤波器满足FIR线性相位的条件，该滤波器是线性相位的FIR滤波器。

图12 所设计的带通滤波器的幅度和相位响应

由滤波器的相位特性也可以看出该滤波器是线性相位的FIR滤波器。接下来我们来看看此题所设计的滤波器的滤波效果：S为含有3个频率成分的信号（归一化频率（w/2 ）分别为0.05、0.2、0.45），用所设计的滤波器滤除归一化频率为0.05和0.45的成分。

s1=sin(0.1*pi*l);s2=sin(0.4*pi*l);s3=sin(pi*0.9*l);

s=s1+s2+s3;

sf=filter(h,1,s);

subplot(311)

stem(l,s);title('滤波前的信号');

subplot(312);

stem(l,sf);title('滤波后的信号');

subplot(313);

stem(l,s2);title('想要保留的信号');

图13 采用filter函数进行数字滤波前后信号比较示意图

由上图可知滤波后的信号和想要保留的信号的幅度和频率基本一致（滤波后的信号相对于想要保的信号有一个相位延迟，这是线性相位FIR滤波器的群延迟引起的，此滤波器留的群延迟r=(M-1)/2=37），所以我们可以说该滤波器基本满足了以上所提出的滤波要求。

② FIR数字滤波器的频率采样法

此题的MATLAB程序为：

M=40;

%取滤波器的阶数为40

al=(M-1)/2;

%群时延

n=0:M-1;

T2=0.59417456;

T1=0.109021;

Hrs=[zeros(1,5),T1,T2,ones(1,7),T2,T1,zeros(1,9),T1,T2,ones(1,7),T2,T1,zeros(1,4)];

%采样值的幅值

k1=0:floor((M-1)/2);k2=floor((M-1)/2)+1:M-1;

angH=[-al*(2*pi)/M*k1,al*(2*pi)/M*(M-k2)];

%采样值的相位

H=Hrs.*exp(j*angH);

h=real(ifft(H,M));

%长度为M的单位脉冲响应

stem(n,h);

title('滤波器的实际单位脉冲响应');

freqz(h,1,512);

title('幅度响应和相位响应');

图14 所设计的滤波器的实际单位脉冲响应

由图14可知滤波器的实际脉冲响应h是偶对称的，即h(n)=h(M-1-n),故该滤波器满足FIR线性相位的条件，该滤波器是线性相位的FIR滤波器。

图15 所设计的带通滤波器的幅度和相位响应

由滤波器的相位特性也可以看出该滤波器是线性相位的FIR滤波器。此滤波器的群延时为al=(M-1)/2=19.5。接下来我们来看看此题所设计的滤波器的滤波效果：S为含有3个频率成分的信号（归一化频率（w/2 ）分别为0.05、0.2、0.45），用所设计的滤波器滤除归一化频率为0.05和0.45的成分。

l=0:100;

s1=sin(0.1*pi*l);s2=sin(0.4*pi*l);s3=sin(pi*0.9*l);

s=s1+s2+s3;

sf=filter(h,1,s);

subplot(311)

stem(l,s);title('滤波前的信号');

subplot(312);

stem(l,sf);title('滤波后的信号');

subplot(313);

stem(l,s2);title('想要保留的信号');

图16 采用filter函数进行数字滤波前后信号比较示意图

同上面分析相似，滤波后的信号和想要保留的信号的幅度和频率基本一致（滤波后的信号相对于想要保的信号有一个相位延迟，这是线性相位FIR滤波器的群延迟引起的，此滤波器留的群延迟（r=(M-1)/2=19.5），所以我们可以说该滤波器基本满足了以上所提出的滤波要求。

③ FIR数字滤波器的最优设计法

此题的MATLAB程序为：

%设计指标

ws1=0.2*pi;

wp1=0.35*pi;

wp2=0.65*pi;

ws2=0.8*pi;

rp=1;

as=60;

%设置边界频率和幅度要求

F=[ws1/pi,wp1/pi,wp1/pi,ws2/pi];

A=[0,1,0];

%设置各频带的波纹要求

devp=(10^(rp/20)-1)/(10^(rp/20)+1);

devs=10^(-as/20);

dev=[devs,devp,devs];

%确定remez参数，其中滤波器的阶数为（N+1），程序运行后得到N=26

[N,Fo,Ao,W]=remezord(F,A,dev);

%调用remez函数进行设计

h=remez(N,Fo,Ao,W);

n=0:N;

stem(n,h);

title('滤波器的单位冲激响应');

freqz(h,1,512);

title('幅度响应和相位响应');

图17 所设计的滤波器的实际单位脉冲响应

由图17可知滤波器的实际脉冲响应h是偶对称的，即h(n)=h(N-n),故该滤波器满足FIR线性相位的条件，该滤波器是线性相位的FIR滤波器。

图18 所设计的带通滤波器的幅度和相位响应

由滤波器的相位特性也可以看出该滤波器是线性相位的FIR滤波器。此滤波器的群延时为al=(N)/2=13。接下来我们来看看此题所设计的滤波器的滤波效果：S为含有3个频率成分的信号（归一化频率（w/2 ）分别为0.05、0.2、0.45），用所设计的滤波器滤除归一化频率为0.05和0.45的成分。

l=0:100;

s1=sin(0.1*pi*l);s2=sin(0.4*pi*l);s3=sin(pi*0.9*l);

s=s1+s2+s3;

sf=filter(h,1,s);

subplot(311)

stem(l,s);title('滤波前的信号');

subplot(312);

stem(l,sf);title('滤波后的信号');

subplot(313);

stem(l,s2);title('想要保留的信号');

图19采用filter函数进行数字滤波前后信号比较示意图

同上面分析相似，滤波后的信号和想要保留的信号的幅度和频率基本一致（滤波后的信号相对于想要保的信号有一个相位延迟，这是线性相位FIR滤波器的群延迟引起的，此滤波器留的群延迟r=(N)/2=13），所以我们可以说该滤波器基本满足了以上所提出的滤波要求。

3.2 利用信号处理工具箱SPTOOL设计法

SPTOOL是信号处理工具箱中一个具有交互式图形用户界面的信号处理工具，专门用于完成常用的数字信号处理任务。通过这个工具，只需要鼠标简单的操纵鼠标（点击或拖动），就可以完成载入、观察、分析、实现和设计数字滤波器并进行谱分析等甚至十分复杂的数字信号处理任务，而不需要用户对数字滤波器的设计原理非常熟悉。下面是用SPTOOL工具设计的IIR滤波器对上述信号进行滤波。SPTOOL使用步骤：导入信号、滤波器、频谱，设计滤波器，对信号滤波，分析输入输出信号的谱。

1 IIR滤波器的设计

创建并导入信号源，

在MATLAB命令窗口输入命令：

fs=4000（这里取采样频率为4000Hz）;

n=0:100;

s1=sin(pi*0.2*n);

s2=sin(pi*0.6*n);

s3=sin(pi*0.9*n);

s4=s1+s3;

s=s1+s2+s3;

此时，变量fs、n、s、s1、s2、s3、s4将显示在MATLAB的workspace列表中。在MTLAB命令窗口中输入sptool，即可弹出SPTOOL的主界面,如图所示：

图20 SPTOOL的主界面

这里我们首先来设计前面的IIR带阻滤波器，前面的数字指标转化为模拟指标为：通带频率分别为500Hz和1600Hz，阻带频率为800Hz和1400Hz，通带波动为1db，最小阻带衰减为40db。以上三个信号转换为模拟域的频率分别为400Hz，1200Hz，1800Hz。

1）导入信号。使用【file/import】可以导入信号、滤波器和谱。信号的来源可以是MATLAB工作空间变量或MAT数据文件，也可以是在DATA前面的文本框中直接输入信号的数据。这里信号源仍然利用前面程序设计法中的混频信号，将信号s和采样频率fs=4kHz导入并命名为sig1，将信号s2导入命名为sig2,将信号s4导入命名为sig4。

2）滤波器的设计。单击Filter栏中的New按钮打开滤波器设计工具。在界面的最上面制定滤波器的名字、采样频率和设计方法。这里命名滤波器名称为filt1，采样频率fs=4kHz和设计方法为chebyshev typeII IIR；左侧的Specification面板用来指定设计指标；minimum order,type为bandstop，通带频率分别为500Hz和1600Hz，阻带频率为800Hz和1400Hz，通带波动为1db，最小阻带衰减为40db。完成以上输入之后单击Apply按钮，则完成了滤波器的设计；中间的Frequency Response面板则根据设计指标自动绘出滤波器的幅频响应，而右侧的Measurements面板显示了滤波器设计完成后的实测参数。此时在Filters栏中，选中filt1按下View按钮即弹出界面，我们还可以观察到滤波器的相位响应、幅频相位响应、脉冲响应、阶跃响应等等（只须点击界面上的相应快捷按钮）。

图21 所设计的滤波器的幅值响应

3）将滤波器应用到sig1信号序列。分别在Signals、Filters栏中选择sig2、filt1单击Filters 栏列表下的Apply按钮，在弹出的Apply Filter对话框中将输出信号命名为Sig3（滤波后信号）。同时按下Alt与Shift键选中Sig3与Sig2点击Signals栏下的View可观察到它们的时域波形(其中红色是信号

sig3，蓝色是信号sig2)，如图所示。由图可观察出信号s2通过该带阻滤波器后得到很大的衰减，几乎

衰减为0。

图22 信号s2滤波前后的比较

分别在Signals、Filters栏中选择sig1、filt1单击Filters栏列表下的Apply按钮，在弹出的Apply Filter对话框中将输出信号命名为Sig5（滤波后信号）。同时按下Alt与Shift键选中Sig4与Sig5点击Signals栏下的View可观察到它们的时域波形(其中红色是信号sig5，蓝色是信号sig4)，如图所示。由图可观察出信号s通过滤波器后所得的信号与s4基本吻合，只是相位上有点偏差，即该带阻滤波器基本滤除了信号s2，滤波效果良好。

图23 滤波后的信号与s4的比较

4）进行频谱分析。在Signals中选择sig1，单击Spectra栏下的Create按钮，在弹出的Spectra Viewer 界面中选择method为FFT，nfft=512，单击Apply按钮生成s的频谱spect1同理得到sig4的频谱spect4，sig5的频谱spect5,同时按下Alt与Shift键选中spect1与spect5，点击Spectra栏下的View可观察到它们的频谱,下图是sig1和sig5的频谱，其中红色代表滤波后信号的频谱（sig5），蓝色代表滤波前信号的频谱（sig1）。

图24 滤波前后信号的频谱的比较

由频谱可以看出，混频信号s中频率为1200Hz的部分基本被滤除了，频率为400Hz和1800Hz的部分基本被保留下来，未受滤波的影响。同时按下Alt与Shift键选中spect4与spect5，点击Spectra栏下的View 可观察到它们的频谱,下图是sig4和sig5的频谱，其中红色代表滤波后信号的频谱（sig5），蓝色代表滤波前信号的频谱（sig4），如下图。

图25 滤波后的信号与信号s4频谱的比较

由上图可以看出，s4与滤波后的信号频谱基本一致，所以我们可以得出结论：该带阻滤波器的滤波效果较好。

2 FIR滤波器的设计

与1的设计过程类似，这里就简述过程。

创建并导入信号源，

在MATLAB命令窗口输入命令：

fs=4000Hz;

l=0:100;

s1=sin(0.1*pi*l);

s2=sin(0.4*pi*l);

s3=sin(pi*0.9*l);

s=s1+s2+s3;

s4=s1+s3;

前面的数字指标转化为模拟指标为：通带频率分别为700Hz和1300Hz，阻带频率为400Hz和1600Hz，通带波动为1db，最小阻带衰减为60db，以上三个信号转换为模拟域的频率分别为200Hz，800Hz，1800Hz。

1）导入信号。将信号s和采样频率fs=4kHz导入并命名为sig1，将信号s2导入命名为sig2,将信号s4导入命名为sig4。

2）滤波器的设计。

图26 所设计的滤波器的幅值响应

3）将滤波器应用到sig1信号、sig4信号，生成的信号分别命名为sig5和sig3。

图27 信号s4滤波前后的比较

上图中sig4信号与sig3信号比较（其中红色代表sig3，蓝色代表sig4）可知信号s4通过该带通滤波器后得到大大衰减，几乎衰减为0。

图28 滤波后的信号与信号s2的比较

上图中sig2和sig5比较（其中红色代表sig5，蓝色代表sig2），可知信号s通过该带通滤波器后的信号与s2基本具有相同的频率和幅值（只是相位上有一点误差），所以我们认为该滤波器较好的滤掉了s4，滤波性能优良。

4）进行频谱分析。分别生成sig1的频谱spect1，sig2的频谱spect2，sig5的频谱spect5。下图是sig1和sig5的频谱，其中红色代表滤波后信号的频谱（sig5），蓝色代表滤波前信号的频谱（sig1）。由图可看出该滤波器基本滤除了频率为200Hz和1800Hz的信号，保留了频率为800Hz的信号。