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).