FFTLog-Test

This example is a translation of fftlogtest.f from the Fortran package FFTLog, which was presented in Appendix B of [Hami00] and published at http://casa.colorado.edu/~ajsh/FFTLog. It serves as an example for the python package pyfftlog (which is a Python version of FFTLog), in the same manner as the original file fftlogtest.f serves as an example for Fortran package FFTLog.

What follows is the original documentation from the file `fftlogtest.f`:

This is fftlogtest.f

This is a simple test program to illustrate how FFTLog works. The test transform is:

(1)\[\int^\infty_0 r^{\mu+1} \exp\left(-\frac{r^2}{2} \right)\ J_\mu(k, r)\ k\ {\rm d}r = k^{\mu+1} \exp\left(-\frac{k^2}{2} \right)\]

Disclaimer

FFTLog does NOT claim to provide the most accurate possible solution of the continuous transform (which is the stated aim of some other codes). Rather, FFTLog claims to solve the exact discrete transform of a logarithmically-spaced periodic sequence. If the periodic interval is wide enough, the resolution high enough, and the function well enough behaved outside the periodic interval, then FFTLog may yield a satisfactory approximation to the continuous transform.

Observe:

  1. How the result improves as the periodic interval is enlarged. With the normal FFT, one is not used to ranges orders of magnitude wide, but this is how FFTLog prefers it.
  2. How the result improves as the resolution is increased. Because the function is rather smooth, modest resolution actually works quite well here.
  3. That the central part of the transform is more reliable than the outer parts. Experience suggests that a good general strategy is to double the periodic interval over which the input function is defined, and then to discard the outer half of the transform.
  4. That the best bias exponent seems to be \(q = 0\).
  5. That for the critical index \(\mu = -1\), the result seems to be offset by a constant from the ‘correct’ answer.
  6. That the result grows progressively worse as mu decreases below -1.

The analytic integral above fails for \(\mu \le -1\), but FFTLog still returns answers. Namely, FFTLog returns the analytic continuation of the discrete transform. Because of ambiguity in the path of integration around poles, this analytic continuation is liable to differ, for \(\mu \le -1\), by a constant from the ‘correct’ continuation given by the above equation.

FFTLog begins to have serious difficulties with aliasing as \(\mu\) decreases below \(-1\), because then \(r^{\mu+1} \exp(-r^2/2)\) is far from resembling a periodic function. You might have thought that it would help to introduce a bias exponent \(q = \mu\), or perhaps \(q = \mu+1\), or more, to make the function \(a(r) = A(r) r^{-q}\) input to fhtq more nearly periodic. In practice a nonzero \(q\) makes things worse.

A symmetry argument lends support to the notion that the best exponent here should be \(q = 0,\) as empirically appears to be true. The symmetry argument is that the function \(r^{\mu+1} \exp(-r^2/2)\) happens to be the same as its transform \(k^{\mu+1} \exp(-k^2/2)\). If the best bias exponent were q in the forward transform, then the best exponent would be \(-q\) that in the backward transform; but the two transforms happen to be the same in this case, suggesting \(q = -q\), hence \(q = 0\).

This example illustrates that you cannot always tell just by looking at a function what the best bias exponent \(q\) should be. You also have to look at its transform. The best exponent \(q\) is, in a sense, the one that makes both the function and its transform look most nearly periodic.

import pyfftlog
import numpy as np
import matplotlib.pyplot as plt

Define the parameters you wish to use

The presets are the Reasonable choices of parameters from fftlogtest.f.

# Range of periodic interval
logrmin = -4
logrmax = 4

# Number of points (Max 4096)
n = 64

# Order mu of Bessel function
mu = 0

# Bias exponent: q = 0 is unbiased
q = 0

# Sensible approximate choice of k_c r_c
kr = 1

# Tell fhti to change kr to low-ringing value
# WARNING: kropt = 3 will fail, as interaction is not supported
kropt = 1

# Forward transform (changed from dir to tdir, as dir is a python fct)
tdir = 1

Calculate input function: \(r^{\mu+1}\exp\left(-\frac{r^2}{2}\right)\)

r = 10**(logrc + (np.arange(1, n+1) - nc)*dlogr)
ar = r**(mu + 1)*np.exp(-r**2/2.0)

Initialize FFTLog transform - note fhti resets kr

kr, xsave = pyfftlog.fhti(n, mu, dlnr, q, kr, kropt)
print(f"pyfftlog.fhti: new kr = {kr}")

Out:

pyfftlog.fhti: new kr = 0.9535389675791917

Call pyfftlog.fht (or pyfftlog.fhtl)

logkc = np.log10(kr) - logrc
print(f"Central point in k-space at log10(k_c) = {logkc}")

# rk = r_c/k_c
rk = 10**(logrc - logkc)

# Transform
# ak = pyfftlog.fftl(ar.copy(), xsave, rk, tdir)
ak = pyfftlog.fht(ar.copy(), xsave, tdir)

Out:

Central point in k-space at log10(k_c) = -0.020661554260541743

Calculate Output function: \(k^{\mu+1}\exp\left(-\frac{k^2}{2}\right)\)

k = 10**(logkc + (np.arange(1, n+1) - nc)*dlogr)
theo = k**(mu + 1)*np.exp(-k**2/2.0)

Plot result

plt.figure()

# Input
ax1 = plt.subplot(121)
plt.title(r'$r^{\mu+1}\ \exp(-r^2/2)$')
plt.xlabel('r')

plt.loglog(r, ar, 'k', lw=2)

plt.grid(axis='y', c='0.9')

# Transformed result
ax2 = plt.subplot(122, sharey=ax1)
plt.title(r'$k^{\mu+1} \exp(-k^2/2)$')
plt.xlabel('k')

plt.loglog(k, theo, 'k', lw=2, label='Theoretical')
plt.loglog(k, ak, 'r--', lw=2, label='FFTLog')

plt.legend()
plt.ylim([1e-8, 1e1])

ax2.yaxis.tick_right()
ax2.yaxis.set_label_position("right")
plt.grid(axis='y', c='0.9')

# Switch off spines
ax1.spines['top'].set_visible(False)
ax1.spines['right'].set_visible(False)
ax2.spines['top'].set_visible(False)
ax2.spines['left'].set_visible(False)
plt.tight_layout(rect=[0, 0, 1, .9])

# Main title
plt.suptitle(r"$\int_0^\infty r^{\mu+1}\ \exp(-r^2/2)\ J_\mu(k,r)\ " +
             r"k\ {\rm d}r = k^{\mu+1} \exp(-k^2/2)$", y=0.98)

plt.show()
../_images/sphx_glr_fftlogtest_001.png