Incorrect behaviour of FFT for random density field generation

[zhaotingchen]

When writing some tests I noticed the power spectrum amplitude is incorrect for odd pb.N but correct for even. A code snippet to reproduce the wrong behaviour:

import numpy as np
import matplotlib.pyplot as plt
from powerbox.powerbox import PowerBox,LogNormalPowerBox
from powerbox.tools import get_power
import powerbox.dft as dft

pkfunc = lambda k: k**-2
boxtype = PowerBox
#boxtype = LogNormalPowerBox
for Nbox in (100,101):
    p = [0] * 40
    for i in range(40):
        pb = boxtype(Nbox, dim=2, pk=pkfunc, boxlength=1.0,ensure_physical=False,
                    )

        p[i], k = get_power(dft.fftshift(pb.delta_x()), pb.boxlength,
                           )
    plt.figure()
    plt.errorbar(k,np.array(p).mean(axis=0),yerr=np.array(p).std(axis=0))
    plt.plot(k,pkfunc(k))
    plt.xscale('log')
    plt.yscale('log')
    plt.title('ratio between in/out = %.2f' %np.mean(np.array(p).mean(axis=0)/pkfunc(k)))

You will see that for even N, the power is correct, whereas for odd N the power is exactly half of the input.

The issue seems to be that the Fourier conventions between Hermitian random array delta_k and the dft.ifft are not matching. To see that:

for Nbox in (100,101):
    pb = boxtype(Nbox, dim=2, pk=pkfunc, boxlength=1.0,ensure_physical=False,
                     a=0,b=2*np.pi,seed=41,
                    )
    delta_k = pb.delta_k()
    delta_x = pb.delta_x()
    kmode = pb.k()
    delta_x_test = dft.ifft(
                dft.ifftshift(delta_k),
                L=pb.boxlength,
                a=pb.fourier_a,
                b=pb.fourier_b,
                backend=pb.fftbackend,
            )[0].real
    print((np.abs(dft.fft(delta_x,L=1)[0])**2*kmode**2).mean())
    print((dft.fftshift(np.abs(dft.fft(delta_x_test,L=1)[0]))**2*kmode**2).mean())

When Nbox is 100, both output gives 1 which is desired. When Nbox is 101 on the other hand, the first output is 0.5. To see why, we can check the ifft of delta_k:

for Nbox in (100,101):
    pb = boxtype(Nbox, dim=2, pk=pkfunc, boxlength=1.0,ensure_physical=False,
                     a=0,b=2*np.pi,seed=42,
                    )
    delta_k = pb.delta_k()
    delta_x_test = dft.ifft(
                delta_k,
                L=pb.boxlength,
                a=pb.fourier_a,
                b=pb.fourier_b,
                backend=pb.fftbackend,
            )[0]
    print(np.abs(delta_x_test.imag).mean())
    delta_x_test = dft.ifft(
                dft.ifftshift(delta_k),
                L=pb.boxlength,
                a=pb.fourier_a,
                b=pb.fourier_b,
                backend=pb.fftbackend,
            )[0]
    print(np.abs(delta_x_test.imag).mean())

For Nbox=100, the imaginary component is small (but non-zero, maybe due to clipping the last element in delta_k when even?). For Nbox=101, the delta_x in the code is actually not real-valued but simply set to real later. Half of the amplitude is probably discarded here.

This does not give me a good idea how to fix it though. I do not understand why the even case works. But somewhere the Fourier convention must be messed up I think. I did check the code for dft and can't spot anything wrong, although I do not understand what is the purpose of _adjust_phase.

[zhaotingchen]

If you change the simulation to lognormal then the behaviour becomes harder to track but the idea is the same. I suspect it also relates to #12

[steven-murray]

@zhaotingchen thanks for pointing this out. I think this is a duplicate of #10, and I had a branch that attempted to fix that, but it got lost amongst other work I had to do in the meantime. Does that issue look like the same issue you're having? And do the arguments there make sense?

[steven-murray]

This is finally fixed in #134! Will be included in v1