提高matlab运算速度的几种方法.docx

由于matlab是一种解释性语言，所以在matlab程序中最忌讳直接使用循环语句，如果不得已要使用for循环，可以采用以下方法提高速度。1、使用6.5以上版本，对循环已作优化；2、尽可能转化为矩阵运算；3、转化为二进制执行文件运算，如使用matlab内带的编译系统或matcom以及com组件技术。其中com组件技术最方便的就是利用com builder来实现，这里重点介绍。com builder是matlab6.5才有的，也是mathworks公司推荐使用于混合编程的，这些日子进行了全方位的摸索，感觉是爽呆了，下面我们一起来揭开它的神秘面纱。此系列分为以下几块：1.matlab下做com组件2.vb,c#.net实现调用3.vc实现调用4.打包5.优缺点评注其中2,3部分可以选择一个看matlab下做com组件com是component object module的简称，它是一种通用的对象接口，任何语言只要按照这种接口标准，就可以实现调用它。matlab6.5新推出来的combuilder就是把matlab下的程序做成com组件，供其他语言调用。我们先准备两个测试文件，并copy一个图片到c盘下，起名叫1.jpg(这些你都可以改的，我这儿是为了程序方便)：第一个叫im_test.m如下:function im_test %这个文件不带输入与输出I=imread(c:1.jpg); %因为以前带有imshow的程序用mcc编成dll后用不%了，所以想试combuilder是否imshow(I); %能支持这些函数第二个叫split2rgb.m，就是前些日子Zosco发给我的那个程序，因为它用mcc编成dll后有问题，所以我在这儿继续将它进行测试，而且它也带有多个输入及输出参数，所以也正好拿来测试在matlab的workspace下打comtool，就打开了matlab com builder，点击file-new project，新建一个工程，在component name里填上comtest，Class name里填上一个sgltest(并将自动生成classes里的comtest remove掉),complie code in选c或c+都无所谓，将Complier options里的Use Handle Graphics library的复选框画上，点击ok就行了。然后点击project-Add files,将im_test.m和split2rgb.m添加入工程，然后点Build-Com Object，就会在comtestdistrib文件夹下生成一个comtest_1_0.dll(它就是做好的com组件)，Build时matlab已经自动将此dll在注册表中注册,为了下面能用其他编译器调用，我们还需做一个准备工作，开一个dos窗口，进入/bin/win32目录下(matlabroot为你机器上matlab安装的路径),打regsvr32 mwcomutil.dll，即对mwcomutil.dll进行注册(这个dll是matlab下作的任何com组件都要用到的dll),下面我们在其他编译器下调用时就可以用了。是不是觉得做起com组件来很简单呢，matlab下可以给com组件中的类添加成员、事件、方法等，我这儿其实是给类sgltest添加了两个方法，怎么添加成员和方法可以参看matlab的com builder的帮助。附录：split2rgb.m的源代码%/ 测试文件function m_nHeight,m_nWidth,mOrigR,mOrigG,mOrigB=Split2RGB(FileName)%/ 读入一个Jpg文件，mOrigData=imread(FileName);%mDestData=imresize(mOrigData,0.5);imwrite(mOrigData,c:2.jpg);%/ 用三个变量保存其R,G,B分量mOrigR=mOrigData(:,:,1);mOrigG=mOrigData(:,:,2);mOrigB=mOrigData(:,:,3);%/ 获得图象的高度，宽度m_nHeight,m_nWidth=size(mOrigR);figure(1);set(gcf,MenuBar,none)imshow(mOrigData);title(Orginal Image:,FileName,Color,b,FontSize,14);xlabel(Height: ,num2str(m_nHeight), Width :,num2str(m_nWidth),Color,b,FontSize,12);%/ 写param文件paraFName=FileName(1:length(FileName)-4),.param;fid=fopen(paraFName,w);fwrite(fid,m_nHeight,uint32);fwrite(fid,m_nWidth,uint32);fclose(fid);%/ 写 R 分量文件RFName=FileName(1:length(FileName)-4),_R.rot;fid=fopen(RFName,w);fwrite(fid,mOrigR,uint8);fclose(fid);%/ 写 G 分量文件GFName=FileName(1:length(FileName)-4),_G.rot;fid=fopen(GFName,w);fwrite(fid,mOrigG,uint8);fclose(fid);%/ 写 B 分量文件BFName=FileName(1:length(FileName)-4),_B.rot;fid=fopen(BFName,w);fwrite(fid,mOrigB,uint8);fclose(fid);这一部分讲vb,c#.net下怎么实现调用上一部分我们生成的comtest_1_0.dll(matlab下做的com组件),记得一定先要对mwcomutil.dll进行注册(怎么注册参看上一部分)1.vb下实现调用打开或新建一个vb工程，点project-Reference,在弹出来的窗口中找到comtest 1.0 Type Library,将前面的复选框选上，点击ok，此时便将此com组件添加到工程里面去了，此时你可以用对象浏览器看到comtest下有个sgltest类，这个类里面有两个方法im_test,split2rgb,还有个MWFlags成员(这个成员是自动生成的),vb下测试代码如下:测试im_test方法的代码:Dim st As sgltestSet st = New sgltestCall st.im_test测试split2rgb方法的代码:Dim st As sgltestSet st = New sgltestDim h As Variant, w As Variant, r As Variant, g As Variant, b As Variant, filename As Variantfilename = c:1.jpgCall st.split2rgb(5, h, w, r, g, b, filename) 可见matlab下函数的输入输出参数在com组件里全是variant型的变量，测试大获成功，结果就跟matlab下运行的一摸一样，爽2.c#.net下实现调用打开或新建一个c#项目(我采用的是编辑器)，选中右边的解决方案资源管理器中的引用，点鼠标右键，选添加引用，在弹出来的窗口中选com，然后也找到comtest_1_0.dll，点选择，然后确定就可，此时此com组件也添加到工程里面去了，同样我们可以选择对象浏览器来看comtest及下面的sgltest类，c#测试项目如下:测试im_test方法的代码:comtest.sgltestClass st=new comtest.sgltestClass(); st.im_test();测试split2rgb方法的代码:comtest.sgltestClass st=new comtest.sgltestClass(); object h=null,w=null,r=null,g=null,b=null;object filename=c:1.jpg;st.split2rgb(5,ref h,ref w,ref r,ref g,ref b,filename);可见输入参数是ref object型的，而输出参数是object型的，测试同样大获成功，与matlab下运行的结果一摸一样，爽呆了。这一部分讲vc下实现调用第一部分我们生成的comtest_1_0.dll，同样要记得先对mwcomutil.dll进行注册(怎么注册参看第一部分)，vb和.net下实现对com组件的调用很简单，而vc下实现这一步我可是摸索了好几天(主要是vc用的不好)1.先做一些准备工作，用ole viewer工具开始-程序-Microsoft visual studio6.0-Microsoft visual studio6.0 Tools-OLE viewer(这个工具也可以通过在vc下点Tools-OLE/COM Object Viewer来打开，在Oleviewer工具里，在右边选择Type libraries,将他展开，找到comtest 1.0 Type Library，选中它，点鼠标右键，选view，便又弹出一窗口，点工具栏上的save按钮，分别将他保存为comtest_1_0.h,comtest_1_0.c(也可以存为comtest_1_0.Idl接口文件),我们就可以通过这两个文件实现对comtest_1_0.dll调用2.vc下调用新建或打开一个vc工程，注意，此时不用对编译器进行任何设置(而用mcc生成的dll我们么设置一大堆,我这儿只设置了Precomplied headers,选Automatic use of precompliedheaders,写上 stdafx.h),将comtest_1_0.h和comtest_1_0.c加入工程，并复制一个comtest_1_0.dll到工程目录下，由于comtest_1_0.dll还要用到mwcomutil.dll，所以将/extern/include/下的mwcomutil.h和mwcomtypes.h也加入工程(将这两个文件拷贝入工程路径下，如果设置了library path,可以不加)此时可以通过classView看到多出了_IID、IMWUtil,Isgltest类，Isgltest就是我们在matlab下建起来的sgltest类vc下代码如下/这几个是引入dll和头文件#import mwcomutil.dll#import comtest_1_0.DLL#include mwcomutil.h#include comtest_1_0.h#include comutil.h /此文件是用来处理由char *向VARIANT类型的转换测试im_test方法的代码:if(FAILED(CoInitialize(NULL) /初始化调用comAfxMessageBox(unable to initialize COM);Isgltest *st=NULL; /创建一个com组件，CLSID_sgltest和IID_Isgltest可以从comtest_1_0.h和comtest_1_0.c里找到HRESULT hr=CoCreateInstance(CLSID_sgltest,NULL,CLSCTX_ALL,IID_Isgltest,(void *)&st);if(SUCCEEDED(hr)st-im_test(); AfxMessageBox(succeed);st-Release();elseAfxMessageBox(unsucceed); 如果你的vc工程是console project的话，上述的AfxMessageBox可改成printf或cout，测试split2rgb方法的代码(类型转换我参照visual c的精华区也转换成功了)if(FAILED(CoInitialize(NULL)AfxMessageBox(unable to initialize COM);Isgltest *st=NULL; HRESULT hr=CoCreateInstance(CLSID_sgltest,NULL,CLSCTX_ALL,IID_Isgltest,(void *)&st);VARIANT m,n,r,g,b,filename;VariantInit(&m);VariantInit(&n);VariantInit(&r);VariantInit(&g);VariantInit(&b);VariantInit(&filename);filename.vt=VT_BSTR;filename.bstrVal=_com_util:ConvertStringToBSTR(C:1.jpg);if(SUCCEEDED(hr)st-split2rgb(5,&m,&n,&r,&g,&b,filename);st-Release();AfxMessageBox(succeed);elseAfxMessageBox(unsucceed); 同样，运行结果与matlab下的结果一摸一样，记得我们的im_test里面可是使用了imshow阿，以前用mcc生成的程序中用它可是有错哦，爽呆了。关于VC下用com组件及其类型的转变请参看msdn及其Visual C的精华区combuilder系列可以结尾了一.打包:在matlab下的workspace里打comtool，点file-open project将我们先前建好的comtest.cbl工程文件打开，再点component-package component就实现了打包，此时到comtestdistrib文件夹里看，生成的comtest.exe就是打包后的解压程序，双击它会解压出一些文件，再点击解压出来的_install.bat就可以实现安装(按他的步骤去做就行了)二.优缺点评注这几天用这个combuilder可把我给爽死了，特别是在vc下调用成功时，记得精华区里曾讲combuilder没有什么实质性的突破，我可不这么认为，它的突破可大了1.做出来的是com，通用的，任何编译器只要支持com,就可以实现调用，想c+ builder,Delphi等的，我想只要按照调用com组件去做，也能成功的2.支持imshow等一些原来混编用不了的函数，对图形库的支持也比以前强(这些还需各位大侠共同测试)3.实现方法简单，没有像以前混编还要设置一大堆东东4.能够在matlab下写自己的类，并为自己的类编写成员、方法和事件，管理工程也方便(这个有点像vc、vb下管理工程一样)用的太爽了，一下子还不知道怎么写缺点了，0，我想缺点还需大家一起用来找出我这儿说一个缺点，感觉它的参数全是VARIANT型的，不怎么方便。一、 基本技术-1)MATLAB索引或引用(MATLAB Indexing or Referencing)在MATLAB中有三种基本方法可以选取一个矩阵的子阵。它们分别是下标法，线性法和逻辑法(subscripted, linear, and logical)。如果你已经熟悉这个内容，请跳过本节1.1)下标法非常简单，看几个例子就好。A = 6:12;A(3,5)ans =8 10A(3:2:end)ans =8 10 12A = 11 14 17; .12 15 18; .13 16 19;A(2:3,2)ans =15161.2)线性法二维矩阵以列优先顺序可以线性展开，可以通过现行展开后的元素序号来访问元素。A = 11 14 17; .12 15 18; .13 16 19;A(6)ans =16A(3,1,8)ans =13 11 18A(3;1;8)ans =1311181.3)逻辑法用一个和原矩阵具有相同尺寸的0-1矩阵，可以索引元素。在某个位置上为1表示选取元素，否则不选。得到的结果是一个向量。A = 6:10;A(logical(0 0 1 0 1)ans =8 10A = 1 23 4;B = 1 0 0 1;A(logical(B)ans =1 4-2)数组操作和矩阵操作(Array Operations vs. Matrix Operations)对矩阵的元素一个一个孤立进行的操作称作数组操作；而把矩阵视为一个整体进行的运算则成为矩阵操作。MATLAB运算符*,/,都是矩阵运算，而相应的数组操作则是.*, ./, ., .A=1 0 ;0 1;B=0 1 ;1 0;A*B % 矩阵乘法ans =0 11 0A.*B % A和B对应项相乘ans =0 00 0-3)布朗数组操作(Boolean Array Operations)对矩阵的比较运算是数组操作，也就是说，是对每个元素孤立进行的。因此其结果就不是一个“真”或者“假”，而是一堆“真假”。这个结果就是布朗数组。D = -0.2 1.0 1.5 3.0 -1.0 4.2 3.14;D = 0ans =0 1 1 1 0 1 1如果想选出D中的正元素：D = D(D0)D =1.0000 1.5000 3.0000 4.2000 3.1400除此之外，MATLAB运算中会出现NaN,Inf,-Inf。对它们的比较参见下例Inf=Inf返回真Inf用于计算数组的双重for循环。使用的工具：数组乘法优化前：A = magic(100);B = pascal(100);for j = 1:100for k = 1:100;X(j,k) = sqrt(A(j,k) * (B(j,k) - 1);endend优化后：A = magic(100);B = pascal(100);X = sqrt(A).*(B-1);用一个循环建立一个向量，其元素依赖于前一个元素使用的工具：FILTER, CUMSUM, CUMPROD优化前：A = 1;L = 1000;for i = 1A(i+1) = 2*A(i)+1;end优化后：L = 1000;A = filter(1,1 -2,ones(1,L+1);如果你的向量构造只使用加法或乘法，你可使用cumsum或cumprod函数。优化前：n=10000;V_B=100*ones(1,n);V_B2=100*ones(1,n);ScaleFactor=rand(1,n-1);for i = 2:nV_B(i) = V_B(i-1)*(1+ScaleFactor(i-1);endfor i=2:nV_B2(i) = V_B2(i-1)+3;end优化后：n=10000;V_A=100*ones(1,n);V_A2 = 100*ones(1,n);ScaleFactor=rand(1,n-1);V_A=cumprod(100 1+ScaleFactor);V_A2=cumsum(100 3*ones(1,n-1);向量累加，每5个元素进行一次：工具：CUMSUM , 向量索引优化前：% Use an arbitrary vector, xx = 1:10000;y = ;for n = 5:5:length(x)y = y sum(x(1:n);end优化后(使用预分配)：x = 1:10000;ylength = (length(x) - mod(length(x),5)/5;% Avoid using ZEROS command during preallocationy(1:ylength) = 0;for n = 5:5:length(x)y(n/5) = sum(x(1:n);end优化后(使用矢量化,不再需要预分配)：x = 1:10000;cums = cumsum(x);y = cums(5:5:length(x);操作一个向量，当某个元素的后继元素为0时，重复这个元素：工具：FIND, CUMSUM, DIFF任务：我们要操作一个由非零数值和零组成的向量，要求把零替换成为它前面的非零数值。例如，我们要转换下面的向量：a=2; b=3; c=5; d=15; e=11;x = a 0 0 0 b 0 0 c 0 0 0 0 d 0 e 0 0 0 0 0;为：x = a a a a b b b c c c c c d d e e e e e e;解(diff和cumsum是反函数)：valind = find(x);x(valind(2:end) = diff(x(valind);x = cumsum(x);将向量的元素累加到特定位置上工具：SPARSE优化前：% The values we are summing at designated indicesvalues = 20 15 45 50 75 10 15 15 35 40 10;% The indices associated with the values are summed cumulativelyindices = 2 4 4 1 3 4 2 1 3 3 1;totals = zeros(max(indices),1);for n = 1:length(indices)totals(indices(n) = totals(indices(n) + values(n);end优化后：indices = 2 4 4 1 3 4 2 1 3 3 1;totals = full(sparse(indices,1,values);注意：这一方法开辟了稀疏矩阵的新用途。在使用sparse命令创建稀疏矩阵时，它是对分配到同一个索引的所有值求和，而不是替代已有的数值。这称为向量累加，是MATLAB处理稀疏矩阵的方式。=三、更多资源-10)矩阵索引和运算下面的MATLAB文摘讨论矩阵索引。它比本技术手册的第1节提供了更加详细的信息MATLAB Digest: Matrix Indexing in MATLAB/company/digest/sept01/matrix.shtml下面的说明链接将指导你如何使用MATLAB中的矩阵操作。含括矩阵创建、索引、操作、数组运算、矩阵运算以及其它主题。MATLAB Documentation: Getting Started/access/helpdesk/help/techdoc/learn_MATLAB/learn_MATLAB.shtmlPeter Acklam是我们的一个power user，他建立了一个网页提供MATLAB数组处理中的一些技巧。Peter Acklams MATLAB array manipulation tips and trickshttp:/home.online.no/pjacklam/MATLAB/doc/mtt/index.html-11)矩阵存储管理(Matrix Memory Management)有关预分配技术如何加快计算速度的更多的信息，参见下面的联机解决方案：26623: How do I pre-allocate memory when using MATLAB?/support/solutions/data/26623.shtml下面的技术手册是关于MATLAB存储管理方面的一个包罗广泛的指南：The Technical Support Guide to Memory Management/support/tech-notes/1100/1106.shtml-12)发挥MATLAB的最高性能(Maximizing MATLAB Performance)这个技术手册对代码矢量化进行一般的讨论，而许多时候，我们的目的是加快代码的速度。因此，本节提供一些附加的资源，它们涉及如何是代码达到最高性能。加速代码的常用技巧：3257: How do I increase the speed or performance of MATLAB?/support/solutions/data/3257.shtml编译你的代码并从MATLAB中调用它22512: What are the advantages of using the MATLAB Compiler tocompile my M-files?/support/solutions/data/22512.shtml=结论正如本文所介绍的，针对代码矢量化和代码提速，有一些常规方法，同时，还有性能