
The URL of this page is http://casa.colorado.edu/~ajsh/FFTLog/ .
Questions? Comments? Email me at Andrew.Hamilton@Colorado.EDU.
This page contains mathematical symbols and Greek letters, and some letters from the extended Latin ISO88591 character set. If these do not show correctly, then
Since Netscape and other browsers do not yet implement symbols correctly in printed copies, if you want a printed version of this document you will be better off printing this printable gzip'd PostScript version (129K; version of 13 March 2000). The HTML document you are reading was adapted from the Appendix of Hamilton (2000), and formal references to FFTLog should refer to that paper.
The FFTLog code is built on top of the NCAR suite of FFT routines, and on a modified version of an implementation of the complex Gammafunction from Takuya Ooura's gamerf package. The required routines are included in fftlog.tar.gz .
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 (§5), so also FFTLog gives the exact Fourier or Hankel transform, of arbitrary order m, of a logarithmically spaced periodic sequence (§6).
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 (§10).
The FFTLog algorithm was originally proposed by Talman (1978).
In cosmology, fluctuations in the matter density of the Universe are thought to have been laid down during an inflationary epoch in the first few moments following the Big Bang (Turner 1997). Vacuum fluctuations in the field that drives inflation should produce a Gaussian distribution of density fluctuations with a near scaleinvariant power spectrum P(k) µ k. That primordial spectrum was processed prior to Recombination by the action of gravity modulated by the pressure of radiation. Following Recombination, when the Universe was about 300,000 years old, the matter power spectrum was further processed by nonlinear gravitational clustering, up to the present time.
The cosmological power spectrum P(k),
a function of wavenumber k,
is the 3dimensional Fourier transform
of the cosmological correlation function
x(r),
a function of spatial separation r.
With the conventional normalization used by cosmologists,
 (1) 
Figure 1: Cosmological power spectrum P(k) of matter fluctuations predicted by the socalled LCDM model, a flat (W = 1) Universe dominated by a cosmological constant (W_{L} = 0.7), and Cold Dark Matter (W_{m} = 0.3), including a sprinkling of baryons (W_{b} = 0.05). The LCDM power spectrum was computed from the formulae of Eisenstein & Hu (1998), nonlinearly evolved according to the formula of Peacock & Dodds (1996). The spectrum is normalized to the amplitude of fluctuations observed by the COBE satellite. 
Figure 2:
Cosmological correlation function x(r)
corresponding to the
LCDM power spectrum
shown in Figure 1.
The top panel of Figure 2 shows the correlation function x(r) computed with FFTLog at two different resolutions, plotted on top of each other: (red, low resolution) with 96 points over the range r = 10^{3} to 10^{3} h^{1}Mpc, and (blue, high resolution) with 768 points over the range r = 10^{6} to 10^{6} h^{1}Mpc. The lines are dashed where the correlation function is negative, at separations r > 119 h^{1}Mpc. The low and high resolution curves are almost indistinguishable except at r > 200 h^{1}Mpc, where the low resolution curve goes to a constant, while the high resolution curve declines as a power law ~ r^{4}. The disagreement is caused by aliasing (see §10) of small and large separations in the low resolution case. Aliasing is almost eliminated in the high resolution case because the range r = 10^{6} to 10^{6} h^{1}Mpc over which the transform was computed is much broader than the range plotted. The straight dashed line shows a power law (r/5 h^{1}Mpc)^{1.8} for reference. Both low and high resolution cases used an unbiased (q = 0) transform, §9, and a lowringing value of k_{0} r_{0}, §8 (actually the choice of k_{0} r_{0} made little difference here). The middle panel of Figure 2 shows the ratio x_{low}/x_{high} of the low to high resolution correlation functions. The bottom panel of Figure 2 shows the ratio x_{FFT}/x_{FFTLog} of the correlation function x_{FFT} computed with a normal FFT (sine transform) with 1023 points over the range r = 0.125 to 128 h^{1}Mpc, to the high resolution correlation function x_{FFTLog} computed with FFTLog. The FFT'd correlation function x_{FFT} rings at the ±5% level. 
In this particular instance, FFTLog outperforms the normal FFT on all counts: it is more accurate, with fewer points, over a larger range, and it shows no signs of ringing. This does not mean that FFTLog is always better than FFT. Rather, FFTLog is well matched to the problem at hand: the cosmological power spectrum extends over many orders of magnitude in wavenumber k, and varies smoothly in lnk.
Consider the continuous Hankel (= FourierBessel) transform pair
 (2) 
 (3) 
 (4) 
Fourier sine and cosine transforms can be regarded as special
cases of Hankel transforms with
m = ±1/2,
since
 (5) 
As first noted by Siegman (1977), if the product kr in the Hankel transform is written as e^{lnk + lnr}, then the transform becomes a convolution integral in the integration variable lnr or lnk. Convolution is equivalent to multiplication in the corresponding Fourier transform space. Thus the Hankel transform can be computed numerically by the algorithm:
An advantage of FFTLog, emphasized by Talman (1978), is that the order m of the Bessel function may be any arbitrary real number. In particular, FFTLog works for 1/2integral m, so includes the cases of Fourier sine and cosine transforms, and spherical Hankel transforms involving the spherical Bessel functions j_{l}(x) º (p/2 x)^{1/2} J_{l+1/2}(x).
 (6) 
 (7) 
 (8) 
The sampling theorem (e.g. Press et al. 1986 §12.1)
asserts that,
given a function a(r)
satisfying equation (7),
the Fourier coefficients c_{m}
can be expressed in terms of the values
a_{n} º a(r_{n})
of the function a(r)
at the N discrete points
r_{n} = n R / N
for
n = 0, ±1, ..., ±[N/2].
For even N, the periodicity of a(r) ensures that
a_{N/2} = a_{N/2}.
Specifically, the sampling theorem asserts that
the Fourier coefficients in the expansion (7) satisfy
 (9) 
 (10) 
Equations (9) and (10) constitute a discrete Fourier transform pair relating two periodic, linearly spaced sequences a_{n} and c_{m} of length N. The standard FFT evaluates the discrete Fourier transform exactly (that is, to machine precision).
 (11) 
 (12) 
 (13) 
 (14) 
The continuous Hankel transform ã(k),
equation (2),
of a function a(r) of the form (12) is
 (15) 
 (16) 
 (17) 
 (18) 
The sampling theorem requires that
u_{N/2} = u_{N/2}
for even N,
which is not necessarily satisfied by equation (18).
However,
at the discrete points
k_{n} = k_{0} e^{nL/N}
considered by the sampling theorem, the contributions at
m = ±N/2
to the sum on the right hand side of equation (17)
are
()^{n} c_{N/2} (u_{N/2} + u_{N/2}*)/2 = ()^{n} c_{N/2} Re u_{N/2} .
Thus the equality (17) remains true
at the discrete points k_{n}
if u_{±N/2}
are replaced by their real parts,
 (19) 
 (20) 
 (21) 
Putting together
equations (13), (14),
(20) and (21)
yields the discrete Hankel transform pair
 (22) 
 (23) 
 (24) 
 (25) 
 (26) 
FFTLog evaluates the forward and inverse discrete Hankel transforms given by equations (22), (23), exactly (to machine precision).
The reciprocal
1/u_{m}(m,q)
in equation (25)
is equal to
u_{m}(m,q),
according to
equations (16) and (18),
 (27) 
In the continuous case, the inverse Hankel transform is equal to the forward transform with q ® q, equations (2). In the discrete case this remains true for odd N, but it is not generally true for even N (the usual choice) except in the important special case discussed in §8.
In the general discrete case (i.e. if the condition [28] in §8 is not satisfied), the inverse discrete Hankel mode v^{}_{n}(m,q), equation (25), differs from the forward Hankel mode v^{+}_{n}(m,q), equation (24), only for even N and only in the coefficient of the highest frequency Fourier component, 1/u_{m}(m,q) versus u_{m}(m,q) for m = ±N/2. To the extent that the highest frequency Fourier coefficient c_{±N/2} of a sequence a_{n} is small, the difference between its inverse discrete Hankel transform and its forward transform with q ® q should be small.
It is possible for the inverse discrete Hankel transform to be singular, if u_{±N/2} is purely imaginary, so that its real part vanishes, making v^{}_{n}(m,q) singular. As discussed in §8, this singularity can be avoided by choosing a lowringing value of k_{0} r_{0}, equation (30).
The forward (inverse) discrete Hankel transforms are also singular at special values of m and q, namely where m+1+q (or m+1q in the inverse case) is zero or an even negative integer, because u_{0}(m,q) = U_{m}(q) is singular at these points. This singularity reflects a real singularity in the corresponding continuous Hankel transform (unlike the singularity of the previous paragraph, which is an avoidable artefact of discreteness). The singularity in u_{0} leads to an additive infinite constant in the discrete Hankel transform. In physical problems this additive infinite constant may somehow cancel out (for example, in the difference between two Hankel transforms). FFTLog's strategy in these singular cases is to evaluate the discrete Hankel transform with the infinite constant set to zero, and to issue a warning.
multiply by u_{m} given by equations (18) and (19) to obtain c_{m} u_{m};
FFT c_{m} u_{m} back to obtain the discrete Hankel transform ã_{n}, equation (21).
A variant of the algorithm is to sandwich the above operations with power law biasing and unbiasing operations. For example, one way to take the unbiased continuous Hankel transform Ã(k) of a function A(r), equation (4), is to bias A(r) and Ã(k) with power laws, equation (3), and take a biased Hankel transform, equation (2). The discrete equivalent of this is:
FFT a_{n} to obtain the Fourier coefficients c_{m}, equation (13);
multiply by u_{m} given by equations (18) and (19) to obtain c_{m} u_{m};
FFT c_{m} u_{m} back to obtain the discrete Hankel transform ã_{n}, equation (21);
Unbias ã_{n} with a power law to obtain Ã_{n} = ã_{n}k_{n}^{q}, equation (3).
Although in the continuous limit the result would be identical to an unbiased Hankel transform, in the discrete case the result differs. With a simple unbiased discrete Hankel transform, it is the sequence A_{n} that is taken to be periodic, whereas in the algorithm above it is not A_{n} but rather a_{n} that is periodic.
The inverse discrete Hankel transform is accomplished by the same series of steps, except that c_{m} is divided instead of multiplied by u_{m}.
The FFTLog code is built on top of the NCAR suite of FFT routines (Swarztrauber 1979), and a modified version of an implementation of the complex Gammafunction from the gamerf package by Ooura (1996).
FFTLog includes driver routines for the specific cases of the Fourier sine and cosine transforms.
 (28) 
In addition to reducing ringing, the condition (28)
means that equation (27) remains true also at
m = ±N/2,
so is true for all m.
In this case the inverse Hankel mode
v^{}_{n}(m,q),
equation (25), is equal to the forward Hankel mode
v^{+}_{n}(m,q)
with q of the opposite sign
 (29) 
In other words, if condition (28) is satisfied, then the inverse discrete Hankel transform equals the forward discrete Hankel transform with q ® q. This is like the continuous Hankel transform, equations (2), where the inverse transform equals the forward transform with q ® q.
The periodicity condition (28) on
u_{±N/2}
translates,
for real
m and q,
into a condition on k_{0} r_{0}
 (30) 
The lowringing condition (30) is a condition on the phasing of the discrete points r_{n} and k_{n} at which the discrete Hankel transform is specified. The condition is analogous to, albeit more complicated than, the condition on the usual FFT that discrete frequencies be phased so that their wavenumbers are integers, equation (7).
FFTLog can be set to use automatically the lowringing value of k_{0} r_{0} nearest to any input value of k_{0} r_{0}.
Note that the lowringing value of k_{0} r_{0} from condition (30) differs for different m, q, and L/N. For example, the sine transform (m = 1/2) and cosine transform (m = 1/2) have different lowringing values of k_{0} r_{0}.
How else does the choice of k_{0} r_{0} affect the Hankel transform? Increasing the value of ln(k_{0} r_{0}) by one notch L/N cyclically shifts the discrete Hankel transform ã_{n}, equation (21), by one notch to the left, ã_{n} ® ã_{n1}. In other words, changing ln(k_{0} r_{0}) by an integral number of notches shifts the origin of the transform, but leaves the transform otherwise unchanged, as might have been expected.
In practice, since in most cases one is probably using the discrete Hankel transform as an approximation to the continuous transform, one would probably want to use k_{0} r_{0} » 1 (or 2, or p, according to taste).
The discrete Hankel modes
v_{n}(m,0) = v^{+}_{n}(m,0) = v^{}_{n}(m,0)
in the lowringing unbiased (q = 0) case
are periodic, orthonormal, and selfsimilar,
equation (26),
 (31) 
Like any orthogonal transformation, the lowringing unbiased (q = 0) Hankel transform commutes with the operations of matrix multiplication, inversion, and diagonalization (for nonlowringing or biased Hankel transforms, q ¹ 0, the operations do not commute). That is, the Hankel transform of the product of two matrices is equal to the product of their Hankel transforms, and so on.
All else being equal (which it may not be), given a choice between applying an unbiased (q = 0) or biased (q ¹ 0) Hankel transform, and between a lowringing k_{0} r_{0}, equation (30), or otherwise, one would be inclined to choose the lowringing unbiased transform, because of its orthogonality property.
Usually one is interested in the discrete Fourier or Hankel transform not for its own sake, but rather as an approximation to the continuous transform. The usual procedure would be to apply the discrete transform to a finite segment of the function a(r) to be transformed. For FFTLog, the procedure can be regarded as involving two steps: truncating the function to a finite logarithmic interval, which causes ringing of the transform; followed by periodic replication of the function in logarithmic space, which causes aliasing.
Figure 3:
Illustrates the ringing and aliasing that occurs
when the continuous Hankel transform of a function
is approximated by the discrete Hankel transform of a finite segment
of the function.
The Figure shows the result of taking an
unbiased (q = 0) Hankel transform,
equation (2),
of order
m = 1/2
(cosine transform) of a function that is Gaussian in the log
The top panels of Figure 3 show on the left the original function a(r), and on the right its Hankel transform ã(k). The middle panels of Figure 3 show on the left the truncated function a(r), and on the right its Hankel transform ã(k). Truncation of a(r) leads to ringing of its transform ã(k) at large wavenumbers k. The oscillations at large k are actually uniformly spaced in k, but appear bunched up because of the logarithmic plotting. The bottom panels of Figure 3 show on the left the truncated, periodically replicated function a(r), and on the right its corresponding periodically replicated Hankel transform ã(k), which is aliased. Vertical blue lines demarcate periodic intervals. Periodic replication means taking a sum of copies shifted by integral periods. From the definition (2) of the continuous Hankel transform, it can be seen that periodically replicating a function a(r) in logarithmic space lnr and then taking its continuous Hankel transform is equivalent to Hankel transforming the function a(r) and then periodically replicating the transform ã(k) in lnk. But truncating a function does not truncate its transform. So whereas a truncated, periodically replicated function a(r) contains contributions from only one period at each point r, the periodically replicated transform contains overlapping contributions from many periods at each point k. This is aliasing. The aliasing is visible in the bottom right panel of Figure 3 as an enhancement of the periodically replicated transform ã(k) on the high k side of the periodic interval. 
Ringing and aliasing can be reduced by taking suitable precautions.
The ringing that results from taking the discrete transform of a finite segment of a function can be reduced by arranging that the function folds smoothly from large to small scales. It may help to bias the function with a power law before transforming it, as in the second algorithm in §7. It may also help to use a lowringing value of k_{0} r_{0}, §8.
Aliasing can be reduced by enlarging the periodic interval. Aliasing can be eliminated (to machine precision) if the interval can be enlarged to the point where the transform ã(k) goes sensibly to zero at the boundaries of the period. Note that it is not sufficient to enlarge the interval to the point where a(r) is sensibly zero at the period boundaries: what is important is that the transform ã(k) goes to zero at the boundaries.
Last updated 9 Dec 2005.