3

我对使用 Python 计算 3D 空间中粒子系统（~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是一个模拟数组，最终的形状是 100,000 x 3。代码的输出给了我类似的东西： 在此处输入图像描述

如您所见，我认为这给了我三个 1D 功率谱（我的数据的每列 1 个），但实际上我想要一个功率谱作为半径的函数。

有没有人知道计算功率谱的任何建议或替代方法/包（我什至会满足于两点自相关函数）。

4

1 回答 1

4

它并不完全按照你设置的方式工作......

你需要一个函数，我们称之为它f(x, y, z)，它描述了空间中的质量密度。在您的情况下，您可以将星系视为点质量，因此您将有一个以每个星系位置为中心的 delta 函数。正是对于这个函数，您可以计算三维自相关，从中可以计算出功率谱。

如果你想使用 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.correlatewithmode='constant'为您完成所有肮脏的工作。

于 2013-05-29T20:52:42.017 回答