0

I am trying to generate a 3D gaussian random field with a power spectrum of P(k) = 1/k^2, and then measure the power spectrum of the generated field as a consistency check (the measured power spectrum should of course match the analytic one, P(k) = 1/k^2). To generate the 3D field, I used the code https://github.com/cphyc/FyeldGenerator/blob/master/FyeldGenerator/core.py which I found in Creating a 2D Gaussian random field from a given 2D variance. To measure the power spectrum, I used the code in https://github.com/nualamccullagh/zeldovich-bao/blob/master/spatial_stats.py.

Here is the code I used:

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

The problem is that the measured power spectrum does not agree with the analytic (i.e. input) power spectrum, as can be seen in this plot: plot The shape of the measured power spectrum (magenta line) is correct, but it should lie right on top of the analytic one (cyan line). I need the measured power spectrum to ALWAYS match the analytic one (that is, even if I change ngrid values).

Any help is greatly appreciated.

Sana Elgamal
  • 1
  • 1
  • Could you please copy-paste the code without line numbers? This numbering makes it harder for anyone wanting to replicate your problem by running your code. – Cris Luengo Jul 06 '21 at 19:38
  • It looks by eye as if the two curves differ by a scale factor somewhat larger than 10**4. Could it be `ngrid**2` (16384)? I haven't done much with FFT's for over 30 years, and never dealt with multidimensional spectra, but are you sure about the divisor on line 111? – pjs Jul 06 '21 at 20:17
  • @CrisLuengo thank you for pointing this out, now the code is without the line numbers. – Sana Elgamal Jul 07 '21 at 10:50
  • Hey @SanaElgamal! Have you found a solution that works? Because I would really need that as well :) Thanks – Valentina May 06 '22 at 12:45

1 Answers1

0

The reason this isn't working lies in the choice of your distribution:

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

The mean of this distribution lies around zero. This means that on average the frequencies will have lower amplitudes (near zero) than described by your power spectral density. You can fix this by choosing a uniformly distributed phase on the domain [0, 2\pi) as in this paper by Ruan & McLaughlin, see Eq. 15 or step 4 on the same page.

So to solve your issue change the distribution to:

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)

The reason we want unity norm of the random variable is that the spatial frequency should be present at the strength dictated by the power spectral density.

Mio
  • 101
  • 3