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.
 Version of 13 Mar 2000.
 For more information about FFTLog, see http://casa.colorado.edu/~ajsh/FFTLog.
 Andrew J S Hamilton March 1999.
 Refs: [Talm78].
FFTLog computes a discrete version of the Hankel Transform (= FourierBessel Transform) with a power law bias \((k r)^q\)
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\)
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\)
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
to be about 1 (or 2, or pi, to taste).
The FFTLog algorithm is (see [Hami00]):
 FFT the input array \(a_j\) to obtain the Fourier coefficients \(c_m\) ;
 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+1x)/2]\) to obtain \(c_m u_m\);
 FFT \(c_m u_m\) back to obtain the discrete Hankel transform \(\tilde{a}_j\).
The Fourier sine and cosine transforms¶
may be regarded as special cases of the Hankel transform with \(\mu = 1/2\) and \(1/2\) since
The Fourier transforms may be done by making the substitutions
and Hankel transforming \(a(r)\) with a power law bias \((k r)^q\)
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¶
may be done by making the substitutions
and Hankel transforming \(a(r)\) with a power law bias \((k r)^q\)
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.
 subroutine `fhti(n,mu,q,dlnr,kr,kropt,wsave,ok)` is an initialization routine.
 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.
 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.
 subroutine `fhtq(n,a,dir,wsave)` computes the biased discrete Hankel transform of a logarithmically spaced periodic sequence. This is the basic FFTLog routine.
 real*8 function `krgood(mu,q,dlnr,kr)` takes an input kr and returns the nearest lowringing 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. selfinverse, provided also that, for n even, kr is lowringing. Nonzero 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+1j) = k_(n+1j) 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 lowringing kr, quietly;
 2 to change kr to nearest lowringing 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^{q1/2} \quad \text{and} \quad \tilde{A}(k) = \tilde{a}(k) k^{q1/2}\]and applying a biased Hankel transform to \(a(r)\).
The steps are: 1. \(a(r) = A(r) r^[dir (q0.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[(jjc) 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 lowringing kr if n is even.
Returns:  a : array
Transformed array Ã(k): a(j) is Ã(k_j) at k_j = k_c exp[(jjc) 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[(jjc) 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 lowringing kr if n is even.
Returns:  a : array
Transformed array Ã(k): a(j) is Ã(k_j) at k_j = k_c exp[(jjc) 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[(jjc) 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 lowringing kr if n is even.
Returns:  a : array
Transformed periodic array ã(k): a(j) is ã(k_j) at k_j = k_c exp[(jjc) 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. selfinverse, provided also that, for n even, kr is lowringing. Nonzero 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+1j) = k_(n+1j) r_j . Normally one would choose kr to be about 1 (default) (or 2, or pi, to taste).
Returns:  krgood : float
lowringing value of kr nearest to input kr. ln(krgood) is always within dlnr/2 of ln(kr).