How to do this more efficiently
显示 更早的评论
Hello,
I have this piece of code in MATLAB:
for kk=0:N-1
for mm=0:N-1
for pp=1:Np
for qq=1:Np
if ((kk*Ts+tau(pp)))<=(mm*Ts+tau(qq))) && ((kk+1)*Ts+tau(pp))>(mm*Ts+tau(qq)))
thetaSR=(((kk+1)*Ts+tau(pp)))-((mm*Ts+tau(qq))));
F_SR_MR(kk+1,mm+1)=F_SR_MR(kk+1,mm+1)+conj(H(pp))*H(qq)*(thetaSR*exp(1i*pi*fc**thetaSR)*sinc(fc*thetaSR));
end
which obviously is not very efficient. How can I re-write it more efficiently?
Thanks
6 个评论
the cyclist
2014-5-24
It's easier to help you if you include a self-contained piece of code (including values of N, tau, etc) so that we can run it without guessing what these values might be.
Jan
2014-5-24
The code does not run. Some end 's are missing an is not a valid Matlab operator.
S. David
2014-5-24
the cyclist
2014-5-25
Do you have a reference (e.g. an image you can post) that shows the mathematical formula you are trying to replicate? It might be easier to build the code directly from that rather than optimizing yours.
S. David
2014-5-25
采纳的回答
Roger Stafford
2014-5-26
It isn't necessary to do the inequality test for every possible pair of kk and mm values. Since kk and mm are integers and Ts must surely be a positive number, your pair of inequalities is logically equivalent to
kk-mm == d
where d = floor((tau(qq)-tau(pp))/Ts). Therefore you can simply add the appropriate vectors along corresponding diagonals of F. For large N, doing it this way should save quite a bit of computation time.
F=zeros(N);
for pp=1:Np
for qq=1:Np
d = floor((tau(qq)-tau(pp))/Ts);
kk = max(d,0):min(N-1,N-1+d);
mm = kk-d;
theta=(d+1)*Ts-(tau(qq)-tau(pp));
ix = d+1+(N+1)*mm;
F(ix) = F(ix)+conj(H(pp))*H(qq)*theta.*exp(1i*pi*fc*theta).*sinc(fc*theta);
end
end
12 个评论
Roger Stafford
2014-5-26
It just occurred to me that since F is constant along its diagonals, the code can be revised to simply compute the first column and first row and then use the 'toeplitz' function to fill out the remainder of matrix F. That could be a bit faster.
S. David
2014-5-26
Roger Stafford
2014-5-26
I don't see any attachment here.
the cyclist
2014-5-26
Eq.png is attached in the comment area of the question.
S. David
2014-5-27
Roger Stafford
2014-5-27
S. David, the code you are showing today (May 27) is not the same as that you gave three days ago (May 24) where you said "This one is working", and it will give different results. This apparently is the "second condition" you referred to. In particular, the diagonals will no longer be a constant value because of the introduction of 'theta2'. The code I gave you was based on that earlier version and of course will not give the correct results.
However, you can still make the new version far more efficient than the four-nested-for-loops method you have shown. Your two conditions
kk*Ts+tau(pp)<=mm*Ts+tau(qq) & (kk+1)*Ts+tau(pp)>mm*Ts+tau(qq)
and
mm*Ts+tau(qq)<=kk*Ts+tau(pp) & (mm+1)*Ts+tau(qq)>kk*Ts+tau(pp)
are logically equivalent to the two respective conditions
kk-mm == floor((tau(pp)-tau(qq))/Ts)
and
kk-mm == ceil((tau(pp)-tau(qq))/Ts)
This means that if (tau(pp)-tau(qq))/Ts is not an integer, you will be altering two, rather than one, diagonals within the N-by-N matrix F. It would therefore be far more efficient to simply write the code to do that directly, rather than blindly searching over all N^2 possible pairings of kk and mm values. It amounts to a rather minor modification of the same code I gave you.
One curious consequence of your new version is that if (tau(pp)-tau(qq))/Ts is an exact integer, then only one diagonal of F is to be modified. However, due to computer imprecision from round-off error that will rarely be the case even though in theory it might be true. Nevertheless, I would have thought that in such a case the "correct" thing to do would be to apply the same modification twice to the corresponding single diagonal. In other words, if both of your two conditions are simultaneously true, wouldn't it make more sense to apply both additions to the common diagonal rather than only the first one, as you have done with the if-elseif construct?
S. David
2014-5-28
Roger Stafford
2014-5-29
As to the equivalence of your two conditions with those I gave yesterday, first, the following is easily seen to be true: if x is any real number and n any integer, then 1) n<=x<n+1 is equivalent to n == floor(x) and 2) n-1< x<=n is equivalent to n == ceil(x). For your first condition, since Ts must be positive, the conditions
kk*Ts+tau(pp) <= mm*Ts+tau(qq) & (kk+1)*Ts+tau(pp) > mm*Ts+tau(qq)
can readily be seen by the appropriate manipulation to be equivalent to
kk-mm <= (tau(qq)-tau(pp))/Ts < kk-mm+1
and in turn, given that kk and mm are integers, by the above statement this is equivalent to
kk-mm == floor((tau(qq)-tau(pp))/Ts)
The second condition
mm*Ts+tau(qq) <= kk*Ts+tau(pp) & (mm+1)*Ts+tau(qq) > kk*Ts+tau(pp)
is similarly equivalent to
kk-mm-1 < (tau(qq)-tau(pp))/Ts <= kk-mm
which in turn is equivalent to
kk-mm == ceil((tau(qq)-tau(pp))/Ts)
This means that each of your two conditions will always be true on one and only one of the diagonals of F. If (tau(qq)-tau(pp))/Ts is not an integer then two adjacent diagonals will be altered since the 'ceil' value will be one greater than the 'floor' value. If it is an integer, only one diagonal is altered since the 'ceil' and 'floor' value would be the same.
Thus your four nested for-loops can be replaced by just two for-loops:
F=zeros(N);
for pp=1:Np
for qq=1:Np
d1 = floor((tau(qq)-tau(pp))/Ts);
kk = max(d1,0):min(N-1,N-1+d1);
mm = kk-d1;
theta1 = (kk+1)*Ts+tau(pp)-mm*Ts+tau(qq);
theta2 = (kk+1)*Ts+tau(pp)+mm*Ts+tau(qq);
F(kk+1+N*mm) = F(kk+1+N*mm) + conj(H(pp))*H(qq)* ...
theta1.*exp(1i*pi*fc*theta2).*sinc(fc*theta1);
d2 = ceil((tau(qq)-tau(pp))/Ts);
if d1~=d2
kk = max(d2,0):min(N-1,N-1+d2);
mm = kk-d2;
theta1 = (mm+1)*Ts+tau(qq)-kk*Ts+tau(pp);
theta2 = (mm+1)*Ts+tau(qq)+kk*Ts+tau(pp);
F(kk+1+N*mm) = F(kk+1+N*mm) + conj(H(pp))*H(qq)* ...
theta1.*exp(1i*pi*fc*theta2).*sinc(fc*theta1);
end
end
end
(Note: The kk+1+N*mm expression in F(kk+1+N*mm) acts as a linear index which accesses the corresponding diagonal in F defined by the kk and mm vectors as kk-mm is held constant.)
It should be noted that small round-off errors by the computer can very easily change the quantity (tau(qq)-tau(pp))/Ts from being an integer to a non-integer or visa versa, which in turn would affect the result in the F matrix. In one case only one diagonal would receive an addition and in the other case two diagonals would be involved.
Roger Stafford
2014-5-29
@S. David, I have a few comments about the computation you have presented here in this thread.
You have stated in an earlier answer to Cyclist's question that your computation was in accordance with the formula shown in
http://www.mathworks.com/matlabcentral/answers/uploaded_files/13251/Eq.png
That formula would imply that PHI(k,m), which I assume is your F(k,m), is constant along diagonal lines where k-m is held constant. This is because an equal change in k and m will induce an equal shift in the two 'g' quantities, and therefore the integral of their product will remain unchanged. However, with 'theta2' as you have defined it, F is definitely not constant along its diagonals. How do you explain that?
A second comment is that it looks as though there is an error in the Eq.png formula itself. The two summations shown there both contain tau(p), and this would cause the two exponential expressions to cancel one another. I suspect the second summation, which is taken with respect to q, should have had tau(q) in its exponential expression. Do you agree with that?
Finally, could you please explain how the product
theta1*exp(1i*pi*fc*theta2)*sinc(fc*theta1)
relates to the formula in Eq.png. I don't follow that at all and it doesn't look right to me.
S. David
2014-5-29
S. David
2014-6-7
更多回答（1 个）
the cyclist
2014-5-24
Especially if N is large, you might get a huge speedup if you preallocate the memory for F_SR_MR. Put the line
F_SR_MR = zeros(N,N);
ahead of the loops.
类别
在 帮助中心 和 File Exchange 中查找有关 Type Checking 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
