How to vectorize this for loop
1 visualización (últimos 30 días)
Mostrar comentarios más antiguos
for i = 1:no_frame
frame(:,i) = y(start:stop).*wind_sam;
start = start+h_size;
stop = start+w_sam-1;
end
3 comentarios
Jan
el 23 de En. de 2018
Editada: Jan
el 23 de En. de 2018
What is an "array function"? Is wind_sam a vector, a matrix or is it really a "function"? Please provide some example data, e.g. produced by rand.
Did you pre-allocate the output frame? If not, this might be more efficient than a vectorization, if creating the required temporary array to store all y(start:stop) is time consuming.
Jan
el 23 de En. de 2018
@serena: Please post a working example, which e.g. explains the initial value of start and whether the subvectors of y overlap or not. Can hsize be smaller than w_sam - 1? Seeing a part of the code only impedes the suggestion of improvements.
Respuestas (3)
Walter Roberson
el 23 de En. de 2018
You have an initial value for start. Remove the first start-1 values from y. Reshape into rows of h_size. Take the first w_sam rows of the result.
1 comentario
Jan
el 23 de En. de 2018
hsize could be smaller than w_sam-1. Then the overlapping subvectors of y cannot be represented by reshaping the vector.
David Goodmanson
el 23 de En. de 2018
Editada: David Goodmanson
el 23 de En. de 2018
Hi serena,
I believe this works, assuming that wind_sam is a column vector of length w_sam. The idea is to make an index matrix and use it to address the contents of y.
I did not use variable names start and stop because they are existing Matlab function names.
st0 = 5; % whatever the first start index is
st = st0 + h_size*(0:no_frame-1); % vector of starting indices
ind = st + (0:w_sam-1)'; % matrix of start-to-stop indices, in columns
frame1 = y(ind).*wind_sam;
If you have an older Matlab version and line for ind does not work, you could use
ind = repmat(st,w_sam,1) + repmat((0:w_sam-1)',1,no_frame);
11 comentarios
Walter Roberson
el 28 de En. de 2018
Jan, your description of the benefits of y(a:b) over y([a:b]) do not agree with the example documentation at https://www.mathworks.com/help/matlab/ref/subsref.html#br5htww-6 which show that subsref is called with expanded indices.
Jan
el 28 de En. de 2018
Editada: Jan
el 28 de En. de 2018
@Walter: You are right, there is a discrepancy. Then I dare to consider the explanations at the link as simplification, which does not take the JIT acceleration into account. Usually I trust the documentation for good reasons, but here I am suspicious because of the timings:
function Untitled
x = rand(1, 1e6);
y = zeros(1000, 1000);
% Warm up - not used for measuring:
for k = 1:1000
a = (k - 1) * 1000 + 1;
y(:, k) = x([a:a+999]);
end
tic
for rep = 1:100
for k = 1:1000
a = (k - 1) * 1000 + 1;
y(:, k) = x([a:a+999]);
end
end
toc
tic
for rep = 1:100
for k = 1:1000
a = (k - 1) * 1000 + 1;
v = a:a+999;
y(:, k) = x(v);
end
end
toc
tic
for rep = 1:100
for k = 1:1000
a = (k - 1) * 1000 + 1;
y(:, k) = x(a:a+999);
end
end
toc
Elapsed time is 1.640906 seconds. % x([a:b])
Elapsed time is 1.636042 seconds. % v=a:b; x(v)
Elapsed time is 0.737971 seconds. % x(a:b)
I can imagine 2 explanations for the speed up of x(a:b):
- a:b is not created explicitly and range checks are omitted for the inner elements.
- a faster memcpy method instead of an elementwise copy (in addition, therefore 1. must be assumed also).
The 2. idea can be excluded by using a:2:b for indexing, which shows a comparable advantage for x(a:2:b) compared to x([a:2:b]).
Therefore I claim boldly, that 1. is applied.
What a pity that the JIT is not documented. But, wait, even blaming the JIT is a too cheeky speculation:
feature jit off
feature accel off
and enabling the profiler does not influence the timings remarkably. Anyway, something smarter than calling subref must happen here. How would you explain the speed difference?
Maybe we should ask TMW for an explanation. Their argument for not documenting the JIT was: "If we publish the methods, the users will optimize their code for the JIT. But we want to optimize the JIT for the real user code." But this is not mutual exclusive.
Walter Roberson
el 26 de En. de 2018
If you have the Communications toolbox, I suggest calling buffer()
0 comentarios
Ver también
Categorías
Más información sobre Logical en Help Center y File Exchange.
Productos
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!