matlab语言编程self.ppt

1、2.5MATLAB程序MATLAB是解释性语言，输入一行语句后回车，就会立即执行得出结果。如果要实现比较复杂的功能，单靠一条一条地在命令窗口中输入指令执行，效率是很低的。如何解决这个问题呢？为此MATLAB提供了扩展名为“.m”的文本文件，在文件中事先写入一行行的MATLAB命令，存盘后从MATLAB的命令窗口调入执行（类似于DOS下的批处理），这种文件称为底稿文件或MATLAB脚本文件。,在文本编辑器中输入语句，然后将文件存盘，例如命名文件名为“my1stprg.m”，保存于MATLAB的默认工作路径中。然后回到MATLAB命令窗口，在提示符“”下键入文件名(可以省略扩展名)后回车，即可运行

2、程序得到结果。,2.5.1M文件简介MATLAB除了如前所述的在命令窗口进行的直接交互的指令操作方式外，另外一种更为重要的工作方式就是M文件的编程工作方式。,M文件有两种形式：一种是脚本文件，另一种是函数文件。,脚本文件（命令文件）实际上是一串指令的集合，与在命令窗口逐行执行文件中的所有指令，其结果是一样的。如果M文件的第一行是以关键字“function”开头的，则就是函数文件。,M文件,二者在内存使用上存在重要的区别：,函数文件内部所定义的变量仅仅在该函数文件内部有效，函数返回后这些内部变量将自动被清除，也就是说，函数内部所定义的变量仅在函数内部起作用，是局部变量。,而脚本文件中所定义的变量

3、，在使用“clear”命令清除之前，始终存在于工作空间当中，是全局变量。,注意:不要和MATLAB已经定义的系统函数和其它自定义的函数同名。,另外，函数文件的文件名必须和函数名相同.,（1）函数定义行（关键字function）functionout1,out2,.=filename(in1,in2,.)输入和输出（返回）的参数个数分别由nargin和nargout两个MATLAB保留的变量来给出。（2）第一行帮助行，即H1行以（%）开头，作为lookfor指令搜索的行（3）函数体说明及有关注解以（%）开头，用以说明函数的作用及有关内容（4）函数体语句函数体内使用的除返回和输入变量这些在func

4、tion语句中直接引用的变量以外的所有变量都是局部变量，即在该函数返回之后，这些变量会自动在MATLAB的工作空间中清除掉。如果希望这些中间变量成为在整个程序中都起作用的变量，则可以将它们设置为全局变量。,函数M文件,1.M文件的建立M文件是一个文本文件，它可以用任何编辑程序来建立和编辑，而一般常用且最为方便的是使用MATLAB提供的文本编辑器,有3种方法：(1)菜单操作。File菜单中选择New菜单项，再选择M-file命令。(2)命令操作。在MATLAB命令窗口输入命令edit。(3)命令按钮操作。单击MATLAB主窗口工具栏上的NewM-File命令按钮。,只是一个简单的ASCII码文本

5、文件，执行程序时逐行解释运行程序，matlab是解释性的编程语言。,M文件的打开(1)菜单操作。从File菜单中选择Open命令，在Open对话框中选中所需打开的M文件。(2)命令操作。在MATLAB命令窗口输入命令：edit文件名。(3)命令按钮操作。单击MATLAB主窗口工具栏上的OpenFile命令按钮，再从弹出的对话框中选择所需打开的M文件。,3.数据和文件的输入输出1)inputinput指令提示用户从键盘输入数据、字符串或表达式，并接受该输入。指令的调用格式有两种：A=input(prompt)%输入数据或表达式A=input(prompt,s)%输入字符串例如：xm=input(

6、Whatsyourname?,s);,2)数据的输出MATLAB提供的命令窗口输出函数主要有disp函数，其调用格式为disp(输出项)其中输出项既可以为字符串，也可以为矩阵。,例输入x,y的值，并将它们的值互换后输出。程序如下：x=input(Inputxplease.);y=input(Inputyplease.);z=x;x=y;y=z;disp(x)disp(y),数据文件.matmat文件是matlab以标准二进制格式保存的数据文件，可将工作空间中有用的数据变量保存下来。mat文件的生成和调用是由函数save和load完成的。,(3)matlab文件的类型,m文件,m文件有两类:独立

7、的m文件称命令文件(脚本文件)可调用m文件称函数文件,4)利用文件输入输出数据指令“save”和“load”用于MATLAB与磁盘的数据交换。“save”指令将MATLAB工作空间的数据存入磁盘；“load”指令则将数据从磁盘送入到MATLAB工作空间。举例说明如下：A=1,2;3,4;b=81;save(Cmydate.mat,A,b)%将变量A和b保存在文件C:mydate.mat,load(Cmydate.mat,A,b);%调入数据文件whos%显示调入的变量NameSizeBytesClassA2x232doublearrayb11x18doublearrayGrandtotalis

8、5elementsusing40bytes,2.5.2程序控制任何计算机语言，只要利用顺序结构、循环结构以及分支结构，就可以完成任何程序功能。由于MATLAB语言矩阵计算功能十分强大，常常仅使用顺序结构借以矩阵的逻辑运算就可以完成计算任务。由于循环结构和分支结构在MATLAB语言中的运行速度相对较慢，因此在算法优化的编程中应当尽可能避免使用，而代之以矩阵运算，从而提高程序运行速度，简化程序代码，并使得程序代码更加接近于数学表达。,1）.for循环语句基本格式for循环变量起始值：步长：终止值循环体end步长缺省值为1，可以在正实数或负实数范围内任意指定。对于正数，循环变量的值大于终止值时，循环

9、结束；对于负数，循环变量的值小于终止值时，循环结束。循环结构可以嵌套使用。,1.循环结构,例如，最简单的forend循环：fork=1:10x(k)=k.2;end事实上，采用矩阵思想也可以获得相同的结果，但编程更加简单明了：k=1:10;x=k.2,绝大部分循环都可以遵从这样的方法变成向量化(矩阵化)的算法，避免采用循环语句，从而大大提高程序的执行效率。,2）.while循环语句基本格式while表达式循环体end若表达式为真，则执行循环体的内容，执行后再判断表达式是否为真，若不为真，则跳出循环体，向下继续执行。,While循环和for循环的区别在于，while循环结构的循环体被执行的次数不

10、是确定的，而for结构中循环体的执行次数是确定的。,与循环结构相关的语句还有break语句和continue语句,它们一般与if语句配合使用。break语句用于终止循环的执行。当在循环体内执行到该语句时，程序将跳出循环，继续执行循环语句的下一语句。continue语句控制跳过循环体中的某些语句。当在循环体内执行到该语句时，程序将跳过循环体中所有剩下的语句，继续下一次循环。,2.条件分支结构if，switch,if表达式语句1else语句2end,if表达式执行语句end,(1)ifelseend语句有3种形式,if表达式1多分支语句1elseif表达式2语句2elseif表达式3语句3else

11、语句nend,这种结构实现了一种多路选择，比较复杂。可替代一般高级语言中的switch-case-end语句,在条件表达式中，通常都是由关系操作符,=,=,=等逻辑操作符%清除MATLAB内存空间,这一命令常用于脚本文件的首句%以避免前面命令在内存空间形成的变量的意外影响F=inline(t,t);%用inline函数建立积分的被积函数forx=1,3,5;f_x(k)=x3+x+log(x)*sin(x)+quadl(F,0,x);%计算表达式值k=k+1;%数组下标加1end,f-x%显示计算结果文件编辑存盘后,在MATLAB命令窗口执行:,helpcalcfx这是开头的注释行,可以用he

12、lp命令显示程序calcfx.m的功能是计算表达式calcfxf-x=2.500034.6550140.9567,下面是函数编程的例子，函数的文件名与函数名相同，为stat.m。程序2-2functionmean,stdev=stat(x)n=length(x);%求输入向量x的长度mean=sum(x)/n;%求平均值stdev=sqrt(sum(x-mean).2/n);%求均方根值,m-函数实例:利用级数展开近似计算sin,sin的级数展开为,1:functions=powersin(x);2:%POWERSIN.Powerseriesforsin(x).3:%POWERSIN(x)tr

13、iestocomputesin(x)fromapowerseries,4:s=0;5:t=x;6:n=1:,7:whileabs(t)eps8:s=s+t;9:t=-x.2/(n+1)*(n+2).*t;10:n=n+2;11:end,利用“tic”和“toc”指令可以对程序段的执行时间进行测定，从而估计出程序执行效率，并找出改进程序、提高效率的方法。“tic”用于计时开始，而“toc”用于计时结束并显示计时结果。,6.测定程序执行时间和时间分配,7.提高程序执行速度的原则要提高编程的执行效率，则一定要根据MATLAB的特点来编写程序，具体有以下原则：,其次，在使用大型数组或矩阵之前对其进行初

14、始化，即采用指令“zeros”或“ones”对矩阵定维，这样可以减少MATLAB在内存分配过程中的耗时，大大提高速度。,首先，在编程中要尽量避免采用循环语句。利用向量化语句来代替循环语句可以大大提高程序运行速度。如果不得不采用多重循环，那么内循环的次数应该尽可能多于外循环的次数。,第三，应该优先考虑使用MATLAB的内在函数。MATLAB的内在函数是采用C语言优化构造的，并固化在MATLAB的内核中，其运行速度可以和C语言的速度等价，而用户自行编制的M文件则是在MATLAB中解释执行的。另外，采用更先进更有效的算法也可以提供计算效率。快速傅立叶变换就是一个典型的例子。,2.5.4MATLAB编

15、程仿真,仿真也称为模拟，在本质上，系统的计算机仿真就是根据物理系统的运行原理建立相应的数学描述并进行计算机数值求解的过程。,建模和仿真过程,建立系统数学模型或仿真模型。编写系统的计算机仿真模型。求解。分析。,求解算法算法的改进简化系统模型,关键问题：,1.系统仿真的数学基础,2.计算机仿真的一般过程,建立计算机程序（层次化的）。执行仿真。对仿真模型和仿真结果的检验。,对仿真数学模型有效性的验证。对计算机仿真模型（程序）的验证。对仿真算法的验证。仿真结果置信度分析。,仿真验证包含以下方面内容：,例1：试对空气中在重力作用下不同质量物体的下落过程进行建模和仿真。已知重力加速度g=9.8m/s2，在

16、初始时刻t0=0s时物体由静止开始坠落。空气对落体的影响可以忽略不计。,质量为m的物体在自由坠落过程中受到竖直向下的恒定重力的作用，由牛顿第二定律，我们知道，重力F，加速度a以及物体质量m之间的关系是：,F=ma其中a=g,初始时刻为t0=0，物体的初始速度为v(t0)=0，并设物体下落的瞬时速度为v(t)。,（1）建立数学模型,设物体在t时刻的位移为s(t),并设初始位移为零，即s(t0)=0,根据加速度、速度、位移三者之间的微积分关系，我们得到一组数学方程：,a=dt/dvv=dt/dsF=ma,以及初始条件（也称为方程的边界条件）,v(t0)=0s(t0)=0,（2）数学模型的解析分析,

17、将方程转换为一种在自变量（时间）上的“递推”表达式,（3）根据数学模型建立计算机仿真模型（编程）,v(t+dt)=v(t)+dv=v(t)+adts(t+dt)=s(t)+ds=s(t)+v(t)dt,（4）执行仿真和结果分析,%ch1example1prg1.mg=9.8;%重力加速度v=0;%设定初始速度条件s=0;%设定初始位移条件t=0;%设定起始时间dt=0.1;%设置计算步长N=20;%设置仿真递推次数.仿真时间等于N与dt的乘积fork=1:Nv=v+g*dt;%计算新时刻的速度s(k+1)=s(k)+v*dt;%新位移t(k+1)=t(k)+dt;%时间更新end,%理论计算,

18、以便与仿真结果对照t_theory=0:0.01:N*dt;%设置解析计算的时间点v_theory=g*t_theory;%解析计算的瞬时速度s_theory=1/2*g*t_theory.2;%解析计算的瞬时位移%作图:仿真结果与解析结果对比t=0:dt:N*dt;plot(t,s,o,t_theory,s_theory,-);xlabel(时间t);ylabel(位移s);legend(仿真结果,理论结果);,5)仿真程序的功能扩展-以动态方式来观察物体坠落的过程,将作图语句放在递推计算循环内，并设置即时作图刷新方式，从而得到这种“动画”仿真的效果。,%ch1example1prg2.mg

19、=9.8;%重力加速度forL=1:5%仿真重复5次以便于观察v=0;%初始速度s=0;%初始位置t=0;dt=0.01;%计算步长fork=1:200v=v+g*dt;%速度s=s+v*dt;%位移t=t+dt;%时间,plot(0,-s,o);axis(-22-200);%坐标范围固定text(0.5,-1,当前时间：t=,num2str(t);text(0.5,-2,当前速度：v=,num2str(v);text(0.5,-3,当前位置：s=,num2str(s);set(gcf,DoubleBuffer,on);%双缓冲避免作图闪烁drawnow;%立即作图endend,例2对乒乓球的

20、弹跳过程进行仿真。忽略空气对球的影响，乒乓球垂直下落，落点为光滑的水平面，乒乓球接触落点立即反弹。如果不考虑弹跳中的能量损耗，则反弹前后的瞬时速率不变，但方向相反。如果考虑撞击损耗，则反弹速率有所降低。我们希望通过仿真得出乒乓球位移随时间变化的关系曲线，并进行弹跳过程的“实时”动画显示。,设球是刚性的，质量为m，垂直下落。碰击面为水平光滑平面。在理想情况下碰击无能量损耗。如果考虑碰击面损耗，则碰击前后速度方向相反，大小按比例系数K(0K1)下降。在t时刻的速度设为v=v(t),位移设为y=y(t),并以碰击点为坐标原点，水平方向为坐标横轴建立直角坐标系。球体的速度以竖直向上方向为正方向。重力加

21、速度为g=9.8m/s2。,首先对乒乓球弹跳过程进行一些理想化假设。,设初始时刻t0=0球体的初始速度为v0=v(t0),初始位移为y0=y(t0)。,v(t+dt)=v(t)-gdt在t+dt时刻小球的位移为y(t+dt)=y(t)+v(t)dt,在小球撞击水平面的瞬间，即y(t)=0的时刻，它的速度方向改变，大小按比例K衰减。当K=1时，就是无损耗弹跳情况。因此，小球反弹瞬间（t+dt时刻）的速度为,v(t+dt)=-Kv(t)-gdt;0K1,初始条件假设：,在空中时小球受重力F=mg作用，其中，g=-(dt/dv)。则在t+dt时刻小球的速度为:,受力分析：,反弹瞬间的位移为,y(t+

22、dt)=y(t)-Kv(t)dt=-Kv(t)dt,(2)仿真模型设计（程序）从数学模型中可见，小球在空中自由运动时刻与撞击时刻的动力方程不同。通过小球所处位置（位移）是否为零可判定小球处于何种状态。,%ch1example2prg1.mg=9.8;%重力加速度v0=0;%初始速度y0=1;%初始位置m=1;%小球质量t0=0;%起始时间K=0.85;%弹跳的损耗系数N=5000;%仿真的总步进数dt=0.001;%仿真步长v=v0;%初状态y=y0;,fork=1:Nif(y0)|(v0)%小球在空中的(含刚刚弹起瞬间)动力方程计算v=v-g*dt;y=y+v*dt;else%碰击瞬间的计算

23、y=y-K.*v*dt;v=-K.*v-g*dt;ends(k)=y;%将当前位移记录到s数组中以便作图end,t=t0:dt:dt*(N-1);%仿真时间长度plot(t,s);xlabel(时间t);ylabel(位移y(t);axis(0501.1);,%ch1example2prg2.mg=9.8;%重力加速度v0=0;%初始速度y0=1;%初始位置m=1;%小球质量t0=0;%起始时间K=0.85;%弹跳的损耗系数N=5000;%仿真的总步进数dt=0.005;%仿真步长v=v0;%初状态y=y0;,5)仿真程序的功能扩展-以动态方式来观察物体坠落的过程,fork=1:Nify0%小

24、球在空中的动力方程计算v=v-g*dt;y=y+v*dt;else%碰击瞬间的计算y=-K.*v*dt;v=-K.*v-g*dt;endplot(0,y,o);axis(-2201);set(gcf,doubleBuffer,on)%双缓冲避免作图闪烁drawnow;end,实例3试仿真得出一个幅度调制系统的输入输出波形。设输入被调制信号是一个幅度为2v，频率为1000Hz的余弦波，调制度为0.5,调制载波信号是一个幅度为5v，频率为10KHz的余弦波。所有余弦波的初相位为0。,（1）数学模型,该调幅系统的输入输出关系表达式为：,y(t)=(M+maMcos2fmt)cos2fct,其中，M=

25、2是被调信号的振幅，fm=1000是其频率。A=5是载波信号的幅度，fc=104是其频率。ma=0.5是调制度。,（2）编程实现在这个系统中的信号最高工作频率为(fm+fc)=11KHz，根据取样定理，只要离散取样率高于该频率的2倍即可无失真。在计算量和数据存储量许可的条件下，取样率可以设置更高，以使仿真计算的结果波形图显示更加光滑。我们设置取样率为105，即在一个载波周期上取样10次，相应的取样间隔为t=10-5s。本例中，取样间隔也作为仿真步进。程序代码如下:,%ch2example1prg1.mdt=1e-5;%仿真采样间隔T=3*1e-3;%仿真终止时间t=0:dt:T;input=2

26、*cos(2*pi*1000*t);%输入被调信号carrier=5*cos(2*pi*1e4*t);%载波output=(2+0.5*input).*carrier;%调制输出,%作图:观察输入信号,载波,以及调制输出subplot(3,1,1);plot(t,input);xlabel(时间t);ylabel(被调信号);subplot(3,1,2);plot(t,carrier);xlabel(时间t);ylabel(载波);subplot(3,1,3);plot(t,output);xlabel(时间t);ylabel(调幅输出);,在实际中用双踪示波器来测量调幅的基带输入波形和调制输

27、出波形时，只要基带输入余弦波的频率与调幅载波频率不是整数倍关系，示波器就只能够对其中一个信号进行同步显示，而另外一个信号在示波器上是滑动的。,如果调整示波器的扫描周期使之与两个测试信号的周期都不呈整数倍关系，那么示波器屏幕上显示的两个信号都是滑动。能否通过仿真程序来再现这种波形显示“滑动”现象呢？,（3）让输出波形“动”起来,考虑一个更接近真实物理系统的调幅模型。考虑了相位噪声n(t)的载波表达为：,v(t)=Acos(2fct+n(t)）,设信道是无衰减的，其加性噪声为n(t)，那么在接收机所收到的调幅信号r(t)为：,r(t)=y(t)+n(t),程序中，我们将仿真时间区间划分为若干段，每

28、段称为一帧，示波器的扫描周期等于帧周期，即我们每仿真得出一帧的数据就让显示刷新一次。程序中故意将输入被调信号的频率设置为1005Hz，这样其信号周期就与仿真的帧周期不是整数倍关系，运行后将看到“不断”滑动的被调信号，载波显示出相位抖动现象，而接收信号则沾染了噪声。,%ch2example1prg3.mdt=1e-6;%仿真采样间隔T=2*1e-3;%仿真的帧周期forN=0:500%总共仿真的帧数t=N*T+(0:dt:T);%帧中的取样时刻input=2*cos(2*pi*1005*t);%输入被调信号carrier=5*cos(2*pi*(1e4)*t+0.1*randn);%载波outp

29、ut=(2+0.5*input).*carrier;%调制输出noise=randn(size(t);%噪声r=output+noise;%调制信号通过加性噪声信道,%作图:观察输入信号,载波,以及调制输出subplot(3,1,1);plot(0:dt:T,input);xlabel(时间t);ylabel(被调信号);text(T*2/3,1.5,当前帧数：N=,num2str(N);subplot(3,1,2);plot(0:dt:T,carrier);xlabel(时间t);ylabel(载波);subplot(3,1,3);plot(0:dt:T,r);xlabel(时间t);yla

30、bel(调幅输出);set(gcf,DoubleBuffer,on);%双缓冲避免作图闪烁drawnow;end,实例4仿真电容的充电过程。,一个电压源通过电阻与电容串联的网络对电容充电。设t=0为初始时刻（初始时刻之前电路断开，不工作），电压源输出电压x(t)为单位阶跃函数，电容两端的电压为y(t)，回路电流为i(t)，并将电压源视为系统输入，电容上的电压视为系统输出。电路的初始状态为y(0)。如图：,电容的充电电路以及等价系统,（1）数学分析首先根据网络拓扑和元件伏安特性建立该电路方程组,并化简得：,该方程也称为系统的状态方程。在方程中，变量y代表电容两端的电压，是电容储能的函数。本例中它既是系统的状态变量，又是系统的输出变量。,（2）数值求解最直接的求解方法是将上式转换为以时间向前递进的计算递推式，并以微小仿真时