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

1、实验三 用双线性变换法设计IIR数字滤波器实验项目名称：用双线性变换法设计IIR数字滤波器实验项目性质：验证性实验所属课程名称：数字信号处理实验计划学时：2一. 实验目的(1) 熟悉用双线性变换法设计IIR数字滤波器的原理与方法。(2) 掌握数字滤波器的计算机仿真方法。(3) 通过观察对实际心电图信号的滤波作用，获得数字滤波的感性知识。二. 实验内容和要求(1) 用双线性变换法设计一个巴特沃斯低通IIR数字滤波器。设计指标参数为：在通带内频率低于0.2时，最大衰减小于1dB；在阻带内0.3，频率区间上，最小衰减大与15dB。(2) 以0.02为采样间隔，打印出数字滤波器在频率区间0，/2的幅频

2、响应特性曲线。(3) 用所设计的滤波器对实际心电图信号采样序列（在本实验后面给出）进行仿真滤波处理，并打印出滤波前后的心电图信号波形图，观察总结滤波作用与效果。三. 实验主要仪器设备和材料计算机，MATLAB6.5或以上版本四. 实验方法、步骤及结果测试(1) 复习有关巴特沃斯模拟滤波器设计和用双线性变换法设计IIR数字滤波器的内容，用双线性变换法设计数字滤波器系统函数。其中满足本实验要求的数字滤波器系统函数为： （3.1）式中： （3.2）根据设计指标，调用MATLAB信号处理工具箱buttord和butter，也可以得到。图3-1 滤波器的组成由公式（3.1）和（3.2）可见，滤波器由三个

3、二阶滤波器、和级联而成，如图3-1所示。(2) 编写滤波器仿真程序，计算对心电图采样序列x(n)的响应序列y(n)。设yk(n)为第k级二阶滤波器Hk(z)的输出序列，yk-1(n)为输入序列，如图3-1所示。由（3.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) 在通用计算机上运行仿真滤波程序，并调用通用绘图

4、子程序，完成实验内容（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

5、， 0， 0， 0， 0， -2， -4， 0， 0， 0， -2， -2， 0，0， -2， -2， -2， -2， 0MATLAB程序清单：%实验三，用双线性变换法设计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; %更

6、正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

7、singal x x=y; if k=2 A=A1; end if k=3 A=A2; end k=k+1;endsubplot(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,

8、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('滤波器的幅频响应曲线');程序运行结果：直接运行程序，结果输出滤波器幅频特性曲线图，有噪声的心电图采集信号波形图和经过三级二阶滤波器滤波后的心电

9、图信号波形图，见图3-2，对比分析就可以看出低通滤波器滤除信号中高频噪声的滤波效果。图3-2 心电图信号波形图五. 实验报告要求(1) 简述实验原理及目的。(2) 由所打印的| H(ejw) |特性曲线及设计过程简述双线性变换法的特点。(3) 对比滤波前后的心电图信号波形，说明数字滤波器的滤波过程与滤波作用。(4) 简要回答思考题。六. 思考题(1) 用双线性变换法设计数字滤波器过程中，变换公式中T的取值对设计结果有无影响？为什么？实验四 用窗函数法设计FIR数字滤波器实验项目名称：用窗函数法设计FIR数字滤波器实验项目性质：验证性实验所属课程名称：数字信号处理实验计划学时：2一. 实验目的(

10、1) 掌握用窗函数法设计FIR数字滤波器的原理与方法。(2) 熟悉线性相位FIR数字滤波器的特性。(3) 了解各种窗函数对滤波特性的影响。二. 实验内容和要求(1) 复习用窗函数法设计FIR数字滤波器一节内容，阅读本实验原理，掌握设计步骤。(2) 用升余弦窗设计一线性相位低通FIR数字滤波器，截止频率。窗口长度N =15，33。要求在两种窗口长度情况下，分别求出，打印出相应的幅频特性和相频特性曲线，观察3dB带宽和20dB带宽。总结窗口长度N 对滤波器特性的影响。设计低通FIR数字滤波器时，一般以理想低通滤波特性为逼近函数，即其中(3) ，用四种窗函数设计线性相位低通滤波器，绘制相应的幅频特性

11、曲线，观察3dB带宽和20dB带宽以及阻带最小衰减，比较四种窗函数对滤波器特性的影响。三. 实验主要仪器设备和材料计算机，MATLAB6.5或以上版本四. 实验方法、步骤及结果测试如果所希望的滤波器的理想的频率响应函数为，则其对应的单位脉冲响应为 （4.1）窗函数设计法的基本原理是用有限长单位脉冲响应序列逼近。由于往往是无限长序列，而且是非因果的，所以用窗函数将截断，并进行加权处理，得到：（4.2）就作为实际设计的FIR数字滤波器的单位脉冲响应序列，其频率响应函数为（4.3）式中，N为所选窗函数的长度。我们知道，用窗函数法设计的滤波器性能取决于窗函数的类型及窗口长度N的取值。设计过程中，要根据

12、对阻带最小衰减和过渡带宽度的要求选择合适的窗函数类型和窗口长度N 。各种类型的窗函数可达到的阻带最小衰减和过渡带宽度见表4.1。表4.1 各种窗函数的基本参数窗函数旁瓣峰值幅度/dB过渡带宽阻带最小衰减/dB矩形窗-134/N-12三角形窗-258/N-25汉宁窗-318/N-44哈明窗-418/N-53不莱克曼窗-5712/N-74凯塞窗(=7.865)-5710/N-80这样选定窗函数类型和长度N之后，求出单位脉冲响应，并按照式（4.3）求出。是否满足要求，要进行演算。一般在尾部加零使长度满足2的整数次幂，以便用FFT计算。如果要观察细节，补零点数增多即可。如果不满足要求，则要重新选择窗函

13、数类型和长度N ，再次验算，直至满足要求。如果要求线性相位特性，则还必须满足根据上式中的正、负号和长度N的奇偶性又将线性相位FIR滤波器分成四类。要根据所设计的滤波特性正确选择其中一类，例如，要设计线性相位低通特性，可以选择这一类，而不能选择这一类。主程序框图如图4.1所示。其中幅度特性要求用dB表示。开始读入窗口长度N计算hd(n)调用窗函数子程序求w(n)调用子程序（函数）计算H(k)=DFTh(n)调用绘图子程序（函数）绘制H(k)幅度相位曲线结束图4-1 主程序框图计算h(n)= hd(n) w(n)设画图时，用打印幅度特性。第k点对应的频率。为使曲线包络更接近的幅度特性曲线，DFT变

14、换区间要选大些。例如窗口长度N=33时，可通过在末尾补零的方法，使长度变为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

15、=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','

16、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); els

17、e 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&

18、#39;); %求其频率响应 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)');

19、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

20、('相位(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; endendtemp=menu('Close All F

21、igure ?','Yes','No');if temp=1 close allend%实验四，用窗函数法设计FIR数字滤波器%实验四中的子函数：产生理想低通滤波器单位脉冲响应hd(n)新建.m文件，把下列绿色区域语句复制，保存为ideal.mfunction 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;

22、 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:

23、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; subpl

24、ot(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.

25、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) 简述实验原理及目

26、的。(2) 按照实验步骤以及要求，比较各种情况下的滤波性能，说明窗口长度N和窗函数类型对滤波特性的影响。(3) 总结用窗函数法设计FIR滤波器的主要特点。(4) 简要回答思考题。六. 思考题(1) 如果给定通带截止频率和阻带截止频率以及阻带最小衰减，如何用窗函数法设计线性相位低通滤波器，写出设计步骤。(2) 如果要求用窗函数法设计带通滤波器，而且给定上、下边带截止频率为和，试求理想带通的单位脉冲响应。附录一 MATLAB信号处理工具箱函数MATLAB信号处理工具箱包含了许多信号处理函数。这些函数的格式和使用说明很容易通过HELP命令查到。调用这些函数很容易验证数字信号处理习题答案，使数字滤波器

27、分析与设计的程序非常简单。例如，已知系统函数H(z)，手工画其幅频响应曲线和相频曲线是相当困难的。但只要调用函数freqz(b，a)就可以画出来(b和a分别为H(z)的分子和分母多项式系数)。为了便于读者分类快速查询，下面按功能列出MATLAB信号处理工具箱函数。一. 表附1-1 波形产生函数名功能sawtooth产生锯齿波或三角波square产生方波sinc产生sinc或者函数diric产生Dirichlet或者周期sinc函数二. 表附1-2 滤波器分析和实现函数名功能abs求绝对值（幅值）angle求相角conv求卷积fftfilt重叠相加法FFT滤波器实现filter直接滤波器实现fi

28、ltfilt零相位数字滤波filticfilter函数初始条件选择freqs模拟滤波器频率响应freqspace频率响应中的频率间隔freqz数字滤波器频率响应grpdelay平均滤波器延迟（群延迟）impz数字滤波器的冲激响应zplane离散系统零极点图三. 表附1-3 线性系统变换函数名功能convmtx卷积矩阵poly2rc从多项式系数中计算反射系数rc2poly从反射系数中计算多项式系数residuezZ变换部分分式展开或留数计算sos2ss变系统二阶分割形式为状态空间形式sos2tf变系统二阶分割形式为传递函数形式sos2zp变系统二阶分割形式为零极点增益形式ss2sos变系统状态空

29、间形式为二阶分割形式ss2tf变系统状态空间形式为传递函数形式ss2zp变系统状态空间形式为零极点增益形式tf2ss变系统传递函数形式为状态空间形式tf2zp变系统传递函数形式为零极点增益形式zp2sos变系统零极点增益形式为二阶分割形式zp2ss变系统零极点增益形式为状态空间形式zp2tf变系统零极点增益形式为传递函数形式四. 表附1-4 IIR滤波器设计函数名功能besselfBessel（贝塞尔）模拟滤波器设计butterButterworth（比特沃思）滤波器设计cheby1Chebyshev （切比雪夫）I型滤波器设计cheby2Chebyshev （切比雪夫）II型滤波器设计ell

30、ip椭圆滤波器设计yulewalk递归数字滤波器设计五. 表附1-5 IIR滤波器阶的选择函数名功能buttordButterworth（比特沃思）滤波器阶的选择cheb1ordChebyshev （切比雪夫）I型滤波器阶的选择cheb2ordChebyshev （切比雪夫）II型滤波器阶的选择ellipord椭圆滤波器阶的选择六. 表附1-6 FIR滤波器设计函数名功能fir1基于窗函数的FIR滤波器设计标准响应fir2基于窗函数的FIR滤波器设计任意响应firls最小二乘FIR滤波器设计intfilt内插FIR滤波器设计remezParks-McCellan最优FIR滤波器设计remezo

31、rdParks-McCellan最优FIR滤波器阶估计七. 表附1-7 变换函数名功能czt线性调频Z变换dct离散余弦变换（DCT）idct逆离散余弦变换dftmtx离散傅立叶变换矩阵fft一维快速傅立叶变换ifft一维逆快速傅立叶变换fftshift重新排列FFT的输出hilbertHilbert（希尔伯特）变换八. 表附1-8 统计信号处理函数名功能cov协方差矩阵xcov互协方差函数估计corrcoef相关系数矩阵xcorr互相关函数估计cohere相关函数平方幅值估计csd互谱密度（CSD）估计psd信号功率谱密度（PSD）估计tfe从输入输出中估计传递函数九. 表附1-9 窗函数函

32、数名功能boxcar矩形窗triang三角窗bartlettBartlett（巴特利特）窗hammingHamming（哈明）窗hanningHanning（汉宁）窗blackmanBlackman（布莱克曼）窗chebwinChebyshev （切比雪夫）窗kaiserKaiser（凯泽）窗十. 表附1-10 参数化建模函数名功能invfreqs模拟滤波器拟合频率响应invfreqz离散滤波器拟合频率响应prony利用Prony法的离散滤波器拟合时间响应stmcb利用Steiglitz-Mcbride迭代方法求线性模型levinsonLevinson-Durbin递归算法lpc线性预测系数十一. 表附1-11 特殊操作函数名功能rceps实倒谱和最小相位重构cceps倒谱分析和最小相位重构decimate降低序列的取样速率interp提高取样速率（内插）resample改变取样速率medifiltl一维中值滤波deconv反卷积和多项式除法modulate通讯仿真中的调制domod通讯仿真中的解调vco电压控制振荡器specgram频谱分析十二. 表附1-12 模拟原形滤波器设计函数名功能besselapBessel（贝塞尔）模