Different shape arrays operations

Question

A bit of background: I want to calculate the array factor of a MxN antenna array, which is given by the following equation:

$AF = \sum_{i}{w_i e^{k (x_i\cos{\phi}\sin{\theta}+y_i\sin{\theta}\sin{\phi}+z_i \cos{\theta}) }}$ Where w_i are the complex weight of the i-th element, (x_i,y_i,z_i) is the position of the i-th element, k is the wave number, theta and phi are the elevation and azimuth respectively, and i ranges from 0 to MxN-1.

In the code I have:

-theta and phi are np.mgrid with shape (200,200) each,

-w_i, and (x,y,z)_i are np.array with shape (NxM,) each

so AF is a np.array with shape (200,200) (sum over i).There is no problem so far, and I can get AF easily doing:

af = zeros([theta.shape[0],phi.shape[0]])
for i in range(self.size[0]*self.size[1]):
            af = af + ( w[i]*e**(-1j*(k * x_pos[i]*sin(theta)*cos(phi) + k * y_pos[i]* sin(theta)*sin(phi)+ k * z_pos[i] * cos(theta))) )

Now, each w_i depends on frequency, so AF too, and now I have w_i with shape (NxM,1000) (I have 1000 samples of each w_i in frequency). I tried to use the above code changing

af = zeros([1000,theta.shape[0],phi.shape[0]])

but I get 'operands could not be broadcast together'. I can solve this by using a for loop through the 1000 values, but it is slow and is a bit ugly. So, what is the correct way to do the summation, or the correct way to properly define w_i and AF ?

Any help would be appreciated. Thanks.

edit The code with the new dimension I'm trying to add is the next:

from numpy import *

class AntennaArray:
    def __init__(self,f,asize=None,tipo=None,dx=None,dy=None):
        self.Lambda = 299792458 / f
        self.k = 2*pi/self.Lambda
        self.size = asize
        self.type = tipo
        self._AF_DATA_SIZE = 200
        self.theta,self.phi =  mgrid[0 : pi : self._AF_DATA_SIZE*1j,0 : 2*pi : self._AF_DATA_SIZE*1j]

        self.element_pos = None
        self.element_amp = None
        self.element_pha = None

        if dx == None:
            self.dx = self.Lambda/2
        else:
            self.dx = dx

        if dy == None:
            self.dy = self.Lambda/2
        else:
            self.dy = dy


        self.generate_array()


    def generate_array(self):


        M = self.size[0]
        N = self.size[1]
        dx = self.dx
        dy = self.dy

        x_pos = arange(0,dx*N,dx)
        y_pos = arange(0,dy*M,dy)
        z_pos = 0

        ele = zeros([N*M,3])

        for i in range(M):
            ele[i*N:(i+1)*N,0] = x_pos[:]


        for i in range(M):
            ele[i*N:(i+1)*N,1] = y_pos[i]


        self.element_pos = ele


        #self.array_factor = self.calculate_array_factor()

    def calculate_array_factor(self):

        theta,phi = self.theta,self.phi

        k = self.k

        x_pos = self.element_pos[:,0]
        y_pos = self.element_pos[:,1]
        z_pos = self.element_pos[:,2]

        w = self.element_amp*exp(1j*self.element_pha)



        if len(self.element_pha.shape) > 1:
            #I have f_size samples of w_i(f)
            f_size = self.element_pha.shape[1]
            af = zeros([f_size,theta.shape[0],phi.shape[0]])
        else:
            #I only have w_i
            af = zeros([theta.shape[0],phi.shape[0]])


        for i in range(self.size[0]*self.size[1]):
        **strong text**#This for loop does the summation over i
            af = af + ( w[i]*e**(-1j*(k * x_pos[i]*sin(theta)*cos(phi) + k * y_pos[i]* sin(theta)*sin(phi)+ k * z_pos[i] * cos(theta))) )
        return af

I tried to test it with the next main

from numpy import *
f_points = 10

M = 2
N = 2

a = AntennaArray(5.8e9,[M,N])

a.element_amp = ones([M*N,f_points])
a.element_pha = zeros([M*N,f_points])

af = a.calculate_array_factor()

But I get

ValueError: 'operands could not be broadcast together with shapes (10,) (200,200) '

Note that if I set

a.element_amp = ones([M*N])
a.element_pha = zeros([M*N])

This works well.

Thanks.

Why the for loop to implement the equation? You said you have numpy arrays, you don't need the for loop. And on what? Seems you have flattened everything. This is totally not needed. What is self? Is that inside a method in a class? If you give us more background on your code, is better. — Valentino
– Valentino, Commented Oct 29, 2019 at 23:02
I have updated the post with the complete code. The for is for the summation only. The others arrays are only 1D because I am summing over them, and the final shape is the (200,200) of the numpy.mgrid. Thanks. — rdm
– rdm, Commented Oct 30, 2019 at 13:48

Community · Accepted Answer · 2020-06-20 09:12:55Z

I had a look at the code, and I think this for loop:

af = zeros([theta.shape[0],phi.shape[0]])
for i in range(self.size[0]*self.size[1]):
            af = af + ( w[i]*e**(-1j*(k * x_pos[i]*sin(theta)*cos(phi) + k * y_pos[i]* sin(theta)*sin(phi)+ k * z_pos[i] * cos(theta))) )

is wrong in many ways. You are mixing dimensions, you cannot loop that way.
And by the way, to make full use of numpy efficiency, never loop over the arrays. It slows down the execution significantly.

I tried to rework that part.

First, I advice you to not use from numpy import *, it's bad practice (see here). Use import numpy as np. I reintroduced the np abbreviation, so you can understand what comes from numpy.

Frequency independent case

This first snippet assumes that w is a 1D array of length 4: I am neglecting the frequency dependency of w, to show you how you can get what you already obtained without the for loop and using instead the power of numpy.

af_points = w[:,np.newaxis,np.newaxis]*np.e**(-1j*
           (k * x_pos[:,np.newaxis,np.newaxis]*np.sin(theta)*np.cos(phi) + 
            k * y_pos[:,np.newaxis,np.newaxis]*np.sin(theta)*np.sin(phi) + 
            k * z_pos[:,np.newaxis,np.newaxis]*np.cos(theta)
           ))
af = np.sum(af_points, axis=0)

I am using numpy broadcasting to obtain a 3D array named af_points, whose shape is (4, 200, 200). To do it, I use np.newaxis to extend the number of axis of an array in order to use broadcasting correctly. More here on np.newaxis.
So, w[:,np.newaxis,np.newaxis] is an array of shape (4, 1, 1). Similarly for x_pos[:,np.newaxis,np.newaxis], y_pos[:,np.newaxis,np.newaxis] and z_pos[:,np.newaxis,np.newaxis]. Since the angles have shape (200, 200), broadcasting can be done, and af_points has shape (4, 200, 200). Finally the sum is done by np.sum, summing over the first axis to obtain a (200, 200) array.

Frequency dependent case

Now w has shape (4, 10), where 10 are the frequency points. The idea is the same, just consider that the frequency is an additional dimension in your numpy arrays: now af_points will be an array of shape (4, 10, 200, 200) where 10 are the f_points you have defined.

To keep it understandable, I've split the calculation:

#exp_point is only the exponent, frequency independent. Will be a (4, 200, 200) array.
exp_points = np.e**(-1j*
            (k * x_pos[:,np.newaxis,np.newaxis]*np.sin(theta)*np.cos(phi) + 
             k * y_pos[:,np.newaxis,np.newaxis]*np.sin(theta)*np.sin(phi) +
             k * z_pos[:,np.newaxis,np.newaxis]*np.cos(theta)
            ))
af_points = w[:,:,np.newaxis,np.newaxis] * exp_points[:,np.newaxis,:,:]
af = np.sum(af_points, axis=0)

And now af has shape (10, 200, 200).

Now I understand what you mean with flattened. I did not understand the concept of broadcasting until your explanation. Thanks a lot.

Collectives™ on Stack Overflow

Different shape arrays operations

1 Answer 1

Frequency independent case

Frequency dependent case

1 Comment

Your Answer

Linked

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

Frequency independent case

Frequency dependent case

1 Comment

Your Answer

Sign up or log in

Post as a guest

Linked

Related