数字信号处理实验指导书2016

1、数字信号处理实验指导书王宏志 吕洪武 张李梅 孙慧然计算机科学与工程学院2016年3月目录网络资源-基于WEB的数字信号处理实验教学系统- 1 -实验一 熟悉MATLAB语言环境- 2 -实验二 离散系统分析- 7 -实验三 用FFT进行信号的频谱分析- 9 -实验四 用窗函数设计FIR数字滤波器- 12 -实验五 设计IIR数字滤波器- 13 -实验六 随机功率谱估计及MATLAB实现- 18 -附录： MATLAB简介- 19 - 27 -网络资源 基于WEB的数字信号处理实验教学系统单击”软件下载”即可下载实验系统,如下图所示:软件安装界面如下所示:软件运行界面如下所示:实验一 熟悉MA

2、TLAB语言环境一、实验目的1 熟悉Matlab的基本使用方法，重点掌握常用于数字信号处理的相关指令。2 利用Matlab实现序列的显示，运算等，加深对信号处理原理课程所学内容的理解。二、实验内容及步骤1. 双击Matlab 图标，进入主窗口，如图1所示。其中右侧的是命令窗口(Command Window)。Current Directory显示的是当前的目录，如果要在命令窗口中直接调用一个M文件，则必须保证其在这个目录下。不同的机器上，Current Directory可能不同。图1 进入Matlab时的画面2. 下面的操作都在命令窗口(Command Window)中进行。1） 键入hel

3、p subplot 并回车，将会看到命令subplot的使用方法的帮助。以后遇到不会的命令可以用键入“ help <某命令>”的方式来查看其使用方法。subplot的用法 ：subplot(m,n,p)是将一个窗口分成m*n个小窗口，p是小窗口的编号，方向是从左至右。2） 键入help stem 并回车，将会看到命令stem的使用方法的帮助。stem用来画离散序列的柄状图。3. 下面我们新建立一个M文件。首先点击图1所示的主菜单File下的“新建”图标，将弹出图2 所示的窗口。接着在这个弹出的窗口点击“保存”图标，建立一个文件impseq.m，将其保存在Current Direct

4、ory下（当前的目录，见图1）。图2 新建一个文件impseq.m4. 我们在这个impseq.m文件中输入如下语句并保存，注意以 %开头来写注释。functionx,n=impseq(np,ns,nf)% 生成x(n)=delta(n-np);nsnnf，即单位冲激序列(n)% np代表脉冲的位置, ns为序列起始位置，nf为序列终止位置。% 调用方式x,n=impseq(np,ns,nf)if ns>np|ns>nf|np>nferror('输入位置参数不满足ns=<n<=nf')else n=ns:nf;x=(n-np)=0;%if判断语句到

5、这里结束%上面的这条语句是关键语句，作用是如果满足条件(n-np)=0，则x=1。%在n=ns:nf的一串值中，只有一个值会满足这个逻辑式，因而只在这个n=np%处,x为1，其余的n值 处，x均为0.这样就构成了延迟np的单位冲激序列end 5.下面我们再回到图1所示的窗口1)打算显示一个序列x(n)=1.5*(n+1)(n3),下面的操作都在命令窗口(Command Window)中进行。注意：如果一个语句后面有分号，就不显示中间结果。n1=-4:5; x1=1.5*impseq(-1,-4,5)-impseq(3,-4,5); % 列出x1序列subplot(2,2,1); stem(n1

6、,x1,'.'); stem(n1,x1,'.'); title('x(n) 的序列图')ylabel('x1(n)'); axis(-5,5,-2,3);text(5.5,-2,'n');显示的结果如图3所示。图3 序列x(n)=1.5*(n+1)(n3)2)打算显示另外一个序列x(n)=nu(n)u(n8)10exp(0.3(n10)*u (n10)u(n16)，0n20新建一个文件stepseq.m,用于生成延迟的单位阶跃序列，其输入参数为序列起始位置ns,序列终止位置nf，及阶跃位置np.将刚才编辑的imp

7、seq.m另存为stepseq.m，并改动两条语句：functionx,n=stepseq(np,ns,nf)x=(n-np)>=0;存盘。下面的操作都在命令窗口(Command Window)中进行。n2=0:20;x21 = n2.*(stepseq(0,0,20)-stepseq(8,0,20); % 列出x21序列x22 = 10*exp(-0.3*(n2-10).*(stepseq(10,0,20)-stepseq(16,0,20);% 列出x22序列x2 = x21-x22; % x2序列是x21和x22之和subplot(2,2,2);stem(n2,x2,'.&#

8、39;);title('x(n) 的序列图')ylabel('x2(n)');参见附录 Matlab简介实验二 离散系统分析一、实验目的：1熟悉和掌握差分方程的求解方法。2熟悉和掌握离散付立叶变换和快速付立叶的算法和实现。3熟悉和掌握卷积运算。二、实验内容：在实验二中，主要有差分方程的求解，DFT及FFT的实现，信号卷积。差分方程可表示为：在设计中用到的方程是： 初始条件是： 在实验过程中可以改变系数，在这里只是为了简单，同学们可以在程序中修改一些参数来观察输出的变化。DFT的实现更简单，在MATLAB语言中有专门的函数，在这里我们选择的输入信号是最简单的余弦函

9、数：，至于序列的长度的选择，同学们可以在实验中自己选择，建议不要选择太长的序列长度，以免在'处理过程中影响速度。建议选择1264之间。'FFT是DFT的快速算法，它的执行速度远远快于DFT，所以可以选择较大的序列长度。信号卷积为，在MATLAB中，提供了卷积函数conv，调用十分方便。在实验中选定的两个信号为：系统单位冲击响应： 输入信号： 系统的序列长度可以选择的范围是50到100，信号的序列长度选择范围是20到50，采样频率一般选择在1000到3000Hz，至于参数选定，同学们可以自己来决定，实验中用到是 当然在以上设计中信号的选择可以是随意的，并不是固定的。三、实验报告要

10、求：1选择适当数据完成上述实验。2给出实验结果并进行分析。四、参考程序： 参见Untitled.m文件。实验三 用FFT进行信号的频谱分析一、实验目的：1在理论学习的基础上，通过实验，加深对快速傅立叶变换的理解，熟悉FFT算法；2熟悉应用FFT对典型信号进行频谱分析的方法；3了解应用FFT进行信号频谱分析过程中可能出现的问题，以便在实际中正确应用FFT。二、实验内容：在实验中，我们分析三个序列的频谱，分别是高斯序列，衰减正弦序列，还有一般序列。高斯序列表达式如下： 其中，,可以在实验中自己选定，一般选择规则如下，选择范围6-10，为2-8，序列长度选择为16-32，因为太长的话，处理速度会相对

11、比较慢。衰减正弦序列表示如下： 其中，可以选定的范围为0-1之间，序列长度选择范围为16-32，采样选择为1000-2000之间。实验中一般序列采用的是： 其中，可选的范围为1000-2000,序列长度为16-32。三、实验报告要求：1选择适当数据完成上述实验。2给出实验结果并进行分析。四、参考程序： 参见Untitled.m文件。1)编写程序，实现基二DIT(时域抽取) FFT。新建一个文件myditfft.m,输入如下语句：function y=myditfft(x)% 用MATLAB语言编写的基2 DIT FFT子程序% % y=myditfft(x)% -% 本程序对输入序列 x 实现

12、DIT-FFT基2算法，% 点数取大于等于x长度的2的幂次% x为给定时间序列% y为x的离散傅立叶变换% m=nextpow2(x);N=2m; % 求x的长度对应的2的最低幂次mif length(x)<N x=x,zeros(1,N-length(x); % 若x的长度不是2的幂，补零到2的整数幂end nxd=bin2dec(fliplr(dec2bin(1:N-1,m)+1; % 求1:2m数列的倒序y=x(nxd); % 将x倒序排列作为y的初始值for mm=1:m% 将DFT作m次基2分解,从左到右，对每次分解作DFT运算Nmr=2mm;u=1; % 旋转因子u初始化为W

13、N0=1WN=exp(-i*2*pi/Nmr); % 本次分解的基本DFT因子WN=exp(-i*2*pi/Nmr) for j=1:Nmr/2 % 本次跨越间隔内的各次蝶形运算 for k=j:Nmr:N % 本次蝶形运算的跨越间隔为Nmr=2mm kp=k+Nmr/2; % 确定蝶形运算的对应单元下标 t=y(kp)*u; % 蝶形运算的乘积项 y(kp)=y(k)-t; % 蝶形运算 y(k)=y(k)+t; % 蝶形运算 end u=u*WN; % 修改旋转因子,多乘一个基本DFT因子WN endend对这个程序的说明（1） y=nextpow2(x):如果x是一个数，用来求最靠近x并

14、大于x 的二的幂次。例如nextpow2(9.5)=4，因为最靠近9.5，并比它大的2的乘幂是24=16;如果x是一个向量，则nextpow2(x)= nextpow2(length(x)。即求x长度最近的二的幂次。（2） Matlab信号处理工具箱中也有倒序函数bitrevorder,为弄清意义，可以自编。（3） dec2bin(x)把十进制数转换为二进制。bin2dec(x) 把二进制数转换为十进制。fliplr(x)把x数组排列左右反转。（4） 变量mm 的含义相当于FFT分段的序数，最大相当于log2N2) 验证myditfft的正确性，下面的操作都在命令窗口(Command Wind

15、ow)中进行。K=input('K= ');%设定数据长度的二的幂次Kx=randn(1,2K);%q先生成一个x向量tic ,X=myditfft(x),toc %验证myditfft的正确性，并测试所需运行时间。实验四 用窗函数设计FIR数字滤波器一、实验目的：1熟悉FIR滤波器设计的基本方法；掌握用窗函数设计FIR数字滤波器的原理及方法；熟悉线性相位FIR滤波器的幅频特性和相位特性；2了解各种不同窗函数对滤波器性能的影响。二、实验内容：这次实验是为了了解各种窗函数的特性，所以实验中用到的窗函数比较多，分别如下：矩形窗，三角窗，汉宁窗，海明窗，布拉克曼窗，凯泽窗。关于这几个

16、窗函数的表达形式同学们可以在教材中查到，这里就不一一列举出来了，实验中信号的采样频率可以选择为1000-3000Hz，数字滤波器的截至频率为200Hz，阶数选择范围为41-100。窗函数选择原则（1）具有较低的旁瓣幅度，尤其是第一旁瓣的幅度；（2）旁瓣的幅度下降的速率要快，以利于增加阻带的衰减；（3）主瓣的宽度要窄，这样可以得到比较窄的过渡带。实验中还用到了一个特例，具体指标如下： 其中阶数的选择范围为32-100之间三、实验报告要求：1选择适当数据完成上述实验。2给出实验结果并进行分析。 3简单说明用不同窗函数设计的滤波器之间的性能差异；四、参考程序： 参见Untitled.m文件。实验五

17、设计IIR数字滤波器一、实验目的1、了解两种工程上最常用的变换方法：脉冲响应不变法和双线性变换法。2、掌握双线性变换法设计IIR滤波器的原理及具体设计方法，熟悉用双线性设计法设计低通、带通和高通IIR数字滤波器的计算机程序。3、观察用双线性变换法设计的滤波器的频域特性，并与脉冲响应不变法相比较，了解双线性变换法的特点。二、实验内容及步骤(i)在各种信号序列中，以巴特渥斯低通数字滤波器设计为例，可以将双线性变换法设计数字滤波器的步骤归纳如下：1、确定数字滤波器的性能指标。这些指标包括：通带、阻带临界频率；通带内的最大衰减；阻带内的最小衰减；采样周期T。2、确定相应的数字频率，；3、计算经过频率预

18、畸的相应参考模拟低通原型的频率，。4、计算低通原型阶数N；计算3dB归一化频率c，从而求得低通原型的传递函数Ha(s)。5、用变换公式 ，代入Ha(s)，求得数字滤波器传递函数.6、分析滤波器频域特性，检查其指标是否满足要求。(ii)自编写实现双线性变换的文件bilinear0.m，内容如下：function bd,ad = bilinear0(ba,aa,Fs)% % 双线性系数变换子程序% % bd,ad = bilinear0(ba,aa,Fs)% -% 从s域到z域的双线性频带变换% 实现:% b(z) b(s)|% - = -| 1-z(-1) Nz(z)% a(z) a(s)|s

19、= 2Fs- = -% 1+z(-1) Dz(z)% ba,aa按s的正幂降序排列至常数项结束% Nz,Dz按z的负幂降序排列，从常数项开始，双线性变换：Nz2*Fs*1,-1;Dz=1,1;% bd,ad按z的负幂降序排列，从常数项开始% 从常数项开始，把ba,aa按s的负幂降序排列。即将ba,aa中的短者左边补零成同长lba=length(ba);laa=length(aa);ld=laa-lba;if ld>=0 ba=zeros(1,ld),ba; % 若aa长度大于ba,給ba前补ld个零else aa=zeros(1,-ld),aa; % 若ba长度大于aa,給aa前补ld个

20、零endNz = 2*Fs*1,-1; Dz=1,1; % 双线性变换分子分母系数向量N = max(lba,laa)-1; % 模拟系统阶次Nbd = 0; ad = 0; % bd,ad系数向量初始化for k = 0:N pld = 1;pln = 1; % 双线性变换分子分母系数乘积项初始化 for l = 0:k-1 pld = conv(pld,Dz); % 求双线性变换分母系数k次幂Dzk end for l = 0:N-k-1 pln = conv(pln,Nz); % 求双线性变换分子系数(N-k)次幂Nz(N-k) end bd = bd+ba(k+1)*conv(pln,

21、pld); % 分子系数多项式向量求和 ad = ad+aa(k+1)*conv(pln,pld); % 分母系数多项式向量求和endad1 = ad(1); ad = ad/ad1; bd = bd/ad1; % 用分母系数多项式首项使分子分母系数归一化(iii)用巴特沃斯滤波器原型设计一个低通滤波器，满足：p =0.2，R p =1dB; s =0.3，A s =15dB;滤波器的采样频率为1000Hz自编写实现此功能的文件hc835.m，内容如下：% % 双线性变换法由巴特沃思变数字滤波器% % 数字滤波器指标:wp = 0.2*pi; % 数字通带频率(弧度)ws = 0.3*pi;

22、% 数字阻带频率(弧度)Rp = 1; % 通带波动(dB)As = 15; % 阻带衰减(dB) % 模拟原型指标的频率逆映射T = 0.001; Fs = 1/T; % 置 T=0.001OmegaP = (2/T)*tan(wp/2); % 原型通带频率预修正OmegaS = (2/T)*tan(ws/2); % 原型阻带频率预修正ep = sqrt(10(Rp/10)-1); % 通带波动参数Ripple = sqrt(1/(1+ep*ep); % 通带波动Attn = 1/(10(As/20); % 阻带衰减% 模拟巴特沃思原型滤波器计算:N,OmegaC =buttord(Omeg

23、aP,OmegaS,Rp,As,'s'); % 原型的阶数和截止频率计算%*巴特沃思滤波器阶次 = 6 z0,p0,k0 = buttap(N); % 归一化巴特沃思原型设计函数p = p0*OmegaC; z = z0*OmegaC; % 将零极点乘以Omegac，得到非归一化零极点k = k0*OmegaCN; % 将k0乘以OmegacN,得到非归一化kba0 = real(poly(z0);ba0 = k0*ba0 % 由零点计算分子系数向量aa0 = real(poly(p0) % 由极点计算分母系数向量ba = real(poly(z);ba = k*ba % 由零

24、点计算分子系数向量aa = real(poly(p) % 由极点计算分母系数向量bd,ad = bilinear(ba,aa,Fs); % 双线性变换:bd1,ad1 = bilinear(ba0,aa0,Fs/OmegaC); % 双线性变换:sos,G = tf2sos(bd,ad) % 变为二阶环节级联结构 % 绘图figure(1); subplot(1,1,1)db,mag,pha,grd,w = myfreqz(bd1,ad1); % 检验频率响应subplot(2,2,1); plot(w/pi,mag); title('幅度响应')xlabel('

25、9;); ylabel('|H|');axis(0,1,0,1.1);set(gca,'XTickMode','manual','XTick',0,0.2,0.3,1); grid % 画刻度线%set(gca,'YTickmode','manual','YTick',0,Attn,Ripple,1)subplot(2,2,3); plot(w/pi,db); title('模值(dB)');xlabel('频率：（单位：pi）'); ylabel(&

26、#39;分贝'); axis(0,1,-40,5);set(gca,'XTickMode','manual','XTick',0,0.2,0.3,1); % 画刻度线set(gca,'YTickmode','manual','YTick',-50,-15,-1,0); gridset(gca,'YTickLabelMode','manual','YTickLabels','50''15'' 1'&

27、#39; 0')subplot(2,2,2); plot(w/pi,pha/pi); title('相位响应')xlabel(''); ylabel('单位:pi'); axis(0,1,-1,1); set(gca,'XTickMode','manual','XTick',0,0.2,0.3,1); % 画刻度线set(gca,'YTickmode','manual','YTick',-1,0,1); gridsubplot(2,2,4);

28、plot(w/pi,grd); title('群延迟')xlabel('频率（单位：pi）'); ylabel('样本'); axis(0,1,0,10)set(gca,'XTickMode','manual','XTick',0,0.2,0.3,1); % 画刻度线set(gca,'YTickmode','manual','YTick',0:2:10); gridset(gcf,'color','w') % 置图形背景

29、色为白还需要编写函数myfreqz，新建一个文件myfreqz.m，输入如下内容：function db,mag,pha,grd,w = myfreqz(b,a);% % freqz 子程序的扩展版本，可得出分贝幅特性和群迟延特性% % db,mag,pha,grd,w = myfreqz(b,a);% -% db = 0 到pi弧度区间内的相对振幅(db)% mag = 0 到pi弧度区间内的绝对振幅% pha = 0 到pi弧度区间内的相位响应% grd = 0 到pi弧度区间内的群延迟% w = 0 到pi弧度区间内的501个频率样本向量% b = Ha(z)的分子多项式系数(对FIR

30、b=h)% a = Ha(z)的分母多项式系数(对 FIR: a=1)%H,w = freqz(b,a,1000,'whole'); % 用MATLAB中的freqz函数计算数字系统频率响应H = (H(1:1:501)' w = (w(1:1:501)' % 取其前一半，并化为列向量mag = abs(H); % 求其幅特性db = 20*log10(mag+eps)/max(mag); % 化为分贝值pha = angle(H); % 求其相特性grd = grpdelay(b,a,w); % 求其群迟延特性程序运行之后，可以观察到频率响应曲线如图所示。图:

31、频率响应曲线实验六 随机功率谱估计及MATLAB实现一、实验目的：1了解功率谱估计的几种算法，比较它们之间的优缺点；2熟悉谱估计的算法，以便在实际中应用。二、实验内容：实验中用到的方法有周期性图法，Bartlett平均周期图法，Weleh法及一般方法。在周期图法中，采样频率选择范围为1000-3000Hz，FFT算法取样点选择范围为512-1024。用到的含有噪声的随机序列为： 其中，表示随机序列是Matlab语言自动产生的。在Bartllett平均周期图法中，采样频率选择范围为1000-3000 Hz，FFT算法取样点选择范围为512-1024。用到的含有噪声的随机序列为： 其中，表示随机序

32、列是Matlab语言自动产生的。窗函数阶数选择为32-64。置信区间小于1大于0就行。Welch法中，采样频率选择范围为1000-3000 Hz，FFT算法取样点选择范围为512-1024。用到的含有噪声的随机序列为： 其中，表示随机序列是Matlab语言自动产生的。窗函数阶数选择为32-64。置信区间小于1大于0就行。三、实验报告要求：1选择适当数据完成上述实验。2给出实验结果并进行分析。 3比较几种方优缺点。四、参考程序：参见Untitled.m文件。附录： MATLAB简介一、 MATLAB搜索路径MATLAB是通过搜索路径来查找M文件的，因此，MATLAB系统文件、Toolboxes工

33、具箱函数、用户自己编写的M文件等都应在搜索路径之内。MATLAB 6.1提供了更强的有关搜索路径的管理功能。当用户输入一个标识符比如Value时，MATLAB按下列步骤处理：检查Value是否为变量；检查Value是否为内部函数；在当前目录下是否存在Value.m文件；在MATLAB搜索路径上是否存在Value.m文件。如果在搜索路径上存在多个Value.m文件，则只执行所找到的第一个Value.m文件；如果找不到这一文件，则给出出错信息。利用MATLAB提供的path,addpath和rmpath函数，可修改搜索路径：·path可显示出当前的搜索路径；·path(s)表示

34、将路径设置到字符串s；·addpath/matlab/design和path(path,/matlab/design)表示将目录C:matlabdesign加入到搜索路径（设当前盘为C：）；·rmpath/matlab/design表示从搜索路径中删去目录C:matlabdesign。搜索路径保存在local子目录的pathdef.m文件中，用户可利用普通编辑器，直接修改搜索路径。MATLAB还专门提供了管理搜索路径的路径浏览器，这可从File菜单中选择Set Path（设置路径）命令或在工具栏中点击Path Browser（路径浏览器）图标。在搜索路径浏览窗口中，可很方便

35、地加入或删去子目录，也可以在左边目录区中通过拖拉方式改变搜索顺序。搜索路径修改后，按Refresh按钮可使修改的路径在本次任务中起作用，按Save Settings可使修改结果在以后启动时生效。另外，what命令可显示出搜索路径上的文件名，如：whatwhat matlab/design可分别显示出当前目录和matlabdesign目录中的文件名。要显示指定文件（如value.m）的内容，可采用type命令：type value要对文件value.m进行编辑，可输入：edit value二、 MATLAB集成环境MATLAB既是一种语言，又是一种编程环境。在这一环境中，系统提供了许多编写、调试

36、和执行MATLAB系统命令窗口。1 MATLAB命令窗口MATLAB的工作空间是输入命令和输出结果的窗口，在这里输入的命令会立即得到执行，并输出结果，这非常适用于编写短小的程序。对编写大型、复杂程序应用采用M文件编程方法。MATLAB命令窗口除了提供File（文件）、Edit（编辑）、Options（选项）、Window（窗口）、Help（帮助）这些菜单命令外，还提供了快捷操作的工具栏。通过工具栏的图标，可快捷地进行各种操作。在MATLAB工作空间，输入的命令行保存在缓冲区里，因此可很容易调出以前输入的命令并加以修改。MATLAB提供的窗口命令编辑键如下表所示。利用这些键可方便地修改以前的命令

37、。按键复合健功 能Ctrl-p调出前一行Ctrl-n调出下一行Ctrl-b光标向后移一个字符Ctrl-f光标向前移一个字符Ctrl-Ctrl-r光标向右移一个字在MATLAB执行程序过程中，可按Ctrl-c键以中断程序的执行。MATLAB程序结果的显示，可利用format命令加以控制。下面以变量x为例，给出各种格式及显示结果。x=4/3 1.2345e-6format short %短格式（缺省情况）1.3333 0.0000 format short e 1.3333e+000 1.2345e-006 format short g1.3333 1.2345-006 format long %

38、长格式 1.33333333333333 0.00000123450000 format long e 1.33333333333333+000 1.234500000000000-006 format long g 1.33333333333333 1.2345e-006 format bank %银行格式1.33 0.00format rat %比率格式 4/3 1/810045format hex %十六进制格式 3ff5555555555555 3eb4b6231abfd271 除了上述这些格式命令之外，MATLAB缺省显示为隔行显示（即format loose格式），采用 forma

39、t compact%紧凑格式可改变成逐行显示。2 编辑M文件将MATLAB语句按特定的顺序组合在一起来，就得到了MATLAB程序，其文件名的后缀为M，故也称为M文件。MATLAB6.1提供了M文件的专用编辑调试器，在这一编辑器中，会以不同的颜色表示不同的内容，这分成五种：命令、关键字、不完整字符串、完整字符串及其它文本。这样很容易发现输入错误，缩短调试时间。启动编辑器的方法有两种：（） 在工作空间中键入：edit fname这时可启动编辑器，并打开fname.m文件。（） 在命令窗口的File菜单或工具栏上选择New命令或New.File图标。为方便操作，编辑器除了提供菜单命令外，也提供了工具

40、栏。MATLAB编辑调试器提供了编辑M文件和调试M文件这两大功能。这里先介绍编辑功能。MATLAB编辑器与其它Windows编辑程序类似，这里不再赘述。只对下列几点作特别说明：（） 在编辑M文件时，可直接转到指定的行，这可从Edit菜单中选择Go To Line命令。（） 可直接计算M文件中表达式的值，结果显示在命令窗口中。这可通过选择表达式，然后在View菜单中选择Selection命令来实现。（） 可根据MATLAB的句法自动缩排，以增加M文件的可读性。先选择文本块，然后在View菜单中选择Auto Indent Selection命令。（） 可按要求设置自动缩排的格式。在View菜单中选

41、择Options命令。三、 MATLAB程序设计初步在MATLAB工作环境下，我们很容易输入各种命令，以便完成指定的功能。然而直接在MATLAB环境下输入命令，边解释边运行，这多少给人带来不便之处：输入等待、修改不便、程序保存和检查困难等等。幸好MATLAB提供了更方便的方法来进行程序设计。1 脚本文件和函数文件定义MATLAB的M文件有两类：脚本文件和函数文件。我们将原本要在MATLAB环境下直接输入的语句，放在一个以.m为后缀的文件中，这一文件就称为脚本文件。有了脚本文件，可直接在MATLAB中输入脚本文件名（不含后缀），这时MATLAB会打开这一脚本文件，并依次执行脚本文件中的每一条语句

42、，这与在MATLAB中直接输入语句的结果完全一致。另一类M文件是函数文件，它的第一行必须是函数定义行。函数文件由五部分构成：·函数定义行；·H1行；·函数帮助文本；·函数体；·注释。例如，函数文件mean.m为function y=mean(x) 函数定义行%MEAN Average or mean value. H1行% For vectors, MEAN(X) is the mean of X. % For matrices, MEAN(X) is a row vector 函数帮助文本 % containing the mean valu

43、e of each column. % m,n=size(x); if m=1 m=n; 函数体 end y=sum(x)/m; 我们以这个函数为例来说明函数的各个部分。（1）函数定义行function y=mean(x)其中，function为函数定义的关键字，mean为函数名，y为输出变量，x为输入变量。当函数具有多个输出变量时，则以方括号括起；当函数具有多个变量时，则直接用圆括号起。例如：function x,y,z=sphere(theta, phi, rho) 当函数不含输出变量，则直接略去输出部分或采用空方括号表示，如：function printresults(x) functi

44、on =printresults(x) 所有在函数中使用和生成的变量为局部变量（除非利用global语句定义），这些变量值只能通过输入和输出进行传递。因此，在调用函数时应通过输入变量将参数传递给函数；函数调用返回时也应通过输出变量将运算结果传递给函数调用者；其它在函数中产生的变量在返回时被全部清除。（2）H1行在脚本和函数文件中，以开头的行称为注释行，也就是说，之后的字符不被MATLAB执行。在函数文件中，其第二行一般是注释行，这一行称为H1行，实际上它是帮助文本中的第一行。H1行不仅可以由help function_name命令显示，而且，lookfor命令只在H1行内搜索，因此这一行内容提

45、供了这个函数的重要信息。（3）函数帮助文本这部分内容是以开头的帮助文本，它用来比较详细地说明这一函数。当在MATLAB输入help function_name时，可显示出H1行和函数帮助文本。这部分文本从H1行开始，到第一个非开头的行结束。（4）函数体函数体是完成指定功能的语句实体，它可采用任何可用的MATLAB命令，包括MATLAB提供的函数和用户设计的M函数。（5）注释在函数文件中，除了函数定义行和函数体之外，其它部分都是可以省略的，不是必须有的。但作为一个函数，为了提高函数的可用性，应加H1行和函数帮助文本；为了提高函数的可读性，应加上适当的注释。2脚本文件和函数文件比较脚本文件和函数文

46、件之间有一些本质上的差异。函数文件去掉其第一行的定义行，就转变成了脚本文件。但这样一来，原先在函数内使用的局部变量就成了基本工作空间中的变量，这会带来几个问题：·基本工作空间中与脚本文件中同名的变量会引导起冲突；·使基本工作空间中变量数急剧增加，造成内存紧张；·编程时要细心考虑各个脚本文件所用到的变量。这些问题在函数文件中不复存在，MATLAB通过实参与形参一一对应的方式来实现函数的调用，这极大地方便了程序设计。例如，分别编写出求取平均值与标准差的脚本文件stat1.m和函数文件stat2.m：stat1.m %脚本文件求阵列x的平均值和标准差m,n=size(x

47、); if m=1 m=n; end s1=sum(x); s2=sum(x.2); mean1=s1/m; stdev=sqrt(s2/m-mean1.2); stat2.m function mean1,stdev=stat2(x) %函数文件求阵列x的平均值和标准差调用格式为mean,stdev=stat2(x) %m,n=size(x); if m=1 m=n; end s1=sum(x); s2=sum(x.2); mean1=s1/m; stdev=sqrt(s2/m-mean1.2); 编写好这两个文件后，我们在MATLAB下来执行这两个文件，从而对脚本文件和函数文件有一个基本的

48、了解。在MATLAB中输入：clear x=rand(4,4)+2; stat1执行后我们可检查一下基本工作空间中的变量情况：whos Name Size Bytes Class m 11 8 double array mean1 14 32 double array n 11 8 double array s1 14 32 double array s2 14 32 double array stdev 14 32 double array x 44 128 double array这说明，在脚本文件中产生的所有变量都返回到了基本工作空间。再检查执行结果：disp(mean1;stdev)

49、2.5685 2.5321 2.6684 2.5605 0.2787 0.3359 0.1513 0.2817 另一方面，我们可通过函数文件来进行同样的操作，这时输入：clear m n s1 s2 mean1 stdev m1,st=stat2(x);执行后，同样检查基本工作空间的变量情况：whos Name Size Bytes Class m 1 14 32 double array st 14 32 double array x 14 128 double array这说明，在基本工作空间中，除了原本产生的x矩阵，调用函数stat2.m后，只增加了由函数返回的结果。通过disp(m1;

50、st)可检查到与stat1.m相同的结果。四、 图形绘制1简单图形绘制这里以产生一个简单的正弦函数曲线为例来说明图形的绘制，这一过程在MATLAB中是很简单的。设要产生从0到之间的正弦函数，则可按下列步骤进行：(1) 产生x轴、y轴数据：x=0:pi/20:2*piy=sin(x);(2) 打开一个新的图形窗口：figure(1)(3) 绘制出正弦曲线：plot(x,y,r-) r-表示以红色实线绘制出正弦曲线。(4) 给图形加上栅格线：gird on这样就可得到正弦曲线。从这一过程可以看出，在MATLAB中建立曲线图是很方便的。我们还可以将图形窗口进行分割，从而可在每个窗口各绘制一条曲线。例如，将图形窗口分割成的窗格，分别绘制出