pyfftlog

Version: 0.2.0 ~ Date: 16 May 2020

pyfftlog - A python version of FFTLog

This is a python version of the logarithmic FFT code FFTLog as presented in Appendix B of Hamilton (2000) and published at casa.colorado.edu/~ajsh/FFTLog.

A simple f2py-wrapper (fftlog) can be found on github.com/prisae/fftlog. Tests have shown that fftlog is a bit faster than pyfftlog, but pyfftlog is easier to implement, as you only need NumPy and SciPy, without the need to compile anything.

I hope that FFTLog will make it into SciPy in the future, which will make this project redundant. (If you have the bandwidth and are willing to chip in have a look at SciPy PR #7310.)

Be aware that pyfftlog has not been tested extensively. It works fine for the test from the original code, and my use case, which is pyfftlog.fftl with mu=0.5 (sine-transform), q=0 (unbiased), k=1, kropt=1, and tdir=1 (forward). Please let me know if you encounter any issues.

Description of FFTLog from the FFTLog-Website

FFTLog is a set of fortran subroutines that compute the fast Fourier or Hankel (= Fourier-Bessel) transform of a periodic sequence of logarithmically spaced points.

FFTLog can be regarded as a natural analogue to the standard Fast Fourier Transform (FFT), in the sense that, just as the normal FFT gives the exact (to machine precision) Fourier transform of a linearly spaced periodic sequence, so also FFTLog gives the exact Fourier or Hankel transform, of arbitrary order m, of a logarithmically spaced periodic sequence.

FFTLog shares with the normal FFT the problems of ringing (response to sudden steps) and aliasing (periodic folding of frequencies), but under appropriate circumstances FFTLog may approximate the results of a continuous Fourier or Hankel transform.

The FFTLog algorithm was originally proposed by Talman (1978).

For the full documentation, see casa.colorado.edu/~ajsh/FFTLog.

Installation

You can install pyfftlog either via conda:

conda install -c conda-forge pyfftlog

or via pip:

pip install pyfftlog

License, Citation, and Credits

Released to the public domain under the CC0 1.0 License.

All releases have a Zenodo-DOI, which can be found on 10.5281/zenodo.3830364.

Be kind and give credits by citing Hamilton (2000). See the references-section in the manual for full references.

Manual and API

pyfftlog – Python version of FFTLog

This is a Python version of the FFTLog Fortran code by Andrew Hamilton, [Hami00].

The function scipy.special.loggamma replaces the file cdgamma.f in the original code, and the functions scipy.fftpack.rfft() and scipy.fftpack.irfft() replace the files drffti.f, drfftf.f, and drfftb.f in the original code.

The original documentation has just been adjusted where necessary, and put into a more pythonic format (e.g. using Parameters and Returns in the documentation’).

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

THE FFTLog CODE

FFTLog computes the discrete Fast Fourier Transform or Fast Hankel Transform (of arbitrary real index) of a periodic logarithmic sequence.

FFTLog computes a discrete version of the Hankel Transform (= Fourier-Bessel Transform) with a power law bias \((k r)^q\)

(1)\[\tilde{a}(k) = \int^\infty_0 a(r) (k r)^q J_{\mu} (k r) k \,dr \, ,\]
(2)\[a(r) = \int^\infty_0 \tilde{a}(k) (k r)^{-q} J_{\mu} (k r) r \,dk \, ,\]

where \(J_{\mu}\) is the Bessel function of order \(\mu\). The index \(\mu\) may be any real number, positive or negative.

The input array \(a_j\) is a periodic sequence of length \(n\), uniformly logarithmically spaced with spacing \(dlnr\)

(3)\[a_j = a(r_j) \quad \text{at} \quad r_j = r_c \exp[(j-j_c) dlnr]\]

centred about the point \(r_c\). The central index \(j_c = (n+1)/2\) is 1/2 integral if \(n\) is even. Similarly, the output array \(\tilde{a}_j\) is a periodic sequence of length \(n\), also uniformly logarithmically spaced with spacing \(dlnr\)

(4)\[\tilde{a}_j = \tilde{a}(k_j) \quad \text{at} \quad k_j = k_c \exp[(j-j_c) dlnr]\]

centred about the point \(k_c\).

The centre points \(r_c\) and \(k_c\) of the periodic intervals may be chosen arbitrarily; but it would be normal to choose the product

(5)\[kr = k_c r_c = k_j r_{(n+1-j)} = k_{(n+1-j)} r_j\]

to be about 1 (or 2, or pi, to taste).

The FFTLog algorithm is (see [Hami00]):

  1. FFT the input array \(a_j\) to obtain the Fourier coefficients \(c_m\) ;
  2. Multiply \(c_m\) by \(u_m = (kr)^{- i 2 m \pi/(n dlnr)} U_{\mu}[q + i 2 m \pi/(n dlnr)]\) where \(U_{\mu}(x) = 2^x \Gamma[(\mu+1+x)/2] / \Gamma[(\mu+1-x)/2]\) to obtain \(c_m u_m\);
  3. FFT \(c_m u_m\) back to obtain the discrete Hankel transform \(\tilde{a}_j\).
The Fourier sine and cosine transforms
(6)\[\tilde{A}(k) = \sqrt{2/\pi} \int^\infty_0 A(r) \sin(k r) \,dr \, ,\]
(7)\[\tilde{A}(k) = \sqrt{2/\pi} \int^\infty_0 A(r) \cos(k r) \,dr \, ,\]

may be regarded as special cases of the Hankel transform with \(\mu = 1/2\) and \(-1/2\) since

(8)\[\sqrt{2/\pi} \sin(x) = \sqrt(x) J_{1/2} (x) \, ,\]
(9)\[\sqrt{2/\pi} \cos(x) = \sqrt(x) J_{-1/2} (x) \, .\]

The Fourier transforms may be done by making the substitutions

(10)\[A(r) = a(r) r^{q-1/2} \quad \text{and} \quad \tilde{A}(k) = \tilde{a}(k) k^{-q-1/2}\]

and Hankel transforming \(a(r)\) with a power law bias \((k r)^q\)

(11)\[\tilde{a}(k) = \int^\infty_0 a(r) (k r)^q J_{\pm 1/2} (k r) k \,dr \, .\]

Different choices of power law bias \(q\) lead to different discrete Fourier transforms of \(A(r)\), because the assumption of periodicity of \(a(r) = A(r) r^{-q+(1/2)}\) is different for different \(q\).

If \(A(r)\) is a power law, \(A(r)\) proportional to \(r^{q-(1/2)}\), then applying a bias \(q\) yields a discrete Fourier transform \(\tilde{A}(k)\) that is exactly equal to the continuous Fourier transform, because then \(a(r)\) is a constant, which is a periodic function.

The Hankel transform
(12)\[\tilde{A}(k) = \int^\infty_0 A(r) J_{\mu} (k r) k \,dr\]

may be done by making the substitutions

(13)\[A(r) = a(r) r^q \quad \text{and} \quad \tilde{A}(k) = \tilde{a}(k) k^{-q}\]

and Hankel transforming \(a(r)\) with a power law bias \((k r)^q\)

(14)\[\tilde{a}(k) = \int^\infty_0 a(r) (k r)^q J_{\mu} (k r) k \,dr \, .\]

Different choices of power law bias \(q\) lead to different discrete Hankel transforms of \(A(r)\), because the assumption of periodicity of \(a(r) = A(r) r^{-q}\) is different for different \(q\).

If \(A(r)\) is a power law, \(A(r)\) proportional to \(r^q\), then applying a bias \(q\) yields a discrete Hankel transform \(\tilde{A}(k)\) that is exactly equal to the continuous Hankel transform, because then \(a(r)\) is a constant, which is a periodic function.

There are five routines:

Comments in the subroutines contain further details.

  1. subroutine `fhti(n,mu,q,dlnr,kr,kropt,wsave,ok)` is an initialization routine.
  2. subroutine `fftl(n,a,rk,dir,wsave)` computes the discrete Fourier sine or cosine transform of a logarithmically spaced periodic sequence. This is a driver routine that calls fhtq.
  3. subroutine `fht(n,a,dir,wsave)` computes the discrete Hankel transform of a logarithmically spaced periodic sequence. This is a driver routine that calls fhtq.
  4. subroutine `fhtq(n,a,dir,wsave)` computes the biased discrete Hankel transform of a logarithmically spaced periodic sequence. This is the basic FFTLog routine.
  5. real*8 function `krgood(mu,q,dlnr,kr)` takes an input kr and returns the nearest low-ringing kr. This is an optional routine called by fhti.

END of the original documentation from the file `fftlog.f`

pyfftlog.pyfftlog.fhti(n, mu, dlnr, q=0, kr=1, kropt=0)[source]

Initialize the working array xsave used by fftl, fht, and fhtq.

fhti initializes the working array xsave used by fftl, fht, and fhtq. fhti need be called once, whereafter fftl, fht, or fhtq may be called many times, as long as n, mu, q, dlnr, and kr remain unchanged. fhti should be called each time n, mu, q, dlnr, or kr is changed. The work array xsave should not be changed between calls to fftl, fht, or fhtq.

Parameters:
n : int

Number of points in the array to be transformed; n may be any positive integer, but the FFT routines run fastest if n is a product of small primes 2, 3, 5.

mu : float

Index of J_mu in Hankel transform; mu may be any real number, positive or negative.

dlnr : float

Separation between natural log of points; dlnr may be positive or negative.

q : float, optional

Exponent of power law bias; q may be any real number, positive or negative. If in doubt, use q = 0, for which case the Hankel transform is orthogonal, i.e. self-inverse, provided also that, for n even, kr is low-ringing. Non-zero q may yield better approximations to the continuous Hankel transform for some functions. Defaults to 0 (unbiased).

kr : float, optional

k_c r_c where c is central point of array = k_j r_(n+1-j) = k_(n+1-j) r_j . Normally one would choose kr to be about 1 (default) (or 2, or pi, to taste).

kropt : int, optional; {0, 1, 2, 3}
  • 0 to use input kr as is (default);
  • 1 to change kr to nearest low-ringing kr, quietly;
  • 2 to change kr to nearest low-ringing kr, verbosely;
  • 3 for option to change kr interactively.
Returns:
kr : float, optional

kr, adjusted depending on kropt.

xsave : array

Working array used by fftl, fht, and fhtq. Dimension: - for q = 0 (unbiased transform): n+3 - for q != 0 (biased transform): 1.5*n+4 If odd, last element is not needed.
pyfftlog.pyfftlog.fftl(a, xsave, rk=1, tdir=1)[source]

Logarithmic fast Fourier transform FFTLog.

This is a driver routine that calls fhtq().

fftl computes a discrete version of the Fourier sine (if mu = 1/2) or cosine (if mu = -1/2) transform

\[\tilde{A}(k) = \sqrt{2/\pi} \int^\infty_0 A(r) \sin(k r) \,dr \, ,\]
\[\tilde{A}(k) = \sqrt{2/\pi} \int^\infty_0 A(r) \cos(k r) \,dr \, ,\]

by making the substitutions

\[A(r) = a(r) r^{q-1/2} \quad \text{and} \quad \tilde{A}(k) = \tilde{a}(k) k^{-q-1/2}\]

and applying a biased Hankel transform to \(a(r)\).

The steps are: 1. \(a(r) = A(r) r^[-dir (q-0.5)]\) 2. call fhtq to transform \(a(r) \rightarrow \tilde{a}(k)\) 3. \(\tilde{A}(k) = \tilde{a}(k) k^[-dir (q+0.5)]\)

fhti must be called before the first call to fftl, with mu=1/2 for a sine transform, or mu=-1/2 for a cosine transform.

A call to fftl with dir=1 followed by a call to fftl with dir=-1 (and rk unchanged), or vice versa, leaves the array a unchanged.

Parameters:
a : array

Array A(r) to transform: a(j) is A(r_j) at r_j = r_c exp[(j-jc) dlnr], where jc = (n+1)/2 = central index of array.

xsave : array

Working array set up by fhti.

rk : float, optional

r_c/k_c = r_j/k_j (a constant, the same constant for any j); rk is not (necessarily) the same quantity as kr. rk is used only to multiply the output array by sqrt(rk)^dir, so if you want to do the normalization later, or you don’t care about the normalization, you can set rk = 1. Defaults to 1.

tdir : int, optional; {1, -1}
  • 1 for forward transform (default),
  • -1 for backward transform.

A backward transform (dir = -1) is the same as a forward transform with q -> -q and rk -> 1/rk, for any kr if n is odd, for low-ringing kr if n is even.
Returns:
a : array

Transformed array Ã(k): a(j) is Ã(k_j) at k_j = k_c exp[(j-jc) dlnr].
pyfftlog.pyfftlog.fht(a, xsave, tdir=1)[source]

Fast Hankel transform FHT.

This is a driver routine that calls fhtq().

fht computes a discrete version of the Hankel transform

\[\tilde{A}(k) = \int^\infty_0 A(r) J_{\mu} (k r) k \,dr \,\]

by making the substitutions

\[A(r) = a(r) r^q \quad \text{and} \quad \tilde{A}(k) = \tilde{a}(k) k^{-q}\]

and applying a biased Hankel transform to \(a(r)\).

The steps are: 1. \(a(r) = A(r) r^{-dir q}\) 2. call fhtq to transform \(a(r) \rightarrow \tilde{a}(k)\) 3. \(\tilde{A}(k) = \tilde{a}(k) k^{-dir q}\)

fhti must be called before the first call to fht.

A call to fht with dir=1 followed by a call to fht with dir=-1, or vice versa, leaves the array a unchanged.

Parameters:
a : array

Array A(r) to transform: a(j) is A(r_j) at r_j = r_c exp[(j-jc) dlnr], where jc = (n+1)/2 = central index of array.

xsave : array

Working array set up by fhti.

tdir : int, optional; {1, -1}
  • 1 for forward transform (default),
  • -1 for backward transform.

A backward transform (dir = -1) is the same as a forward transform with q -> -q, for any kr if n is odd, for low-ringing kr if n is even.
Returns:
a : array

Transformed array Ã(k): a(j) is Ã(k_j) at k_j = k_c exp[(j-jc) dlnr].
pyfftlog.pyfftlog.fhtq(a, xsave, tdir=1)[source]

Kernel routine of FFTLog.

This is the basic FFTLog routine.

fhtq computes a discrete version of the biased Hankel transform

\[\tilde{a}(k) = \int^\infty_0 a(r) (k r)^q J_{\mu} (k r) k \,dr \, .\]

fhti must be called before the first call to fhtq.

A call to fhtq with dir=1 followed by a call to fhtq with dir=-1, or vice versa, leaves the array a unchanged.

Parameters:
a : array

Periodic array a(r) to transform: a(j) is a(r_j) at r_j = r_c exp[(j-jc) dlnr] where jc = (n+1)/2 = central index of array.

xsave : array

Working array set up by fhti.

tdir : int, optional; {1, -1}
  • 1 for forward transform (default),
  • -1 for backward transform.

A backward transform (dir = -1) is the same as a forward transform with q -> -q, for any kr if n is odd, for low-ringing kr if n is even.
Returns:
a : array

Transformed periodic array ã(k): a(j) is ã(k_j) at k_j = k_c exp[(j-jc) dlnr].
pyfftlog.pyfftlog.krgood(mu, q, dlnr, kr)[source]

Return optimal kr.

Use of this routine is optional.

Choosing kr so that

\[(k r)^{- i pi/dlnr} U_{\mu}(q + i pi/dlnr)\]

is real may reduce ringing of the discrete Hankel transform, because it makes the transition of this function across the period boundary smoother.

Parameters:
mu : float

index of J_mu in Hankel transform; mu may be any real number, positive or negative.

q : float

exponent of power law bias; q may be any real number, positive or negative. If in doubt, use q = 0, for which case the Hankel transform is orthogonal, i.e. self-inverse, provided also that, for n even, kr is low-ringing. Non-zero q may yield better approximations to the continuous Hankel transform for some functions.

dlnr : float

separation between natural log of points; dlnr may be positive or negative.

kr : float, optional

k_c r_c where c is central point of array = k_j r_(n+1-j) = k_(n+1-j) r_j . Normally one would choose kr to be about 1 (default) (or 2, or pi, to taste).
Returns:
krgood : float

low-ringing value of kr nearest to input kr. ln(krgood) is always within dlnr/2 of ln(kr).

Examples

Note

Click here to download the full example code

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
Download all examples in Python source code: examples_python.zip
Download all examples in Jupyter notebooks: examples_jupyter.zip

Gallery generated by Sphinx-Gallery

References

[Hami00]Hamilton, A. J. S., 2000, Uncorrelated modes of the non-linear power spectrum: Monthly Notices of the Royal Astronomical Society, 312, pages 257–284; DOI: 10.1046/j.1365-8711.2000.03071.x; Website of FFTLog: casa.colorado.edu/~ajsh/FFTLog.
[Talm78]Talman, J. D., 1978, Numerical Fourier and Bessel transforms in logarithmic variables: Journal of Computational Physics, 29, pages 35–48; DOI: 10.1016/0021-9991(78)90107-9.

Changelog

v0.2.0 : First packaged release

2020-05-16

First packaged release on PyPi and conda-forge. This includes:

v0.1.1 : Bugfix uneven values

2019-08-16

  • Small bugfix for uneven values.
v0.1.0 : Initial upload to GitHub

2016-12-09

  • Initially working version uploaded to GitHub.