Numpy中数据的功率谱和自相关

Question

我对用Python计算三维空间中粒子系统(~100,000)的功率谱很感兴趣。到目前为止，我发现的是Numpy中的一组函数(fft，fftn，.)计算离散傅里叶变换，其中绝对值的平方是功率谱。我的问题是我的数据是如何表示的--而且说实话，答案可能相当简单。

我所拥有的数据结构是一个形状为(n,2)的数组，n是我所拥有的粒子数，每一列代表n个粒子的x、y和z坐标。我认为我应该使用的函数是fftn()函数，它采用了n维数组的离散傅里叶变换，但它并没有提到格式。如何将数据表示为要输入fftn的数据结构

下面是我到目前为止尝试测试该函数的内容：

import numpy as np
import random
import matplotlib.pyplot as plt

DATA = np.zeros((100,3))

for i in range(len(DATA)):
    DATA[i,0] = random.uniform(-1,1)
    DATA[i,1] = random.uniform(-1,1)
    DATA[i,2] = random.uniform(-1,1)

FFT = np.fft.fftn(DATA)
PS = abs(FFT)**2

plt.plot(PS)
plt.show()

名为DATA的数组是一个模拟数组，最终的形状是100000×3。代码的输出给我提供了如下内容：

​

正如你所看到的，我认为这是给我三个一维功率谱(我的数据的每一列1)，但是我真的想要一个功率谱作为半径的函数。

有没有人有任何建议或替代的方法/包，他们知道计算功率谱(我甚至满足于两点自相关函数)。

Answer 1

它不太像你设定的方式.

你需要一个函数，我们把它叫做f(x, y, z)，它描述了空间中质量的密度。在你的例子中，你可以把星系看作点质量，所以在每个星系的位置都有一个以δ函数为中心的星系。对于这个函数，你可以计算三维自相关，你可以用它来计算功率谱。

如果您想要使用numpy为您这样做，您将首先必须离散您的函数。一个可能的模拟例子是：

import numpy as np
import matplotlib.pyplot as plt

space = np.zeros((100, 100, 100), dtype=np.uint8)

x, y, z = np.random.randint(100, size=(3, 1000))
space[x, y, z] += 1

space_ps = np.abs(np.fft.fftn(space))
space_ps *= space_ps

space_ac = np.fft.ifftn(space_ps).real.round()
space_ac /= space_ac[0, 0, 0]

现在，space_ac拥有数据集的三维自相关函数.这并不是你想要的，要得到一维的相关函数，你必须把球壳的数值平均在原点附近：

dist = np.minimum(np.arange(100), np.arange(100, 0, -1))
dist *= dist
dist_3d = np.sqrt(dist[:, None, None] + dist[:, None] + dist)
distances, _ = np.unique(dist_3d, return_inverse=True)
values = np.bincount(_, weights=space_ac.ravel()) / np.bincount(_)

plt.plot(distances[1:], values[1:])

这样做还有另一个问题:当你像上面那样计算功率谱时，数学上就好像你的三维阵列包裹在边界上，即点[999, y, z]是[0, y, z]的邻居。所以你的自相关可以显示出两个非常遥远的星系是近邻。处理这一问题的最简单方法是，使数组沿每个维度的大小增加一倍，填充额外的零，然后丢弃额外的数据。

或者，您可以使用scipy.ndimage.filters.correlate和mode='constant'为您完成所有的脏工作。