Normalisation of 1D FFT in Numpy

Question

I am trying to measure the power spectrum of Gaussian realisations, but I do not know how to fix the normalisation of the FFTs. In the current version, the measured power spectrum is dependent on the Nmesh I use. I want the normalisation to be valid not only for the Gaussian field but also for even powers of the field. Does anyone know how to fix it?

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

norm_FFT = "ortho"

def MeasurePS(grid_abs_k, field, k_min, k_max, dk):
    """
    Function to measure the power spectrum
    """
    
    global norm_FFT
    
    #FFT of the input field
    _field_Fourier = np.fft.fft(field, norm = norm_FFT)*Normk
    
    #Compute the absolute value of the field in Fourier space
    _field_Fourier_abs = np.abs(_field_Fourier)**2
    
    #Generate the array that will be used to compute the power spectrum
    k_bins = np.arange(k_min, k_max + dk, dk)
    k_ctrs = 0.5*(k_bins[1:] + k_bins[:-1])
    
    #Generate the array to storage the power spectrum
    ps_output = np.zeros_like(k_ctrs)
    
    #Estimate the power spectrum
    for i in range(len(ps_output)):
        _ = _field_Fourier_abs
        mask_modes = (grid_abs_k > k_bins[i]) & (grid_abs_k <= k_bins[i+1])
        _ = _[mask_modes]
        ps_output[i] = np.average(_)
    return k_ctrs, ps_output


fig, ax = plt.subplots(1,3, figsize = (15,4))
for factor in [2**n for n in range(1,8)]:
    np.random.seed(1)

    #Input specs.
    Nmesh = 2**10*factor
    H = L/Nmesh #cell size (resolution)
    kf = 2*np.pi/L
    k_Ny = Nmesh*kf/2
    
    Normx = np.sqrt(Nmesh)
    Normk = 1.0/Normx

    #Compute the 1D momentum array
    k_array = np.fft.fftfreq(Nmesh)*2*np.pi/H
    grid_k_abs = np.sqrt(np.abs(k_array)**2)

    np.random.seed(1)
    gaussian_seed = np.random.normal(0,1,size=[Nmesh])
    
    #Take the fourier transform of the Gaussian seed
    gaussian_seed_fourier = np.fft.fft(gaussian_seed, norm = norm_FFT)*Normk

    #Compute the power Spectrum
    ps = linear_power(grid_k_abs,A,R)

    #Get the desired Gaussian field by multiplying the white-noise by sqrt(P(k))
    G_k = gaussian_seed_fourier*ps**0.5

    #Take the inverse transform
    G = np.fft.ifft(G_k, norm =  norm_FFT).real
    G *= Normx
    G2 = G**2 - np.mean(G**2)
    G4 = G**4 - np.mean(G**4)
    
    k, ps = MeasurePS(grid_k_abs,G, kf,k_Ny,70*kf)
    k, ps2 = MeasurePS(grid_k_abs,G2, kf,k_Ny,70*kf)
    k, ps4 = MeasurePS(grid_k_abs,G4, kf,k_Ny,70*kf)
    ax[0].plot(k,k*ps)
    ax[1].plot(k,k*ps2)
    ax[2].plot(k,k*ps4)
    
ax[0].set_xlim((0,15))
ax[0].set_xlabel("waveumber")
ax[0].set_title("Power spectrum of G")
ax[1].set_xlim((0,15))
ax[1].set_xlabel("waveumber")
ax[1].set_title("Power spectrum of $G^2$")
ax[2].set_xlim((0,15))
ax[2].set_title("Power spectrum of $G^4$")
ax[2].set_xlabel("waveumber")

I need to fix the normalisation in such a way that there is no dependence on the measured power spectrum on the Nmesh I use.