How to do this more efficiently

Question

S. David 2014-5-24

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/130980-how-to-do-this-more-efficiently

评论： S. David 2014-6-7

采纳的回答： Roger Stafford

在 MATLAB Online 中打开

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 个评论
显示 4更早的评论隐藏 4更早的评论

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

Eq.png

I attached the formula. g(t) in it is a rectangular pulse of magnitude one over the period [0,Ts).

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Roger Stafford 2014-5-26

2
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/130980-how-to-do-this-more-efficiently#answer_138260

在 MATLAB Online 中打开

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 个评论
显示 10更早的评论隐藏 10更早的评论

S. David 2014-5-27

在 MATLAB Online 中打开

Thanks "the cyclist" for mentioning where the attachment is.

Roger Stafford, the complete code is as follows (I included yours as well):

clear all;
clc
N=512; 
fc=12*10^3; 
B=8000; 
df=B/N;
T=1/df;
Ts=T/N;
tau=[0 5 7 20].*10^-3;
h=[1 0.8 0.7 0.5];
Np=length(tau);
H=h.*exp(-1i*2*pi*fc*tau);
F1=zeros(N);
tic
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)
                   theta1=(kk+1)*Ts+tau(pp)-mm*Ts+tau(qq); 
                   theta2=(kk+1)*Ts+tau(pp)+mm*Ts+tau(qq);
                   F1(kk+1,mm+1)=F1(kk+1,mm+1)+conj(H(pp))*H(qq)*(theta1*exp(1i*pi*fc*theta2)*sinc(fc*theta1));
                elseif mm*Ts+tau(qq)<=kk*Ts+tau(pp) && (mm+1)*Ts+tau(qq)>kk*Ts+tau(pp)
                   theta1=(mm+1)*Ts+tau(qq)-kk*Ts+tau(pp); 
                   theta2=(mm+1)*Ts+tau(qq)+kk*Ts+tau(pp);
                   F1(kk+1,mm+1)=F1(kk+1,mm+1)+conj(H(pp))*H(qq)*(theta1*exp(1i*pi*fc*theta2)*sinc(fc*theta1));
                end
            end
        end
    end
end
toc
tic
F2=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;
    F2(ix) = F2(ix)+conj(H(pp))*H(qq)*theta.*exp(1i*pi*fc*theta).*sinc(fc*theta);
  end
end
toc

It took 1.251529 s for my way while 0.004646 s for yours, which is 250 times faster!! But they are not equivalent, because the second condition is not met I guess. I appreciate if you explain to me why you did it this way, and how to include the second condition as well.

Thanks in advance

Roger Stafford 2014-5-27

在 MATLAB Online 中打开

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?

Roger Stafford 2014-5-29

在 MATLAB Online 中打开

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

在 MATLAB Online 中打开

@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-5-29

在 MATLAB Online 中打开

EqModified.png

I didn't see your last response when I replied the last time.

I am glad you asked about this, because I was going to mention it next.

Actually, your are right, according to the equation attached earlier, the matrix should be circular. However, my code is a bit more complicated, and I though writing it this way would make things easier to understand. So, yes, the thetas here are meaningless.

The (k,m)th element of F is given by the equation attached below. You can notice the presence of the exponential inside the integral changes the situation, and the matrix is no longer circular. I think this answers all your questions.

In this case the code would be:

clear all;
clc
N=512; 
fc=12*10^3; 
B=8000; 
df=B/N;
T=1/df;
Ts=T/N;
tau=[0 5 7 20].*10^-3;
h=[1 0.8 0.7 0.5];
Np=length(tau);
H=h.*exp(-1i*2*pi*fc*tau);
F=zeros(N);
a=[0.001 0.0012 0.0024 0.003];
for kk=0:N-1
    for mm=0:N-1
        for pp=1:Np
            for qq=1:Np
                if ((kk*Ts+tau(pp))/(1+a(pp)))<=(mm*Ts+tau(qq))/(1+a(qq)) && ((kk+1)*Ts+tau(pp))/(1+a(pp))>(mm*Ts+tau(qq))/(1+a(qq))
                   theta1=(((kk+1)*Ts+tau(pp))/(1+a(pp)))-((mm*Ts+tau(qq))/(1+a(qq))); 
                   theta2=(((kk+1)*Ts+tau(pp))/(1+a(pp)))+((mm*Ts+tau(qq))/(1+a(qq)));
                   F(kk+1,mm+1)=F(kk+1,mm+1)+conj(H(pp))*H(qq)*(theta1*exp(1i*pi*fc*(a(qq)-a(pp))*theta2)*sinc(fc*(a(qq)-a(pp))*theta1));
                elseif (mm*Ts+tau(qq))/(1+a(qq))<=(kk*Ts+tau(pp))/(1+a(pp)) && (((mm+1)*Ts+tau(qq))/(1+a(qq)))>((kk*Ts+tau(pp))/(1+aSR))
                   theta1=(((mm+1)*Ts+tau(qq))/(1+a(qq)))-((kk*Ts+tau(pp))/(1+a(pp))); 
                   theta2=(((mm+1)*Ts+tau(qq))/(1+a(qq)))+((kk*Ts+tau(pp))/(1+a(pp)));
                   F(kk+1,mm+1)=F(kk+1,mm+1)+conj(H(pp))*H(qq)*(theta1*exp(1i*pi*fc*(a(qq)-a(pp))*theta2)*sinc(fc*(a(qq)-a(pp))*theta1));
                end
            end
        end
    end
end

Hope it makes more sense now.

Can we still follow the same approach you did previously in this case?

Thanks

S. David 2014-6-7

Any hint?

请先登录，再进行评论。

Answer 2

the cyclist 2014-5-24

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/130980-how-to-do-this-more-efficiently#answer_138116

在 MATLAB Online 中打开

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.

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

S. David 2014-5-24

It is already there.

请先登录，再进行评论。

How to do this more efficiently

6 个评论
显示 4更早的评论隐藏 4更早的评论

采纳的回答

12 个评论
显示 10更早的评论隐藏 10更早的评论

更多回答（1 个）

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

另请参阅

类别

标签

产品

Community Treasure Hunt

How to do this more efficiently

6 个评论 显示 4更早的评论隐藏 4更早的评论

采纳的回答

12 个评论 显示 10更早的评论隐藏 10更早的评论

更多回答（1 个）

1 个评论 显示 -1更早的评论隐藏 -1更早的评论

另请参阅

类别

标签

产品

Community Treasure Hunt

6 个评论
显示 4更早的评论隐藏 4更早的评论

12 个评论
显示 10更早的评论隐藏 10更早的评论

1 个评论
显示 -1更早的评论隐藏 -1更早的评论