用双线性变换法设计IIR数字滤波器

实验三 用双线性变换法设计IIR数字滤波器

实验项目名称：用双线性变换法设计IIR数字滤波器 实验项目性质：验证性实验 所属课程名称：数字信号处理 实验计划学时：2 一. 实验目的

(1)熟悉用双线性变换法设计IIR数字滤波器的原理与方法。 (2)掌握数字滤波器的计算机仿真方法。

(3)通过观察对实际心电图信号的滤波作用，获得数字滤波的感性知识。 二. 实验内容和要求

(1)用双线性变换法设计一个巴特沃斯低通IIR数字滤波器。设计指标参数为：在通带内频率低于0.2π时，最大衰减小于1dB；在阻带内[0.3π，π]频率区间上，最小衰减大与15dB。 (2)以0.02π为采样间隔，打印出数字滤波器在频率区间[0，π/2]的幅频响应特性曲线。 (3)用所设计的滤波器对实际心电图信号采样序列（在本实验后面给出）进行仿真滤波处理，并打印出滤波前后的心电图信号波形图，观察总结滤波作用与效果。 三. 实验主要仪器设备和材料 计算机，MATLAB6.5或以上版本 四. 实验方法、步骤及结果测试

(1)复习有关巴特沃斯模拟滤波器设计和用双线性变换法设计IIR数字滤波器的内容，用双线性变换法设计数字滤波器系统函数

H?z?。其中满足本实验要求的数字滤波器系统函数为：

H?z???1?1.2686z3k?1?1?0.7051z?2??0.00073781?z?1 ?1?2?1?21?1.0106z?0.3583z1?0.9044z?0.2155z??6?????Hk?z?

（3.1）

A1?2z?1?z?2式中： Hk?z??，k?1，2，3 ?1?21?Bkz?CkzA?0.09036B1?1.2686，C1??0.7051?? （3.2） B2?1.0106，C2??0.3583B3?0.9044，C3??0.2155根据设计指标，调用MATLAB信号处理工具箱buttord和butter，也可以得到H?z?。

由公式（3.1）和（3.2）可见，滤波器H?z?由三个二阶滤波器H1?z?、

H2?z?和H3?z?级联而成，如图3-1所示。

x?n?

H1?z? y1?n? H2?z? y2?n? H3?z? y3?n??y?n?

图3-1 滤波器H?z?的组成

(2)编写滤波器仿真程序，计算H?z?对心电图采样序列x(n)的响应序列y(n)。

设yk(n)为第k级二阶滤波器Hk(z)的输出序列，yk-1(n)为输入序列，如

图3-1所示。由（3.2）式可得到差分方程为：

yk?n??Ayk?1?n??2Ayk?1?n?1??Ayk?1?n?2??Bkyk?n?1??Ckyk?n?2?（3.3）

当k=1时，yk-1(n)=x(n)。所以H(z)对x(n)的总响应序列y(n)可以用顺序迭代算法得到。即依次对k=1，2，3，求解差分方程（3.3），最后得到y3(n)=y(n)。仿真程序就是实现上述求解差分方程和顺序迭代算法的通用程序。也可以直接调用MATLAB filter函数实现仿真。 (3)在通用计算机上运行仿真滤波程序，并调用通用绘图子程序，完成实验内容（2）和（3）。 (4)本实验中用到的心电图信号采用序列x(n)

人体心电图信号在测量过程往往受到工业高频干扰，所以必须经过低

通滤波处理后，才能作为判断心脏功能的有用信息。下面给出一实际心电图信号采样序列样本x(n)，其中存在高频干扰。在实验中，以x(n)作为输入序列，滤除其中的干扰成分。

{x(n)}={-4， -2， 0， -4， -6， -4， -2， -4， -6， -6，

-4， -4， -6， -6， -2， 6， 12， 8， 0，-16， -38，-60，-84，-90，-66，-32， -4， -2， -4， 8， 12， 12， 10， 6， 6， 6， 4， 0， 0， 0， 0， 0， -2， -4， 0， 0， 0， -2， -2， 0， 0， -2， -2， -2， -2， 0}

MATLAB程序清单：

%实验三，用双线性变换法设计IIR数字滤波器

x=[-4,-2,0,-4,-6,-4,-2,-4,-6,-6,-4,-4,-6,-6,-2,6,12,8,0,-16,-38,-60,-84,-90,-66,-32,-4,-2,-4,8,12,12,10,6,6,6,4,0,0,0,0,0,-2,-4,0,0,0,-2,-2,0,0,-2,-2,-2,-2,0]; k=1; close all; figure(1); subplot(2,2,1); n=0:55; %更正 stem(n,x,'.');

axis([0,56,-100,50]); %更正

hold on; n=0:60; m=zeros(61); plot(n,m); xlabel('n'); ylabel('x(n)');

title('心电图信号采集序列x(n)'); B=[0.09036,2*0.09036,0.09036]; A=[1.2686,-0.7051]; A1=[1.0106,-0.3583]; A2=[0.9044,-0.2155]; while(k<=3)

y=filter(B,A,x); %The function is to filte the singal x x=y; if k==2 A=A1; end if k==3 A=A2; end k=k+1; end

subplot(2,2,3) n=0:55; %更正 stem(n,y,'.'); axis([0,56,-15,5]); hold on; n=0:60; m=zeros(61);

plot(n,m); xlabel('n'); ylabel('y(n)');

title('三级滤波后的心电图信号'); %求数字滤波器的幅频特性 A=[0.09036,0.1872,0.09036]; B1=[1,-1.2686,0.7051]; B2=[1,-1.0106,0.3583]; B3=[1,-0.9044,0.2155]; [H1,w]=freqz(A,B1,100); [H2,w]=freqz(A,B2,100); [H3,w]=freqz(A,B3,100); H4=H1.*(H2); H=H4.*(H3); mag=abs(H);

db=20*log10((mag+eps)/max(mag)); subplot(2,2,2) plot(w/pi,db); axis([0,0.5,-50,10]); grid on; %更正

title('滤波器的幅频响应曲线');

程序运行结果：

直接运行程序，结果输出滤波器幅频特性曲线图，有噪声的心电图采集信号波形图和经过三级二阶滤波器滤波后的心电图信号波形图，见图3-2，对比分析就可以看出低通滤波器滤除信号中高频噪声的滤波效果。

心电图信号采集序列x(n)5000滤波器的幅频响应曲线x(n)-20-50-40-1000204000.10.20.30.4n三级滤波后的心电图信号50y(n)-5-10-15020n40

五. 实验报告要求 (1) 简述实验原理及目的。

(2) 由所打印的| H(ejw) |特性曲线及设计过程简述双线性变换法的特点。 (3) 对比滤波前后的心电图信号波形，说明数字滤波器的滤波过程与滤波作用。 (4) 简要回答思考题。 六. 思考题

(1) 用双线性变换法设计数字滤波器过程中，变换公式

21?z?1s? ?1T1?z中T的取值对设计结果有无影响？为什么？

实验四 用窗函数法设计FIR数字滤波器

实验项目名称：用窗函数法设计FIR数字滤波器 实验项目性质：验证性实验 所属课程名称：数字信号处理 实验计划学时：2 一. 实验目的

(1)掌握用窗函数法设计FIR数字滤波器的原理与方法。 (2)熟悉线性相位FIR数字滤波器的特性。 (3)了解各种窗函数对滤波特性的影响。 二. 实验内容和要求

(1) 复习用窗函数法设计FIR数字滤波器一节内容，阅读本实验原理，掌握设计步骤。

(2) 用升余弦窗设计一线性相位低通FIR数字滤波器，截止频率??c?rad。窗口长度N =15，33。要求在两种窗口长度情况下，分别4求出h?n?，打印出相应的幅频特性和相频特性曲线，观察3dB带宽和20dB带宽。总结窗口长度N 对滤波器特性的影响。

设计低通FIR数字滤波器时，一般以理想低通滤波特性为逼近函数Hej?，即

Hdej?N?1 2?????j???e，???c??? ??0，?c????其中??1hd?n??2?????Hde??j?1ed??2?j?????ce?j??ej??d??csin??c?n?a?? ??n?a?(3) N?33，?c??4，用四种窗函数设计线性相位低通滤波器，绘制相应的幅频特性曲线，观察3dB带宽和20dB带宽以及阻带最小衰减，比较四种窗函数对滤波器特性的影响。 三. 实验主要仪器设备和材料 计算机，MATLAB6.5或以上版本 四. 实验方法、步骤及结果测试

如果所希望的滤波器的理想的频率响应函数为Hdej?，则其对应的单位脉冲响应为

1hd?n??2?????H?e?ed? （4.1） ?jj?d??窗函数设计法的基本原理是用有限长单位脉冲响应序列h?n?逼近由于hd?n?往往是无限长序列，而且是非因果的，所以用窗函数??n?hd?n?。

将hd?n?截断，并进行加权处理，得到：

h?n??hd?n???n?

（4.2）

h?n?就作为实际设计的FIR数字滤波器的单位脉冲响应序列，其频率响应函数Hej?为

??He????h?n?ej?n?0N?1j?n （4.3）

式中，N为所选窗函数??n?的长度。

我们知道，用窗函数法设计的滤波器性能取决于窗函数??n?的类型及窗口长度N的取值。设计过程中，要根据对阻带最小衰减和过渡带宽度的

要求选择合适的窗函数类型和窗口长度N 。各种类型的窗函数可达到的阻带最小衰减和过渡带宽度见表4.1。

表4.1 各种窗函数的基本参数

窗函数 矩形窗 三角形窗 汉宁窗 哈明窗 不莱克曼窗 凯塞窗(α=7.865)

旁瓣峰值幅度/dB 过渡带宽 阻带最小衰减/dB

-13 -25 -31 -41 -57 -57

4π/N 8π/N 8π/N 8π/N 12π/N 10π/N -12 -25 -44 -53 -74 -80

这样选定窗函数类型和长度N之后，求出单位脉冲响应

h?n??hd?n????n?，并按照式（4.3）求出Hej?。Hej?是否满足要求，要进行演算。一般在h?n?尾部加零使长度满足2的整数次幂，以便用FFT计算Hej?。如果要观察细节，补零点数增多即可。如果Hej?不满足要求，则要重新选择窗函数类型和长度N ，再次验算，直至满足要求。

如果要求线性相位特性，则h?n?还必须满足

????????h?n???h?N?1?n?

根据上式中的正、负号和长度N的奇偶性又将线性相位FIR滤波器分成四类。要根据所设计的滤波特性正确选择其中一类，例如，要设计线性相位低通特性，可以选择h?n??h?N?1?n?这一类，而不能选择

h?n???h?N?1?n?这一类。

主程序框图如图4.1所示。其中幅度特性要求用dB表示。

开始 读入窗口长度N 计算hd(n) 调用窗函数子程序求w(n) 计算h(n)= hd(n) w(n) 调用子程序（函数）计算H(k)=DFT[h(n)] 调用绘图子程序（函数）绘制H(k)幅度相位曲线 结束 图4-1 主程序框图

设

H(k)?DFT[h(n)]H(k)?HR(k)?jHI(k)

2H(k)?HR(k)?HI2(k)画图时，用20lgH(k)打印幅度特性。第k点对应的频率?k?2?k。N为使曲线包络更接近Hej?的幅度特性曲线，DFT变换区间要选大些。例如窗口长度N=33时，可通过在h?n?末尾补零的方法，使长度变为64，再进行64点DFT，则可以得到更精确的幅度衰减特性曲线。

??

下面给出MATLAB主程序：

%实验四，用窗函数法设计FIR数字滤波器 b=1; close all; i=0; while(b);

temp=menu('选择窗函数长度

N','N=10','N=15','N=20','N=25','N=30','N=33','N=35','N=40','N=45','N=50','N=55','N=60','N=64');

menu1=[10,15,20,25,30,33,35,40,45,50,55,60,64]; N=menu1(temp);

temp=menu('选择逼近理想低通滤波器截止频率

Wc','Wc=pi/4','Wc=pi/2','Wc=3*pi/4','Wc=pi','Wc=0.5','Wc=1.0','Wc=1.5','Wc=2.0','Wc=2.5','Wc=3.0');

menu2=[pi/4,pi/2,3*pi/4,pi,0.5,1,1.5,2,2.5,3]; w=menu2(temp); n=[0:(N-1)];

hd=ideal(w,N); %得到理想低通滤波器 k=menu('请选择窗口类

型:','boxcar','hamming','hanning','blackman'); if k==1

B=boxcar(N);

string=['Boxcar','N=',num2str(N)]; else if k==2

B=hamming(N);

string=['Hamming','N=',num2str(N)]; else if k==3

B=hanning(N);

string=['Hanning','N=',num2str(N)]; else if k==4

B=blackman(N);

string=['Blackman','N=',num2str(N)]; end end end end

h=hd.*(B)'; %得到FIR数字滤波器

[H,m]=freqz(h,[1],1024,'whole'); %求其频率响应 mag=abs(H); %得到幅值 db=20*log10((mag+eps)/max(mag)); pha=angle(H); %得到相位 i=i+1; figure(i) subplot(2,2,1); n=0:N-1; stem(n,h,'.');

axis([0,N-1,-0.1,0.3]); hold on; n=0:N-1; x=zeros(N); plot(n,x,'-'); xlabel('n'); ylabel('h(n)');

title('实际低通滤波器的h(n)'); text((0.3*N),0.27,string); hold off; subplot(2,2,2); plot(m/pi,db); axis([0,1,-100,0]);

xlabel('w/pi'); ylabel('dB');

title('衰减特性(dB)'); grid;

subplot(2,2,3); plot(m,pha); hold on; n=0:7; x=zeros(8); plot(n,x,'-'); title('相频特性'); xlabel('频率(rad)'); ylabel('相位(rad)'); axis([0,3.15,-4,4]); subplot(2,2,4); plot(m,mag); title('频率特性'); xlabel('频率W(rad)'); ylabel('幅值'); axis([0,3.15,0,1.5]); text(0.9,1.2,string);

b=menu('Do You want To Continue ?','Yes','No'); if b==2 b=0; end end

temp=menu('Close All Figure ?','Yes','No'); if temp==1

close all end

%实验四，用窗函数法设计FIR数字滤波器

%实验四中的子函数：产生理想低通滤波器单位脉冲响应hd(n) 新建.m文件，把下列绿色区域语句复制，保存为ideal.m function hd=ideal(w,N); alpha=(N-1)/2; n=[0:(N-1)]; m=n-alpha+eps; hd=sin(w*m)./(pi*m);

新建.m文件，把下列红色区域语句复制，保存为main.m 注意：ideal.m 与main.m 一定要在同一文件夹内。

close all; i=0; N=10; w=pi/4; n=[0:(N-1)];

hd=ideal(w,N); %得到理想低通滤波器 B=boxcar(N);

h=hd.*(B)'; %得到FIR数字滤波器

[H,m]=freqz(h,[1],1024,'whole'); %求其频率响应 mag=abs(H); %得到幅值 db=20*log10((mag+eps)/max(mag)); pha=angle(H); %得到相位 i=i+1; figure(i) subplot(2,2,1); n=0:N-1;

stem(n,h,'.');

axis([0,N-1,-0.1,0.3]); hold on; n=0:N-1; x=zeros(N); plot(n,x,'-'); xlabel('n'); ylabel('h(n)');

title('实际低通滤波器的h(n)'); text((0.3*N),0.27,'boxcar,N=10'); hold off; subplot(2,2,2); plot(m/pi,db); axis([0,1,-100,0]); xlabel('w/pi'); ylabel('dB');

title('衰减特性(dB)'); grid;

subplot(2,2,3); plot(m,pha); hold on; n=0:7; x=zeros(8); plot(n,x,'-'); title('相频特性'); xlabel('频率(rad)'); ylabel('相位(rad)'); axis([0,3.15,-4,4]); subplot(2,2,4);

plot(m,mag); title('频率特性'); xlabel('频率W(rad)'); ylabel('幅值'); axis([0,3.15,0,1.5]); text(0.9,1.2,'boxcar,N=10'); 程序运行结果：

运行程序，根据实验内容要求和程序提示选择你要进行的实验参数。三个实验参数选定后，程序运行输出用所选窗函数设计的实际FIR低通数字滤波器的单位脉冲响应h(n)、幅频衰减特性(20lg|H(ejw)|)、相频特性及幅频特性|H(ejw)|的波形，h(n)和|H(ejw)|图中标出了所选窗函数类型及其长度N值。对四种窗函数（N=15和N=33）的程序运行结果如图4-2到图4-9所示，由图可以看出用各种窗函数设计的FIR滤波器的阻带最小衰减及过渡带均与教材中一致。在通带内均为严格相位特性。

五. 实验报告要求 (1) 简述实验原理及目的。

(2) 按照实验步骤以及要求，比较各种情况下的滤波性能，说明窗口长

度N和窗函数类型对滤波特性的影响。 (3) 总结用窗函数法设计FIR滤波器的主要特点。 (4) 简要回答思考题。 六. 思考题

(1) 如果给定通带截止频率和阻带截止频率以及阻带最小衰减，如何用窗函数法设计线性相位低通滤波器，写出设计步骤。 (2) 如果要求用窗函数法设计带通滤波器，而且给定上、下边带截止频率为?1和?2，试求理想带通的单位脉冲响应hd?n?。

附录一 MATLAB信号处理工具箱函数

MATLAB信号处理工具箱包含了许多信号处理函数。这些函数的格式和使用说明很容易通过HELP命令查到。调用这些函数很容易验证数字信号处理习题答案，使数字滤波器分析与设计的程序非常简单。例如，已知系统函数H(z)，手工画其幅频响应曲线和相频曲线是相当困难的。但只要调用函数freqz(b，a)就可以画出来(b和a分别为H(z)的分子和分母多项式系数)。为了便于读者分类快速查询，下面按功能列出MATLAB信号处理工具箱函数。

一. 表附1-1 波形产生

函数名 sawtooth square sinc diric

二. 表附1-2 滤波器分析和实现

产生锯齿波或三角波 产生方波 产生sinc或者

sin??t?函数 ?t功能

产生Dirichlet或者周期sinc函数

函数名 abs angle conv fftfilt filter

求绝对值（幅值） 求相角 求卷积

功能

重叠相加法FFT滤波器实现 直接滤波器实现

filtfilt filtic freqs freqspace freqz grpdelay impz zplane

零相位数字滤波 filter函数初始条件选择 模拟滤波器频率响应 频率响应中的频率间隔 数字滤波器频率响应 平均滤波器延迟（群延迟） 数字滤波器的冲激响应 离散系统零极点图

三. 表附1-3 线性系统变换

函数名 convmtx poly2rc rc2poly residuez sos2ss sos2tf sos2zp ss2sos ss2tf ss2zp tf2ss

卷积矩阵

功能

从多项式系数中计算反射系数 从反射系数中计算多项式系数 Z变换部分分式展开或留数计算 变系统二阶分割形式为状态空间形式 变系统二阶分割形式为传递函数形式 变系统二阶分割形式为零极点增益形式 变系统状态空间形式为二阶分割形式 变系统状态空间形式为传递函数形式 变系统状态空间形式为零极点增益形式 变系统传递函数形式为状态空间形式

tf2zp zp2sos zp2ss zp2tf

变系统传递函数形式为零极点增益形式 变系统零极点增益形式为二阶分割形式 变系统零极点增益形式为状态空间形式 变系统零极点增益形式为传递函数形式

四. 表附1-4 IIR滤波器设计

函数名 besself butter cheby1 cheby2 ellip yulewalk

功能

Bessel（贝塞尔）模拟滤波器设计 Butterworth（比特沃思）滤波器设计 Chebyshev （切比雪夫）I型滤波器设计 Chebyshev （切比雪夫）II型滤波器设计 椭圆滤波器设计 递归数字滤波器设计

五. 表附1-5 IIR滤波器阶的选择

函数名 buttord cheb1ord cheb2ord ellipord

功能

Butterworth（比特沃思）滤波器阶的选择 Chebyshev （切比雪夫）I型滤波器阶的选择 Chebyshev （切比雪夫）II型滤波器阶的选择 椭圆滤波器阶的选择

附录二 实验中用到的一些子程序

%实验二中的子函数:通用画图程序 function printf(j,k,N,b) subplot(2,2,1); if(b~=9)

if b==1|b==2|b==3

M=3.2; n=0:7;

stem(n,j,'.'); n=0:7; m=zeros(8);

%确定图中标注的X轴坐标

else

M=N*0.4; n=0:N-1; stem(n,j,'.'); n=0:N-1; m=zeros(N);

end

hold on; plot(n,m); t=max(j); xlabel('n');

string=['x',num2str(b),'(n)的波形']; ylabel('x(n)'); text(M,(t*0.8),string); end

subplot(2,2,3); n=0:N-1; stem(n,k,'.'); t=max(k); xlabel('K');

string=['x',num2str(b),'(n)的N=',num2str(N),'点FFT']; text((N*0.4),(t*0.8),string); ylabel('|X(k)|');

%实验四中的子函数：产生理想低通滤波器单位脉冲响应hd(n) function hd=ideal(w,N); alpha=(N-1)/2; n=[0:(N-1)]; m=n-alpha+eps; hd=sin(w*m)./(pi*m);

本文来源：https://www.bwwdw.com/article/hae2.html