《数字信号处理》实验指导书学生版

1、数字信号处理实验指导书赣南师范学院物理与电子信息学院实验指导书 20152016学年 第 1学期学院物理与电子信息学院教研室信息与通信工程教研室课程名称数字信号处理授课班级13电信本、13电信卓越主讲教师黄隆胜职称副教授2014年03月目 录前 言1实验一 熟悉MATLAB环境4实验二 用MATLAB进行离散系统的Z域分析6实验三 傅立叶变换10实验四 IIR及FIR滤波器的MATLAB实现13前 言MATLAB是由美国Math Works公司推出的软件产品。MATLAB是“Matrix Laboratory”的缩写，意及“矩阵实验室”。MATLAB是一完整的并可扩展的计算机环境，是一种进行科

2、学和工程计算的交互式程序语言。它的基本数据单元是不需要指定维数的矩阵，它可直接用于表达数学的算式和技术概念，而普通的高级语言只能对一个个具体的数据单元进行操作。因此，解决同样的数值计算问题，使用MATLAB要比使用Basic、Fortran和C语言等提高效率许多倍。许多人赞誉它为万能的数学“演算纸”。MATLAB采用开放式的环境，你可以读到它的算法，并能改变当前的函数或增添你自己编写的函数。在欧美的大学和研究机构中，MATLAB是一种非常流行的计算机语言，许多重要的学术刊物上发表的论文均是用MATLAB来分析计算以及绘制出各种图形。它还是一种有利的教学工具，它在大学的线性代数课程以及其它领域的

3、高一级课程的教学中，已成为标准的教学工具。最初的MATLAB是用FORTRAN编写的，在DOS环境下运行。新版的MATLAB 是C语言编写的高度集成系统。它在几乎所有流行的计算机机种，诸如PC、MACINTOSH、SUN、VAX上都有相应的MATLAB版本。新版的MATLAB增强了图形处理功能，并在WINDOWS环境下运行。现今，MATLAB的发展已大大超出了“矩阵实验室”的范围，在许多国际一流专家学者的支持下，Maths Works公司还为MATLAB配备了涉及到自动控制、信息处理、计算机仿真等种类繁多的工具箱（Tool  Box），这些工具箱有数理统计、信号处理、系统辨识、最优化

4、、稳健等等。近年来一些新兴的学科方向，Maths Works公司也很快地开发了相应的工具箱，例如：神经网络、模糊逻辑等。本课程实验要求学生运用MATLAB编程完成一些数字信号处理的基本功能，加深对教学内容的理解。实验一 熟悉MATLAB环境实验目的1熟悉MATLAB主界面，并学会简单的菜单操作。2学会简单的矩阵输入与信号输入。3掌握部分绘图函数。 实验内容    1数组的加减乘除和乘方运算输入A=1 2 3，B=4 5 6，求C=A+B，D=A-B，E=A.*B，F=A./B，G=A.B,并用stem画出A，B，C，D，E，F，G。再输入一些数组，进行类似运算。S

5、tem函数格式可通过help stem命令查询2给出实现单位脉冲序列、单位阶跃序列的MATLAB函数，并用stem函数画出图形。注：(1) 单位脉冲序列可以通过借助MATLAB中的零矩阵函数zeros表示。全零矩阵zeros(1,N)产生一个由N个零组成的列向量，对于有限区间的可以通过以下MATLAB程序表示% 单位抽样序列实现程序k=-30:30;delta=zeros(1,30),1,zeros(1,30);stem(k,delta)(2) 单位阶跃序列可以通过借助MATLAB中的单位矩阵函数ones表示。单位矩阵ones(1,N)产生一个由N个1组成的列向量，对于有限区间的可以通过以下M

6、ATLAB程序表示% 单位阶跃序列实现程序k=-30:30;uk=zeros(1,30),ones(1,31);stem(k,uk)3已知两离散序列分别为和，试用MATLAB绘出它们的波形及的波形。注：对于离散序列来说，序列相加是将信号对应时间序号的值逐项相加，在这里不能象连续时间信号那样用符号运算来实现，而必须用向量表示的方法，即在MATLAB中离散序列的相加需表示成两个向量的相加，因而参加运算的两序列向量必须具有相同的维数。实现离散序列相加的MATLAB实用子程序如下function f,k=lsxj(f1,f2,k1,k2)%实现f(k)=f1(k)+f2(k),f1,f2,k1,k2是

7、参加运算的二离散序列及其对应的时间序列向量,f和k为返回的和序列及其对应的时间序列向量k=min(min(k1),min(k2):max(max(k1),max(k2);%构造和序列长度s1=zeros(1,length(k);s2=s1; %初始化新向量s1(find(k>=min(k1)&(k<=max(k1)=1)=f1;%将f1中在和序列范围内但又无定义的点赋值为零s2(find(k>=min(k2)&(k<=max(k2)=1)=f2;%将f2中在和序列范围内但又无定义的点赋值为零f=s1+s2; %两长度相等序列求和stem(k,f,'

8、;filled')axis(min(min(k1),min(k2)-1),(max(max(k1),max(k2)+1),(min(f)-0.5),(max(f)+0.5)%坐标轴显示范围该实验首先建立一个lsxj.m文件，即上段程序，然后存储在默认路径下。该实验参考程序：%求两离散序列之和实现程序f1=-2:2;k1=-2:2;f2=1 1 1;k2=-1:1;subplot 221;stem(k1,f1),axis(-3 3 -2.5 2.5);title('f1k');subplot 222;stem(k2,f2),axis(-3 3 -2.5 2.5);titl

9、e('f2k');subplot 223;f,k=lsxj(f1,f2,k1,k2);%调用自定义的lsxj函数stem(k,f);title('fk=f1k+f2(k)');4试用MATLAB绘出3中两离散序列乘法的波形。注：与离散序列加法相似，这里参加运算的两序列向量必须具有相同的维数。实现离散时间信号相乘的MATLAB实用子程序如下function f,k=lsxc(f1,f2,k1,k2)%实现f(k)=f1(k)+f2(k),f1,f2,k1,k2是参加运算的二离散序列及其对应的时间序列向量,f和k为返回的和序列及其对应的时间序列向量k=min(min

10、(k1),min(k2):max(max(k1),max(k2);%构造和序列长度s1=zeros(1,length(k);s2=s1; %初始化新向量s1(find(k>=min(k1)&(k<=max(k1)=1)=f1;%将f1中在和序列范围内但又无定义的点赋值为零s2(find(k>=min(k2)&(k<=max(k2)=1)=f2;%将f2中在和序列范围内但又无定义的点赋值为零f=s1.*s2; %两长度相等序列求和stem(k,f,'filled')axis(min(min(k1),min(k2)-1),(max(max(k

11、1),max(k2)+1),(min(f)-0.5),(max(f)+0.5)%坐标轴显示范围同样，在该实验中现建立lsxj.m文件存放在默认路径下，然后再调用该文件实现两序列相乘。该实验参考程序如下%求两离散序列之积实现程序f1=-2:2;k1=-2:2;f2=1 1 1;k2=-1:1;subplot 221;stem(k1,f1),axis(-3 3 -2.5 2.5);title('f1k');subplot 222;stem(k2,f2),axis(-3 3 -2.5 2.5);title('f2k');subplot 223;f,k=lsxc(f1,

12、f2,k1,k2);title('fk=f1k*f2(k)'); 实验报告要求1 简述实验目的。2 按实验步骤附上实验程序及结果。实验二 用MATLAB进行离散系统的Z域分析一、实验目的掌握应用MATLAB实现IIR、FIR数字滤波器结构的方法，即熟悉应用MATLAB实现IIR、FIR数字滤波器结构的相关函数。二、实验仪器PC机（装有MATLAB软件）三、MATLAB几个信号处理工具箱函数1impz：求H(z)的Z反变换h(n)h,T=impz(B,A,N)h为存放h(n)的列向量，时间变量N存放在列向量T中，当N为标量时，表示T=0：N-1，计算h(n)，n=0，1，2，N-

13、1；当N为向量时，T=N，仅计算N指定的整数点上的h(n)。2freqz：求数字滤波器H(z)的频率响应函数H=freqzB,A,w计算由向量w指定的数字频率点上数字滤波器H(z)的频率响应，结果存在H向量中。H,w=freqzB,A,M计算出M个频率点上的频率响应，存放在H向量中，M个频率存放在w向量中，freqz函数自动将这M个频率点均匀设置在频率范围0，之间。若缺省w和M时，函数自动选取512个频率点计算。不带输出向量的freqz函数将自动绘制幅频和相频曲线。3zplane:绘制H(z)的零极点图zplane(z,p)绘制出列向量z中的零点（以符号“”表示）和列向量p中的极点（以符号“&

14、#215;”表示）以及参考单位圆，并在多阶零点和极点的右上角标出其阶数。如果z和p为矩阵，则会以不同颜色绘出z和p各列中的零点和极点。zplane(B,A)绘制出系统函数H(z)的零极点图。四、实验内容1利用MATLAB绘制离散系统的零极图对于离散系统其系统函数可由式表示，则系统函数的零点和极点可以用MATLAB的多项式求根函数roots()来实现，调用该函数的命令格式为：p=roots(A)其中A为待求根的多项式的系数构成的行向量，返回向量p则包含该多项式所有的根位置列向量。值得注意的是：求系统函数零极点时，离散系统的系统函数可能有两种形式，一种是分子和分母多项式按Z的降幂次序排列，另一种是

15、分子和分母多项式按的升幂次序排列。若是以Z的降幂次序排列，则系数向量一定要由多项式的最高幂次开始，一直到常数项，缺项要用0补齐; 若是以按的升幂次序排列，则分子和分母多项式系数向量的维数一定要相同，不足的要用0补齐，不则的零点或极点就可能被漏掉。下面是求系统函数零极点，并绘制其零极点图的MATLAB实用函数ljdt(),该函数在绘制系统零极点图的同时，还绘出了Z平面的单位圆。function ljdt(A,B)% The function to draw the pole-zero diagram for discrete systemp=roots(A); %求系统极点q=roots(B);

16、 %求系统零点p=p' %将极点列向量转置为行向量q=q' %将零点列向量转置为行向量x=max(abs(p q 1); %确定纵坐标范围x=x+0.1;y=x; %确定横坐标范围clfhold onaxis(-x x -y y) %确定坐标轴显示范围w=0:pi/300:2*pi;t=exp(i*w);plot(t) %画单位园axis('square')plot(-x x,0 0) %画横坐标轴plot(0 0,-y y) %画纵坐标轴text(0.1,x,'jImz')text(y,1/10,'Rez')plot(real(

17、p),imag(p),'x') %画极点plot(real(q),imag(q),'o') %画零点title('pole-zero diagram for discrete system') %标注标题hold off完成下列实验内容：（1）已知某离散系统的系统函数为试用MATLAB求出该系统的零极点，并画出零极点分布图，判断系统是否稳定。（2） 已知某离散系统的系统函数为试用MATLAB求出该系统的零极点，并画出零极点分布图，求系统的单位冲激响应和幅频响应，并判断系统的是否稳定。2. 利用MATLAB分析离散系统的零极图分布与系统单位响应时域

18、特性的关系已知离散系统的零极分布分别如图2-1所示，其中虚线表示单位圆，试用MATLAB分析系统单位响应的时域特性。图2-1 离散系统的零极点分布图因系统的零极点分布图已知，则系统的系统函数就可知，故可以利用MATLAB函数impz()函数可求出系统单位响应。对于图2-1所示的系统，其系统函数分别为、实验报告要求 1 简述实验目的。2 按实验步骤编写相应的程序，有实验报告中附上实验程序及结果。实验三 傅立叶变换实验目的1 熟悉傅立叶变换的各种性质2 熟悉基本信号的频域转换3 熟悉应用FFT对典型信号进行频谱分析的方法4 了解应用FFT进行信号频谱分析过程中可能出现的问题，以便在实际中

19、正确应用FFT实验原理时域信号处理我们已经比较熟悉，信号的频谱函数对于我们却是一个全新的概念。一个信号的时域转换可以通过傅立叶变换（DFT）来完成。有限长序列的DFT是其Z变换在单位圆上的等距采样，或者说是序列傅立叶变换的等距采样，因此可以用于序列的频谱分析。FFT并不是与DFT不同的另一种变换，而是为了减少DFT运算次数的一种快速算法。它是对变换式进行一次次分解，使其成为若干小数点的组合，从而减少运算量。常用的FFT是以2为基数的，其长度N=2L。从频率采样定理知道，N点序列x(n)的N个离散时间傅立叶变换等间隔样本能唯一地重构x(n)。这些单位圆上的N个样本叫做离散傅立叶系数。设，为一周期

20、（具有无限持续时间）序列，则它的主周期为具有有限持续时间的离散傅立叶变换，N点序列的离散傅立叶变换由下式给出： N点序列的离散傅立叶反变换为： 在MATLAB中用fft函数实现快速傅立叶变换，如fft(x)实现2点fft变换，fft(x,N) 实现N点fft变换，N必须为2n。实验内容1、MATLAB实现基本序列的离散傅里叶变换计算已知复正弦序列，余弦序列，分别对序列求当和时的DFT，并绘出幅频特性曲线，并分析两种N值下DFT是否有差别，及产生原因。程序清单如下：N=16;N1=8;n=0:N-1;k=0:N1-1;x1n=exp(j*pi*n/8); %产生x1(n)X1k=fft(x1n,

21、N); %计算N点DFTx1(n)X2k=fft(x1n,N1); %计算N1点DFTx1(n)x2n=cos(pi*n/8); %产生x2(n)X3k=fft(x2n,N); %计算N点DFTx2(n)X4k=fft(x2n,N1); %计算N1点DFTx2(n)Subplot(2,2,1);stem(n,abs(X1k),'.');axis(0,20,0,20);ylabel('|X1(k)|')title('16点的DFTx1(n)')subplot(2,2,2);stem(n,abs(X3k),'.');axis(0,20

22、,0,20);ylabel('|X2(k)|')title('16点的DFTx2(n)')subplot(2,2,3);stem(k,abs(X2k),'.');axis(0,20,0,20);ylabel('|X1(k)|')title('8点的DFTx1(n)')subplot(2,2,4);stem(k,abs(X4k),'.');axis(0,20,0,20);ylabel('|X2(k)|')title('8点的DFTx2(n)')N点离散傅里叶变换的一种物

23、理解释就是是以N为周期的周期延拓序列的离散傅立叶级数系数的主值区间序列，即。当N=16时， 和正好分别是、的一个周期，所以和的周期延拓序列就是这两个单一频率的正弦序列，其离散傅里叶级数的系数分别如图3.14中的16点的DFTx1(n)和16点的DFTx2(n)所示。而当N=8时， 和正好分别是、的半个周期，所以和的周期延拓序列就不再是单一频率的正弦序列，而是含有丰富的谐波成分，其离散傅立叶级数的系数与N=16时的差别很大，因此对信号进行谱分析的时候，一定要截取整个周期，否则得到错误的频谱。2、MATLAB验证N点DFT的物理意义已知，绘制相应的幅频和相频曲线，并计算图示N=8和N=16时的DF

24、T。程序清单如下：N1=8;N2=16; % 两种FFT的变换长度n=0:N1-1;k1=0:N1-1; k2=0:N2-1;w=2*pi*(0:2047)/2048;Xw=(1-exp(-j*4*w)./(1-exp(-j*w); %对x(n)的频谱函数采样2048个点可以近似的看作是连续的频谱xn=(n>=0)&(n<4); %产生x(n)X1k=fft(xn,N1); %计算N1=8点的X1(k)X2k=fft(xn,N2); %计算N2=16点的X2(k)subplot(3,2,1);plot(w/pi,abs(Xw);xlabel('w/')sub

25、plot(3,2,2);plot(w/pi,angle(Xw);axis(0,2,-pi,pi);line(0,2,0,0);xlabel('w/')subplot(3,2,3);stem(k1,abs(X1k),'.');xlabel('k(=2k/N1)');ylabel('|X1(k)|');hold onplot(N1/2*w/pi,abs(Xw) %图形上叠加连续频谱的幅度曲线subplot(3,2,4);stem(k1,angle(X1k);axis(0,N1,-pi,pi);line(0,N1,0,0);xlabel

26、('k(=2k/N1)') ;ylabel('ArgX1(k)');hold onplot(N1/2*w/pi,angle(Xw) %图形上叠加连续频谱的相位曲线subplot(3,2,5);stem(k2,abs(X2k),'.');xlabel('k(=2k/N2)');ylabel('|X2(k)|');hold onplot(N2/2*w/pi,abs(Xw)subplot(3,2,6);stem(k2,angle(X2k),'.');axis(0,N2,-pi,pi);line(0,N2,

27、0,0);xlabel('k(=2k/N2)') ;ylabel('ArgX2(k)');hold onplot(N2/2*w/pi,angle(Xw)前面的理论已经知道，序列x(n)的N点DFT的物理意义是对在0，2上进行N点的等间隔采样。由结果可以直观的看到X(k)与之间的采样关系。3、MATLAB实现快速卷积用快速卷积法计算下面两个序列的卷积，给出实验结果并进行分析。快速卷积计算框图如图所示：图 快速卷积框图FFT的变换长度L必须满足LN+M-1，输出y(n)才等于x(n)和h(n)的线性卷积。计算两个序列的卷积时，也可以直接调用函数conv来计算，因为M

28、ATLAB中的计时比较粗糙，所以只有当N和M较大的时候，才能比较两种方法的执行时间快慢。实验报告要求 1. 简述实验目的。2. 按实验步骤附上实验程序及结果，并给出相应的分析。实验四 IIR及FIR滤波器的MATLAB实现 一、实验目的掌握应用MATLAB实现IIR、FIR数字滤波器结构的方法，即熟悉应用MATLAB实现IIR、FIR数字滤波器结构的相关函数。二、实验仪器PC机（装有MATLAB软件）三、实验内容编写MATLAB 程序，实现以下题目：（一）IIR滤波器结构实现1在MATLAB中，用filter函数实现IIR的直接形式用直接型实现系统函数为的IIR数字滤波器，求单位脉冲

29、响应和单位阶跃信号的输出。2在MATLAB中给定级联型系统函数，由扩展函数casfiltr实现IIR的级联形式系统函数，用级联型结构实现。程序清单如下：b0=4;B=1,1,0;1,-1.4142136,1;A=1,-0.5,0;1,0.9,0.81;N=60;delta=impseq(0,0,N);h=casfiltr(b0,B,A,delta);x=ones(1,5),zeros(1,N-5);y=casfiltr(b0,B,A,x);subplot(1,2,1);stem(h);title('直接型h(n)');subplot(1,2,2);stem(y);title(&

30、#39;直接型y(n)');扩展函数casfiltr程序：function y=casfiltr(b0,B,A,x);% IIR各FIR滤波器的级联型的实现% % y=casfiltr(b0,B,A,x);% y=输出序列% b0=级联型的增益系数% B=包含各bk和K乘3维实系数矩阵% A=包含各ak的K乘3维实系数矩阵% x=输入序列%K,L=size(B);N=length(x);w=zeros(K+1,N);w(1,:)=x;for i=1:1:K w(i+1,:)=filter（B（i,:）,A(i,:),w(i,:)）;endy=b0 * w(K+1,:);单位抽样序列的生

31、成函数impseq.mfunction x,n=impseq(n0,ns,nf)n=ns:nf;x=(n-n0)=0; %序列的起点为ns，终点为nf，在n=n0点处生成一个单位脉冲。3在MATLAB中用给定并联系统函数，由扩展函数parfiltr实现IIR的并联形式。用并联型实现系统函数为的滤波器。程度清单如下：C=0;B=-14.75,-12.90;24.50,26.82;A=1,-7/8,3/32;1,-1,0.5;N=60;delta=impseq(0,0,N);h=parfiltr(C,B,A,delta);x=ones(1,5),zeros(1,N-5);y=parfiltr(C,

32、B,A,x);subplot(1,2,1);stem(h);title('并联型h(n)');subplot(1,2,2);stem(y);title('并联型y(n)');扩展函数parfiltr:function y=parfiltr(C,B,A,x);% IIR滤波器的并联型实现% % y=parfiltr(C,B,A,x);% y=输出序列%C=当M>=N时(FIR)的多项式部分%B=包含各bk的K乘2维实系数矩阵%A=包含各ak的K乘3维实系数矩阵% x=输入序列%K,L=size(B);N=length(x);w=zeros(K+1,N);w(1,:)=filter(C,1,x);for i=1:1:K w(i+1,:)=filter(B(i,:),A(i,:),x);endy=sum(w);（二）FIR滤波器结构实现1直接形式FIR直接型系统函数表示为，即FIR结构的直接形式由包含系数的行向量b描述。在MATLAB中，把分母矢量a置为1，这种结构用filter函数实现。2级联形式FIR级联型系统函数表示为，即级联型FIR滤波器可以通过casfiltr函数实现。但这种形式与IIR形式类似，也可以使用dir2cas函数，把分母矢量a置为1，用cas2dir函数从级联形式转换为直接形式而获得。FIR滤波器的系统函数为