Matlab中文教程培训讲学.doc

《Matlab中文教程培训讲学.doc》由会员分享，可在线阅读，更多相关《Matlab中文教程培训讲学.doc（81页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、Good is good, but better carries it.精益求精，善益求善。Matlab中文教程-MatLab简介MATLAB是什么？典型的使用包括：数学和计算算术发展模型，模拟，和原型数据分析，开发，和可视化科学和工程图学应用发展包括图形用户界面设计MATLAB表示矩阵实验室。MATLAB系统MATLAB系统由5主要的部分构成：1.MATLAB语言。这是高阶的矩阵/数组语言,带控制流动陈述，函数，数据结构，输入/输出，而且面向目标的编程特点。Ops操作符和特殊字符。Lang程序设计语言作。strfun字符串。iofun输入/输出。timefun时期和标有日期。datatype

2、s数据类型和结构。2.MATLAB工作环境。这是你作为MATLAB用户或程序编制员的一套工具和设施。3.制图这是MATLAB制图系统。它为2维上，而且三维的数据可视化，图象处理，动画片制作和表示图形包括高阶的指令在内。它也为包括低阶的指令在内,允许你建造完整的图形用户界面（GUIs），MATLAB应用。制图法功能在MATLAB工具箱中被组织成5文件夹：graph2d2-的维数上的图表。graph3d三维的图表。specgraph专业化图表。graphics制图法。uitools图形用户界面工具。4.MATLAB的数学的函数库。数学和分析的功能在MATLAB工具箱中被组织成8文件夹。elmat初

3、步矩阵，和矩阵操作。elfun初步的数学函数。specfun专门的数学函数。matfun矩阵函数用数字表示的线性的代数。datafun数据分析和傅立叶变换。polyfun插入物，并且多项式。funfun功能函数。sparfun稀少矩阵。5.MATLAB应用程序接口（API）。这是允许你写C、Fortran语言与MATLAB交互。关于SimulinkSimulink？MATLAB为做非线性的动态的系统的模拟实验的交互式的系统。它是允许你通过把方框图拉到屏幕，灵活地窜改它制作系统的模型的用图表示的鼠标驱动的程序。实时工作室？允许你产生来自你的图表块的C代码，使之能用于各种实时系统。关于工具箱工具箱

4、是为了解答特别种类的问题扩展MATLAB环境的MATLAB函数的综合的（M-文件）收集MatLab工作环境命令窗口若输入A=123;456;7810按下回车键后显示如下A=1234567810清除命令窗口clc这并不清除工作间，只是清除了显示，仍可按上箭头看到以前发出的命令数据格式命令x=4/31.2345e6formatshort1.33330.0000formatshorte1.3333e+0001.2345e006formatshortg1.33331.2345e006formatlong1.333333333333330.00000123450000formatlonge1.33333

5、3333333333e+0001.234500000000000e006formatlongg1.333333333333331.2345e006formatbank1.330.00format+formatrat4/31/810045formathex3ff55555555555553eb4b6231abfd271若最大的元素大于1000或小于0.001,则显示short或long格式时时会加上一个比例还有两个格式:formatcompactformatloose禁止结果的显示在命令后加上分号,则屏幕上不会立即显示出结果,这在运算大的数据量时十分有用,如下命令产生100*100的幻方:A=m

6、agic(100);长命令行如想另起一行输入命令,在末尾加上.即可,如:s=11/2+1/31/4+1/51/6+1/7.1/8+1/91/10+1/111/12;MatLab工作间你可用who或whos来察看当前工作间中有哪些变量,如:whosNameSizeBytesClassA4x4128doublearrayD3x5120doublearrayM10x140cellarrayS1x3628structarrayh1x1122chararrayn1x18doublearrays1x510chararrayv1x1428chararrayGrandtotalis93elementsusin

7、g984bytes若要从工作间中删除所有的变量,用clear保存、重载工作间你可以将工作间保存为一个二进制的M文件，以后还可以恢复回来：savejune10也可只保存工作间中的部分变量值savejune10xyz重载时只需输入loadjune10文件名保存在字符串中这样可以像调用函数一样调用工作间save(myfile,VAR1,VAR2)A=myfile;load(A)与下面的命令相同savemyfileVAR1VAR2loadmyfile下面的命令把1至10的平方值分别存放在data1至data10中：file=data;fori=1:10j=i.2;save(fileint2str(i)

8、,j);end查找路径当你输入“yourpig时发生了什么呢?1:察看是否是变量;2:察看是否是内建函数;3:察看当前目录下是否有文件:yourpig.m;4:察看查找目录下是否有文件:yourpig.m;对于查找路径中的文件,what显示当前目录下的文件,加上路径后可显示输入的路径下所有的MatLab文件.如:whatmatlab/elfun以下二命令分别显示、编辑m文件typerankeditrank图像窗口下面的命令产生一个与命令窗口隔离的图形窗口,figureplot函数则会在新的窗口中绘制图形,如t=0:pi/100:2*pi;y=sin(t);plot(t,y)则有如下图形:寻求帮

9、助下面的函数在寻求帮助时十分有用:help列出你所寻求帮助的函数的功能描述；lookfor列出所有函数的功能描述中含有你所输入的内容的函数的简介如:helpinverse显示inverse.mnotfound.但如输入lookforinverse则显示INVHILBInverseHilbertmatrixACOSHInversehyperboliccosineERFINVInverseoftheerrorfunctionINVMatrixinversePINVPseudoinverseIFFTInversediscreteFouriertransformIFFT2Twodimensionali

10、nversediscreteFouriertransformICCEPSInversecomplexcepstrumIDCTInversediscretecosinetransform数据分析和统计面向列的数据集这年头似乎十分风行”面向”这个词，这儿故也套用，其英文为Column-OrientedDataSets,可理解为MatLab按列的存储方式来分析数据,下面是一个例子:TimeLocation1Location2Location301h001111902h007131103h0014172004h001113905h0043516906h0038467607h006113218608h0

11、07513518009h00388811510h0028365511h0012121412h0018273013h0018192914h0017151815h0019364816h0032471017h0042659218h00576615119h0044559020h0011414525721h0035586822h0011121523h001391524h001097以上数据被保存在一个称为count.dat的文件中1111971311141720111394351693846766113218675135180388811528365512121418273018192917151819

12、36483247104265925766151445590114145257355868111215139151097下面，我们调入此文件，并看看文件的一些参数loadcount.datn,p=size(count)n=24p=3创建一个时间轴后，我们可以把图画出来：t=1:n;set(0,defaultaxeslinestyleorder,-|-|-.)set(0,defaultaxescolororder,000)plot(t,count),legend(Location1,Location2,Location3,0)xlabel(Time),ylabel(VehicleCount),gr

13、idon足以证明，以上是对个对象的次观测基本数据分析函数（一定注意是面向列的）继续用上面的数据，其每列最大值均值及偏差分别为：mx=max(count)mu=mean(count)sigma=std(count)mx=114145257mu=32.000046.541765.5833sigma=25.370341.405768.0281重载函数，还可以定位出最大最小值的位置mx,indx=min(count)mx=797indx=22324试试看，你能看懂下面的命令是干什么的吗？n,p=size(count)e=ones(n,1)x=counte*mu点这看看答案！下面这句命令则找出了整个矩阵

14、的最小值:min(count(:)ans=7协方差及相关系数下面,我们来看看第一列的方差:cov(count(:,1)ans=643.6522cov()函数作用于矩阵,则会计算其协方差矩阵.corrcoef()用于计算相关系数,如:corrcoef(count)ans=1.00000.93310.95990.93311.00000.95530.95990.95531.0000数据的预处理未知数据NaN(NotaNumber-不是一个数)被定义为未经定义的算式的结果,如0/0.在处理数据中,NaN常用来表示未知数据或未能获得的数据.所有与NaN有关的运算其结果都是NaN.a=magic(3);a

15、(2,2)=NaNa=8163NaN7492sum(a)ans=15NaN15在做统计时，常需要将NaN转化为可计算的数字或去掉,以下是几种方法:注:判断一个值是否为NaN,只能用isnan(),而不可用x=NaN;i=find(isnan(x);x=x(i)先找出值不是NaN的项的下标,将这些元素保留x=x(find(isnan(x)同上,去掉NaNx=x(isnan(x);更快的做法x(isnan(x)=;消掉NaNX(any(isnan(X),:)=;把含有NaN的行都去掉用此法可以从数据中去掉不相关的数据,看看下面的命令是干什么用的:mu=mean(count);sigma=std(c

16、ount);n,p=size(count)outliers=abs(countmu(ones(n,1),:)3*sigma(ones(n,1),:);nout=sum(outliers)nout=100count(any(outliers),:)=;点这看看答案回归与曲线拟合我们经常需要把观测到的数据表达为函数，假如有如下的对时间的观测：t=0.3.81.11.62.3;y=0.50.821.141.251.351.40;plot(t,y,o),gridon多项式回归由图可以看出应该可以用多项式来表达：y=a0+a1*t+a2*t2系数a0,a1,a2可以由最小平方拟合来确定，这一步可由反除号

17、来完成解下面的三元方程组可得：X=ones(size(t)tt.2X=1.0000001.00000.30000.09001.00000.80000.64001.00001.10001.21001.00001.60002.56001.00002.30005.2900a=Xya=0.53180.91910.2387a即为待求的系数，画图比较可得T=(0:0.1:2.5);Y=ones(size(T)TT.2*a;plot(T,Y,t,y,o,),gridon结果令人失望，但我们可以增加阶数来提高精确度，但更明智的选择是用别的方法线性参数回归形如：y=a0+a1*exp(-t)+a2*t*exp(

18、-t)计算方法同上：X=ones(size(t)exp(t)t.*exp(t);a=Xya=1.39740.89880.4097T=(0:0.1:2.5);Y=ones(size(T)exp(T)T.exp(T)*a;plot(T,Y,t,y,o),gridon看起来是不是好多了！例子研究：曲线拟合下面我们以美国人口普查的数据来研究一下有关曲线拟合的问题（MatLab是别人的,教学文档是别人的,例子也是别人的，我只是一个翻译而已）loadcensus这样我们得到了两个变量，cdate是1790至1990年的时间列向量(10年一次),pop是相应人口数列向量.上一小节所讲的多项式拟合可以用函数p

19、olyfit()来完成,数字指明了阶数p=polyfit(cdate,pop,4)Warning:Matrixisclosetosingularorbadlyscaled.Resultsmaybeinaccurate.RCOND=5.429790e20p=1.0e+05*0.00000.00000.00000.01266.0020产生警告的原因是计算中的cdata值太大，在计算中的Vandermonde行列式使变换产生了问题，解决的方法之一是使数据标准化预处理：标准化数据数据的标准化是对数据进行缩放，以使以后的计算能更加精确，一种方法是使之成为0均值：sdate=(cdatemean(cdat

20、e)./std(cdate)现在再进行曲线拟合就没事了！p=polyfit(sdate,pop,4)p=0.70470.921023.470673.859862.2285pop4=polyval(p,sdate);plot(cdate,pop4,cdate,pop,+),gridon在上面的数据标准化中，也可以有别的算法，如令年的人口数为余量分析p1=polyfit(sdate,pop,1);pop1=polyval(p1,sdate);plot(cdate,pop1,cdate,pop,+)res1=poppop1;figure,plot(cdate,res1,+)p=polyfit(sda

21、te,pop,2);pop2=polyval(p,sdate);plot(cdate,pop2,cdate,pop,+)res2=poppop2;figure,plot(cdate,res2,+)p=polyfit(sdate,pop,4);pop4=polyval(p,sdate);plot(cdate,pop4,cdate,pop,+)res4=poppop4;figure,plot(cdate,res4,+)可以看出，多项式拟合即使提高了阶次也无法达到令人满意的结果指数拟合从人口增长图可以发现人数的增长基本是呈指数增加的，因此我们可以用年份的对数来进行拟合，这儿，年数是标准化后的！log

22、p1=polyfit(sdate,log10(pop),1);logpred1=10.polyval(logp1,sdate);semilogy(cdate,logpred1,cdate,pop,+);gridonlogres2=log10(pop)polyval(logp2,sdate);plot(cdate,logres2,+)上面的图不令人满意，下面，我们用二阶的对数分析：logp2=polyfit(sdate,log10(pop),2);logpred2=10.polyval(logp2,sdate);semilogy(cdate,logpred2,cdate,pop,+);grido

23、nr=pop10.(polyval(logp2,sdate);plot(cdate,r,+)这种余量分析比多项式拟合的余量分析图案要随机的多（没有很强的规律性），可以预见，随着人数的增加，余粮所反映的不确定度也在增加，但总的说来，这种拟合方式要强好多！误差边界误差边界常用来反映你所用的拟合方式是否适用于数据，为得到误差边界，只需在polyfit()中传递第二个参数，并将其送入polyval().下面是一个二阶多项式拟合模型,年份已被标准化,下面的代码用了2,对应于95%的可置信度:p2,S2=polyfit(sdate,pop,2);pop2,del2=polyval(p2,sdate,S2)

24、;plot(cdate,pop,+,cdate,pop2,g,cdate,pop2+2*del2,r:,.cdate,pop22*del2,r:),gridon差分方程和滤波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/41/41/41/4;

25、我们载入数据,取其第一列,并计算有:loadcount.datx=count(:,1);y=filter(b,a,x);t=1:length(x);plot(t,x,.,t,y,),gridonlegend(OriginalData,SmoothedData,2)实现所表示的就是滤波后的数据,它代表了4小时的平均车流量MatLab的信号处理工具箱中提供了很多用来滤波的函数,可用来处理实际问题!快速傅立叶变换(FFT)傅立叶变换能把信号按正弦展开成不同的频率值,对于取样信号,用的是离散傅立叶变换.FFT是离散傅立叶变换的一种高速算法,在信号和图像处理中有极大的用处!fft离散傅立叶变换fft2二

26、维离散傅立叶变换fftnn维离散傅立叶变换ifft离散傅立叶反变换ifft2二维离散傅立叶反变换ifftnn维离散傅立叶反变换abs幅度angle相角unwrap相位按弧度展开,大于的变换为2的补角fftshift把零队列移至功率谱中央cplxpair把数据排成复数对nextpow2下两个更高的功率向量x的FFT可以这样求:x=43791000y=fft(x)y=6.000011.48532.7574i2.000012.0000i5.4853+11.2426i18.00005.485311.2426i2.0000+12.0000i11.4853+2.7574ix虽然是实数,但y是复数,其中,第

27、一个是因为它是常数相加的结果,第五个则对应于奈奎斯特频率,后三个数是由于负频率的影响,它们是前面三个数的共轭值!下面,让我们来验证一下太阳黑子活动周期是11年!Wolfer数记录了300年太阳黑子的数量及大小:loadsunspot.datyear=sunspot(:,1);wolfer=sunspot(:,2);plot(year,wolfer)title(SunspotData)现在来看看其FFT:Y=fft(wolfer);Y的幅度是功率谱,画出功率谱和频率的对应关系就得出了周期图,去掉第一点,因为他只是所有数据的和,画图有:N=length(Y);Y(1)=;power=abs(Y(1

28、:N/2).2;nyquist=1/2;freq=(1:N/2)/(N/2)*nyquist;plot(freq,power),gridonxlabel(cycles/year)title(Periodogram)上面的图看起来不大方便,下面我们画出频谱-周期图period=1./freq;plot(period,power),axis(04002e7),gridonylabel(Power)xlabel(Period(Years/Cycle)为了得出精确一点的解,如下:mpindex=max(power);period(index)ans=11.0769变换后的幅度和相位abs()和angl

29、e()是用来计算幅度和相位的先创建一信号,再进行分析,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),gridonsubplot(2,1,2),plot(f,p*180/pi)ylabel(PhaseDegrees),gridonxlabel(FrequencyHertz)可以发现幅度曲线关

30、于奈奎斯特频率对称,只有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=111123136B

31、=magic(3)B=816357492还有一个3*2的随机矩阵:C=fix(10*rand(3,2)C=942867看看列矩阵，行矩阵，以及常数的表达：u=3;1;4v=201s=7产生的矩阵是：u=314v=201s=7加减法X=A+BX=92747105128Y=XAY=816357492若二矩阵维数不统一，则会出错！X=A+CErrorusing=+Matrixdimensionsmustagree.向量的乘积与转置x=v*ux=2X=u*vX=603201804X=BX=834159672如与均是列向量，则x*y无解,但下二表达式却可以:x*yy*x称内积或点积.下面的语句产生单位矩

32、阵eye(m,n)若用eye(n)则产生n*n的方阵解线性方程情况一:x=Aux=10125又如：X=ABX=193117413606情况二；y是不同时刻t时的观测值:t=0.3.81.11.62.3;y=.82.72.63.60.55.50;若函数形式是:y(t)=c1+c2*exp(t);构造矩阵:E=ones(size(t)exp(t)E=1.00001.00001.00000.74081.00000.44931.00000.33291.00000.20191.00000.1003则可求得系数c1及c2c=Eyc=0.47600.3413表明：y(t)=0.4760+0.3413*exp

33、(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=331352121若c不是方阵,则用pinv:X=pinv(C)X=0.11590.07290.01710.05340.11520.0418那么我们可以发现，下面个命令具有同样的功效（是*的矩阵，）：x=Abx=pinv(A)*bx=inv(A*A)*A*bLU.RQ.及Cholesky分解MatLab求解线性方程建立在以下三个分解之上:Cholesky分解Guass(高斯)分解正交分解Cho

34、lesky分解A=p*p让我们临时把A变一变:A=pascal(6)A=111111123456136101521141020355615153570126162156126252A是二项式系数,每一项是其左方与上方系数之和,求其Cholesky分解系数有:R=chol(A)R=11111101234500136100001410000015000001认识二项式系数这样对于线性方程便可化简：A*x=bR*R*x=bx=R(Rb)复杂度由O(n3)变为O(n2);LU分解A=LU其中,L时下三角阵,U是上三角阵,如:L,U=lu(B)L=1.0000000.37500.54411.00000.

35、50001.00000U=8.00001.00006.000008.50001.0000005.2941同样：A*x=b可以解为x=U(Lb)QR分解正交阵有如下性质:QQ=I正交阵的好处在于,她保持了原阵的长度,角度,并且在计算的过程中不会扩大误差.RQ分解如下:A=QR或AP=QR其中,Q是正交阵,R是上三角阵.矩阵的幂与指数若A是方阵,p是正数,则X=A2X=361061425102546若是方阵，且是非奇异的，则X=A(-P)将inv(A)P次方,如:Y=B(3)Y=0.00530.00680.00180.00340.00010.00360.00160.00700.0051分数词幂将由

36、的特征值决定若是对矩阵的每个元素进行幂，用.,如X=A.2A=1111491936sqrtm(A)计算A(1/2)，但要更精确，而sqrt(A)则计算A.(1/2),是一个元素一个元素的算dx/dt=Ax,可以表示为x(t)=exp(tA)*x(0);下面来看看如何计算:-expm(A)A=061621652010x0=111计算如下：X=;fort=0:.01:1X=Xexpm(t*A)*x0;end作图有：plot3(X(1,:),X(2,:),X(3,:),o)特征值Av=v若L是特阵值矩阵,则特征向量是V:AV=VL;如下:A=061621652010lambda=eig(A)lamb

37、da=3.07102.4645+17.6008i2.4645-17.6008i由exp(t)可以看出exp(At)(见上小节)若用二参数调用函数eig(),则返回特征向量及特征值矩阵:V,D=eig(A)V=0.83260.1203+0.2123i0.12030.2123i0.35530.4691+0.4901i0.46910.4901i0.42480.62490.2997i0.6249+0.2997iD=3.07100002.4645+17.6008i0002.464517.6008i对于下面的矩阵：A=61219920334915V=0.47410.40820.40820.81270.81

38、650.81650.33860.40820.4082D=1.00000001.00000001.0000可以看出，有二特征值是一样的，其特征向量仅差一个符号，在SymbolicMathToolbox中提供了Jordan标准型的函数,如下:X,J=jordan(A)X=1.75001.50002.75003.00003.00003.00001.25001.50001.2500J=100011001常微分方程常微分方程(OdinaryDifferentialEquations-ODE)类别函数描述解常微分方程ode45ode23ode113ode15sode23sode23tode23tb常微分方

39、程选项odesetodeset常微分方程输出选项odeplotodephas2odephas3odeprint如何表述问题多项式与插值多项式多项式的表达MatLab中用按降幂排列的多项式系数组成的行向量表示多项式,如:p(x)=x3-2x-5被表示为:p=1025;多项式的根r=roots(p)r=2.09461.0473+1.1359i1.04731.1359i根被储存为列向量若要由方程的根构造多项式，则p2=poly(r)p2=18.8818e-1625多项式估计可以用多项式估计出多项式在某一点的值：polyval(p,5)ans=110同样也可以估计矩阵多项式的值p(X)=X32X5I,

40、X=245;103;715;Y=polyvalm(p,X)Y=37717943911181136490253639卷积多项式相乘是一个卷积的过程,conv()a=123;b=456;c=conv(a,b)c=413282718多项式相除是其逆过程，用deconv():q,r=deconv(c,a)q=456r=00000多项式曲线逼近polyfit(x,y,n)能用多项式逼近由x,y向量提供的数据,n是其阶数,如:x=12345;y=5.543.1128290.7498.4;p=polyfit(x,y,3)p=0.191731.582160.326235.3400将图画出x2=1:.1:5;y2=polyval(p,x2);plot(x,y,o,x2,y2)gridon分式多项式分解residue()可将分式多项式分解如下:对于下式分解为:b=48