Matlab中文教程共63页
- 格式:doc
- 大小:275.00 KB
- 文档页数:61
Matlab中文手册目录 (1)第1章MATLAB 6.5环境 (11)1.1 MATLAB简介 (11)1.1.1 MATLAB工具箱 (11)1.1.2 MATLAB功能和特点 (12)1.2 MATLAB 6.5环境设置 (13)1.2.1 菜单栏 (13)1.2.2 工具栏 (18)1.2.3 通用操作界面窗口 (18)1.3 MATLAB 6.5帮助 (32)1.4 MATLAB 6.5其他管理 (35)1.4.1 MATLAB用户文件格式 (35)1.4.2设置搜索路径 (36)1.4.3文件管理命令 (37)1.4.4 退出MATLAB (39)1.5 一个实例 (39)第2章MATLAB数值计算 (43)2.1 变量和数据 (43)2.1.1数据类型 (43)2.1.2数据 (43)2.1.3变量 (46)2.2 矩阵和数组 (47)2.2.1矩阵输入 (47)2.2.2矩阵元素和操作 (54)2.2.3字符串 (65)2.2.4矩阵和数组运算 (72)2.2.5多维数组 (92)2.3稀疏矩阵 (98)2.3.1稀疏矩阵的建立 (98)2.3.2稀疏矩阵的存储空间 (103)2.3.3稀疏矩阵的运算 (105)2.4多项式 (105)2.4.1多项式的求值、求根和部分分式展开 (105)2.4.2多项式的乘除法和微积分 (109)2.4.3多项式拟合和插值 (112)2.5元胞数组和结构数组 (115)2.5.1元胞数组 (116)2.5.2结构数组 (121)2.6数据分析 (127)2.6.1数据统计和相关分析 (127)2.6.2差分和积分 (129)2.6.3卷积和快速傅里叶变换 (134)2.6.4向量函数 (137)第3章MATLAB符号计算 (138)3.1 符号表达式的建立 (138)3.1.1 创建符号常量 (138)3.1.2 创建符号变量和表达式 (141)3.1.3 符号矩阵 (143)3.2符号表达式的代数运算 (145)3.2.1符号表达式的代数运算 (146)3.2.2 符号数值任意精度控制和运算 (149)3.2.3 符号对象与数值对象的转换 (152)3.3符号表达式的操作和转换 (154)3.3.1符号表达式中自由变量的确定 (154)3.3.2符号表达式的化简 (156)3.3.3符号表达式的替换 (161)3.3.4求反函数和复合函数 (164)3.3.5 符号表达式的转换 (166)3.4 符号极限、微积分和级数求和 (169)3.4.1符号极限 (169)3.4.2符号微分 (171)3.4.3符号积分 (174)3.4.4符号级数 (176)3.5 符号积分变换 (178)3.5.1傅里叶(Fourier)变换及其反变换 (178)3.5.2拉普拉斯(Laplace)变换及其反变换 (180)3.5.3 Z变换及其反变换 (181)3.6符号方程的求解 (183)3.6.1代数方程 (183)3.6.2符号常微分方程 (185)3.7符号函数的可视化 (186)3.7.1符号函数的绘图命令 (186)3.7.2图形化的符号函数计算器 (189)3.8 Maple函数的使用 (189)3.8.1访问Maple函数 (190)3.8.2 获得Maple的帮助 (191)第4章MATLAB计算的可视化和GUI设计 (192)4.1二维曲线的绘制 (192)4.1.1基本绘图命令plot (192)4.1.2绘制曲线的一般步骤 (198)4.1.3多个图形绘制的方法 (199)4.1.4曲线的线型、颜色和数据点形 (202)4.1.5设置坐标轴和文字标注 (204)4.1.6交互式图形命令 (210)4.2 MATLAB的三维图形绘制 (211)4.2.1绘制三维线图命令plot3 (211)4.2.2绘制三维网线图和曲面图 (212)4.2.3立体图形与图轴的控制 (217)4.2.4色彩的控制 (220)4.3 MATLAB的特殊图形绘制 (225)4.3.1条形图 (225)4.3.2面积图和实心图 (227)4.3.3直方图 (228)4.3.4饼图 (230)4.3.5离散数据图 (231)4.3.6对数坐标和极坐标图 (232)4.3.7等高线图 (234)4.3.8复向量图 (235)4.4图形窗口的功能 (236)4.5对话框 (238)4.6句柄图形 (242)4.6.1句柄图形体系 (242)4.6.2图形对象的操作 (242)4.6.3图形对象属性的获取和设置 (248)4.7图形用户界面(GUI)设计 (251)4.7.1可视化的界面环境 (252)4.7.2菜单 (253)4.7.5回调函数 (258)4.7.6 GUI应用举例 (258)4.8动画 (261)4.8.1以电影方式产生动画 (261)4.8.2以对象方式产生动画 (262)第5章MATLAB程序设计 (264)5.1脚本文件和函数文件 (264)5.1.1 M文本编辑器 (265)5.1.2 M文件的基本格式 (265)5.1.3 M脚本文件 (266)5.1.4 M函数文件 (268)5.2程序流程控制 (269)5.2.1 for ... end循环结构.. (269)5.2.2 while ... end循环结构. (271)5.2.3 If...else...end条件转移结构 (272)5.2.4 switch...case开关结构.. (273)5.2.5 try... catch... end试探结构 .. (275)5.2.6流程控制语句 (276)5.3函数调用和参数传递 (280)5.3.1子函数和私有函数 (280)5.3.2局部变量和全局变量 (281)5.3.4程序举例 (287)5.4 M文件性能的优化和加速 (289)5.4.1 P码文件 (289)5.4.2 M文件性能优化 (290)5.4.3 JIT和加速器 (292)5.5内联函数 (295)5.6利用函数句柄执行函数 (298)5.6.1函数句柄的创建 (298)5.6.2用feval命令执行函数 (299)5.7利用泛函命令进行数值分析 (301)5.7.1求极小值 (302)5.7.2求过零点 (304)5.7.3数值积分 (305)5.7.4微分方程的数值解 (306)第6章线性控制系统分析与设计 (309)6.1线性系统的描述 (309)6.1.1状态空间描述法 (309)6.1.2传递函数描述法 (311)6.1.3零极点描述法 (312)6.1.4离散系统的数学描述 (313)6.2线性系统模型之间的转换 (317)6.2.1连续系统模型之间的转换 (317)6.2.2连续系统与离散系统之间的转换 (325)6.2.3模型对象的属性 (330)6.3结构框图的模型表示 (334)6.4线性系统的时域分析 (348)6.4.1零输入响应分析 (348)6.4.2脉冲响应分析 (350)6.4.3阶跃响应分析 (351)6.4.4任意输入的响应 (353)6.4.5系统的结构参数 (356)6.5线性系统的频域分析 (358)6.5.1频域特性 (358)6.5.2连续系统频域特性 (359)6.5.3幅值裕度和相角裕度 (366)6.5.4离散系统频域分析 (367)6.6线性系统的根轨迹分析 (367)6.6.1绘制根轨迹 (367)6.6.2根轨迹的其它工具 (370)6.7线性系统的状态空间设计 (373)6.7.1极点配置法 (373)6.7.2最优二次型设计 (374)第7章Simulink仿真环境 (376)7.1演示一个Simulink的简单程序 (376)7.2 Simulink的文件操作和模型窗口 (379)7.2.1 Simulink的文件操作 (379)7.2.2 Simulink的模型窗口 (380)7.3 模型的创建 (382)7.3.1模块的操作 (382)7.3.2信号线的操作 (385)7.3.3给模型添加文本注释 (387)7.4 Simulink的基本模块 (387)7.4.1基本模块 (387)7.4.2常用模块的参数和属性设置 (390)7.5复杂系统的仿真与分析 (395)7.5.1仿真的设置 (395)7.5.2连续系统仿真 (397)7.5.3离散系统仿真 (400)7.5.4仿真结构参数化 (402)7.6子系统与封装 (403)7.6.1建立子系统 (403)7.6.2条件执行子系统 (405)7.6.3子系统的封装 (407)7.7用MATLAB命令创建和运行Simulink模型 (413)7.7.1用MATLAB命令创建Simulink模型 (413)7.7.2用MATLAB命令运行Simulink模块 (417)7.8以Simulink为基础的模块工具箱简介 (418)第8章MATLAB高级应用 (419)8.1 MATLAB应用接口 (419)8.1.1 MEX文件 (419)8.1.2 使用MATLAB编译器生成MEX和EXE文件 (424)8.2 低级文件的输入输出 (426)8.2.1打开和关闭文件 (426)8.2.2读写格式化文件 (428)8.2.3读写二进制数据 (433)8.2.4文件定位 (435)8.3 图形文件的转储 (437)8.4 Notebook (438)8.4.1 Notebook的安装 (438)8.4.2 Notebook的启动 (438)8.4.3 Notebook的使用 (441)8.4.4 Notebook中MATLAB的使用 (445)第1章MATLAB 6.5环境1.1MATLAB简介●MATLAB(Matrix Laborator)是MathWorks公司开发科学与工程计算软件;●广泛应用于自动控制、数学运算、信号分析、计算机技术、图像信号处理、财务分析、航天工业、汽车工业、生物医学工程、语音处理和雷达工程等行业;●国内外高校和研究部门科学研究的重要工具;●MATLIB 已成为数学计算工具方面事实上的标准,MATLIB 6.5是最新版本。
MatLab简介MATLAB是什么?典型的使用包括:数学和计算算术发展模型,模拟,和原型数据分析,开发,和可视化科学和工程图学应用发展包括图形用户界面设计MATLAB表示矩阵实验室。
MATLAB系统MATLAB系统由5主要的部分构成:1. MATLAB语言。
这是高阶的矩阵/数组语言,带控制流动陈述,函数,数据结构,输入/输出,而且面向目标的编程特点。
Ops 操作符和特殊字符。
Lang 程序设计语言作。
strfun 字符串。
iofun 输入/输出。
timefun 时期和标有日期。
datatypes数据类型和结构。
2. MATLAB工作环境。
这是你作为MATLAB用户或程序编制员的一套工具和设施。
3. 制图这是MATLAB制图系统。
它为2维上,而且三维的数据可视化,图象处理,动画片制作和表示图形包括高阶的指令在内。
它也为包括低阶的指令在内,允许你建造完整的图形用户界面(GUIs),MATLAB应用。
制图法功能在MATLAB工具箱中被组织成5文件夹:graph2d 2-的维数上的图表。
graph3d 三维的图表。
specgraph 专业化图表。
graphics 制图法。
uitools 图形用户界面工具。
4. MATLAB的数学的函数库。
数学和分析的功能在MATLAB工具箱中被组织成8文件夹。
elmat 初步矩阵,和矩阵操作。
elfun 初步的数学函数。
specfun 专门的数学函数。
matfun 矩阵函数-用数字表示的线性的代数。
datafun 数据分析和傅立叶变换。
polyfun 插入物,并且多项式。
funfun 功能函数。
sparfun 稀少矩阵。
5. MATLAB应用程序接口(API)。
这是允许你写C、Fortran语言与MATLAB 交互。
关于 SimulinkSimulink ? MATLAB为做非线性的动态的系统的模拟实验的交互式的系统。
它是允许你通过把方框图拉到屏幕,灵活地窜改它制作系统的模型的用图表示的鼠标驱动的程序。
实时工作室?允许你产生来自你的图表块的C代码,使之能用于各种实时系统。
关于工具箱工具箱是为了解答特别种类的问题扩展MATLAB环境的MATLAB函数的综合的(M-文件)收集MatLab工作环境命令窗口若输入A = [1 2 3; 4 5 6; 7 8 10]按下回车键后显示如下A =1 2 34 5 67 8 10清除命令窗口clc这并不清除工作间,只是清除了显示,仍可按上箭头看到以前发出的命令数据格式命令x = [4/3 1.2345e–6]format short1.3333 0.0000format short e1.3333e+000 1.2345e–006format short g1.3333 1.2345e–006format long1.33333333333333 0.00000123450000format long e1.333333333333333e+000 1.234500000000000e–006format long g1.33333333333333 1.2345e–006format bank1.33 0.00format +++format rat4/3 1/810045format hex3ff5555555555555 3eb4b6231abfd271若最大的元素大于1000或小于0.001,则显示short或long格式时时会加上一个比例还有两个格式:format compactformat loose禁止结果的显示在命令后加上分号,则屏幕上不会立即显示出结果,这在运算大的数据量时十分有用,如下命令产生100*100的幻方:A = magic(100);长命令行如想另起一行输入命令,在末尾加上"..."即可,如:s = 1 – 1/2 + 1/3 – 1/4 + 1/5 – 1/6 + 1/7 ... – 1/8 + 1/9 – 1/10 + 1/11 – 1/12;MatLab工作间你可用who或whos来察看当前工作间中有哪些变量,如: whosName Size Bytes ClassA 4x4 128 double arrayD 3x5 120 double arrayM 10x1 40 cell arrayS 1x3 628 struct arrayh 1x11 22 char arrayn 1x1 8 double arrays 1x5 10 char arrayv 1x14 28 char arrayGrand total is 93 elements using 984 bytes若要从工作间中删除所有的变量,用clear保存、重载工作间你可以将工作间保存为一个二进制的M文件,以后还可以恢复回来:save june10也可只保存工作间中的部分变量值save june10 x y z重载时只需输入load june10文件名保存在字符串中这样可以像调用函数一样调用工作间save(’myfile’,’VAR1’,’VAR2’)A = ’myfile’;load(A)与下面的命令相同save myfile VAR1 VAR2load myfile下面的命令把1至10的平方值分别存放在data1至data10中:file = ’data’;for i = 1:10 j = i.^2;save([file int2str(i)],’j’);end查找路径当你输入“yourpig"时发生了什么呢?1:察看是否是变量;2:察看是否是内建函数;3:察看当前目录下是否有文件:yourpig.m;4:察看查找目录下是否有文件:yourpig.m;对于查找路径中的文件,what显示当前目录下的文件,加上路径后可显示输入的路径下所有的MatLab文件.如:what matlab/elfun以下二命令分别显示、编辑m文件type rankedit rank图像窗口下面的命令产生一个与命令窗口隔离的图形窗口,figureplot函数则会在新的窗口中绘制图形,如t = 0:pi/100:2*pi;y = sin(t);plot(t,y)则有如下图形:寻求帮助下面的函数在寻求帮助时十分有用:help 列出你所寻求帮助的函数的功能描述;lookfor 列出所有函数的功能描述中含有你所输入的内容的函数的简介如:help inverse显示inverse.m not found.但如输入lookfor inverse则显示INVHILB Inverse Hilbert matrixACOSH Inverse hyperbolic cosineERFINV Inverse of the error functionINV Matrix inversePINV PseudoinverseIFFT Inverse discrete Fourier transformIFFT2 Two–dimensional inverse discrete Fourier transform ICCEPS Inverse complex cepstrumIDCT Inverse discrete cosine transform数据分析和统计面向列的数据集这年头似乎十分风行”面向”这个词,这儿故也套用,其英文为"Column-Oriented Data Sets",可理解为MatLab按列的存储方式来分析数据,下面是一个例子:Time Location 1 Location 2 Location 301h00 11 11 902h00 7 13 1103h00 14 17 2004h00 11 13 905h00 43 51 6906h00 38 46 76 07h00 61 132 186 08h00 75 135 180 09h00 38 88 115 10h00 28 36 55 11h00 12 12 14 12h00 18 27 30 13h00 18 19 29 14h00 17 15 18 15h00 19 36 48 16h00 32 47 10 17h00 42 65 92 18h00 57 66 151 19h00 44 55 90 20h00 114 145 257 21h00 35 58 68 22h00 11 12 15 23h00 13 9 15 24h00 10 9 7 以上数据被保存在一个称为count.dat的文件中.7 13 1114 17 2011 13 943 51 6938 46 7661 132 186 75 135 180 38 88 115 28 36 5512 12 1418 27 3018 19 2917 15 1819 36 4832 47 1042 65 9257 66 151 44 55 90 114 145 257 35 58 6811 12 1510 9 7下面,我们调入此文件,并看看文件的一些参数load count.dat[n,p] = size(count)n =24p =3创建一个时间轴后,我们可以把图画出来:t = 1:n;set(0,'defaultaxeslinestyleorder’,’-|--|-.’)set(0,'defaultaxescolororder’,[0 0 0])plot(t,count), legend('Location 1','Location 2','Location 3',0) xlabel('Time'), ylabel('Vehicle Count'), grid on足以证明,以上是对3个对象的24次观测.基本数据分析函数(一定注意是面向列的)继续用上面的数据,其每列最大值.均值.及偏差分别为:mx = max(count)mu = mean(count)sigma = std(count)mx =114 145 257mu =32.0000 46.5417 65.5833 sigma =25.3703 41.4057 68.0281重载函数,还可以定位出最大.最小值的位置[mx,indx] = min(count)mx =7 9 7indx =2 23 24试试看,你能看懂下面的命令是干什么的吗?[n,p] = size(count)e = ones(n,1)x = count – e*mu点这看看答案!下面这句命令则找出了整个矩阵的最小值:min(count(:))ans =7协方差及相关系数下面,我们来看看第一列的方差:cov(count(:,1))ans =643.6522cov()函数作用于矩阵,则会计算其协方差矩阵. corrcoef()用于计算相关系数,如:corrcoef(count)ans =1.0000 0.9331 0.95990.9331 1.0000 0.95530.9599 0.9553 1.0000数据的预处理未知数据NaN(Not a Number--不是一个数)被定义为未经定义的算式的结果,如 0/0.在处理数据中,NaN常用来表示未知数据或未能获得的数据.所有与NaN有关的运算其结果都是NaN.a = magic(3);a(2,2) = NaNa =8 1 63 NaN 74 9 2sum(a)ans =15 NaN 15在做统计时,常需要将NaN转化为可计算的数字或去掉,以下是几种方法: 注:判断一个值是否为NaN,只能用 isnan(),而不可用 x==NaN;用此法可以从数据中去掉不相关的数据,看看下面的命令是干什么用的:mu = mean(count);sigma = std(count);[n,p] = size(count)outliers = abs(count —mu(ones(n, 1),:)) > 3*sigma(ones(n, 1),:); nout = sum(outliers)nout =1 0 0count(any(outliers'),:) = [];点这看看答案回归与曲线拟合我们经常需要把观测到的数据表达为函数,假如有如下的对时间的观测:t = [0 .3 .8 1.1 1.6 2.3]’;y = [0.5 0.82 1.14 1.25 1.35 1.40]’;plot(t,y,’o’),grid on多项式回归由图可以看出应该可以用多项式来表达:y=a0+a1*t+a2*t^2系数a0,a1,a2可以由最小平方拟合来确定,这一步可由反除号"\"来完成解下面的三元方程组可得:X = [ones(size(t)) t t.^2]X =1.0000 0 01.0000 0.3000 0.09001.0000 0.8000 0.64001.0000 1.1000 1.21001.0000 1.60002.56001.00002.3000 5.2900a = X\ya =0.5318 0.9191 –0.2387a即为待求的系数,画图比较可得T = (0:0.1:2.5)’;Y = [ones(size(T)) T T.^2]*a;plot(T,Y,'–',t,y,'o',), grid on结果令人失望,但我们可以增加阶数来提高精确度,但更明智的选择是用别的方法.线性参数回归形如:y=a0+a1*exp(-t)+a2*t*exp(-t)计算方法同上:X = [ones(size(t)) exp(– t) t.*exp(– t)];a = X\ya =1.3974 – 0.8988 0.4097T = (0:0.1:2.5)';Y = [ones(size(T)) exp(– T) T.exp(– T)]*a;plot(T,Y,'–',t,y,'o'), grid on看起来是不是好多了!例子研究:曲线拟合下面我们以美国人口普查的数据来研究一下有关曲线拟合的问题(MatLab是别人的,教学文档是别人的,例子也是别人的,我只是一个翻译而已...)load census这样我们得到了两个变量,cdate是1790至1990年的时间列向量(10年一次),pop是相应人口数列向量.上一小节所讲的多项式拟合可以用函数polyfit()来完成,数字指明了阶数p = polyfit(cdate,pop,4)Warning: Matrix is close to singular or badly scaled.Results may be inaccurate. RCOND = 5.429790e–20p =1.0e+05 *0.0000 –0.0000 0.0000 –0.0126 6.0020产生警告的原因是计算中的cdata值太大,在计算中的Vandermonde行列式使变换产生了问题,解决的方法之一是使数据标准化.预处理:标准化数据数据的标准化是对数据进行缩放,以使以后的计算能更加精确,一种方法是使之成为0均值:sdate = (cdate – mean(cdate))./std(cdate)现在再进行曲线拟合就没事了!p = polyfit(sdate,pop,4)p =0.7047 0.9210 23.4706 73.8598 62.2285pop4 = polyval(p,sdate);plot(cdate,pop4,'–',cdate,pop,'+'), grid on在上面的数据标准化中,也可以有别的算法,如令1790年的人口数为0.余量分析p1 = polyfit(sdate,pop,1);pop1 = polyval(p1,sdate);plot(cdate,pop1,'–',cdate,pop,'+')figure, plot(cdate,res1,'+')p = polyfit(sdate,pop,2);pop2 = polyval(p,sdate);plot(cdate,pop2,'–',cdate,pop,'+')figure, plot(cdate,res2,’+’)p = polyfit(sdate,pop,4);pop4 = polyval(p,sdate);plot(cdate,pop4,'–',cdate,pop,'+')figure, plot(cdate,res4,'+')可以看出,多项式拟合即使提高了阶次也无法达到令人满意的结果指数拟合从人口增长图可以发现人数的增长基本是呈指数增加的,因此我们可以用年份的对数来进行拟合,这儿,年数是标准化后的!logp1 = polyfit(sdate,log10(pop),1);logpred1 = 10.^polyval(logp1,sdate);semilogy(cdate,logpred1,'–',cdate,pop,'+');grid onlogres2 = log10(pop) –polyval(logp2,sdate);plot(cdate,logres2,'+')上面的图不令人满意,下面,我们用二阶的对数分析:logp2 = polyfit(sdate,log10(pop),2);logpred2 = 10.^polyval(logp2,sdate);semilogy(cdate,logpred2,'–',cdate,pop,'+');grid onr = pop – 10.^(polyval(logp2,sdate));plot(cdate,r,'+')这种余量分析比多项式拟合的余量分析图案要随机的多(没有很强的规律性),可以预见,随着人数的增加,余粮所反映的不确定度也在增加,但总的说来,这种拟合方式要强好多!误差边界误差边界常用来反映你所用的拟合方式是否适用于数据,为得到误差边界,只需在polyfit()中传递第二个参数,并将其送入polyval().下面是一个二阶多项式拟合模型,年份已被标准化,下面的代码用了2σ,对应于95%的可置信度:[p2,S2] = polyfit(sdate,pop,2);[pop2,del2] = polyval(p2,sdate,S2);plot(cdate,pop,'+',cdate,pop2,'g–',cdate,pop2+2*del2,'r:',...cdate,pop2–2*del2,'r:'),grid on差分方程和滤波MatLab中的差分和滤波基本都是对向量而言的,向量则是存储取样信号或序列的.函数y = filter(b, a, x)将用a,b描述的滤波器处理向量x,然后将其存储在向量y中,filter()函数可看为是一差分方程a1y(n)=b1*x(1)+b2*x(2)+...-a2*y(2)-...如有以下差分方程:y(n)=1/4*x(n)+1/4*x(n-1)+1/4*x(n-2)+1/4*x(n-3),则a = 1;b = [1/4 1/4 1/4 1/4];我们载入数据,取其第一列,并计算有:load count.datx = count(:,1);y = filter(b,a,x);t = 1:length(x);plot(t,x,'–.',t,y,'–'),grid onlegend('Original Data','Smoothed Data',2)实现所表示的就是滤波后的数据,它代表了4小时的平均车流量MatLab的信号处理工具箱中提供了很多用来滤波的函数,可用来处理实际问题!快速傅立叶变换(FFT)傅立叶变换能把信号按正弦展开成不同的频率值,对于取样信号,用的是离散傅立叶变换.FFT是离散傅立叶变换的一种高速算法,在信号和图像处理中有极大的用处!fft离散傅立叶变换fft2二维离散傅立叶变换fftn n维离散傅立叶变换ifft离散傅立叶反变换ifft2二维离散傅立叶反变换ifftn n维离散傅立叶反变换abs幅度angle相角unwrap相位按弧度展开,大于π的变换为2π的补角fftshift把零队列移至功率谱中央cplxpair把数据排成复数对nextpow2下两个更高的功率向量x的FFT可以这样求:x = [4 3 7 –9 1 0 0 0]’y = fft(x)y =6.000011.4853 –2.7574i–2.0000 –12.0000i–5.4853 +11.2426i18.0000–5.4853 –11.2426i–2.0000 +12.0000i11.4853 +2.7574ix虽然是实数,但y是复数,其中,第一个是因为它是常数相加的结果,第五个则对应于奈奎斯特频率,后三个数是由于负频率的影响,它们是前面三个数的共轭值!下面,让我们来验证一下太阳黑子活动周期是11年!Wolfer数记录了300年太阳黑子的数量及大小:load sunspot.datyear = sunspot(:,1);wolfer = sunspot(:,2);plot(year,wolfer)title(’Sunspot Data’)现在来看看其FFT:Y = fft(wolfer);Y的幅度是功率谱,画出功率谱和频率的对应关系就得出了周期图,去掉第一点,因为他只是所有数据的和,画图有:N = length(Y);Y(1) = [];power = abs(Y(1:N/2)).^2;nyquist = 1/2;freq = (1:N/2)/(N/2)*nyquist;plot(freq,power),grid onxlabel(’cycles/year’)title(’Periodogram’)上面的图看起来不大方便,下面我们画出频谱-周期图period = 1./freq;plot(period,power),axis([0 40 0 2e7]),grid onylabel(’Power’)xlabel(’Period(Years/Cycle)’)为了得出精确一点的解,如下:[mp index] = max(power);period(index)ans =11.0769变换后的幅度和相位abs()和angle()是用来计算幅度和相位的先创建一信号,再进行分析,unwarp()把相位大于π的变换为2π的补角: t = 0:1/99:1;x = sin(2*pi*15*t) + sin(2*pi*40*t);y = fft(x);m = abs(y);p = unwrap(angle(y));f = (0:length(y)–1)'*99/length(y);subplot(2,1,1), plot(f,m),ylabel('Abs. Magnitude'), grid onsubplot(2,1,2), plot(f,p*180/pi)ylabel('Phase [Degrees]'), grid onxlabel('Frequency [Hertz]')可以发现幅度曲线关于奈奎斯特频率对称,只有0-50Hz的信息是有用的! FFT的长度与速度可以为FFT加上第二个参数,告诉MatLab这是n点FFT.如y = fft(x,n),若x长度大于n,软件自动补0,否则截取x.若:1. n为2的幂,软件将执行基2快速傅立叶算法,这时的运算速度是最快的2. n为合数,软件将n分解为素数来算,计算量与n的值有关.n为1013将比1000点的速度慢的多!3. n为素数,软件执行DFT的公式,此时最慢矩阵和线性代数MatLab中的矩阵MatLab中有好多函数可以产生不同的矩阵,下面就让我们产生两个3*3的矩阵,这一章中,我们的学习就靠她们了!!!A = pascal(3)A =1 1 11 2 31 3 6B = magic(3)B =8 1 63 5 74 9 2还有一个3*2的随机矩阵:C = fix(10*rand(3,2))C =9 42 86 7看看列矩阵,行矩阵,以及常数的表达:u = [3; 1; 4]v = [2 0 —1]s = 7产生的矩阵是:u =314v =2 0 —1 s =7加减法X = A + BX =9 2 74 7 105 12 8Y = X –AY =8 1 63 5 74 9 2若二矩阵维数不统一,则会出错!X = A + CError using ==> +Matrix dimensions must agree.向量的乘积与转置x = v*ux =2X = u*vX =6 0 —32 0 —18 0 —4X = B'X =8 3 41 5 96 7 2如x与y均是列向量,则x*y无解,但下二表达式却可以:x'*yy'*x称内积或点积.下面的语句产生单位矩阵eye(m,n)若用eye(n)则产生n*n的方阵解线性方程情况一:x = A\ux =10—125又如:X = A\BX =19 –3 —1—17 4 136 0 —6情况二;y是不同时刻t时的观测值: t = [0 .3 .8 1.1 1.6 2.3]';y = [.82 .72 .63 .60 .55 .50]';若函数形式是:y(t)=c1+c2*exp(t); 构造矩阵:E = [ones(size(t)) exp(–t)]E =1.0000 1.00001.0000 0.74081.0000 0.44931.0000 0.33291.0000 0.20191.0000 0.1003则可求得系数c1及c2c = E\yc =0.4760 0.3413表明:y(t)=0.4760+0.3413*exp(t) 画图如下:T = (0:0.1:2.5)';Y = [ones(size(T)) exp(–T)]*c; plot(T,Y,'–',t,y,'o')转置与行列式若A是方阵,且是非奇异的,则:d = det(A)X = inv(A)d =1X =3 —3 1—3 5 —21 —2 1若c不是方阵,则用pinv:X = pinv(C)X =0.1159 —0.0729 0.0171—0.0534 0.1152 0.0418那么我们可以发现,下面3个命令具有同样的功效(A是m*n的矩阵,m>n):x = A\bx = pinv(A)*bx = inv(A’*A)*A’*bLU.RQ.及Cholesky分解MatLab求解线性方程建立在以下三个分解之上:Cholesky分解Guass(高斯)分解正交分解Cholesky分解A=p*p'让我们临时把A变一变:A = pascal(6)A =1 1 1 1 1 11 2 3 4 5 61 3 6 10 15 211 4 10 20 35 561 5 15 35 70 1261 6 21 56 126 252A是二项式系数,每一项是其左方与上方系数之和,求其Cholesky分解系数有:R = chol(A)R =1 1 1 1 1 10 1 2 3 4 50 0 1 3 6 100 0 0 1 4 100 0 0 0 1 50 0 0 0 0 1R认识二项式系数.这样对于线性方程便可化简:A*x = bR'*R*x = bx = R\(R'\b)复杂度由O(n^3)变为O(n^2);LU分解A = L U其中,L时下三角阵,U是上三角阵,如: [L,U] = lu(B)L =1.0000 0 00.3750 0.5441 1.00000.5000 1.0000 0U =8.0000 1.0000 6.00000 8.5000 —1.00000 0 5.2941同样:A*x = b可以解为x = U\(L\b)QR分解正交阵有如下性质:Q'Q = I正交阵的好处在于,她保持了原阵的长度,角度,并且在计算的过程中不会扩大误差.RQ分解如下:A = Q R或 A P = Q R其中,Q是正交阵,R是上三角阵.矩阵的幂与指数若A是方阵,p是正数,则X = A^2X =3 6 106 14 2510 25 46若A是方阵,且是非奇异的,则X=A^(-P)将inv(A) P次方,如:Y = B^(–3)Y =0.0053 —0.0068 0.0018—0.0034 0.0001 0.0036—0.0016 0.0070 —0.0051分数词幂将由A的特征值决定.若是对矩阵的每个元素进行幂,用.^,如X = A.^2A =1 1 11 4 91 9 36sqrtm(A)计算A^(1/2),但要更精确,而sqrt(A)则计算A.^(1/2),是一个元素一个元素的算.dx/dt=Ax,可以表示为x(t)=exp(tA)*x(0);下面来看看如何计算:--expm(A)A =0 —6 —16 2 —16—5 20 —10x0 =111计算如下:X = [];for t = 0:.01:1X = [X expm(t*A)*x0];end作图有:plot3(X(1,:),X(2,:),X(3,:),'–o')特征值Av=λv若L是特阵值矩阵,则特征向量是V: AV=VL;如下:A =0 —6 —16 2 —16—5 20 —10lambda = eig(A)lambda =—3.0710—2.4645+17.6008i—2.4645-17.6008i由exp(λt)可以看出exp(At)(见上小节)若用二参数调用函数eig(),则返回特征向量及特征值矩阵: [V,D] = eig(A)V =—0.8326 —0.1203+ 0.2123i —0.1203– 0.2123i —0.3553 0.4691+ 0.4901i 0.4691– 0.4901i—0.4248 0.6249– 0.2997i 0.6249+ 0.2997iD =—3.0710 0 00 —2.4645+17.6008i 00 0 —2.4645—17.6008i对于下面的矩阵:A =6 12 19—9 —20 —334 9 15V =0.4741 0.4082 —0.4082—0.8127 —0.8165 0.81650.3386 0.4082 —0.4082D =—1.0000 0 00 1.0000 00 0 1.0000可以看出,有二特征值是一样的,其特征向量仅差一个符号,在Symbolic Math Toolbox中提供了Jordan标准型的函数,如下:[X,J] = jordan(A)X =—1.7500 1.5000 2.75003.0000 —3.0000 —3.0000—1.2500 1.5000 1.2500J =—1 0 00 1 10 0 1常微分方程常微分方程(Odinary Differential Equations---ODE)类别函数描述解常微分方程ode45 ode23 ode113 ode15s ode23s ode23t ode23tb常微分方程选项odeset odeset常微分方程输出选项odeplot odephas2 odephas3 odeprint如何表述问题多项式与插值多项式多项式的表达MatLab中用按降幂排列的多项式系数组成的行向量表示多项式,如: p(x)=x^3-2x-5被表示为:p = [1 0 –2 –5];多项式的根r = roots(p)r =2.0946–1.0473 + 1.1359i–1.0473 – 1.1359i根被储存为列向量.若要由方程的根构造多项式,则p2 = poly(r)p2 =1 8.8818e-16 –2 –5多项式估计可以用多项式估计出多项式在某一点的值:polyval(p,5)ans =110同样也可以估计矩阵多项式的值p(X) = X^3 – 2X – 5I,X = [2 4 5; –1 0 3; 7 1 5];Y = polyvalm(p,X)Y =377 179 439111 81 136490 253 639卷积多项式相乘是一个卷积的过程,conv()a = [1 2 3];b = [4 5 6];c = conv(a,b)c =4 13 28 27 18多项式相除是其逆过程,用deconv(): [q,r] = deconv(c,a)q =4 5 6r =0 0 0 0 0多项式曲线逼近polyfit(x,y,n)能用多项式逼近由x,y向量提供的数据,n是其阶数,如: x = [1 2 3 4 5];y = [5.5 43.1 128 290.7 498.4];p = polyfit(x,y,3)p =–0.1917 31.5821 –60.3262 35.3400将图画出x2 = 1:.1:5;y2 = polyval(p,x2);plot(x,y,’o’,x2,y2)grid on分式多项式分解residue()可将分式多项式分解如下:对于下式分解为:b = [–4 8];a = [1 6 8];[r,p,k] = residue(b,a)r =–12 8p =–4 –2k =[]重载此函数可以完成分式多项式相加:[b2,a2] = residue(r,p,k)b2 =–4 8a2 =1 6 8插值插值是在已知的数据列中,估计别点的函数值.一维插值一维插值在MatLab中有两种方法:@ 多项式插值@ 建立在FFT上的插值。