西安邮电数字信号处理实验报告

1、西安邮电大学通信与信息工程学院实 验 报 告（2016/2017学年 第 2学期） 实验名称： 数字信号处理实验 学生姓名： 石龙飞 专 业： 通信工程 学 号： 03151307 指导教师： 常虹 2017年06月01日西安邮电大学通院数字信号处理实验独立设课实验过程考核表学生姓名石龙飞班级/学号通工1510/03151307承担任务实验室（单位）三号实验楼423所在部门电科实验室实施时间2017年5月11日 2017年6月1日具体内容第一次熟悉MATLAB的使用第二次信号的表示第三次FFT频谱分析及应用第四次信号的运算卷积第五次线性卷积的快熟处理办法重叠相加法第六次线性卷积的快熟处理办法重

2、叠保留法第七次IIR数字滤波器的设计第八次FIR数字滤波器的设计指导教师姓名常虹职务或职称指导教师对学生的评价学习态度 认真 一般 不认真学习纪律 全勤 偶尔缺勤 经常缺勤实践能力 很强 一般 较差 指导教师评语指导教师（师傅）签字 年 月 日 西安邮电大学通院数字信号处理实验独立设课实验成绩鉴定表学生姓名石龙飞专业班级/学号通工1510/03151307315进行时间2017年5月11日 2017年6月1日考核内容标准 平时成绩（30分）课堂出勤情况、课堂纪律及学习态度（25分）遵守实验室的规章制度情况（5分）实验能力（40分）分析、解决实际问题能力（10分）文件检索或软硬件设计能力（10分

3、）实验完成情况（20分）实习报告（30分）语言及文字表达力（10分）实验报告结构（10分）书写格式规范（10分）总评成绩1 熟悉MATLAB1.1 实验目的1、 了解MATLAB程序设计语言的基本特点，熟悉MATLAB软件的运行环境。2、掌握各种信号的建模方式。3、掌握各种信号的托兴表示方法。4、掌握变量等有关概念，具备初步的将一般数学模型转化为对应的计机模型并进行处理的能力。1.2 实验原理1、图形保持hold on/off命令控制是保持原有图形还是刷新原有图形，不带参数的hold命令在两种状态之间进行切换。2、绘图修饰命令title(图形名称)xlabel(x轴说明)ylabel(y轴说明

4、)text(x,y,图形说明)legend1.3 实验内容以及步骤单位采样序列代码：x=-5:5;y1=zeros(1,5),1,zeros(1,5);stem(x,y1,'r')axis(-5,5,0,2);单位阶跃序列代码：x=-5:5;y1=zeros(1,5),1,ones(1,5);stem(x,y1,'r')axis(-5,5,0,2); 矩形序列代码：x=-5:10;y1=zeros(1,5),ones(1,5),zeros(1,6);plot(x,y1)stem(x,y1,'r')axis(-5,10,0,2);三角波代码：fs=

5、10000;t=-1:1/fs:1;w=0.4;x=tripuls(t,w);plot(t,x)方波代码：t=-2*pi:0.001:2*pi;x=square(t);plot(t,x);xlabel('t'),ylabel('z=square(t)')锯齿波代码：x=0:0.01:5;y=1-mod(x,1);plot(x,y)Sinc函数代码：clearn=-5:0.5:10y=sin(0.5*pi*n);stem(n,y,'g')1.4 实验结果分析学习使用MATLAB7.0软件；学习信号的图形表示方法，掌握各种各种信号的的建模方式；实现单

6、位采样序列、单位阶跃序列、举行序列、三角波、方波、锯齿波、Sinc函数。1.5 实验总结以及体会构成经过本次试验，我对于用MATLAB来作出一些数字信号中的一些基本波形的图像有了深刻的认识，可以熟练而又准确地以MATLAB为工具做出诸如单位阶跃序列，矩形序列等多种序列的图形，不仅如此，在作图中，我对于MATLAB的一些作图的基本指令，以及需要作出的序列有了深刻的认识与见解；本节课做出这些图形并不是重点，重点是了解MATLAB程序设计语言的基本特点，掌握这一个重要的工具，并且可以将数学模型熟练地转化为计算机模型，并进行处理，这才是本次实验真正的作用与意义，也是作为当代大学生的我们所必须具备的基本

7、能力。1.6 指导老师评语2 信号的表示2.1 实验目的1、了解MATLAB 程序设计语言的基本特点，熟悉MATLAB 软件运行环境。2、掌握各种信号的建模方式。3、掌握各种信号的图形表示方法。4、掌握变量等有关概念，具备初步的将一般数学模型转化为对应的计算机模型并进行处理的能力2.2 实验原理常用的 MATLAB 绘图语句有figure、plot、subplot、stem 等，图形修饰语具有title、axis、text 等。2.2.1 figure 语句figure 有两种用法。当只有一句figure 命令时，程序会创建一个新的图形窗口，并返回一个整数型的窗口编号。当采用figure(n)

8、时，表示将第n 个图形窗口作为当前的图形窗口，将其显示在所有窗口的最前面。如果该图形窗口不存在，则新建一个窗口，并赋以编号n。2.2.2 plot 语句线形绘图函数。用法为 plot(x,y,s)。参数x 为横轴变量，y 为纵轴变量，s 用以控制图形的基本特征如颜色、粗细等，通常可以省略，常用方法如表1-1 所示。表 1-1 plot 命令的参数及其含义参数含义参数含义参数含义Y黄色.点-实线M紫色O圆;虚线C青色X打叉-.点划线R红色+加号-破折线G绿色*星号向上三角形B蓝色S正方形<向左三角形W白色D菱形>向右三角形K黑色v向下三角形P五角星形表1-12.2.3 subplot

9、 语句subplot(m,n,i)是分割显示图形窗口命令，它把一个图形窗口分为m 行n 列共m×n 个小窗口，并指定第i 个小窗口为当前窗口。2.2.4 二维统计分析图在MATLAB 中，二维统计分析图形很多，常见的有条形图、阶梯图、杆图和填充图等，所采用的函数分别是：bar(x,y,选项)stairs(x,y,选项)stem(x,y,选项)fill(x1,y1,选项1,x2,y2,选项2,)2.2.5 图形保持hold on/off 命令控制是保持原有图形还是刷新原有图形，不带参数的hold 命令在两种状态之间进行切换。2.2.6 绘图修饰命令title(图形名称)xlabel(x

10、 轴说明)ylabel(y 轴说明)text(x,y,图形说明)legend(图例1,图例2,)2.2.7 MATLAB 常用信号生成函数（1） ZEROS 功能：产生全零阵列 调用格式：X=ZEROS (N) %产生N 行N 列的全零矩阵 X=ZEROS(M,N) %产生M 行N 列的全零矩阵（2） ONES 功能：产生全 1 阵列 调用格式：X=ONES(N) %产生N 行N 列的全1 矩阵 X=ONES(M,N) %产生M 行N 列的全1 矩阵（3） SINC 功能：辛格函数 调用格式：Y= SINC(X) %（4） RECTPULS 功能：产生矩形脉冲信号 调用格式：Y= RECTPU

11、LS (T) %产生高度为1、宽度为1、关于T0 对称的矩形脉冲 Y= RECTPULS (T，W) %产生高度为1、宽度为W、关于T0 对称的矩形脉冲（5） RAND 功能：产生伪随机序列 调用格式：Y= RAND (1,N) %产生0,1上均匀分布的随机序列 Y= RANDN (1,N) %产生均值为0，方差为1 的白噪声序列（6） SAWTOOTH 功能：产生周期锯齿波或三角波 调用格式：Y= SAWTOOTH (T) %产生幅值为1，1，以2 为周期的方波 Y=SAWTOOTH(T,WIDTH)% 产生幅值为1，1，以WIDTH *2 为周期的方波（7） SQUARE 功能：产生方波

12、调用格式：Y= SQUARE (T) %产生幅值为1，1，以2 为周期的锯齿波 Y= SQUARE(T,DUTY) % 产生幅值为1，1，以占空比为DUTY 的方波 例： t = 0:.0001:.0625; y = SQUARE(2*pi*30*t，80); plot(t,y)产生一个占空比为80%的方波（8） FLIPLR 功能： 序列左右翻转 调用格式：Y = FLIPLR(X) % X = 1 2 3 翻转后 3 2 14 5 6 6 5 4（9） CUMSUM 、SUM 功能： 计算序列累加 调用格式：Y = CUMSUM(X) % 向量X 元素累加，记录每一次的累加结果,而SUM只

13、记录最后的结果。2.3：实验内容以及步骤1、分别以条形图、阶梯图、杆图和填充图形式绘制曲线y=2sin(x)。程序如下：x=0:pi/10:2*pi;y=2*sin(x);subplot(2,2,1);bar(x,y,'g');title('bar(x,y,''g'')');axis(0,7,-2,2);subplot(2,2,2);stairs(x,y,'b');title('stairs(x,y,''b'')');axis(0,7,-2,2);subplot(2,

14、2,3);stem(x,y,'k');title('stem(x,y,''k'')');axis(0,7,-2,2);subplot(2,2,4);fill(x,y,'y');title('fill(x,y,''y'')');axis(0,7,-2,2);2、采用图形保持，在同一坐标内绘制曲线y1=0.2e-0.5xcos(4x) 和y2=2e-0.5xcos(x)。程序如下：x=0:pi/100:2*pi;y1=0.2*exp(-0.5*x).*cos(4*pi*x

15、);plot(x,y1)hold ony2=2*exp(-0.5*x).*cos(pi*x);plot(x,y2);hold off3、在0x2 区间内，绘制曲线y1=2e-0.5x 和y2=cos(4x)，并给图形添加图形标注。程序如下：x=0:pi/100:2*pi;y1=2*exp(-0.5*x);y2=cos(4*pi*x);plot(x,y1,x,y2)title('x from 0 to 2pi'); %加图形标题xlabel('Variable X'); %加X 轴说明ylabel('Variable Y'); %加Y 轴说明tex

16、t(0.8,1.5,'曲线y1=2e-0.5x'); %在指定位置添加图形说明text(2.5,1.1,'曲线y2=cos(4pix)');legend(y1, y2) %加图例4、实现单位采样序列 (n)、单位阶跃序列u(n)、矩形序列 ( ) N R n ，并用图形显示。写出程序及输出图形%单位采样序列 (n)n=-5:10;y=zeros(1,5),1,zeros(1,10);stem(n,y);%单位阶跃序列u(n)axis(-5,10,0,2);n=-5:10;y=zeros(1,5),ones(1,11);stem(n,y,'r')a

17、xis(-5,10,0,2);%矩形序列 ( ) N R n n=-5:10;y=zeros(1,5),ones(1,5),zeros(1,6);plot(n,y)stem(n,y)axis(-5,10,0,2);title('矩形序列');5、实现三角波、方波、锯齿波、Sinc 函数，并用图形显示。写出程序及输出图形%三角波fs=10000;t=-1:1/fs:1;w=0.4;x=tripuls(t,w);plot(t,x)title('三角波');%方波t=-2*pi:0.001:2*pi;x=square(t);plot(t,x);xlabel('

18、t'),ylabel('x=square(t)');%锯齿波x= 0:0.01:5;y=1-mod(x,1);plot(x,y)title('锯齿波')%Sinc函数n=-5:0.5:10;y=sin(0.5*pi*n);stem(n,y,'g')title('正弦序列');2.4：实验结果分析1、 分别以条形图、阶梯图、杆图和填充图形式绘制曲线y=2sin(x)。2、 采用图形保持，在同一坐标内绘制曲线y1=0.2e-0.5xcos(4x) 和y2=2e-0.5xcos(x)。3、 在0x2 区间内，绘制曲线y1=2e-

19、0.5x 和y2=cos(4x)，并给图形添加图标注。4、实现单位采样序列 (n)、单位阶跃序列u(n)、矩形序列 ( ) N R n ，并用图形显示，写出程序及输出图形。5、 实现三角波、方波、锯齿波、Sinc 函数，并用图形显示。写出程序及输出图形2.5：实验总结以及体会构成本次实验主要是熟悉了如何使用MATLAB表示信号，通过MATLAB的函数来实现是主要方法，其次也可以自行编写函数，表示自己需要的特定信号，来适应不同的应用场合。在信号表示中学会画图是非常重要的一件事，几乎所有的信号都可以也需要用图像表示出来，这样对分析和实验都是有着好处的，熟练掌握画图函数，这些便都可以迎刃而解，画图函

20、数有很多种，其参数也是形形色色，记住每个参数的位置意义是非常重要的。因为是对数字信号的处理，因此免不了会用到很多的理论知识，这时就是理论和时间结合的时候，MATLAB是一个可以仿真信号的工具，可以模拟实际，需要用理论知识来设计解决问题所需要的模型，然后用MATLAB将其仿真出来，看看实际结果是否是我们预期的，这个结果也将是我们在现实中实施该方案所能看到的结果。2.6：指导老师评语3 FFT 频谱分析及应用3.1 实验目的1、通过实验加深对FFT 的理解；2、熟悉应用FFT 对典型信号进行频谱分析的方法。3.1 实验原理在各种信号序列中，有限长序列占重要地位。对有限长序列可以利用离散傅立叶变换(

21、DFT)进行分析。DFT 不但可以很好的反映序列的频谱特性，而且易于用快速算法(FFT)在计算机上进行分析。有限长序列的 DFT 是其z 变换在单位圆上的等距离采样，或者说是序列傅立叶的等距离采样，因此可以用于序列的谱分析。FFT 是DFT 的一种快速算法，它是对变换式进行一次次分解，使其成为若干小数据点的组合，从而减少运算量。在 MATLAB 信号处理工具箱中的函数fft（x,n）,可以用来实现序列的N 点快速傅立叶变换。经函数fft求得的序列一般是复序列，通常要求出其幅值和相位。MATLAB中提供了求复数的幅值和相位的函数：abs、angle，这些函数一般和fft同时使用。3.3.3 实验

22、内容以及步骤被噪声污染得信号，比较难看出所包含得频率分量，如一个由50Hz 和120Hz正弦信号构成的信号，受零均值随机噪声的干扰，数据采样率为1000Hz，试用FFT 函数来分析其信号频率成分，要求：1 画出时域波形；2 分析信号功率谱密度。N = 40;n = 0:N-1;t = 0.01*n;x = 2*sin(4*pi*t) + 5*cos(8*pi*t);k = 0:N/2;w = 2*pi/N*k;X = fft(x, N);magX = abs(X(1:N/2+1);subplot(2, 1, 1);stem(n, x, '.');title('signa

23、l x(n)');subplot(2, 1, 2);stem(magX);title('FFT N=40');xlabel('f (unit : pi)'); ylabel('|X|');grid3.3.4 实验结果分析3.3.5实验总结以及体会构成本次实验只要是针对有限长序列，我们可以通过MATLAB对其进行离散傅里叶变化，是非常方便的，因为工具内部已经将函数写好，我们只需要将变化的信号交给函数就可以，可以说大大地减小了手工计算量，并且能过直观得将图画出来。其中需要注意的问题也有很多，比如因为操作失误，将函数名称写错，软件机会报错，这个

24、问题还是比较容易查找的，其次是在表示信号的时候，将一些数值写错，到会最后的变换不是我们需要的，画出的图像有明显的问题，这些都需要细心才能避免，当然还有一些比较头疼的错误，比如函数不会使用，所给的参数错误，或者不符合标准。3.3.6指导老师评语4 信号的运算-卷积4.1 实验目的1、掌握信号的线性卷积运算。2、掌握信号的循环卷积运算。3、掌握信号循环卷积计算线性卷积的条件。4.2 实验原理数字信号卷积定理4.3 实验内容以及步骤1、假设卷积下面信号X(n)= 0n<130 elseh(n)=1 0n<120 else选定循环卷积的长度为N = 21。确定 1y (n) = x(n)h

25、(n)的哪些数值与线性卷积 2y (n) = x(n) *h(n)结果中的数值相同。编写程序代码并输出图形，并分析错误数据的原因，怎样才能使两者数据相同。n1=0:1:12;x1=0.9.n1;h=ones(1,12);h=0 h;N=length(x1)+length(h)-1;n=0:N-1;ny=0:20;y1=circonvt(x1,h,21); y2=circonvt(x1,h,N);x1=x1 zeros(1,N-length(x1);h=h zeros(1,N-length(h);X1= fft(x1,N);H=fft(h,N);X=X1.*H;x=ifft(X);x=real(

26、x);subplot(2,2,1);stem(n,x1);title('x1(n)');axis(0,33,0,1);subplot(2,2,2);stem(n,h);title('h(n)');axis(0,33,0,1);subplot(2,2,3);stem(ny,y1,'fill');title('21点循环卷积');axis(0,33,0,8);hold on;subplot(2,2,4);stem(n,x);title('线性卷积');axis(0,33,0,8);subplot(2,2,3);stem

27、(n,x,'r','-');axis(0,33,0,8);hold off2、已知系统响应为h(n) = sin(0.2n) + cos(0.5n) 0 n < 20，输入为x(n) = exp(0.2n) 0 n <10 ，画出用DFT 方法实现求系统输出的系统框图，编写用DFT 实现的程序代码并输出图形function y = circonvt(x1,x2,N)if (length(x1)>N | length(x2)>N) error('N必须大于等于x的长度');endx1=x1 zeros(1,N-length(

28、x1);x2=x2 zeros(1,N-length(x2);X1=fft(x1,N);X2=fft(x2,N);X=X1.*X2;y = ifft(X,N);y=real(y);n1=0:1:19;hn=sin(0.2*n1)+cos(0.5*n1);n2=0:1:9;xn=exp(0.2*n2);N=length(xn)+length(hn)-1;n=0:N-1;ny=0:20;y1=circonvt(xn,hn,21); y2=circonvt(xn,hn,N);xn=xn zeros(1,N-length(xn);hn=hn zeros(1,N-length(hn);X1= fft(x

29、n,N);H=fft(hn,N);X=X1.*H;x=ifft(X);x=real(x);subplot(2,2,1);stem(n,xn);title('xn(n)');axis(0,33,0,5);subplot(2,2,2);stem(n,hn);title('h(n)');axis(0,33,-2,2);subplot(2,2,3);stem(ny,y1,'fill');title('21点循环卷积');axis(0,33,-20,30);hold on;subplot(2,2,4);stem(n,x);title(

30、9;线性卷积');axis(0,33,-20,30);subplot(2,2,3);stem(n,x,'r','-');axis(0,33,-20,30);hold off4.4 实验结果分析1、假设卷积下面信号X(n)= 0n<130 elseh(n)=1 0n<120 else选定循环卷积的长度为N = 21。确定 1y (n) = x(n)h(n)的哪些数值与线性卷积 2y (n) = x(n) *h(n)结果中的数值相同。编写程序代码并输出图形，并分析错误数据的原因，怎样才能使两者数据相同。2、已知系统响应为h(n) = sin(0.

31、2n) + cos(0.5n) 0 n < 20，输入为x(n) = exp(0.2n) 0 n <10 ，画出用DFT 方法实现求系统输出的系统框图，编写用DFT 实现的程序代码并输出图形4.5 实验总结以及体会构成本次试验只要使用了卷积函数来计算卷积，此处的卷积函数是自己来写的，学要用到所学到的知识，从最基本的计算设计卷积函数，需要考虑到使用时的参数，返回值之类的，一个何时何地都可以用的函数。MATLAB的函数编写也还需要注意几个问题，比如函数的命名是要与m文件一致的，不然无法调用，系统会提示此函数没有定义，或者其他更奇怪的错误。这里也学到了一种查错的方法，就是添加断点，可以跟

32、踪每一个变量，了解他在每个地方的值，这样根据理论值，很容易知道是哪里出了问题，对于写代码和代码查错是十分有帮助的。4.6 指导老师评语5 信号的块处理算法-重叠相加法5.1 实验目的1、掌握信号的块运算重叠相加法原理。2、掌握信号的块运算重叠相加法的实现方法。3、通过实验和分析深入了解重叠相加法。5.2 实验原理重叠相加法把长的输入信号分成小的非重叠段。如果这些段的长度是M ，冲激响应的长度是Lh ,那么一个长度为N > M + Lh 1的循环卷积将会避免所有的时间混叠效应（通过使用填补零）。然而，每一段输出的长度现在大于M 。要将输出的各段合起来，必须把每一段的重叠部分加起来。因此，这

33、个方法也可以称为“重叠输出”法。这里讲解较简略，更详细的论述请在相关教科书中查找。5.3 实验内容以及步骤编写一个实现块卷积重叠相加法的M 文件函数。该函数的输入之一应是段的长度M，或者是FFT 的长度。循环卷积应在DFT 域中进行。最终，这将会使程序最快的运行。function y = test4_1(x, h, Nfft)H=fft(h,Nfft);M=Nfft - length(h)+1; for ix = 1:M:length(x) x_seg = x(ix:ix+M-1); X = fft(x_seg,Nfft); Y = X.*H; y_seg = ifft(Y); y(ix:ix

34、+Nfft-1)= y(ix:ix+Nfft-1) +y_seg(1: Nfft);endy = ifft(Y, Nfft);if any(imag(h) && any(imag(x) y=real(y);endn1 = 0:1:19;hn = sin(0.2*n1) + cos(0.5*n1);n2 = 0:1:9;xn = exp(0.2*n2);z = conv(hn, xn);N = length(xn)+length(hn)-1;hn = hn zeros(1, N-length(hn);R = length(xn);M = length(hn);%y = test4

35、_1(xn, hn, 8);y = fftfilt(xn, hn, 8);k1 = 1:length(y);k = 1:N;plot(k1, y, 'b-o', k, z, 'g-x')xlabel('Time index n');ylabel('Amplitude')legend('yn', 'zn');5.4 实验结果分析5.5 实验总结以及体会构成最基本的卷积运算时所有卷积的基础，是他还远不够实际中的需要，在实际生产生活中有时不仅仅是计算出结果就行，设计者还要考虑这个算法效益，如果数据庞大，计

36、算一个卷积需要几分钟甚至以小时为单位，那么这个算法基本也没什么意义，所以我们还要考虑时间的开销是否合理，于是优化算法就显得尤为重要，将复杂的问题以最短的时间计算出来。在这里我们优化了卷积算法，即重叠相加法，可以比之前的算法更快的完成卷积，为我们节省时间。5.6 指导老师评语6 线性卷积的快速处理方法重叠保留法6.1 实验目的1、掌握信号的块运算重叠保留法原理。 2、掌握信号的块运算重叠保留法的实现方法。3、通过实验和分析深入了解重叠保留法。6.2 实验原理学习使用 MATLAB2014 软件进行建模。学习信号的重叠保留法的基本原理。学习重叠保留法的 MA

37、TLAB 实现。实现信号的重叠保留法运算；通过与卷积运算进行对比，比较其浮点运算次数及消耗的时间。6.3 实验内容以及步骤通过conv函数来测试此函数的正确性。考虑使用下面信号xn=1 0n<170 else hn=sinn13 0n<1000 else% circonvt函数function y = circonvt(x1,x2,N)if (length(x1)>N | length(x2)>N) error('N必须大于等于x的长度');endx1=x1 zeros(1,N-length(x1);x2=x2 zeros(1,N-length

38、(x2);X1=fft(x1,N);X2=fft(x2,N);X=X1.*X2;y = ifft(X,N);y = real(y); % ovrlpsav函数function y = ovrlpsav(x,h,N)Lenx=length(x);M=length(h);if N<M N=M+2;endM1=M-1;L=N-M1;h=h,zeros(1,N-M);x=zeros(1,M1),x,zeros(1,N-1);K=floor(Lenx+M-1)/(L);Y=zeros(K+1,N);for k = 0:K xk = x (k*L+1:k*L+N); Y(k+1,:)=circonv

39、t(xk,h,N);endY=Y(:,M:N)'y=(Y(:)'endxn=ones(1,17);n1=0:1:100; hn=sin(n1*pi/13); z=conv(hn,xn); N=length(xn)+length(hn)-1; hn=hn zeros(1,N-length(hn); R=length(xn); M=length(hn); y=ovrlpsav(xn,hn,10); k1=1:length(y);k=1:N; plot(k1,y,'b-o',k,z,'g-x'); title('x(n)和h(n)的循环卷积图像

40、'); xlabel('n');ylabel('y');6.4 实验结果分析6.5 实验总结以及体会构成最基本的卷积运算时所有卷积的基础，是他还远不够实际中的需要，在实际生产生活中有时不仅仅是计算出结果就行，设计者还要考虑这个算法效益，如果数据庞大，计算一个卷积需要几分钟甚至以小时为单位，那么这个算法基本也没什么意义，所以我们还要考虑时间的开销是否合理，于是优化算法就显得尤为重要，将复杂的问题以最短的时间计算出来。在这里我们优化了卷积算法，即重叠相加法，可以比之前的算法更快的完成卷积，为我们节省时间。6.6 指导老师评语7 IIR数字滤波器的设计7.1

41、实验目的1、掌握脉冲响应不变法和双线性变换法设计IIR数字滤波器的原理和方法。2、观察双线性变换法和脉冲响应不变法设计的滤波器的频率特性，了解双线性变换法和脉冲响应不变法的特点和区别。7.2 实验原理7.2.1 脉冲响应不变法所谓脉冲响应不变法就是使数字滤波器的单位脉冲响应序列h(n)等于模拟滤波器的单位冲激响应和(t)的采样值，即：，其中，T为采样周期。在MATLAB中，可用函数impinvar实现从模拟滤波器到数字滤波器的脉冲响应不变映射，调用格式为：b,a=impinvar(c,d,fs)b,a=impinvar(c,d)其中，c、d分别为模拟滤波器的分子和分母多项式系数向量；fs为采样

42、频率（Hz），缺省值fs=1Hz；b、a分别为数字滤波器分子和分母多项式系数向量。7.2.2 双线性变换法：由于s平面和z平面的单值双线性映射关系为s，其中T为采样周期。因此，若已知模拟滤波器的传递函数，将上式代入即可得到数字滤波器的系统函数H(z)。在双线性变换中，模拟角频率和数字角频率的变换关系为： 可见，和w之间的变换关系为非线性的。在MATLAB中，可用函数bilinear实现从模拟滤波器到数字滤波器的双线性变换映射，调用格式为：b,a=bilinear(c,d,fs)7.2.3 设计步骤：(1) 定技术指标转换为模拟低通原型设计性能指标。(2) 估计满足性能指标的模拟低通性能阶数和截

43、止频率。利用MATLAB中buttord、cheb1ord、cheb2ord、ellipord等函数，调用格式如：其中，Wp为带通边界频率，rad/s；Ws为阻带边界频率，rad/s；Rp为带通波动，dB；Rs为阻带衰减，dB；s表示为模拟滤波器；函数返回值n为模拟滤波器的最小阶数；Wn为模拟滤波器的截止频率(-3dB频率)，rad/s。函数适用低通、高通、带通、带阻滤波器。(3) 设计模拟低通原型。利用MATLAB中buttap、cheb1ap、cheb2ap、elliap等函数，调用格式如z,p,k=buttap(n)。采用上述函数所得到原型滤波器的传递函数为零点、极点、增益表达式，需要和

44、函数c,d=zp2tf(z,p,k)配合使用，以转化为多项式形式。(4) 由模拟低通原型经频率变换获得模拟低通、高通、带通或带阻滤波器。利用MATLAB中lp2lp、lp2hp、lp2bp、lp2bs等函数，调用格式如c1,d1=lp2lp(c,d,Wn)。(5) 利用脉冲响应不变法或双线性不变法，实现模拟滤波器到数字滤波器的映射。本实验用到的特殊函数：db,mag,pha,grd,w=freqz_m(b,a)，计算幅频和相频响应。7.3 实验内容以及步骤用双线性变换法设计一个Chebyshev1型数字带通滤波器，设计指标为T=1ms,Rp=1dB,Wp1=0.35，Wp2=0.65,Rs=6

45、0dB,Ws1=0.2Ws2=0.8。按实验步骤附上所设计滤波器的H（z）及相应的幅频特性曲线定性分析得到的图形，判断设计是否满足要求。Ts=0.001;Fs=1/Ts;Rp=1;Rs=60;wp1=0.35*pi;wp2=0.60*pi;ws1=0.2*pi;ws2=0.8*pi;Wp1=(2/Ts)*tan(wp1/2);Wp2=(2/Ts)*tan(wp2/2);Wp=Wp1,Wp2;Ws1=(2/Ts)*tan(ws1/2);Ws2=(2/Ts)*tan(ws2/2);Ws=Ws1,Ws2;BW=Wp2-Wp1;Omegaw0=sqrt(Wp1*Wp2);N,OmegaC=cheb1o

46、rd(Wp,Ws,Rp,Rs,'s');z0,p0,k0=cheb1ap(N,Rp);AnalogB=k0*real(poly(z0);AnalogA=real(poly(p0);BandB,BandA=lp2bp(AnalogB,AnalogA,Omegaw0,BW);DigitalB,DigitalA=bilinear(BandB,BandA,Fs);sos,G=tf2sos(DigitalB,DigitalA);Hz,Wz=freqz(DigitalB,DigitalA,1024,'whole');dbHz=20*log10(abs(Hz)+eps)/max(abs(Hz);grd=grpdelay(DigitalB,DigitalA,Wz);subplot(2,3,1);p