测量生成的 3D 高斯随机场的功率谱（具有指定的功率谱）

Question

我试图生成一个功率谱为 P(k) = 1/k^2 的 3D 高斯随机场，然后测量生成的场的功率谱作为一致性检查（测量的功率谱当然应该匹配解析式 P(k) = 1/k^2)。为了生成 3D 场，我使用了代码 https://github.com/cphyc/FyeldGenerator/blob/master/FyeldGenerator/core.py，我在根据给定的 2D 方差创建 2D 高斯随机场中找到了该代码。为了测量功率谱，我使用了https://github.com/nualamccullagh/zeldovich-bao/blob/master/spatial_stats.py中的代码。

这是我使用的代码：

import matplotlib.pyplot as plt
import numpy as np
import sys
np.set_printoptions(threshold=sys.maxsize)
import six
import scipy.stats as stats

#Define global variables
ndim = 3
ngrid = boxsize = 128
n = 2
A = 1
shape = (ngrid, ngrid, ngrid)

def generate_field(statistic, power_spectrum, shape, unit_length=1,
                   fft=np.fft, fft_args=dict()):
    """
    Generates a field given a stastitic and a power_spectrum.
    """

    fftfreq = np.fft.fftfreq
    rfftfreq = np.fft.rfftfreq

    #Compute the k grid 
    all_k = [fftfreq(s, d=unit_length) for s in shape[:-1]] + \
            [rfftfreq(shape[-1], d=unit_length)]

    kgrid = np.meshgrid(*all_k, indexing='ij')
    knorm = np.sqrt(np.sum(np.power(kgrid, 2), axis=0))
    fourier_shape = knorm.shape

    fftfield = statistic(fourier_shape)

    power_k = np.where(knorm == 0, 0, np.sqrt(power_spectrum(knorm)))
    fftfield *= power_k

    return (fft.irfftn(fftfield), fftfield)

if __name__ == '__main__':
    def Pkgen(n):
        def Pk(k):
            return A*np.power(k, -n)

        return Pk

    def distrib(shape):
        # Build a unit-distribution of complex numbers with random phase
        a = np.random.normal(loc=0, scale=1, size=shape)
        b = np.random.normal(loc=0, scale=1, size=shape)
        return a + 1j * b

# density is the configuration-space density grid (real-valued, shape = ngrid x ngrid x ngrid)
# nkbins is the number of bins to compute the power spectrum in. It computes it in log-spaced bins in the range (2*Pi/L to Pi*ngrid / L)
def getPk(density, nkbins=100):
    #make sure the density has mean 0
    density=density-np.mean(density)
    ngrid=density.shape[0]
    
    #Fourier transform of density
    deltak=np.fft.rfftn(density)
    sk=deltak.shape
    #print('shape of deltak is', sk)
    
    #Square the density in Fourier space to get the 3D power, make sure k=0 mode is 0
    dk2=(deltak*np.conjugate(deltak)).astype(np.float)
    dk2[0,0,0]=0.0
    
    #set up k-grid
    kmin=2*np.pi/boxsize
    kny=np.pi*ngrid/boxsize
    
    a = np.fromfunction(lambda x,y,z:x, sk).astype(np.float)
    a[np.where(a > ngrid/2)] -= ngrid
    b = np.fromfunction(lambda x,y,z:y, sk).astype(np.float)
    b[np.where(b > ngrid/2)] -= ngrid
    c = np.fromfunction(lambda x,y,z:z, sk).astype(np.float)
    c[np.where(c > ngrid/2)] -= ngrid
    kgrid = kmin*np.sqrt(a**2+b**2+c**2).astype(np.float)
        
    #now we want to compute the 1-D power spectrum which involves averaging over shells in k-space 
    #define the k-bins we want to compute the power spectrum in
    binedges=np.logspace(np.log10(kmin), np.log10(kny),nkbins)
    numinbin=np.zeros_like(binedges)
    pk=np.zeros_like(binedges)
    kmean=np.zeros_like(binedges)
    kgridFlatten=kgrid.flatten()
    dk2 = dk2.flatten()
    index = np.argsort(kgridFlatten)
    
    kgridFlatten=kgridFlatten[index]
    dk2=dk2[index]
    c0=0.*c.flatten()+1.
    c0[np.where(c.flatten()==0.)]-=0.5
    c0=c0[index]
    cuts = np.searchsorted(kgridFlatten,binedges)
    
    for i in np.arange(0, nkbins-1):
        if (cuts[i+1]>cuts[i]):
            numinbin[i]=np.sum(c0[cuts[i]:cuts[i+1]])
            pk[i]=np.sum(c0[cuts[i]:cuts[i+1]]*dk2[cuts[i]:cuts[i+1]])
            kmean[i]=np.sum(c0[cuts[i]:cuts[i+1]]*kgridFlatten[cuts[i]:cuts[i+1]])

    wn0=np.where(numinbin>0.)
    pk=pk[wn0]
    kmean=kmean[wn0]
    numinbin=numinbin[wn0]
    
    pk/=numinbin
    kmean/=numinbin
    
    pk*= boxsize**3/ngrid**6
    
    return kmean, pk, kgrid, kmin, a, b, c


#call functions
densityRealField = np.real(generate_field(distrib, Pkgen(n), shape)[0])
km = getPk(densityRealField)[0]
PdensityMeasured = getPk(densityRealField)[1]

P_analytic = np.zeros(len(km))
#Analytic (input) Power Spectrum
for i in range(len(km)):
        P_analytic[i] = A/(km[i]**(n))

#plot both analytic and measured power spectrum
plt.clf()
line1 = plt.plot(km, P_analytic, color = 'cyan', linestyle = 'dashed', label = r'$P(k) \propto 1/k^{2}$')
line2 = plt.plot(km, PdensityMeasured, color = 'magenta', label = r'measured $P(k)$')
plt.legend()
plt.xscale('log')
plt.yscale('log')
plt.xlabel("$k$")
plt.ylabel("$P(k)$")
plt.tight_layout()
plt.savefig("P_measured_n=2_A=1.png", dpi = 300, bbox_inches = "tight")
plt.show()

问题在于测量的功率谱与解析（即输入）功率谱不一致，如下图所示： 测量的功率谱（洋红色线）的形状是正确的，但它应该位于解析功率谱（青色线）的正上方。我需要测量的功率谱始终与分析功率谱匹配（也就是说，即使我更改 ngrid 值）。

非常感谢任何帮助。

Answer 1

这不起作用的原因在于您的发行版的选择：

def distrib(shape):
    # Build a unit-distribution of complex numbers with random phase
    a = np.random.normal(loc=0, scale=1, size=shape)
    b = np.random.normal(loc=0, scale=1, size=shape)
    return a + 1j * b

该分布的均值在零附近。这意味着平均而言，频率的振幅（接近于零）低于功率谱密度所描述的振幅。您可以通过在域 [0, 2\pi) 上选择均匀分布的相位来解决此问题，如 Ruan & McLaughlin 的这篇论文 所示，请参阅方程 1。 15 或步骤 4 在同一页上。

因此，为了解决您的问题，请将发行版更改为：

def distrib(shape):
    # Build a unit-distribution of complex numbers with random phase
    theta = np.random.uniform(0, 2*np.pi, size=shape)
    return np.exp(1j*theta)

我们想要随机变量统一范数的原因是空间频率应该以功率谱密度决定的强度出现。

Answer 2

虽然我没有完全解决理论功率谱与测量功率谱不同的问题，但我确实取得了一些进展。

我们看到频谱（几乎）相对于波数绘制是恒定的。波动只是由于噪音造成的。

我们看到它很好地遵循幂律拟合函数

a x**b

，其中幂几乎恰好为 3。要解决此问题，请将

pk*= boxsize**3/ngrid**6

更改为

pk*= boxsize**6/ngrid**6

。

不幸的是，当理论功率谱的功率发生变化时，还存在另一个问题（对于您的情况，功率为 -2，因为 P(k) 与 1 / k^2 成正比）。 我真的不确定 -1.84 是从哪里来的，而且我看不出数学有什么其他问题。

我希望这有帮助！