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

Question

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

3

我很感兴趣用Python计算一个由100,000个3D空间的粒子系统的功率谱。目前我发现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行和3列的形状。代码的输出给我了如下图所示的内容： enter image description here

从图中可以看出，我认为这给了我三个1D功率谱（每列数据一个），但实际上我想得到一个随半径变化的功率谱。

是否有人有任何建议或其他方法/软件包来计算功率谱（甚至可以接受两点自相关函数）。

- astromax

你是否正在对位置进行FFT？我不知道你为什么要这样做。FFT的假设是你的数据在某个域中被规则采样，并且数据与这些采样对齐。听起来像是你的粒子在各个地方，位置由某个数组定义（可能还有一个幅度数组？）。如果我理解正确，这很难处理。 - Henry Gomersall

好的，这是我困惑的一部分。我知道可以为一组粒子计算两点自相关函数。在星系调查中经常这样做。“对于给定的距离，两点自相关函数是一个关于一个变量（距离）的函数，它描述了两个星系被这个特定距离分开的概率。它可以被看作是一个不均匀因子 - 在某个距离尺度上值越高，宇宙在该距离尺度上就越不均匀。” http://en.wikipedia.org/wiki/Correlation_function_(astronomy) - astromax

如果我能计算自相关函数，那么功率谱就是它的傅里叶变换。因此理论上我相信可以计算已知三维坐标的粒子群的功率谱。 - astromax

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Jaime · Accepted Answer

您设置的方式并不完全有效...

您需要一个函数，我们将其称为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]的相邻点。因此，您的自相关可能会将两个非常远的星系显示为相邻的邻居。处理这种情况的最简单方法是沿着每个维度使您的数组增大一倍，用额外的零填充，然后丢弃额外的数据。

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