推荐星级:
  • 1
  • 2
  • 3
  • 4
  • 5

Writing Fast MATLAB Code:编写快速的MATLAB代码

更新时间:2019-12-06 10:58:54 大小:335K 上传用户:梦留香查看TA发布的资源 浏览次数:99 下载积分:2分 出售积分赚钱 评价赚积分 ( 如何评价?) 标签:matlab 收藏 评论(1) 举报

资料介绍

Writing Fast MATLAB Code:编写快速的MATLAB代码

部分文件列表

文件名 大小
Writing_Fast_MATLAB_Code:编写快速的MATLAB代码.pdf 335K

部分页面预览

(完整内容请下载后查看)
本人精通MATLAB等编程语言,可以提供以下方向的帮助  
1. MATLAB/GUI/SIMULINK/C++/VC++编程问题;  
2. 线性与非线性控制、智能控制、模糊控制;  
3. 数值计算问题、小波分析算法、有限元问题;  
4. 电机控制、电力系统、机器人路径优化、机器人控制;  
5. 粒子群算法、神经网络、模拟退火算法等智能优化算法;  
6. 图像处理、信号处理、语音信号处理、电子通信等方向;  
有问题的朋友,可以将问题直接发到我的邮箱,24小时内给您答复!非  
常欢迎大家加我为QQ好友,欢迎访问我的空间!  
联系方式:  
QQ:626815632  
邮箱:  
QQ空间:http://626815632.qzone.qq.com/  
声明:本资料来源于网络,切勿用做商业用途!请您支持正版图书!  
Writing Fast MATLAB Code  
Pascal Getreuer, June 2006  
Contents  
1
2
3
The Profiler  
2
3
Array Preallocation  
Vectorization  
5
5
6
3.1 Vectorized Computations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  
3.2 Vectorized Logic . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  
4
5
Inlining Simple Functions  
10  
Numerical Integration  
12  
5.1 One-Dimensional Integration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13  
5.2 Multidimensional Integration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14  
6
7
8
Signal Processing  
16  
19  
22  
Referencing Operations  
Miscellaneous Tricks  
8.1 Clip a value without using if statements . . . . . . . . . . . . . . . . . . . . . . . . . . . 22  
8.2 Convert any array into a column vector . . . . . . . . . . . . . . . . . . . . . . . . . . . 22  
8.3 Find the min/max of a matrix or N-d array . . . . . . . . . . . . . . . . . . . . . . . . . 22  
8.4 Vector Normalization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23  
8.5 Flood filling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23  
8.6 Vectorized use of set on GUI objects . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24  
9
Further Reading  
24  
Introduction  
The Matlab programming language is parsed; code is interpreted during runtime. Languages like  
C and Fortran are faster because they are first compiled into the computer’s native language. The  
advantages of parsing in realtime are greater platform independence, greater language flexibility, and  
easier debugging. The disadvantages are slower speed, increased overhead, and limited low-level control.  
This article discusses strategies for improving the speed of Matlab code. Keep in mind that Matlab  
has gone through many versions and that it is available on many platforms. The fastest method on one  
system may not be the fastest on another. This article provides methods that are generally fast, but  
makes no claim on what is fastest.  
Caution!  
Learn the language first. Optimization requires comfort with the syntax and functions of the  
language. This article is not a tutorial on MATLAB.  
Use comments. Optimized code tends to be terse and cryptic. Help others and yourself by  
remembering to comment.  
Don’t optimize code before its time. Is optimization worth the effort? If the code will soon  
be revised or extended, it will be rewritten anyway.  
Only optimize where necessary. Focus your efforts on bottlenecks.  
1
1 The Profiler  
Matlab 5.0 and newer versions include a tool called the “profiler” that helps determine where the  
bottlenecks are in a program. Consider the following function:  
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  
function result = example1(Count)  
for k = 1:Count  
result(k) = sin(k/50);  
if result(k) < 0.9  
result(k) = gammaln(k);  
end  
end  
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  
To analyze the efficiency this function, first enable the profiler and clear any old profiler data:  
>> profile on  
>> profile clear  
Now run the program. Change the input argument higher or lower so that it takes about a second.  
>> example1(5000);  
Then enter  
>> profile report  
The profiler generates an HTML report on the function and launches a browser window. Depending on  
the system, profiler results may be a little different from this example.  
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  
MATLAB Profile Report: Summary  
Report generated 30-Jul-2004 16:57:01  
Total recorded time:  
Number of M-functions:  
Clock precision:  
3.09 s  
4
0.016 s  
Function List  
Name  
Time  
3.09 100.0%  
Time Time/call Self time  
Location  
example1  
gammaln  
profile  
1
3562  
1
3.094000 2.36 76.3% ~/example1.m  
0.73  
0.00  
23.7%  
0.0%  
0.0%  
0.000206 0.73 23.7% ../toolbox/matlab/specfun/gammaln.m  
0.000000 0.00  
0.000000 0.00  
0.0% ../toolbox/matlab/general/profile.m  
0.0% ../toolbox/matlab/general/profreport.m  
profreport 0.00  
1
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  
2
Clicking the “example1” link gives more details:  
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  
100% of the total time in this function was spent on the following lines:  
3: for k = 1:Count  
2.11 68% 4:  
5:  
result(k) = sin(k/50);  
0.14  
5% 6:  
if result(k) < -0.9  
result(k) = gammaln(k);  
end  
0.84 27% 7:  
8:  
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  
The most time-consuming lines are displayed, along with time, time percentage, and line number. The  
most costly lines are the computations on lines 4 and 7.  
2 Array Preallocation  
Matlab’s matrix variables have the ability to dynamically augment rows and columns. For example,  
>> a = 2  
a =  
2
>> a(2,6) = 1  
a =  
2
0
0
0
0
0
0
0
0
0
0
1
Matlab automatically resizes the matrix. Internally, the matrix data memory must be reallocated with  
larger size. If a matrix is resized repeatedly—like within a for loop—this overhead can be significant.  
To avoid frequent reallocations, “preallocate” the matrix with the zeros command. Consider the code:  
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  
a(1) = 1;  
b(1) = 0;  
for k = 2:8000  
a(k) = 0.99803  
b(k) = 0.06279  
a(k 1) 0.06279  
a(k 1) + 0.99803  
b(k 1);  
b(k 1);  
*
*
*
*
end  
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  
The profiler timed this code to take 0.47 seconds. After the for loop, both arrays are row vectors of  
length 8000, thus to preallocate, create empty a and b row vectors each with 8000 elements.  
3
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  
a = zeros(1,8000);  
b = zeros(1,8000);  
a(1) = 1;  
% Preallocation  
b(1) = 0;  
for k = 2:8000  
a(k) = 0.99803  
b(k) = 0.06279  
end  
a(k 1) 0.06279  
a(k 1) + 0.99803  
b(k 1);  
b(k 1);  
*
*
*
*
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  
With this modification, the code takes only 0.14 seconds (over three times faster). Preallocation is often  
easy to do, in this case it was only necessary to determine the right preallocation size and add two lines.  
What if the final array size can vary? Then use the upper bound on the array size and cut the excess  
after the loop:  
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  
a = zeros(1,10000);  
count = 0;  
% Preallocate  
for k = 1:10000  
v = exp(rand(1) rand(1));  
*
if v > 0.5  
count = count + 1;  
a(count) = v;  
% Conditionally add to array  
end  
end  
a = a(1:count);  
% Trim the result  
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  
The average run time of this program is 0.42 seconds without preallocation and 0.18 seconds with it.  
Preallocation can also be done with cell arrays, using the cell command to create the desired size.  
Using preallocation with a frequently resizing cell array is even more beneficial than with double arrays.  
4
3 Vectorization  
To “vectorize” a computation means to replace parallel operations with vector operations. This strategy  
often improves speed ten-fold. Good vectorization is a skill that must be developed; it requires comfort  
with Matlab’s language and creativity.  
3.1 Vectorized Computations  
Many standard Matlab functions are “vectorized,” they can operate on an array as if the function had  
been applied individually to every element.  
>> sqrt([1,4;9,16])  
>> abs([0,1,2,5,6,7])  
ans =  
ans =  
1
3
2
4
0
1
2
5
6
7
Consider the following function:  
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  
function d = minDistance(x,y,z)  
% Find the min distance between a set of points and the origin  
nPoints = length(x);  
d = zeros(nPoints,1);  
% Preallocate  
for k = 1:nPoints  
% Compute distance for every point  
d(k) = sqrt(x(k)ˆ2 + y(k)ˆ2 + z(k)ˆ2);  
end  
d = min(d);  
% Get the minimum distance  
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  
For every point, the distance between it and the origin is computed and stored in d. The minimum  
distance is then found with min. To vectorize the distance computation, replace the for loop with  
vector operations:  
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  
function d = minDistance(x,y,z)  
% Find the min distance between a set of points and the origin  
d = sqrt(x.ˆ2 + y.ˆ2 + z.ˆ2);  
d = min(d);  
% Compute distance for every point  
% Get the minimum distance  
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  
The modified code performs the distance computation with vector operations. The x, y and z arrays  
are first squared using the per-element power operator, .^ (the per-element operators for multiplication  
and division are .* and ./). The squared components are added with vector addition. Finally, the  
square root of the vector sum is computed per element, yielding an array of distances.  
5
The first version of the minDistance program takes 0.73 seconds on 50000 points. The vectorized  
version takes less than 0.04 seconds, more than 18 times faster.  
Some useful functions for vectorizing computations:  
min, max, repmat, meshgrid, sum, cumsum, diff, prod, cumprod, filter  
3.2 Vectorized Logic  
The previous section shows how to vectorize pure computation. Bottleneck code often involves condi-  
tional logic. Like computations, Matlab’s logic operators are vectorized:  
>> [1,5,3] < [2,2,4]  
ans =  
1
0
1
Two arrays are compared per-element. Logic operations return “logical” arrays with binary elements.  
How is this useful? Matlab has a few powerful functions for operating on logical arrays:  
find: Find indices of nonzero elements.  
any: True if any element of a vector is nonzero (or per-column for a matrix).  
all: True if all elements of a vector are nonzero (or per-column for a matrix).  
>> find([1,5,3] < [2,2,4])  
>> find(eye(3) == 1)  
ans =  
ans =  
1
3
1
5
9
The find function returns the indices where the vector logic operation returns true. In the first example,  
1 < 2 is true, 5 < 2 is false, and 3 < 4 is true, so find reports that the first and third comparisons are  
true. In the second example, find returns the indices where the identity matrix is equal to one. The  
indices 1, 5, and 9 correspond to the diagonal of a 3 by 3 matrix.  
The any and all functions are simple but occasionally very useful. For example, any(x(:) < 0)  
returns true if any element of x is negative.  
6
Example 1: Removing elements  
The situation often arises where array elements must be removed on some per-element condition. For  
example, this code removes all NaN and infinite elements from an array x:  
i = find(isnan(x)  
x(i) = [];  
|
isinf(x));  
% Find bad elements in x  
% and delete them  
Alternatively,  
i = find(˜isnan(x) & ˜isinf(x)); % Find elements that are not NaN and not infinite  
x = x(i); % Keep those elements  
Both of these solutions can be further streamlined by using logical indexing:  
x(isnan(x)  
|
isinf(x)) = [];  
% Delete bad elements  
or  
x = x(˜isnan(x) & ˜isinf(x));  
% Keep good elements  
Example 2: Piecewise functions  
The sinc function has a piecewise definition,  
(
sin(x)/x, x = 0  
1, x = 0  
sinc(x) =  
This code uses find with vectorized computation to handle the two cases separately:  
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  
function y = sinc(x)  
% Computes the sinc function perelement for a set of x values.  
y = ones(size(x));  
% Set y to all ones, sinc(0) = 1  
% Find nonzero x values  
i = find(x ˜= 0);  
y(i) = sin(x(i)) ./ x(i);  
% Compute sinc where x ˜= 0  
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  
An interesting alternative is y = (sin(x) + (x == 0))./(x + (x == 0)).  
Example 3: Drawing images with meshgrid  
The meshgrid function takes two input vectors and converts them to matrices by replicating the first  
over the rows and the second over the columns.  
>> [x,y] = meshgrid(1:5,1:3)  
x =  
1
1
1
2
2
2
3
3
3
4
4
4
5
5
5
7
y =  
1
2
3
1
2
3
1
2
3
1
2
3
1
2
3
The matrices above work like a map for a width 5, height 3 image. For each pixel, the x-location can  
be read from x and the y-location from y. This may seem like a gratuitous use of memory as x and y  
simply record the column and row positions, but this is useful. For example, to draw an ellipse,  
% Create x and y for a width 150, height 100 image  
[x,y] = meshgrid(1:150,1:100);  
% Ellipse with origin (60,50) of size 15 x 40  
Img = sqrt(((x60).ˆ2 / 15ˆ2) + ((y50).ˆ2 / 40ˆ2)) > 1;  
% Plot the image  
imagesc(Img); colormap(copper);  
axis image; axis off;  
Drawing lines is almost the same, just a change in the formula.  
[x,y] = meshgrid(1:150,1:100);  
% The line y = x 0.8 + 20  
*
Img = (abs((x 0.8 + 20) y) > 1);  
*
imagesc(Img); colormap(copper);  
axis image; axis off;  
Polar functions can be drawn by first converting x and y variables with the cart2pol function.  
[x,y] = meshgrid(1:150,1:100);  
[th,r] = cart2pol(x 75,y 50); % Convert to polar  
% Spiral centered at (75,50)  
Img = sin(r/3 + th);  
imagesc(Img); colormap(hot);  
axis image; axis off;  
8

推荐下载

全部评论(1)

  • 2019-12-09 10:28:40suxindg

    谢谢分享