《nupt dsp 实验指导书.doc》由会员分享，可在线阅读，更多相关《nupt dsp 实验指导书.doc（16页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、实验预习：Matlab软件的使用Matlab软件由Mathworks公司推出，至今已发展为适合多学科、多种工作平台的功能强大的大型软件。在欧美高校，Matlab已经成为线性代数、自动控制理论、数理统计、数字信号处理、时间序列分析、动态系统仿真等高级课程的基本教学工具，成为攻读学位的大学生、硕士生、博士生必须掌握的基本技能。在设计研究单位和工业部门，Matlab被广泛用于科学研究和解决各种具体问题。在国内，Matlab也日益盛行。1 Matlab 的主要特点1.语言简洁紧凑，使用方便灵活，库函数极其丰富。利用其丰富的库函数避开烦杂的子程序编程任务，压缩了一切不必要的编程工作。2.运算符丰富。Ma
2、tlab提供了和C语言几乎一样多的运算符，灵活使用Matlab的运算符将使程序变得极为简短3.语法限制不严格，程序设计自由度大。如在Matlab里，无需对矩阵预定义就可使用。4.图形功能强大。在FORTRAN和C语言里，绘图很不容易，在Matlab里数据的可视化非常简单，Matlab 还具有较强的编辑图形界面的能力。5.具有功能强大的工具箱。6.程序的可移植性好。7.Matlab 的缺点是：与其它高级程序相比，程序的执行速度慢。2.Matlab的环境Matlab既是一种语言，又是一个编程环境。作为一个编程环境，Matlab提供了很多方便用户管理变量、输入输出数据以及生成和管理M文件的工具。到M
3、atlab 6.x 则把多种开发工具集成为Matlab桌面系统。包含如下8个组成部分：命令窗口（Command Window）、历史命令窗口（Command History）、资源目录本（Launch Pad）、当前路径浏览器（Current Directory Brower）、帮助浏览器（Help Brower）、工作空间浏览器（Workspace Brower）、数组编辑器（Array Editor）以及程序编辑调试器（EditorDebugger）。它们的功能简述如下：1.Matlab 的命令窗口（Command Window）在Windows桌面上，双击Matlab的图标，就可以进入M
4、atlab的工作环境。首先出现Matlab的标志图形，接着出现其默认的桌面系统。其左上视窗为资源目录（Launch Pad），可切换为工作空间浏览器（Workspace Brower）；其左下视窗为历史命令窗口（Command History），可切换为当前路径浏览器（Current Directory Brower）；右半个视窗则为命令窗口（Command Window）。命令窗是用户和Matlab做人机对话的主要环境。是它的提示符，可以在提示符后键入Matlab的各种命令并读出相应的结果，但在输入命令后加分号结束，则不会显示结果。2.Matlab 的工作空间(workspace)工作空间是
5、指运行Matlab 程序或命令所生成的所有变量和Matlab提供的常量构成的空间。每打开一次Matlab会自动建立一个工作空间，工作空间在Matlab运行期间一直存在，关闭Matlab后自动消失。刚打开的工作空间只有Matlab提供的几个常量，如pi(3.)、虚数单位i等。运行Matlab程序时。程序的变量被加入工作空间。工作空间的存在在调试程序时非常有用。3.Matlab 的程序编辑调试器（EditorDebugger）Matlab的程序编制有两种方式：一种称为行命令方式，这就是在命令窗中一行一行地输入程序，计算机每次对一行命令作出反应，像计算器那样，这只能编辑简单的程序。更好的程序编制方式
6、是把程序写成一个有多行语句组成的文件，让Matlab来执行这个文件。编写和修改这种文件程序就要用到程序编辑调试器。进入程序编辑器可以有三种方式：1.选择菜单栏的“File”项中的“New”或“Open”项；2.选择工具栏的“New”或“Open”按钮；3.在命令窗内输入edit命令；4.Matlab 的路径浏览器Matlab管理的文件范围由它的搜索路径来确定。该搜索路径由Matlab启动文件MATLABroottoolboxlocalmatlabbrc.m来规定。用户只有将要运行的程序所在的目录加入到Matlab 的路径搜索范围内，在Matlab命令窗口键入文件名后才能执行。要显示或修改搜索路
7、径，可以用path命令：1.path 列出Matlab的搜索路径；2.path(path 加入的新搜索路径) 在原搜索路径群中加入一新搜索路径用命令行设置路径要键入文件的全路径名，很容易出错，用菜单工具较好。Matlab6.x所设置的修改路径的菜单工具如下：在命令窗中单击路径浏览器按钮或选择file菜单里的“set path” ，即可设置搜索路径。5. Matlab 的帮助系统最常用的就是在命令窗口键入 “help 命令名” ，即可在命令窗口出现该命令的帮助信息。6. 历史命令窗口（Command History）历史命令窗口用于记录并显示本次工作进程中曾键入的全部行命令。7. 资源目录本（L
8、aunch Pad）资源目录本用于把用户在当前系统中安装的所有Matlab产品说明、演示以及帮助信息的目录集成起来，便于用户迅速调用查阅。8. 数组编辑器（Array Editor）用户可以直接在数组编辑器中修改所打开的数据，甚至可以更改该数据的数学结构以及显示方式。3. M文件的编写* 用Matlab 语言编写的可在Matlab中运行的程序，称为M文件。M文件有两种形式 ：脚本文件（Script File）和函数文件（Function File ）。这两种文件的扩展名，均为“ . m” 。3.1. M脚本文件对于一些比较简单的问题 ，在指令窗中直接输入指令计算 。对于复杂计算，采用脚本文件（
9、Script file）最为合适 。MATLAB只是按文件所写的指令执行 。M脚本文件的特点是：脚本文件的构成比较简单，只是一串按用户意图排列而成的（包括控制流向指令在内的）MATLAB指令集合。脚本文件运行后 ，所产生的所有变量都驻留在 MATLAB基本工作空间（Base workspace）中。只要用户不使用清除指令（clear）， MATLAB指令窗不关闭，这些变量将一直保存在基本工作空间中。3.M文件有两种形式 ：脚本(命令)文件（Script File）和函数文件（Function File ）。这两种文件的扩展名，均为“ . m”1.M脚本(命令)文件对于一些比较简单的问题 ，在指
10、令窗中直接输入指令计算 。对于复杂计算，采用脚本文件（Script file）最为合适 。MATLAB只是按文件所写的指令执行 。M脚本文件的特点是：脚本文件的构成比较简单，只是一串按用户意图排列而成的（包括控制流向指令在内的）MATLAB指令集合。脚本文件运行后 ，所产生的所有变量都驻留在 MATLAB基本工作空间（Base workspace）中。只要用户不使用清除指令（clear）， MATLAB指令窗不关闭，这些变量将一直保存在基本工作空间中。2. M函数文件与脚本文件不同 ，函数文件犹如一个“黑箱”，把一些数据送进并经加工处理，再把结果送出来。MATLAB提供的函数指令大部分都是由函
11、数文件定义的。M函数文件的特点是：从形式上看 ，与脚本文件不同，函数文件的笫一行总是以 “function”引导的“函数申明行”。从运行上看 ，与脚本文件运行不同，每当函数文件运行，MATLAB就会专门为它开辟一个临时工作空间，称为函数工作空间（ Function workspace）。当执行文件最后一条指令时，就结束该函数文件的运行，同时该临时函数空间及其所有的中间变量就立即被清除。MATLAB允许使用比“标称数目”较少的输入输出变量，实现对函数的调用。M函数文件的形式：function 输出参数函数名（输入参数）函数体 在命令文件和函数文件中都可以以%开头加注释行。编写函数文件时，文件名可
12、以和函数名不同，但调用时应使用文件名。不过最好把文件名和函数名统一，以免出错。3.3. 总结：所以，两者的区别在于：1.脚本文件没有输入参数，也不返回输出参数。函数文件可以有输入参数，也可返回输出参数。2.脚本文件对工作空间的变量进行操作。函数文件的变量为局部变量，只有其输入、输出变量保留在工作空间。4 . M文件的调试编写 M文件时，错误（Bug）在所难免。错误有两种：语法（Syntax）错误和运行（Runtime）错误。语法错误是指变量名、函数名的误写，标点符号的缺、漏等。对于这类错误，通常能在运行时发现，终止执行，并给出相应的错误原因以及所在行号。运行错误是算法本身引起的，发生在运行过程
13、中。相对语法错误而言，运行错误较难处理 。尤其是M函数文件，它一旦运行停止，其中间变量被删除一空，错误很难查找。有两种调试方法：直接调试法和工具调试法:1. 直接调试法：可以用下面方法发现某些运行错误。1) 在M文件中，将某些语句后面的分号去掉， 迫使M文件输出一些中间计算结果，以便发现可能的错误。2) 在适当的位置，添加显示某些关键变量值的语句（包括使用 disp 在内）。3) 在原M脚本或函数文件的适当位置，增添指令 keyboard 。 keyboard 语句可以设置程序的断点 。M文件运行时，在keyboard命令的地方将暂停，同时在命令窗口中显示字符“K”，这时，可以检查和修改函数工
14、作空间中的变量的值。在命令窗口键入return，再按回车键即可继续运行程序。4) 通过将原M函数文件的函数申明行注释掉，可使一个中间变量难于观察的M函数文件变为一个所有变量都保留在基本工作空间中的M脚本文件。2. GUI 界面调试法：MATLAB提供了一个基于GUI界面的调试。使用它，可以对函数进行调试。Debug菜单的使用：Continue：恢复程序运行至结束或另一个断点 。Single Step：单步执行函数。Step In：深入下层局部工作区 。Quit Debugging：退出调试状态。Set/Clear Breakpoint：设置/清除光标处的断点 。Clear All Breakp
15、oints：清除程序中的所有断点 。Stop if Error：运行至出错或结束。Stop if Warning：运行至警告消息或结束。Stop if NaN of Inf：运行至运算结果出现 NaN 或 Inf。5 . Matlab与数字信号处理简介5.1信号表示离散时间信号只在离散时间点上有值，是时间上不连续的“序列”。Matlab中的主要数据类型是二维或多维的实矩阵或复矩阵。数字信号处理过程中所用到的基本数据对象（ 例如：一维信号或序列，多通道信号，二维信号等等）都可以用矩阵表示。5.2. 波形发生器许多不同的函数都可以产生波形。其中大部分函数都需要一个时间向量作为参数。如果选择1000
16、Hz的抽样频率产生波形，则合适的时间向量应是：t=0:0.001:1这样就构造了一个有1001个元素的行向量。该向量表示时间从0到1秒，以千分之一秒为步长。例如：产生一正弦波 t=0:0.01:2; x=sin(2*pi*t);5.3 信号处理工具箱（toolbox） Matlab提供了许多现成的信号处理函数，可直接调用。详细内容可点主菜单工具条上的help，浏览Signal Processing Toolbox 即可。6. Matlab的循环控制循环结构、选择结构和顺序结构是复杂程序的基本组成单元，因此熟练掌握Matlab的循环结构和选择结构是编程的基本要求。1. 循环结构Matlab的循环
17、结构可用for语句和while语句来实现。1）for语句使用for语句使用较为灵活，一般用于循环次数已经确定的情况。for语句的调用格式为：for 矩阵变量名表达式1：表达式2：表达式3 语句体end其中，表达式1的值为循环的初值，表达式2的值为步长，表达式3的值为循环的终止值。如果省略表达式2，则默认步长为1。for语句允许嵌套，在程序里，每一个“for”关键字必须和一个“end”关键字配对，否则出错。下面是一个for循环的实例，该循环计算出15的乘法表。for n=1:5 for m=1:n r(m,n)=m*n; % r为一个矩阵， r(m,n)表示矩阵第m行第n列的元素 endend2
18、）while语句的使用while语句一般用于事先不能确定循环次数的情况。while语句的调用格式为：while 表达式 语句体end当表达式的值为真时，执行语句体；当表达式的值为假时，终止循环。当表达式的计算对象为矩阵时，只有当矩阵的所有元素均为真时，才执行语句体。2. 选择结构选择结构由if语句、switch语句和break语句来实现，这里我们主要介绍if语句。1）if语句的使用if语句的调用格式为：if 逻辑表达式1 语句体1elseif 逻辑表达式2 语句体2elseif 逻辑表达式3 语句体3else 语句体end 3. 循环程序的优化循环的向量化前面提过，Matlab的一个缺点是当对
19、矩阵的单个元素作循环时运算速度很慢。编程时，把循环向量化，不但能缩短程序的长度，更能提高程序的执行效率。由于Matlab的基本数据为矩阵和向量，所以编程时，应尽量对向量和矩阵编程，而不是像在其他应用程序里，对矩阵的元素编程。下面是一个向量化编程的例子： %一般循环编程 x=1; for k=1:1001 y(k)=log10(x); x=x+0.01;end % 向量化编程 x=1:0.01:10;y=log10(x);附：简单矩阵操作：若要输入矩阵，则必须在每一行结尾加上分号（；），如下例： A = 1 2 3 4; 5 6 7 8; 9 10 11 12; A A = 1 2 3 4 5
20、6 7 8 9 10 11 12 同样地，我们可以对矩阵进行各种处理：A(2,3) = 5 % 改变位于第二行，第三列的元素值 A = 1 2 3 4 5 6 5 8 9 10 11 12 B = A(2,1:3) % 取出部份矩阵B B = 5 6 5 A = A B % 将B转置後以列向量并入A A = 1 2 3 4 5 5 6 5 8 6 9 10 11 12 5 A(:, 2) = % 删除第二列（：代表所有行） A = 1 3 4 5 5 5 8 6 9 11 12 5 A = A; 4 3 2 1 % 加入第四行 A = 1 3 4 5 5 5 8 6 9 11 12 5 4 3
21、 2 1 A(1 4, :) = % 删除第一和第四行（：代表所有列） A = 5 5 8 6 9 11 12 5实验一一、观察采样引起的混叠设模拟信号为，t的单位为毫秒(ms)。1. 设采样频率为3kHz，确定与混叠的采样重建信号。2. 画出和在范围内的连续波形。（因数字计算机无法真正画出连续波形，可用较密的离散点的连线来近似。）3. 分别用和在两信号波形上标记出3kHz采样点。两信号波形是否相同？采样后的两序列是否相同？编程参考：a) 产生一正弦波: t=0:0.01:2; x=sin(2*pi*t);b）信号相乘x(n)=x1(n)x2(n) Matlab实现：x=x1.*x2； （“.
22、* ”命令完成向量各元素对应相乘）c）基本绘图：* plot是最基本的二维图形命令。例：t=0: pi/20 :2*pi ;plot(t,sin(t);若要画出多条曲线，只需将座标对依次放入plot函数即可： plot(t, sin(t), t, cos(t); * 在调用函数plot时可以指定线型、颜色和数据点的图标，基本调用格式为： plot(x,y,colour-linestyle-marker) 例：plot(t, sin(t), r, t, cos(t), m); 若要同时改变颜色及图线型态（Line style），也是在座标对後面加上相关字串即可：plot(t, sin(t), r
23、o, t, cos(t), m*);关于plot命令的更详细的信息可以在命令窗口键入命令“help plot”查看plot的帮助信息。* Matlab 自动将图形画在图形窗口上，图形窗口和命令窗口是独立的窗口。函数figure 可建立新的图形窗口，并把新建窗口指定为当前窗口用于输出图形。例：t=0: pi/20 :pi ; figure(1); plot(t,sin(t);figure(2); plot(t,cos(t);* 命令subplot 可以把多个图形画在一个图形窗口里。例：t=0: pi/20 :pi ; subplot(2,1,1) ; plot(t,sin(t);subplot(
24、2,1,2) ; plot(t,cos(t);d）用Matlab画模拟信号因为Matlab是对向量操作，所以画出的只能是在离散时间点有值的离散时间信号，那么要用Matlab画模拟信号只能是以尽量高的采样频率近似模拟信号。在本题中模拟信号和可以用采样率画出和来做近似。二、判别离散时间系统的时不变性。设输入序列为，系统实现对的抽取。1. 设。取延迟量D（例如D30）。记，画出、的序列波形。2. 编程求出系统对的响应以及对的响应3. 画出、的波形。该系统是否为时不变的？编程参考：利用help命令学习“length”、“zeros”函数的用法延时、抽取要注意的是序号的变化。矢量的下标只能为正整数，为0
25、或为负就会提示出错。三、利用卷积计算出输入信号通过FIR滤波器的输出，并观察输出信号的input-on暂态、input-off暂态和稳态部分。（计算卷积可用conv命令）考虑下面两个滤波器，第一个的单位脉冲响应为，另一个的单位脉冲响应为；输入为周期方波，在一个周期内。1） 分别画出两个滤波器的输出的波形，并与书上p144例4.1.8的两幅图比较是否一致。2） 计算出图中稳态部分的值。编程参考：1）利用help命令学习产生全1向量的函数“ones”和完成卷积操作的函数“conv”2）产生周期波形的方法：例：x=ones(1,25),zeros(1,25);p=6;x=x*ones(1,p);x=
26、x(:); % x(:)命令完成将x的各列连成一个长的列向量x=x;实 验 二 一、观察序列频谱，观察信号通过系统后波形与频谱的变化已知输入信号，其中，N可取5000点。（1）画出的前100点波形编程参考： 生成正余弦序列的方法在实验一已经给出。 画一序列的某部分样点的波形，可以用下面的简单方法： figure(1),plot(n1:n2,x1(n1:n2);对一图形中的几个不同曲线进行名称的标注，以区分不同曲线，可用legend命令。如：plot(n1:n2,x1,n3:n4,x2,g, n5:n6,x,b.-),legend(x1,x2,x)（2）画出的DTFT频谱（）编程参考：由于计算机
27、无法画出连续频谱，所以可在内均匀取足够密的点数，如M5000个频率点，求出这些频率点上的频谱值，并画出随变化的曲线。例如：M=5000; % set samples number in frequency domain w=pi/M*(0:M-1); % discrete frequencyXw=zeros(1,M);for k=1:M, Xw(k)=sum(x*(exp(-j*w(k)*(0:N-1); endfigure(2),plot(w/pi,abs(Xw),% 求复数的模时用函数abs（）xlabel(omega(pi),ylabel(|X(omega)|) % 对横、纵坐标标注以上
28、程序段完全是DTFT定义式的实现。当然，还可利用Matlab中的Signal Processing Toolbox 中的函数很简洁地画出频谱。但初学时最好能用笨办法做一遍。（3）某LTI系统，画出系统的幅度频响编程参考：求频响可以用（2）中的方法，也可以用函数freqz（）。用以下命令学习其用法： help freqz（4）求系统对的响应(可以自己编程也可利用卷积函数)。画出的波形,并与 的波形比较（各画100点）；画出的幅度谱，并与比较。 编程参考：求复数的模用函数abs，求复数幅角用函数angle。（5）问：滤波器是什么频响类型的滤波器？你从以上实验中观察到什么？与课本上的什么结论相吻合？
29、二、系统函数，根据正准型结构（canonical form）编写样本处理算法。内部状态的初始值设为零，输入信号采用逐个样本手动输入的方式(用input命令)，求输出信号。编程参考： for i=1:10, x=input(input x =);自编“sample-by-sample algorithm”yend实 验 三1 观察窗函数的影响。信号为，KHz， KHz， KHz， 采样频率 KHz。a) 写出（）的频谱；b) 分别画出窗长度 ，的矩形窗频谱和Hamming窗频谱。观察L的变化对窗谱的主瓣宽度、旁瓣密集度、相对旁瓣水平的影响。 c) 时域采样点数分别取，分别画出加矩形窗及加Hamm
30、ing窗时DTFT频谱；你能否从频谱上分辨出信号的三个频率分量？若能分辨出，它们的位置和相对大小是否准确？L的大小对频率的物理分辨率（physical frequency resolution）有何影响？2 理解频率的物理分辨率和计算分辨率的区别信号同上，加矩形窗。时域采样点数分别取，。画出以上各种时长情况下，频域采样点数分别为，时的DFT（在同一个图上用虚线画出相应的DTFT频谱，用于比较）。离散频谱DFT和连续频谱DTFT有什么关系？L一定的情况下，能否通过增加N改善频率的物理分辨率？N的作用是什么？编程参考：生成Hamming窗函数：可用定义自己生成，或者用hamming（）函数。对信号
31、x加窗：可用类似 x=x.*w 语句。其中，“.*”为元素对应相乘。加矩形窗也可用简单的截取。求DFT谱：可用函数 fft ()。求DTFT谱：可用定义或加大点数下的fft ()。实 验 四1 窗口法设计FIR数字滤波器(a) 用矩形窗设计长度分别为N=11、41、81、121的低通FIR滤波器，要求截止频率为。画出滤波器的单位冲激响应和幅度频响曲线。问题：理想滤波器的频响是怎样的？当N增大时，FIR滤波器在附近的最大纹波幅度是否降低？其余纹波的幅度是否减小？纹波的密度怎样变化？过渡带宽度怎样变化？（如有必要可增大N值观察）。在N=11时，画出滤波器的相频曲线。它是否是线性的？（b）用汉明窗再
32、次设计同样的滤波器。用汉明窗设计出的滤波器与用矩形窗相比有什么特点？编程参考： 当k0时为0/0型，Matlab无法求值，但易知。所以在给出FIR的h（n）时注意需赋值h(M+1)=wc/pi 。2 以Butterworth 模拟低通滤波器为原型，设计IIR数字滤波器。（a）截止频率。设计11阶IIR数字低通滤波器，画出幅频、相频曲线。（b）截止频率。设计11阶IIR数字高通滤波器，画出幅频、相频曲线。 问题：所设计的IIR滤波器与FIR滤波器的频率特性有何区别？编程参考：可用函数butter（）。大家最好读一下其help内容的第一段，掌握用法。 参考程序实验一参考程序% =% problem
33、 1% =clear% estimate x(t) and xa(t) with a much higher sampling freq. fs1time_period=6; % unit: msfs1=50; % unit: kHzT1=1/fs1; % unit: msn1=0:fix(time_period/T1);x=cos(5*pi*n1*T1)+4*sin(2*pi*n1*T1).*sin(3*pi*n1*T1);xa=cos(pi*n1*T1);% obtain x(nT) and xa(nT) with given sampling freq. fsfs=3; T=1/fs;
34、n=0:fix(time_period/T);x_sample=cos(5*pi*n*T)+4*sin(2*pi*n*T).*sin(3*pi*n*T);xa_sample=cos(pi*n*T);figure,plot(n1*T1,x,r,n1*T1,xa,b,n*T,x_sample,ro),hold on, stem(n*T,xa_sample,b:x)legend(x(t),xa(t),x(nT),xa(nT),xlabel(t(ms)% =% problem 2% =clear% plot x(n) and x(n-D)D=30;N=500;n=1:N;x=sin(2*pi/100*n);for n=1:N+D, if (n-D)=0, xD(n)=0; else xD(n)=x(n-D); endendfigure,subplot(2,1,1),plot(1:N,x,r:,1:length(xD),xD,b),legend(x(n),xD(n),xlabel(n)% plot y(n) and yD(n)for n=1:fix(N/2) y(n)=x(2*n);endfor n=1:length(y)+D, if (n-D)=L, dft=abs(fft(x,N); else dft=zeros(N,1);