Deviations from Parseval's theorem

[egorssed]

Hey, thanks for the package.

Though there is a bug(feature) in your code.
If you consider Parseval's theorem you can theoretically compute the variance of the generated GRF from its power spectrum.
In your case, the variance of the generated field quite significantly deviates from the theoretical one (right plot). For big-scale perturbations (big power slope) deviation can even reach 300% ! Though, if you average generated variance over ~100 phases you will get the theoretical one.

I assume it is related to the way how you generate phases of the Fourier image. One can either sample phases using straightforward random.uniform or he can use more sophisticated, yet robust, Polar Box-Muller transform.

The figure represents how deviation between generated variance and one from Parseval's theorem depends on the Power Slope of the GRF. (figures were computed for 100x100 pixels GRF and 100 random seeds for every Power slope)

If my guess is correct you might want to change these two functions

powerbox/powerbox/powerbox.py

Line 183 in 315e30b

def gauss_hermitian(self):

powerbox/powerbox/dft.py

Line 247 in 315e30b

def _adjust_phase(ft, left_edge, freq, axes, b):

[steven-murray]

Thanks @egorssed, this is very interesting. I'll have to look at it more closely. I am tied up for the next week, but hopefully will get time to look at it after that. In the meantime, I am more than willing to take PRs :-P

[ajouellette]

I'm potentially interested in looking into this. Is there a reference somewhere for the sampling method? I'm not sure I understand the use of the Box-Muller transform. As far as I can tell, the Box-Muller transform generates normally distributed variables from uniformly distributed ones. How does that apply to the Fourier phases?

[steven-murray]

@ajouellette I haven't had the time to look into this as yet, sorry. @egorssed it would be great to have a reference to detail what you outlined above!

[egorssed]

Hey, It's been half a year. I don't remember your problems already.

@ajouellette
When you simulate field from a Fourier image, this Fourier image is somewhat F(k) = sqrt(P(k)/2)*(R+i*I)
Here F(k) is Fourier image in wavenumber k, P is power spectrum. R and I are standard normal variables that actually define the stochastic realisation of your random field. The trick is to sample those two values independently.
Also think that a Fourier transform of a Gaussian is a Gaussian. So if you want a Gaussian random field in the configuration space, then its Fourier image has to be Gaussian as well.

http://www.gaussianprocess.org/gpml/chapters/RW.pdf

I believe you can try to dig into my code if you want.

[ajouellette]

Thanks @egorssed. I realize it's been a while, but do you have any more information about how you generated the plots above?

[egorssed]

Hmm, I believe I sampled fields with Box-Muller and with this package. Computed variances of every field and compared.
Tho, I would advice to carry out experiment where only sampling varies, not the entire code.