实验三IIR滤波器的设计与信号滤波.docx

实验三IIR滤波器的设计与信号滤波.docx

《实验三IIR滤波器的设计与信号滤波.docx》由会员分享，可在线阅读，更多相关《实验三IIR滤波器的设计与信号滤波.docx（12页珍藏版）》请在冰豆网上搜索。

实验三IIR滤波器的设计与信号滤波

实验三IIR滤波器的设计与信号滤波

实验三　IIR滤波器的设计与信号滤波

１、实验目的

（１）熟悉用双线性变换法设计IIR数字滤波器的原理与方法。

（２）掌握数字滤波器的计算机仿真方法。

（３）通过观察对实际心电图信号的滤波作用，获得数字滤波的感性知识。

2、实验仪器：

PC机一台MATLAB软件

3、实验原理

利用双线性变换设计IIR滤波器（只介绍巴特沃斯数字低通滤波器的设计），首先要设计出满足指标要求的模拟滤波器的传递函数

，然后由

通过双线性变换可得所要设计的IIR滤波器的系统函数

。

如果给定的指标为数字滤波器的指标，则首先要转换成模拟滤波器的技术指标，这里主要是边界频率

的转换，对

指标不作变化。

边界频率的转换关系为

。

接着，按照模拟低通滤波器的技术指标根据相应设计公式求出滤波器的阶数

和

截止频率

；根据阶数

查巴特沃斯归一化低通滤波器参数表，得到归一化传输函数

；最后，将

代入

去归一，得到实际的模拟滤波器传输函数

。

之后，通过双线性变换法转换公式

，得到所要设计的IIR滤波器的系统函数

。

利用所设计的数字滤波器对实际的心电图采样信号进行数字滤波器。

4、实验步骤及内容

（１）复习有关巴特沃斯模拟滤波器的设计和用双线性变换法设计IIR数字滤波器的内容，用双线性变换法设计一个巴特沃斯IIR低通数字滤波器。

设计指标参数为：

在通带内频率低于

时，最大衰减小于

；在阻带内

频率区间上，最小衰减大于

。

（２）绘制出数字滤波器的幅频响应特性曲线。

（３）用所设计的滤波器对实际心电图信号采样序列（实验数据在后面给出）进行仿真滤波处理，并分别绘制出滤波前后的心电图信号波形图，观察总结滤波作用与效果。

（４）输入为20Hz正弦和200Hz的正弦的叠加波形，要求用双线性变换法设计一巴特沃斯数字低通滤波器滤除200Hz的正弦，使输出中只保留20Hz的正弦波。

并绘制出滤波前和滤波后的波形。

5、实验用MATLAB函数介绍

buttord（）;butter（）;bilinear（）;freqz（）;freqs（）;filter（）;figure（）;plot（）;stem（）;abs（）;title（）;xlabel（）;ylabel（）;text（）;holdon;axis（）;gridon;subplot（）;等

6、思考题

（１）用双线性变换法设计数字滤波器过程中，变换公式

中T的取值，对设计结果有无影响?

为什么?

（２）如果用脉冲响应不变法设计该IIR数字低通滤波器，程序如何改动？

7、实验报告要求

（１）简述实验目的及实验原理。

（２）编程实现各实验内容，列出实验清单及说明。

（３）由绘制的

特性曲线及设计过程简述双线性变换法的特点。

（４）对比滤波前后的心电图信号波形，说明数字滤波器的滤波过程与滤波作用。

（５）简要回答思考题。

8、心电图信号采样序列

人体心电图信号在测量过程中往往受到工业高频干扰，所以必须经过低通滤波处理后，才能作为判断心脏功能的有用信息。

下面给出的数据是一实际心电图信号采样序列样本

，其中存在高频干扰。

本实验中，以

作为输入序列，滤除其中的干扰成分。

附录

一、常用的MATLAB函数介绍

1．plot

功能：

线型绘图函数。

格式：

plot（v）；plot（x,y）

说明：

plot（v），v是长度为N的数值向量。

作用是在坐标系中顺序地用直线连接顶点，生成一条折（曲）线。

当向量元素充分多时，即可生成一条光滑的曲线。

在实验中，若FFT点数足够多时，用plot打印的幅频特性就很接近

连续曲线。

plot（x,y），x和y都是长度为n的向量。

作用在坐标系中生成顺序连接顶点的折（曲）线。

这种调用可被用来生成参数方程的图形。

２．stem

功能：

绘制离散序列图。

格式：

stem（y）；stem（x,y）；stem（…,’线端符号’）；stem（…,’线型’）；

stem（…,‘线型’，‘线端符号’）

说明：

　　stem（y）和stem（x,y）分别与plot（v）和plot（x,y）的绘图规则相同，只是stem绘制的是离散序列图。

实验中用于绘制时域序列

的波形图和序列的离散傅里叶变换

的幅度图。

３．subplot

功能：

多坐标设置与定位当前坐标系。

格式：

subplot（m,n,k）

说明：

subplot（m,n,k）将图形窗口分成m行n列的m×n块子区域，按行从上到下，从左到右的顺序，在第i块子区定义一个坐标系，使其成为当前坐标系，随后的绘图函数将在该坐标系输出图形。

另外，同一个图形窗口的坐标系可以重叠，这样可以产生前面的坐标系遮住后面坐标系的各种图形效果。

４．figure

功能：

创建新的图形窗口（用于输出图形的窗口）。

格式：

figure；figure（h）

说明：

figure　函数创建一个新的图形窗口，并成为当前图形窗口，所创建的图形窗口的序号（句柄值）是按同一MATLAB程序中创建的顺序号。

使用figure（h）函数，该方法常用在程序设计中，用于控制将各种波形图输出到相应的图形窗口中。

打印输出或存储时，一个图形窗口打印一张图纸或存储一个图形文件。

５．axis

功能：

设置图形的坐标范围。

格式：

axis（[xminxmaxyminymax]）

６．gridon

功能：

画出图形的分格线。

7．title

功能：

书写图名。

格式：

title（‘s’）

8．xlabel

功能：

横坐标名称。

格式：

xlabel（‘s’）

9．ylabel

功能：

纵坐标名称。

格式：

ylabel（‘s’）

10．text

功能：

在图面（x,y）位置处书写字符注释。

格式：

text（x,y,‘s’）

11．holdon

功能：

通过该命令，保持当前图形不变化，准备在当前的图形窗口上绘制新的曲线。

12．abs

功能：

求绝对值（模值）。

格式：

y=abs（x）

说明：

y=abs（x）用于计算x的绝对值，当x为复数时，得到的是复数的模值。

当x为字符串时，abs（x）得到字符串的各个字符的ASCII码，例如，x=’123’,则，abs（x）得到：

495051。

13．angle

功能：

求相角。

格式：

ψ＝angle（h）

说明：

ψ＝angle（h）用于求复矢量或复矩阵的相角（以弧度为单位），相角介于

和

之间。

14．conv

功能：

求卷积。

格式：

c=conv（a,b）

说明：

conv（a,b）用于求矢量a和b的卷积，即

式中Ｎ为矢量a和b的最大长度。

例如，当a=[123],b=[456]时，则

c=conv（a,b）

c=

413282718

此函数可直接用于求两个有限长序列的卷积。

设x（n）和h（n）的长度分别为M和N，则计算二者卷积的MATLAB语句如下：

y=conv（x,h）

y的长度为M+N-1。

15．filter

功能：

利用IIR滤波器或FIR滤波器对数据进行滤波。

格式：

y=filter（b,a,x）;[y,zf]=filter（b,a,x）;y=filter（b,a,x,zi）

说明：

filter利用数字滤波器对数据进行滤波，其实现采用直接Ⅱ型结构，因而适用于IIR和FIR两种滤波器。

滤波器的系统函数为

即滤波器系数a=[a0a1a2…aN],b=[b0b1…bM]。

这里的标准形式为a0=1，如果输入矢量a时，a0≠1，则MATLAB将自动进行归一化系数的操作；如果a0=０，则给出出错信息。

y=filter（b,a,x）利用给定系数矢量a和b对x中的数据进行滤波，结果放入y矢量中。

y=filter（b,a,x,zi）可在zi中指定x的初始状态。

[y,zf]=filter（b,a,x）除得到矢量y外，还得到x的最终状态矢量zf。

16．fftfilt

功能：

利用重叠相加法实现短序列b（n）和长序列x（n）的块之间的卷积。

格式：

y=fftfilt（b,x）

17．freqz

功能：

数字滤波器的频率响应。

格式：

［H,w］=freqz（b,a,Ｎ）；［Ｈ，f］=freqz（b,a,N,Fs）;

H=freqz（b,a,w）;H=freqz（b,a,f,Fs）;freqz（b,a）

说明：

freqz用于计算数字滤波器H（z）的频率响应函数

。

［H,w］=freqz（b,a,Ｎ）可得到数字滤波器的Ｎ点频率响应值，这Ｎ个点均匀地分布在［０，

］上，并将这Ｎ个频点的频率记录在w中，相应的频响值记录在Ｈ中。

要求Ｎ为大于零的整数，最好为２的整数次幂，以便采用FFT计算，以提高速度。

［Ｈ，f］=freqz（b,a,N,Fs）用于对

在［０，

］上等间隔采样Ｎ点，采样点频率及相应频响值分别记录在f和H中。

由用户指定

（以Hz为单位）值。

H=freqz（b,a,w）用于对

在［０，

］上进行采样，采样频率点由矢量w指定。

H=freqz（b,a,f,Fs）用于对

在［０，

］上采样，采样频率点由矢量f给定。

freqz（b,a）用于在当前窗口中绘制出幅频和相频特性曲线。

18．impz

功能：

计算H（z）相应的单位脉冲响应h（n）。

格式：

［h,t］=impz（b,a）;[h,t]=impz（b,a,z）;

19．fft

功能：

一维快速傅里叶变换（FFT）

格式：

y=fft（x）;y=fft（x,N）

说明：

fft函数用于计算矢量或矩阵的离散傅里叶变换。

y=fft（x）利用FFT算法计算矢量x的离散傅里叶变换，当x为矩阵时，y为矩阵x每一列的FFT。

当x的长度为2的整数次幂时，fft采用基2FFT算法，否则采用稍慢的混合基算法。

y=fft（x,N）采用Ｎ点FFT。

当x长度小于Ｎ时，fft函数自动在x尾部补零，以构成Ｎ点数据；当x的长度大于Ｎ时，fft截取x的前面Ｎ点数据进行FFT。

20．ifft

功能：

一维逆快速傅里叶变换（IFFT）

格式：

y=ifft（x）;y=ifft（x,N）

二、IIR滤波器的设计与实现

基于模拟滤波器变换原理IIR滤波器的经典设计：

首先，根据模拟滤波器的指标设计出相应的模拟滤波器，然后，将设计好的模拟滤波器转换成满足给定技术指标的数字滤波器。

通常算法有脉冲响应不变法和双线性变换法。

在MATLAB的数字信号处理工具箱中提供了相应的设计函数。

常用的有：

（１）Butterworth滤波器阶数选择函数

[N,Wn]=buttord（Wp,Ws,Rp,Rs）

输入参数：

Wp通带截止频率，Ws阻带截止频率，Rp通带最大衰减，Rs阻带最小衰减；

输出参数：

N符合要求的滤波器最小阶数，Wn为Butterworth滤波器固有频率（3dB）。

（２）零极点增益模型到传递函数模型的转换

[num,den]=zp2tf（Z,P,K）;

输入参数：

Z,P,K分别表示零极点增益模型的零点、极点和增益；

输出参数：

num,den分别为传递函数分子和分母的多项式系数。

（３）从低通向低通的转换

[b,a]=lp2lp（Bap,Aap,Wn）;

功能：

把模拟滤波器原型转换成截至频率为Wn的低通滤波器。

（４）双线性变换函数

[bz,az]=bilinear（b,a,Fs）;

功能：

把模拟滤波器的零极点模型转换为数字滤波器的零极点模型，其中Fs是采样频率。

例4：

用双线性变换法设计一个Butterworth数字低通滤波器，要求其通带截至频率

，阻带截至频率

，通带衰减

小于0.75dB，阻带衰减

大于20dB。

MATLAB实现程序：

wp=0.2*pi;ws=0.4*pi;

ap=0.75;as=20;

T=1;

Wp=2/T*tan（wp/2）;%求模拟原型滤波器参数

Ws=2/T*tan（ws/2）;

[N,Wc]=buttord（Wp,Ws,ap,as,'s'）;%求滤波器的阶数

[BS,AS]=butter（N,Wc,'s'）;%得到模拟原型滤波器的系统函数

[BZ,AZ]=bilinear（BS,AS,1/T）;%双线性变换法转换到数字滤波器的系统函数

[H,W]=freqs（BS,AS）;

plot（W/pi,20*log10（abs（H）））;%画模拟原型滤波器的频响

holdon;grid;

[H,W]=freqz（BZ,AZ）;

plot（W/pi,20*log10（abs（H））,'r'）;%画数字滤波器的频响

运行结果：

Butterworth低通滤波器的频率响应图