Andrew Hamilton's Homepage
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 ISO-8859-1 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 Gamma-function 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 scale-invariant 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 3-dimensional Fourier transform
of the cosmological correlation function
a function of spatial separation r.
With the conventional normalization used by cosmologists,
Figure 1: Cosmological power spectrum P(k) of matter fluctuations predicted by the so-called LCDM model, a flat (W = 1) Universe dominated by a cosmological constant (WL = 0.7), and Cold Dark Matter (Wm = 0.3), including a sprinkling of baryons (Wb = 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.
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 103 h-1Mpc, and (blue, high resolution) with 768 points over the range r = 10-6 to 106 h-1Mpc. The lines are dashed where the correlation function is negative, at separations r > 119 h-1Mpc. The low and high resolution curves are almost indistinguishable except at r > 200 h-1Mpc, 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 106 h-1Mpc over which the transform was computed is much broader than the range plotted. The straight dashed line shows a power law (r/5 h-1Mpc)-1.8 for reference.
The middle panel of Figure 2 shows the ratio xlow/xhigh of the low to high resolution correlation functions.
The bottom panel of Figure 2 shows the ratio xFFT/xFFTLog of the correlation function xFFT computed with a normal FFT (sine transform) with 1023 points over the range r = 0.125 to 128 h-1Mpc, to the high resolution correlation function xFFTLog computed with FFTLog. The FFT'd correlation function xFFT 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 (= Fourier-Bessel) transform pair
Fourier sine and cosine transforms can be regarded as special
cases of Hankel transforms with
m = ±1/2,
As first noted by Siegman (1977), if the product kr in the Hankel transform is written as elnk + 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/2-integral m, so includes the cases of Fourier sine and cosine transforms, and spherical Hankel transforms involving the spherical Bessel functions jl(x) º (p/2 x)1/2 Jl+1/2(x).
The sampling theorem (e.g. Press et al. 1986 §12.1)
given a function a(r)
satisfying equation (7),
the Fourier coefficients cm
can be expressed in terms of the values
an º a(rn)
of the function a(r)
at the N discrete points
rn = n R / N
n = 0, ±1, ..., ±[N/2].
For even N, the periodicity of a(r) ensures that
a-N/2 = aN/2.
Specifically, the sampling theorem asserts that
the Fourier coefficients in the expansion (7) satisfy
Equations (9) and (10) constitute a discrete Fourier transform pair relating two periodic, linearly spaced sequences an and cm of length N. The standard FFT evaluates the discrete Fourier transform exactly (that is, to machine precision).
The continuous Hankel transform ã(k),
of a function a(r) of the form (12) is
The sampling theorem requires that
u-N/2 = uN/2
for even N,
which is not necessarily satisfied by equation (18).
at the discrete points
kn = k0 enL/N
considered by the sampling theorem, the contributions at
m = ±N/2
to the sum on the right hand side of equation (17)
(-)n cN/2 (uN/2 + uN/2*)/2 = (-)n cN/2 Re uN/2 .
Thus the equality (17) remains true
at the discrete points kn
are replaced by their real parts,
equations (13), (14),
(20) and (21)
yields the discrete Hankel transform pair
FFTLog evaluates the forward and inverse discrete Hankel transforms given by equations (22), (23), exactly (to machine precision).
in equation (25)
is equal to
equations (16) and (18),
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  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 um(m,-q) for m = ±N/2. To the extent that the highest frequency Fourier coefficient c±N/2 of a sequence an 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 low-ringing value of k0 r0, 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+1-q in the inverse case) is zero or an even negative integer, because u0(m,q) = Um(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 u0 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.
FFT cm um 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 an to obtain the Fourier coefficients cm, equation (13);
FFT cm um back to obtain the discrete Hankel transform ãn, equation (21);
Unbias ãn with a power law to obtain Ãn = ãnkn-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 An that is taken to be periodic, whereas in the algorithm above it is not An but rather an that is periodic.
The inverse discrete Hankel transform is accomplished by the same series of steps, except that cm is divided instead of multiplied by um.
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 Gamma-function from the gamerf package by Ooura (1996).
FFTLog includes driver routines for the specific cases of the Fourier sine and cosine transforms.
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
equation (25), is equal to the forward Hankel mode
with q of the opposite sign
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
m and q,
into a condition on k0 r0
The low-ringing condition (30) is a condition on the phasing of the discrete points rn and kn 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 low-ringing value of k0 r0 nearest to any input value of k0 r0.
Note that the low-ringing value of k0 r0 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 low-ringing values of k0 r0.
How else does the choice of k0 r0 affect the Hankel transform? Increasing the value of ln(k0 r0) by one notch L/N cyclically shifts the discrete Hankel transform ãn, equation (21), by one notch to the left, ãn ® ãn-1. In other words, changing ln(k0 r0) 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 k0 r0 » 1 (or 2, or p, according to taste).
The discrete Hankel modes
vn(m,0) = v+n(m,0) = v-n(m,0)
in the low-ringing unbiased (q = 0) case
are periodic, orthonormal, and self-similar,
Like any orthogonal transformation, the low-ringing unbiased (q = 0) Hankel transform commutes with the operations of matrix multiplication, inversion, and diagonalization (for non-low-ringing 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 low-ringing k0 r0, equation (30), or otherwise, one would be inclined to choose the low-ringing 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.
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,
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 low-ringing value of k0 r0, §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.