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

[the cyclist]

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.

[S. David]

You are right. Here you go:

N=128; 
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);
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)))
                 theta=(((kk+1)*Ts+tau(pp)))-((mm*Ts+tau(qq)))); 
                 F(kk+1,mm+1)=F(kk+1,mm+1)+conj(H(pp))*H(qq)*(theta*exp(1i*pi*fc**theta)*sinc(fc*theta));
              end

[Jan]

The code does not run. Some end 's are missing an is not a valid Matlab operator.

[S. David]

This one is working.

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);
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))
                 theta=(((kk+1)*Ts+tau(pp)))-((mm*Ts+tau(qq))); 
                 F(kk+1,mm+1)=F(kk+1,mm+1)+conj(H(pp))*H(qq)*(theta*exp(1i*pi*fc*theta)*sinc(fc*theta));
              end
          end
      end
  end
end

[the cyclist]

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]

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

Answer 1

[Roger Stafford]

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

[Roger Stafford]

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]

Thanks for your replies, but could you please write beside each line in your code why you are doing so mathematically? Also, in my code I didn't include a second condition which can be seen from the equation I have attached, which is basically interchanging mm with kk. Is this second condition included in your code?

[Roger Stafford]

I don't see any attachment here.

[the cyclist]

Eq.png is attached in the comment area of the question.

[S. David]

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]

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]

I am sorry, but I still don't get it. As a start: How my conditions and yours are equivalent?

[Roger Stafford]

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]

@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]

Wow, now I get it, that is really elegant.

Do you think the round-off error would significantly affect the result? Because I see some elements in the original approach which don't appear in yours.

Thanks

[S. David]

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]

Any hint?

Answer 2

[the cyclist]

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.

[S. David]

It is already there.