How to vectorize the antenna arrayfactor expression in matlab

Question

I have the antenna array factor expression here:

I have coded the array factor expression as given below:

lambda = 1;
M = 100;N = 200; %an M x N array
dx = 0.3*lambda; %inter-element spacing in x direction
m = 1:M; 
xm = (m - 0.5*(M+1))*dx; %element positions in x direction

dy = 0.4*lambda;
n = 1:N;
yn = (n - 0.5*(N+1))*dy;

thetaCount = 360; % no of theta values

thetaRes = 2*pi/thetaCount; % theta resolution

thetas = 0:thetaRes:2*pi-thetaRes; % theta values

phiCount = 180;

phiRes = pi/phiCount;

phis = -pi/2:phiRes:pi/2-phiRes;

cmpWeights = rand(N,M); %complex Weights

AF = zeros(phiCount,thetaCount); %Array factor

tic
for i = 1:phiCount
    for j = 1:thetaCount

        for p = 1:M
            for q = 1:N

                AF(i,j) = AF(i,j) + cmpWeights(q,p)*exp((2*pi*1j/lambda)*(xm(p)*sin(thetas(j))*cos(phis(i)) + yn(q)*sin(thetas(j))*sin(phis(i))));

            end
        end
    end
end

How can I vectorize the code for calculating the Array Factor (AF).

I want the line:

AF(i,j) = AF(i,j) + cmpWeights(q,p)*exp((2*pi*1j/lambda)*(xm(p)*sin(thetas(j))*cos(phis(i)) + yn(q)*sin(thetas(j))*sin(phis(i))));

to be written in vectorized form (by modifying the for loop).

Answer 1

Approach #1: Full-throttle

The innermost nested loop generates this every iteration - cmpWeights(q,p)*exp((2*pi*1j/lambda)*(xm(p)*sin(thetas(j))*cos(phis(i)) + yn(q)*sin(thetas(j))*sin(phis(i)))), which are to summed up iteratively to give us the final output in AF.

Let's call the exp(.... part as B. Now, B basically has two parts, one is the scalar (2*pi*1j/lambda) and the other part (xm(p)*sin(thetas(j))*cos(phis(i)) + yn(q)*sin(thetas(j))*sin(phis(i))) that is formed from the variables that are dependent on the four iterators used in the original loopy versions - i,j,p,q. Let's call this other part as C for easy reference later on.

Let's put all that into perspective:

• Loopy version had AF(i,j) = AF(i,j) + cmpWeights(q,p)*exp((2*pi*1j/lambda)*(xm(p)*sin(thetas(j))*cos(phis(i)) + yn(q)*sin(thetas(j))*sin(phis(i)))), which is now equivalent to AF(i,j) = AF(i,j) + cmpWeights(q,p)*B, where B = exp((2*pi*1j/lambda)*(xm(p)*sin(thetas(j))*cos(phis(i)) + yn(q)*sin(thetas(j))*sin(phis(i)))).

• B could be simplified to B = exp((2*pi*1j/lambda)* C), where C = (xm(p)*sin(thetas(j))*cos(phis(i)) + yn(q)*sin(thetas(j))*sin(phis(i))).

• C would depend on the iterators - i,j,p,q.

So, after porting onto a vectorized way, it would end up as this -

%// 1) Define vectors corresponding to iterators used in the loopy version
I = 1:phiCount;
J = 1:thetaCount;
P = 1:M;
Q = 1:N;

%// 2) Create vectorized version of C using all four vector iterators
mult1 = bsxfun(@times,sin(thetas(J)),cos(phis(I)).'); %//'
mult2 = bsxfun(@times,sin(thetas(J)),sin(phis(I)).'); %//'

mult1_xm = bsxfun(@times,mult1(:),permute(xm,[1 3 2]));
mult2_yn = bsxfun(@times,mult2(:),yn);
C_vect = bsxfun(@plus,mult1_xm,mult2_yn);

%// 3) Create vectorized version of B using vectorized C
B_vect = reshape(exp((2*pi*1j/lambda)*C_vect),phiCount*thetaCount,[]);

%// 4) Final output as matrix multiplication between vectorized versions of B and C
AF_vect = reshape(B_vect*cmpWeights(:),phiCount,thetaCount);

---

Approach #2: Less-memory intensive

This second approach would reduce the memory traffic and it uses the distributive property of exponential - exp(A+B) = exp(A)*exp(B).

Now, the original loopy version was this -

AF(i,j) = AF(i,j) + cmpWeights(q,p)*exp((2*pi*1j/lambda)*...
    (xm(p)*sin(thetas(j))*cos(phis(i)) + yn(q)*sin(thetas(j))*sin(phis(i))))

So, after using the distributive property, we would endup with something like this -

K = (2*pi*1j/lambda)
part1 = K*xm(p)*sin(thetas(j))*cos(phis(i));
part2 = K*yn(q)*sin(thetas(j))*sin(phis(i));
AF(i,j) = AF(i,j) + cmpWeights(q,p)*exp(part1)*exp(part2);

Thus, the relevant vectorized approach would become something like this -

%// 1) Define vectors corresponding to iterators used in the loopy version
I = 1:phiCount;
J = 1:thetaCount;
P = 1:M;
Q = 1:N;

%// 2) Define the constant used at the start of EXP() call
K = (2*pi*1j/lambda);

%// 3) Perform the sine-cosine operations part1 & part2 in vectorized manners
mult1 = K*bsxfun(@times,sin(thetas(J)),cos(phis(I)).'); %//'
mult2 = K*bsxfun(@times,sin(thetas(J)),sin(phis(I)).'); %//'

%// Perform exp(part1) & exp(part2) in vectorized manners
part1_vect = exp(bsxfun(@times,mult1(:),xm));
part2_vect = exp(bsxfun(@times,mult2(:),yn));

%// Perform multiplications with cmpWeights for final output
AF = reshape(sum((part1_vect*cmpWeights.').*part2_vect,2),phiCount,[])

---

Quick Benchmarking

Here are the runtimes with the input data listed in the question for the original loopy approach and proposed approach #2 -

---------------------------- With Original Approach
Elapsed time is 358.081507 seconds.

---------------------------- With Proposed Approach #2
Elapsed time is 0.405038 seconds.

The runtimes suggests a crazy performance improvement with Approach #2!