MATLAB例程.doc

《MATLAB例程.doc》由会员分享，可在线阅读，更多相关《MATLAB例程.doc（13页珍藏版）》请在得力文库 - 分享文档赚钱的网站上搜索。

1、【精品文档】如有侵权，请联系网站删除，仅供学习与交流MATLAB例程.精品文档.Matlab Resources： http:/www.people.cornell.edu/pages/jf262/matlab_main.htm 1. Matlab at MathWorks Full documentation by the developer 2. MatLab Tutorial http:/users.ece.gatech.edu/%7Ebonnie/book/TUTORIAL/tutorial.html This tutorial is available as a su

2、pplement to the textbook Fundamentals of Signals and Systems Using Matlab by Edward Kamen and Bonnie Heck, published by Prentice Hall. Also available as pdf. http:/www.people.cornell.edu/pages/jf262/matlab/useful%20pages/tutorial.pdf 3. MatLab graphics http:/www.people.cornell.edu/pages/jf262/matlab

3、/useful%20pages/Matlab%20graphics.htm Tomi Salminen, CSC 1998. An introduction to Matlab graphics commands with examples and figures. 4. Numerical Computing with MATLAB http:/www.people.cornell.edu/pages/jf262/matlab/ebook1/ncm_index.htm 5 . MatLab vector and matrix tips (pdf) http:/www.people.corne

4、ll.edu/pages/jf262/matlab/useful%20pages/mtt.pdf 6. Matlab Programing: Design & Application (Chinese) http:/www.people.cornell.edu/pages/jf262/matlab/matlablearn/index.htm1、 figure 1） h=bar(clrflux_up10(:,2),clrflux_up10(:,3),Edgecolor,k,Facecolor,k) hold on bar(pollutionflux_up10(:,2),pollutionflux

5、_up10(:,3),Edgecolor,k,Facecolor,w) axis(400 700 0 15); box off xlabel(TOA flux ( Wm-2 ) legend(AODleq0.2,AOD0.5) axis(400 700 0 15); set(gca,TickLength,0.015;0.015) set(gca,TickDir,out);set(gca,position,0.15 0.3 0.7 0.4)plot(DATEnum,BJ2005_MISR_AER(:,4),:ks,LineWidth,1,MarkerEdgeColor,m, MarkerFace

6、Color,m,MarkerSize,5)hold on plot(DATEnum,BJ2005_MISR_AER(:,8),:ko,LineWidth,1,MarkerEdgeColor,b, MarkerFaceColor,b,MarkerSize,5) 附加绘图程序：date=070427;station=bj;n = size(date,1);c=strcat(station,date);a=load(D:BJDOASbj2,c,.dat);figure(Position,1,1,800,800);plot(1:288,a(:,6),*-r,LineWidth,2.5);hold on

7、plot(1:288,a(:,8),o-b,LineWidth,2.5);plot(1:288,a(:,10),-k,LineWidth,2.5);hold offYLabel(fontsize16Mixing Ratio(ppb);XLabel(fontsize16Hours);Title(fontsize14,date(1:2),/,date(3:4),/,date(5:6);Legend(fontsize14SO_2,fontsize14NO_2,fontsize14O_3);legend(boxoff);XLim(0.5,288.5);%YLim(0,max(a(:,6),a(:,8)

8、,a(:,10)*1.2);set(gca,xtick,0.5:24:288.5,xticklabel,0:2:24,fontsize,12);saveas(gcf,D:BJDOASbj5,c,.emf);close;2）文本标注 Ht = text(0,0,fontsize12alpha=0.3) set (ht,Color,k,HorizontalAlignment,Center)3) 双坐标轴绘制 I）粗略绘制 set(gca,TickLength,0.015;0.015) set(gca,position,0.15 0.6 0.5 0.2);axis(0,4000,-15,20,0,1

9、00); AX,H1,H2= plotyy(upair2(1:7,3),upair2(1:7,4),upair2(1:7,3),upair2(1:7,6) set(H1,Color,k);%set(H1,XMinorTick,on);axes(Position,get(AX(1),Position),XColor,k,YColor,k)set(H2,Color,k,LineStyle,:); title(10-02)II）精确绘制 tp: 压力变量 tt：温度变量 clf reset , h_ap = axes(position,0.13,0.13,0.7,0.75);set(h_ap,Xco

10、lor,b,Ycolor,b,Xlim,0,5,Ylim,0,15);nx=10 ; ny=6; % 为两套坐标系定义共同的坐标刻度分段数。pxtick = 0:(5-0)/nx:5; pytick = 0(15-0)/ny):15;set(h_ap,Xtick, pxtick, Ytick, pytick,Xgrid, on, Ygrid, on) % 保证坐标系按约定分段数刻度；h_linet = line(tp,yp,Color,b);set(get(h_ap,Xlabel),String,时间rightarrow分)set(get(h_ap,Ylabel),String,压力right

11、arrow（times105pa）)h_at = axes(position,get(h_ap,Position); %保证两套坐标系重合的技巧。set(h_at,Color,none,Xcolor,r, Ycolor,r);保证重合坐标系都可见的技术，使得后面产生的h_at轴背景“无色透明,而不 遮盖前面的图像。set(h_at,Xaxislocation,top) % 把第二套坐标系的横轴设置在轴位框的“顶边”，把纵轴设置在轴位的“右% 边”，而且增值%的方向是自上而下的set(h_at,Yaxislocation,right,Ydir,rev)set(get(h_at,Xlabel),S

12、tring,时间rightarrow分)set(get(h_at,Ylabel),String,压力rightarrow（times105pa）)set(h_at,Ylim,0,210)line(tt,yt,Color,r,Parent,h_at); % 确保温度曲线绘制在h_at坐标中xpm = get(h_at,Xlim); %获取h_at坐标系横轴的取值范围。txtick = xpm(1):(xpm(2)-xpm(1)/nx):xpm(2);tytick = 0:(210-0)/ny):210;set(h_ap,Xtick, txtick, Ytick, tytick)4）绘制多个图形，

13、在一张图上ff=figure(ss) h1 = subplot(3,1,1) plot(modis_spr(:,1),modis_spr(:,3),g-) xlim(2000.5 2006.5) ylim(0 1.5) hold on plot(modis_sum(:,1),modis_sum(:,3),r-o) hold on plot(modis_aut(:,1),modis_aut(:,3),y-) hold on plot(modis_win(:,1),modis_win(:,3),b-*) legend(MAM,JJA,SON,DJF,Orientation,horizontal)

14、legend(boxoff) %legend() title(ZHANDIANss,:) set(h1,position,0.1,0.65,0.6,0.2) set(gca,ytick,0:0.5:1.5,yticklabel,0.0;0.5;1.0;1.5,xticklabel,fontsize,10) ylabel(AOD at 550 nm) set(gca,yminortick,on); set(gca,TickLength,0.02;0.02); h2 = subplot(3,1,2) plot(spr(:,1),spr(:,4),g-) xlim(2000.5 2006.5) ho

15、ld on plot(sum(:,1),sum(:,4),r-o) hold on plot(aut(:,1),aut(:,4),y-) hold on plot(win(:,1),win(:,4),b-*) %set(gca,yticklabel,0.0 0.1 0.3 0.4,xticklabel,) set(h2,position,0.1,0.45,0.6,0.2,fontsize,10) set(gca,xticklabel,) set(gca,yminortick,on); set(gca,TickLength,0.02;0.02); set(gca,TickLength,0.02;

16、0.02,YAxisLocation,right); ax=axes(Position,get(gca,Position),. YAxisLocation,left,. Color,none,xlim,get(gca,xlim),. ylim,get(gca,ylim),yticklabel,); ylabel(PM10,);% 写成两行 set(gca,yminortick,on,fontsize,10); set(gca,TickLength,0.02;0.02); h3=subplot(3,1,3) plot(spr(:,1),spr(:,3),g-) xlim(2000.5 2006.

17、5) ylim(0 200) hold on plot(sum(:,1),sum(:,3),r-o) hold on plot(aut(:,1),aut(:,3),y-) hold on plot(win(:,1),win(:,3),b-*) set(h3,position,0.1,0.25,0.6,0.2) set(gca,xticklabel,2001:1:2006,fontsize,10) set(gca,TickLength,0.01;0.01); xlabel(year) ylabel(API) set(gca,yminortick,on,fontsize,10); set(gca,

18、TickLength,0.02;0.02); %saveas(ff,E:resultseason_figure,sitess,:,season_modis.emf); %saveas(ff,E:resultseason_figure,sitess,:,season_modis.eps);2、相关系数计算r s=corrcoef(x,y)3、气溶胶插值程序wt=log(0.44/0.47)/log(0.659/0.47) taer44=raddata(ii,20)*(raddata(ii,19)/raddata(ii,20)*wt4、%计算大气层顶向上和地面向下辐射通量的反射率for i=1:x

19、h_num % solardistance: 日地距离 xxx=2*3.1416*(xh_total(i,5)-1)/365 ; solardistance(i)=1.000109+0.033494*cos(xxx)+0.001472*sin(xxx)+0.000768*cos(2*xxx)+0.000079*sin(2*xxx); % 计算日地距离的比值albupflux(i) = xh_total(i,18)/(solardistance(i)*solardistance(i)*xh_total(i,25)*1366);% 计算大气层顶的反照率 albdownflux(i) = xh_to

20、tal(i,22)/(solardistance(i)*solardistance(i)*xh_total(i,25)*1366); % 计算透过 End6、插值程序load E:ceresskytoaoutputHuabeitotalflux1002.dat;datafile=totalflux1002toa_new=zeros(81,81);lat=zeros(81,81);lon=zeros(81,81);for n=1:81lat(1:81,n)=42:-0.05:38;lon(n,1:81)=115:0.05:119;endZI=zeros(81,81);for m=1:81 for

21、 n=1:81 %ZI(m,n) = interp2(bb,dd,ff,lon(m,n),lat(m,n); ZI(m,n) = griddata(datafile(:,7),datafile(:,8),datafile(:,11),lon(m,n),lat(m,n);%进行插值 endend7 计算平均值% 利用07年1月1113日的北京站aeronet数据来标定五台光度计的敏感性参数clear;clc;load E:sunphometersunprogramCIMEL_SUN_070113.mat%对于手持光度计，01号：55列为440nm%_给出变量的行数和列数line01,col01=

22、size(sun01_data);line02,col02=size(sun02_data);line10,col10=size(sun10_data);line11,col11=size(sun11_data);lineLOA,colLOA=size(sunLOA_data);lineBJ,colBJ=size(Beijing_data);%01#光度计标定day=11;12;13;sun_num=0; sun_num2=0;SUN01=0;for j=1:line01 %判断所测数据值间隔时间，间隔小于6分钟，求平均，认为是一次观测 if (j=1) day=sun01_data(j,4)

23、; time=sun01_data(j,5)+sun01_data(j,6)/60; end if (day=sun01_data(j,4) difference = (sun01_data(j,5)+sun01_data(j,6)/60) - time; if(difference 0.05 ) sun_num2=sun_num2+1; sun01_mean(sun_num2,:) = SUN01/sun_num; sun_num=1; SUN01= sun01_data(j,:); time=sun01_data(j,5)+sun01_data(j,6)/60; % 由于判断differe

24、nce0.1要多循环一次，所以必须减去一次循环 end elseif(day=sun01_data(j,4) sun_num2=sun_num2+1; sun01_mean(sun_num2,:) = SUN01/sun_num; sun_num=1; SUN01= sun01_data(j,:); time=sun01_data(j,5)+sun01_data(j,6)/60; day=sun01_data(j,4); endend sun_num2=sun_num2+1; sun01_mean(sun_num2,:) = SUN01/sun_num;%sun01_mean为最终输出变量，为

25、每个时间段记录的平均%_9、读文件ddate=01;02;03;04;05;06;07;08;09;10;11;12;13;14;15;16;17;18;19;20;21;22;23;24;25;26;27;28;29;30;31;fid = fopen(E:atmosphereoutputximeng200701shanghai_ACS_M_Z_54102_200701,ddate(dday,:),.out,r); if fid=-1 %判断要读取的文件是否存在 ddd=str2num(ddate(dday,:); data=fscanf(fid,%i %i %f %f %f %f %f %

26、f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f,28 inf); data=data; d_n,bb=size(data);10、hdf文件读取Hdf-EOS GRID DATAMISR: Name: AerosolParameterAverageDimensions: Name: YDimSize: 360Name: XDimSize: 720Name: AlgorithmTypeSize: 8Tile Dimensions: 45 90Upper Left Grid Corner: -180000000 90000

27、000Lower Right Grid Corner: 180000000 -90000000Rows: 360.0Columns: 720.0Projection: geoZone Code: -1.0Sphere: Clarke 1866Projection Parameters:Origin Code: ulPixel Registration Code: center直接读取：Optical_depth_average = hdfread(E:MISRMISR_AM1.hdf,AerosolParameterAverage,Fields,Optical depth average)%

28、按照INDEX读取：Optical_depth_average = hdfread(E:MISRMISR_AM1.hdf,AerosolParameterAverage,Fields,Optical depth average,Index, 1.0 1.0 1.0 ,1.0 1.0 1.0 ,360.0 720.0 8.0 );% 按照地理坐标读取：Optical_depth_average = hdfread(E:MISRMISR_AM1.hdf,AerosolParameterAverage,Fields,Optical depth average,Box,80.0,100.0,70.0,

29、50.0);按Interpolate读取Optical_depth_average = hdfread(E:MISRMISR_AM1.hdf,AerosolParameterAverage,Fields,Optical depth average,Interpolate,80.0,100.0,70.0,50.0);11、hdf swath数据读取mod28L2File = MOD28L2.A2001185.0830.003.2001308112641.hdf;sstInfo = hdfinfo(mod28L2File,eos);sstSwath = sstInfo.Swath;idx = 1

30、1,4 4,; % Subset by 4 in each dimensionsst = hdfread(sstSwath,Fields,sst,Index,idx);% Read flags and construct land monFlags = hdfread(sstSwath,Fields,common_flags,Index,idx);land = (bitget(commonFlags,8) = 1);% Read SST attributes: scale and units.sstInfoHDF = hdfinfo(mod28L2File);datasets = sstInf

31、oHDF.Vgroup.Vgroup(2).SDS;sstDataset = datasets(strmatch(sst,datasets.Name,exact); %exact:准确匹配sstAttr = sstDataset.Attributes;sstSlope = double(sstAttr(strmatch(Slope, sstAttr.Name).Value)sstIntercept = double(sstAttr(strmatch(Intercept,sstAttr.Name).Value)sstUnits = sstAttr(strmatch(Units,sstAttr.N

32、ame).Value% Read MODIS geolocation data fields.mod03File = MOD03.A2001185.0830.003.2001305034045.hdf;geoInfo = hdfinfo(mod03File,eos);sstLat = double(hdfread(geoInfo.Swath,Fields,Latitude, Index,idx);sstLon = double(hdfread(geoInfo.Swath,Fields,Longitude,Index,idx);sstLat(sstLat = -999) = NaN;sstLon

33、(sstLon = -999) = NaN;12、绘制colormap % Load the colormap and display a standard colorbar.load colormapscmap = sstCmap; ndviCmap; % We saved our colormap in two equal partsfigure % first figurecolormap(cmap)set(axes,Visible,off)h = colorbar(horiz);set(get(h,Title),String,Standard Colorbar)13、读文件，去除头文件

34、 fid = fopen(E:Experimental Data试验数据Lidar Data2006-082006-08-03-Runresult,data1(i,1).name,r); if fid=-1 temp=fgetl(fid); data=fscanf(fid,%f %f %f %f %f %f %f %f,8 1000);%将数据倒过来 end end newname2=num2str(k),h;xlswrite(E:Experimental Data试验数据Lidar Data2006-082006-08-03-Runm,newname2,NuuH);14、多项式拟合及图形绘制

35、 X= 0:0.1:1;y=3 4 5 7.8 9 10 11 13.3 15; Dy = 0.15; for i=1:6 a,S = POLYFIT(x,y,n); An = a; da = dy *sqrt(diag(inv(S.R*S.R) %计算个系数的误差 DAn=da; freedom(n)=S.df % 记录自由度 ye, delta = polyval(a,x,S); % 计算拟和多项式的范围 YEn = ye; Dn = delta ; Chi2(n) = sum(y-ye).2/dy/dy % 计算不同阶次的x2量endQ = 1-chicdf(chi,freedom);

36、% 用于判断拟合良好度Subplot(1,2,1),plot(1:6,abs(chi2-freedom),b)xlabel(阶次),title(chi2与自由度)subplot(1,2,2),plot(1:6,Q,r,1:6,ones(1:6)*0.5)xlabel(阶次),title(Q与0.5线)clf，plot(x,y,b+)errorbar(x,YE3,D3,r)text(0.1,5.5,;text(0.1,5,f)15、如何读文本文件fidin=fopen(test2.txt); % 打开test2.txt文件 fidout=fopen(mkmatlab.txt,w); % 创建MK

37、MATLAB.txt文件 while feof(fidin) % 判断是否为文件末尾 tline=fgetl(fidin); % 从文件读行 if double(tline(1)=48&double(tline(1)=57 ）% 判断首字符是否是数值 fprintf(fidout,%snn,tline); % 如果是数字行，把此行数据写入文件MKMATLAB.txt continue % 如果是非数字继续下一次循环 end end fclose(fidout); MK=importdata(MKMATLAB.txt); %16、一些常用字符串命令用fread读入数据，再利用下面的一些函数处理二

38、进制数及其他字符，再把二进制转化为十进制数。deblank :去掉字符串末尾的空格findstr :查找字符串isstr :字符串判断isletter:字母判断isspace :空子符判断lower :转换为小写strcat :字符串连接组合strcmp :字符串比较strcmpi :字符串比较(忽略大小写)strjust :调整字符串排列位置strmatch:寻找符合条件的行strncmp :比较字符串的前n个字符strrep :寻找和替代strtok :寻找字符串中第一个空字符前边的字符串 /标记查找strvcat :字符串竖向连接(组合)symvar :寻找字符串表达式中的变量(除了i,

39、j,pi,inf,nan,eps以及通用函数 的其他标识符),如果没有找到变量则返回空单元数组strings :Matlab字符串函数说明str2mat :字符串转换成文本texlabel:将字符串转换为Tex格式(推荐理科同学掌握Latex)upper :字符串大写17、MATLAB中的colormap colormap(hot) shading flat % 不同填充格式 shading faceted shading interphsv - hue-saturation-value color maphot - black-red-yellow-white color mapgray - linear gray-scale color mapbone - gray-scale with tinge of blue color mapcopper - linear copper-tone color mappink - pastel shades of pink color mapwhite - all white color mapflag - alternating red, white, blue, and black color maplines - color map with the line col