数字信号处理指导书

实验三 用双线性变换法设计IIR 数字滤波器一、实验目的：1、熟悉用双线性变换法设计IIR 数字滤波器的原理与方法。

2、掌握数字滤波器的计算机仿真方法。

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

二、实验内容：1、 用双线性变换法设计一个巴特沃斯低通IIR 数字滤波器。

设计指标参数为：在通带内频率低于0.2π时，最大衰减小于1dB ；在阻带内[0.3π，π]频率区间上，最小衰减大于15dB 。

2、 以0.02π为采样间隔，打印出数字滤波器在频率区间[0，2π]上的幅频响应特性曲线。

3、 用所设计的滤波器对实际心电图信号采样序列(在本实验后面给出)进行仿真滤波处理，并分别打印出滤波前后的心电图信号波形图，观察总结滤波作用与效果。

三、实验步骤：1、 复习有关巴特沃斯模拟滤波器设计和用双线性变换法设计IIR 数字滤波器的内容，用双线性变换法设计数字滤波器系统函数)(z H 。

教材中满足本实验要求的数字滤波器函数()()()()()()51121212310.00737811 1.26860.70511 1.0100.35831 1.90440.2155k k z H z zz z z z z H z -------=+=-+-+-+=∏式中：()()2121121------++=CkzBkz z z A z H k ，k=1，2，3 A=0.09036，B 1=1.2686，C 1=-0.7051 B 2=1.0106， C 2=-0.3583 B 3=0.9044， C 3=-0.2155根据设计指标，调用MATLAB 信号处理工具函数buttord 和butter ，也可得到()z H 。

由滤波器的函数可见，滤波器()z H 由三个二阶滤波器()z H 1、()z H 2、()z H 3级联组成。

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

()n y k 为第k 阶滤波器()z H k 的输出序列，()n y k 1-为输入序列，根据滤波器()z H 的组成可得差分方程()()()()()()21212111-+-+---+=---n y C n y B n Ay n Ay n Ay n y k k k k k k k k 当1=k 时，()()n x n y k =-1。

所以()z H 对()n x 的总响序列()n y 可以用顺序迭代算法得到。

即依次对k=1，2，3求解差分方程，最后得到()()n y n y =3。

仿真程序就是实现上述求解差分方程和顺序迭代算法的通用程序，也可直接调用MATLAB filter 函数实现仿真。

3、 在通用计算机上运行仿真滤波程序，完成实验内容（2）和（3）。

四、实验报告要求：1、 简述实验目的及原理。

2、 由所打印的)(ωj e H 特性曲线及设计过程简述双线性变换法的特点。

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

下面给出一实际心电图信号采样序列样本)(n x ，中存在高频干扰。

在实验中，以)(n x 作为输入序列，滤除其中的干扰成分。

}{⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧---=,0,0,12,38,4,4)(n x ,2,0,12,60,4,2---- ,22,10,84,6,0---- ,2,4,6,90,6,4----- ,2,0,6,66,2,6---}0,0,6,32,6,4--,0,4,4,12,2--,2,0,2,8,4---,2,0,4,0,6---,0,0,8,16,6--实验四 用窗函数法设计FIR 数字滤波器一、实验目的：1、掌握用窗函数法设计FIR 数字滤波器的原理和方法。

2、熟悉线性相位FIR 数字滤波器特性。

3、了解各种窗函数对滤波特性的影响。

二、实验内容及步骤：1、复习用窗函数法设计FIR 数字滤波器一节内容。

一般，设计线性相位FIR 数字滤波器采用窗函数法或频率抽样法，本实验采用窗函数法，采用矩形窗。

下面简要介绍窗函数法设计FIR 理论知识，更详细的相关知识请参考教材。

如果所希望的滤波器的理想频率响应函数为)(jw d e H ，如理想的低通，由信号系统的知识知道，在时域系统的冲击响应h d (n)将是无限长的，如图2、图3所示。

H d (w)-w c w c图2图3若时域响应是无限长的，则不可能实现，因此需要对其截断，即设计一个FIR 滤波器频率响应∑-=-=1)()(N n jwn jwe n h e H 来逼近)(jw d e H ，即用一个窗函数w(n)来截断h d (n)，如式3所示：)()()(n w n h n h d =（式3）。

最简单的截断方法是矩形窗，实际操作中，直接取h d (n)的主要数据即可。

)(n h 作为实际设计的FIR 数字滤波器的单位脉冲响应序列，其频率响应函数为：∑-=-=1)()(N n jwn jwe n h e H （式4）令jwe z =，则∑-=-=10)()(N n n z n h z H （式5），式中，N 为所选窗函数)(n w 的长度。

如果要求线性相位特性，)(n h 还必须满足：)1()(n N h n h --±=（式6），根据式6中的正、负和长度N 的奇偶性又将线性相位FIR滤波器分成四类。

要根据所设计的滤波器特性正确选择其中一类。

例如：要设计线性相位低通特性，可选择)1()(n N h n h --=类。

有关窗函数的种类及特性，请参照教材。

2、编写程序1）编写能产生矩形窗、升余弦窗、改进升余弦窗和二阶升余弦窗的窗函数子程序，也可调用MATLAB 相关函数。

2）编写低通FIR 数字滤波器主程序。

3）上机实验内容。

A 、用升余弦窗设计一线性相位低通FIR 数字滤波器，截止频率rad c 4πω=。

窗口长度N ＝15,33。

要求在两种窗口长度情况下，分别求出h (n)，打印出相应的幅频特性和相频特性曲线，观察3dB 带宽和20dB 带宽。

总结窗口长度N 对滤波特性的影响。

B 、n ＝33,c ω＝4π，用四种窗函数设计线性相位低通滤波器。

绘制相应的幅频特性曲线，观察3dB 和20 dB 带宽以及阻带最小衰减，比较四种窗函数对滤波器特性的影响。

三、思考题：1、 如果给定通带截止频率和阻带截止频率以及阻带最小衰减，如何用窗函数法设计线性相位低通滤波器?写出设计步骤。

2、 如果要求用窗函数法设计带通滤波器，且给定上、下边带截止频率为1ω和2ω，试求理想带通的单位脉冲响应)(n h d 。

四、实验报告要求：1、 简述实验目的及原理。

2、 按照实验步骤及要求，比较各种情况下的滤波性能，说明窗口长度N 和窗函数类型对滤波特性的影响。

3、 总结用窗函数法设计FIR 滤波器的主要特点。

4、 简要回答思考题。

参考：用于数字信号处理的MATLAB 操作或者函数1、向量的生成：利用冒号（：）生成向量。

1）X=j:k (j<k)生成X=[j,j+1,j+2,…,k-1,k]2）X=j:I:k如果I>0 且j<k或I<0且j>k,则生成向量X=[j,j+I,j+2I,…,k-I,k]如果I<0 且j>k,则生成向量X=[j,j+I,j+2I,…,k]例： X1=1：5 X1=[1 2 3 4 5]X2=1：0.5：3 X2=[1.000 1.500 2.000 2.500 3.000] X3=5：-1：1 X3=[5 4 3 2 1]2、矩阵的生成：Zeros 生成全0阵B= Zeros(n) 生成n*n的全0矩阵B= Zeros(m,n) 生成m*n的全0矩阵B= Zeros(Size(A)) 生成与与矩阵A大小相同的全0矩阵注：m,n必须是一个非负数ones生成全1阵rand生成均匀分布的随机阵3、矩阵的算术运算1）加法和减法对于同维矩阵指令为：A+B、A-B；对于矩阵和标量（一个数）的加减运算，指令为：A+3、A-9 2）乘法和除法运算A*B 是数学中的矩阵乘法，遵循矩阵乘法规则；A．*B 是同维矩阵对应位置元素做乘法；B=inv（A）是求矩阵的逆；A/B 是数学中的矩阵除法，遵循矩阵除法规则；A．/B 是同维矩阵对应位置元素相除；'A表示矩阵的转置运算；4、数组函数下面列举一些基本函数，他们的用法和格式都相同。

Sin(A),cos(A),exp(A),log(A)(相当于ln)Sgrt(A)开平方 abs(A)求模 real(A)求实部 imag(A)求虚部式中A可以是标量也可以是矩阵例：利用等差向量产生一个正弦值向量t=0:0.1:10A=sin(t)PLot(A)这时候即可看到一个绘有正弦曲线的窗口弹出。

另：每条语句后面加“；”表示不要显示当前语句的执行结果；不加“；”表示要显示当前语句的执行结果。

5、绘图（二维绘图）利用前例如果要使向量的横纵坐标一一对应，则应写为plot(t,A)。

如不要横坐标对应时，可只写plot(A)。

如果A是一个矩阵不是向量，则在同一窗口中可绘出与矩阵行数相同的曲线且颜色不同。

同时也可指定曲线的颜色和格式。

如plot(t,A,’b.’)则原来的连续曲线就变成了蓝色的点线。

下面给出常用选项。

r红 g绿 b蓝 y 黄 m 洋红 c青 w白 k黑━实线┅┅虚线*用星号标数据点。

用点标 o（字母‘o‘）用圆圈标 x（字母’x‘）用叉标-.点划线下面给出窗口分割语句 subplot(x,y,n)，表示将一个窗口分割成x行，y列，共x*y个窗口，取第n个窗口为当前窗口。

例subplot(3,2,2) 分割3*2个窗口，取第二个窗口，如果下面有绘图语句，就表示要在第二个窗口中绘图。

6、联机帮助1)Help 准确指令，如Help elfunHelp expHelp ftt2)look for —有关指令，如Look for fourier本实验可能用到的一些序列和函数1）单位冲击序列()n N δ-：function [x,n]=impseq(n0,n1,n2) n=[n1:n2];x=[(n-n0)==0];利用逻辑运算实现在n=N 处为1。

如：impseq(5,0,10) 2）单位阶跃序列()u n N -）：function [x,n]=stepseq(n0,n1,n2) n=[n1:n2]x=[(n-n0)>=0]; 利用逻辑运算实现在n >=N 处为1。

tepseq(5,0,10) 3）矩形序列：function [x,n]=stepseq(n0,n1,n2) n=[n1:n2]x=[(n-n0)>=0] stepseq(0,0,9)实现了()10R n也可以用function [x,n]=squseq(M ，N,n1,n2) n=[n1:n2]x=[(n-M)>=0&(n-N)<0] 利用逻辑运算实现在N>n >=M 处为1 4)正弦序列()()sin *x n n ω=x(n)=sin(n*pi) x(t)=sin(2*t)5)也可调用库函数产生一些信号如sinc 6）直接输入序列或者函数表达式 x(t)=sin(2*t)x(n)=[1,1,1,1,1,1]3、掌握二维绘图方法 1）plot----连续曲线 t=0:pi/50pi; y=sin(t); plot(t,y,)grid on 保留前面坐标系中已经存在的图形plot(t,y,s1,s2,s3) s1,s2,s3表示线型、颜色等p133 2）stem 离散图形[n,y]=squseq(4，4,0,10); Stem(n,y,o)3)一个窗口中绘制多个图形 Subplot(2,3,1) Plot(x)Subplot(2,3,2) Plot(y)滤波器设计中的可能要调用的函数见下面滤波器的设计中数字巴特沃兹滤波器的设计方法方法一%获得模拟滤波器的性能指标，边界频率进行预畸变校正，按公式2tan2Tω⎛⎫Ω= ⎪⎝⎭进Wp=0.2*piWs=0.3*pi;Rp=1;%单位DBRs=15; %单位DBFs=0.02*pi; %单位Hz%模拟滤波器设计[n,Wn]=buttord(Wp,Ws,Rp,Rs);%获得阶数和通带截止频率[z,p,k]=buttap(n);%设计模拟低通滤波器原型[b0,a0]=zp2tf(z,p,k);%转化为传递函数表示的滤波器[b,a]=lp2lp(b0,a0,Wn);%将滤波器原形转化为要求设计的低通滤波器%使用双线形变换法得到设计的数字滤波器[bz,az]=bilinear(b,a,Fs);[h,w]=freqz(bz,az);%计算频响绘频率响应曲线plot(h);方法二1，数字滤波器的性能指标归一Wp=0.2Ws=0.3Rp=1;%单位DBRs=15; %单位DB2，获得滤波器的阶数和通带截止频率[n,Wn]=buttord(Wp,Ws,Rp,Rs)设计巴特沃兹低通数字滤波器[b,a]＝butter(n,Wn)Freqz=(b,a)方法三见本实验附录程序A=[0.09036,2*0.09036,0.09036];B1=[1,-1.2686,0.7051];B2=[1,-1.0106,0.3583];B3=[1,-0.9044,0.2155];[H1,w]=frez(A,B1,100);[H2,w]=frez(A,B2,100);[H3,w]=frez(A,B3,100);H=H1.*H2.*H3;mag=abs9h);db=20log((mag+eps)/max(mag));subplot(2,2,1)plot(w/pi,db)DFT的源程序functionX=dft(x)N=length(x);n=0:N-1;k=0:N-1;WN=exp(-j*2*pi/N);%求得WN因子%Nk=n`*k;WNnk=WN.^nk;%求得DFT运算矩阵，注意向量乘为点乘% X=x*WNnk%利用矩阵运算计算DFT%IDFT源程序。