Ceyron / exponax

Efficient Differentiable n-d PDE solvers in JAX.
https://fkoehler.site/exponax/
MIT License
6 stars 1 forks source link

Gaussian Random Field uses wrong exponent #9

Closed Ceyron closed 1 week ago

Ceyron commented 2 months ago

There is a $\dots/2$ factor in the exponent. I think the division by two is not necessary because we already take the wavenumber norm before.

Ceyron commented 1 week ago

Thinking about it further: The exponent seems correct with the factor of 1/2 because we want the power spectrum to follow a power-law.

For example, assume I would draw random white noise in terms of pixel-wise normal variates $u_i \propto \mathcal{N}$ (for all $i$ degrees of freedom) and transformed the corresponding state into Fourier space

$$ \hat{u}_h = \mathcal{F}_h(u_h) $$

Its power spectrum would then be

$$ P = \frac{1}{2} | \hat{u}_h |^2 $$

(important is that the power spectrum is the absolute squared). If I wanted to have this follow a power law with the wavenumber norm $| k |_2$, I'd have

$$ P \sim (| k |_2)^{-n} $$

Now, I want to design a filter on the original transformed white noise $\hat{u}_h$ (not on the power spectrum), I would have to further take the sqrt of the power law, i.e., we get a Gaussian Random Field with power law exponent $n$ if we draw a white noise (also possible directly in Fourier space) $\hat{u}_h$ and multiply it with the filter

$$ W =( | k | )^{-n/2} $$

Ceyron commented 1 week ago

Some more hand-wavy evidence, below is a series of Gaussian Random Fields created with Exponax under different power-law exponents and they (visually) match Fig. 2 of this blog post

image

Corresponding code:

generators = [ex.ic.GaussianRandomField(2, powerlaw_exponent=i, max_one=True) for i in range(0, 4+1)]
ic_s = [g(300, key=jax.random.PRNGKey(1)) for g in generators]
ic = jnp.concatenate(ic_s, axis=0)
ex.viz.plot_state_2d_facet(ic, titles=[f"Powerlaw exponent {i}" for i in range(0, 4+1)])