异形数组操作

Question

一些背景知识：我想计算 MxN 天线阵列的阵列因子，由以下等式给出：

其中w_i是第i个元素的复数权重，(x_i,y_i,z_i)是第i个元素的位置，k是wave number,theta,phi分别为仰角和方位角，i取值范围为0到MxN-1。

在我的代码中：

-theta 和 phi 是 np.mgrid，每个形状为 (200,200)，

-w_i，并且 (x,y,z)_i 是 np.array，每个形状为 (NxM,)

所以 AF 是一个 np.array，形状为 (200,200)（对 i 求和）。到目前为止没有问题，我可以很容易地得到 AF：

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))) )

现在，每个 w_i 都取决于频率，所以 AF 也是如此，现在我有 w_i 形状为 (NxM,1000)（我有每个 w_i 的 1000 个样本频率）。我尝试使用上面的代码更改

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

但我得到 'operands could not be broadcast together'。我可以通过对 1000 个值使用 for 循环来解决这个问题，但它很慢而且有点难看。那么，求和的正确方法是什么，或者正确定义 w_i 和 AF 的正确方法是什么？

如有任何帮助，我们将不胜感激。谢谢。

编辑我要添加的具有新维度的代码是下一个：

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

我试着用下一个主

测试它

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()

但我明白了

值错误：'operands could not be broadcast together with shapes (10,) (200,200) '

请注意，如果我设置

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

这个效果很好。

谢谢。

Answer 1

我看了一下代码，我认为这个 for 循环：

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))) )

在很多方面都是错误的。你在混合维度，你不能那样循环。
顺便说一句，为了充分利用 numpy 效率，从不循环数组。它会显着减慢执行速度。

我试着重做那部分。

首先，我建议您不要使用 from numpy import *，这是不好的做法（参见 here）。使用 import numpy as np。我重新引入了 np 缩写，因此您可以理解 numpy.

的含义

与频率无关的情况

第一个片段假定 w 是长度为 4 的一维数组：我忽略了 w 的频率依赖性，以向您展示如何在没有 for 循环并使用 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)

我在 np.newaxis 上使用 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。
所以，w[:,np.newaxis,np.newaxis] 是一个形状为 (4, 1, 1) 的数组。 x_pos[:,np.newaxis,np.newaxis]、y_pos[:,np.newaxis,np.newaxis] 和 z_pos[:,np.newaxis,np.newaxis] 也类似。由于角度的形状为 (200, 200)，因此可以进行广播，并且 af_points 的形状为 (4, 200, 200)。最后通过 np.sum 求和，对第一个轴求和以获得 (200, 200) 数组。

频率依赖案例

现在 w 的形状为 (4, 10)，其中 10 是频率点。这个想法是一样的，只是考虑频率是你的 numpy 数组中的一个额外维度：现在 af_points 将是一个形状数组 (4, 10, 200, 200) 其中 10 是 f_points你已经定义了。

为了便于理解，我拆分了计算：

#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)

现在 af 的形状为 (10, 200, 200)。

异形数组操作

Different shape arrays operations

numpy

python-3.x

array-broadcasting

与频率无关的情况

频率依赖案例