Home Andrew Hamilton's Homepage

FFTLog
A code to take the fast Fourier or Hankel transform of a discrete periodic sequence of logarithmically spaced points.

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.

Contents

  1. Download
  2. What is FFTLog?
  3. Motivation and Example
  4. Introduction
  5. Normal discrete Fourier transform
  6. Discrete Hankel transform
  7. FFTLog algorithm
  8. Low-ringing condition on k0 r0
  9. Unitary Hankel transform
  10. Ringing and Aliasing
  11. Troubleshooting
  12. References, Links

1. Download

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 .

2. What is FFTLog?

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

3. Motivation and Example

FFTLog emerged from a problem in cosmology (Hamilton 2000). The problem required Fourier transforming a function that extended over many orders of magnitude, and was `smooth' in logarithmic space. Actually, it was necessary to transform whole matrices of such functions, so a fast transform method was desirable.

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 x(r), a function of spatial separation r. With the conventional normalization used by cosmologists,

P(k) = ó¥
õ0		 x(r) sin(kr)
kr
 4pr2 dr  ,    x(r) = ó¥
õ0		 P(k) sin(kr)
kr
 4pk2 dk
(2p)3
  .
(1)

Cosmological power spectrum 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 DMR icon COBE satellite.

Cosmological correlation function 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 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.

Both low and high resolution cases used an unbiased (q = 0) transform, §9, and a low-ringing value of k0 r0, §8 (actually the choice of k0 r0 made little difference here).

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.

4. Introduction

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

Consider the continuous Hankel (= Fourier-Bessel) transform pair

ã(k) = ó¥
õ0		 a(r) (k r)q Jm(k r) k dr  ,    a(r) = ó¥
õ0		 ã(k) (k r)-q Jm(k r) r dk  .
(2)
If the substitution
a(r) = A(r) r-q    and     ã(k) = Ã(k) kq
(3)
is made, then the Hankel transform pair (2) becomes equivalent to the transform pair
Ã(k) = ó¥
õ0		 A(r) Jm(k r) k dr  ,    A(r) = ó¥
õ0		 Ã(k) Jm(k r) r dk  .
(4)
Although the Hankel transform (2) with a power law bias (k r)±q is thus equivalent in the continuous case to the unbiased Hankel transform (4), the transforms are different when they are discretized and made periodic; for if a(r) is periodic, then A(r) = a(r) rq is not periodic. FFTLog evaluates discrete Hankel transforms (2) with arbitrary power law bias.

Fourier sine and cosine transforms can be regarded as special cases of Hankel transforms with m = ±1/2, since

J1/2(x) = (2/px)1/2 sin(x)  ,    J-1/2(x) = (2/px)1/2 cos(x) .
(5)

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:

This is the idea behind a number of Fast Hankel Transform (FHT) algorithms (Sieman 1977; Sheng & Siegman 1980; Candel 1981; Anderson 1982; Hansen 1985; Fanning 1996) including FFTLog (Talman 1978).

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

5. Normal discrete Fourier transform

First, recall the essential properties of the standard discrete Fourier transform of a periodic sequence of linearly spaced points. Suppose that a(r) is a continuous, in general complex-valued, function that is periodic with period R,
a(r + R) = a(r)  .
(6)
Without loss of generality, take the fundamental interval to be [-R/2, R/2] , centred at zero. Since a(r) is periodic, its continuous Fourier transform contains only discrete Fourier modes e2pim r/R with integral wavenumbers m. Suppose further that the function a(r) is `smooth' in the specific sense that it is some linear combination only of the N lowest frequency Fourier modes, m = 0, ±1, ..., ±[N/2], where [N/2] denotes the largest integer less than or equal to N/2,
a(r) = [N/2]
S´
m = -[N/2] 
 cm e2pimr/R
(7)
the outermost Fourier coefficients being equal, c-N/2 = cN/2, in the case of even N. The primed sum in equation (7) signifies a sum over integral m from -[N/2] to [N/2], with the proviso that for even N the outermost elements of the sum receive only half weight:
[N/2]
S´
n = -[N/2] 
 xn º [N/2]
S
n = -[N/2] 
 wn xn
(8)
with wn = 1 except that w-N/2 = wN/2 = 1/2 if N is even.

The sampling theorem (e.g. Press et al. 1986 §12.1) asserts that, 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 for 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

cm = 1
N		 [N/2]
S´
n = -[N/2] 
 an e-2pimn/N
(9)
the discrete points an themselves satisfying
an = [N/2]
S´
m = -[N/2] 
 cm e2pimn/N
(10)
in accordance with equation (7).

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

6. Discrete Hankel transform

Now suppose that the, in general complex-valued, function a(r), instead of being periodic in ordinary space r, is periodic in logarithmic space lnr, with logarithmic period L,
a(r eL) = a(r)  .
(11)
Take the fundamental interval to be [lnr0 - L/2, lnr0 + L/2] , centred at lnr0. As in §5, the periodicity of a(r) implies that its Fourier transform with respect to lnr contains only discrete Fourier modes e2pim ln(r/r0)/L with integral wavenumbers m. Suppose further, as in §5 eq. (7), that a(r) contains only the N lowest frequency Fourier modes
a(r) = [N/2]
S´
m = -[N/2] 
 cm e2pim ln(r/r0)/L
(12)
with c-N/2 = cN/2 for even N. The sampling theorem asserts that the Fourier coefficients cm are given by
cm = 1
N		 [N/2]
S´
n = -[N/2] 
 an e-2pimn/N
(13)
where an º a(rn) are the values of the function a(r) at the N discrete points rn = r0 enL/N for n = 0, ±1, ..., ±[N/2],
an = [N/2]
S´
m = -[N/2] 
 cm e2pimn/N  .
(14)

The continuous Hankel transform ã(k), equation (2), of a function a(r) of the form (12) is

ã(k) = [N/2]
S´
m = -[N/2] 
 cm ó¥
õ0		 e2pim ln(r/r0)/L (k r)q Jm(k r) k dr  .
(15)
The integrals on the right hand side of equation (15) can be done analytically, in terms of
Um(x) º ó¥
õ0		 tx Jm(t) dt = 2x G[(m+1+x)/2]
G[(m+1-x)/2]
(16)
where G(z) is the usual Gamma-function. Thus equation (15) reduces to
ã(k) = [N/2]
S´
m = -[N/2] 
 cm um e-2pim ln(k/k0)/L
(17)
where um is
um(m,q) º ( k0 r0 )-2pim/L  Um æ
ç
è		 q + 2pim
L
 ö
÷
ø		  .
(18)
Notice that um* = u-m, which ensures that ã(k) is real if a(r) is real. Equation (17) gives the (exact) continuous Hankel transform ã(k) of a function a(r) of the form (7). Like a(r), the Hankel transform ã(k) is periodic in logarithmic space lnk, with period L. The fundamental interval is [lnk0 - L/2, lnk0 + L/2] , centred at lnk0, which may be chosen arbitrarily (but see §8 below).

The sampling theorem requires that u-N/2 = uN/2 for even N, which is not necessarily satisfied by equation (18). However, 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) are (-)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 if u±N/2 are replaced by their real parts,

u±N/2 ® Re uN/2  .
(19)
With the replacement (19), the sampling theorem asserts that the coefficients cm um in the sum (17) are determined by the values ãn º ã(kn) of the Hankel transform at the N discrete points kn = k0 en L / N for n = 0, ±1, ..., ±[N/2]
cm um = 1
N		 [N/2]
S´
n = -[N/2] 
 ãn e2pimn/N
(20)

ãn = [N/2]
S´
m = -[N/2] 
 cm um e-2pimn/N  .
(21)

Putting together equations (13), (14), (20) and (21) yields the discrete Hankel transform pair

ãn = [N/2]
S´
m = -[N/2] 
 am v+m+n(m,q)
(22)

am = [N/2]
S´
n = -[N/2] 
 ãn v-m+n(m,q)
(23)
in which the forward discrete Hankel mode v+n(m,q) is the discrete Fourier transform of um(m,q) given by equations (18) and (19),
v+n(m,q) = 1
N		 [N/2]
S´
m = -[N/2] 
 um(m,q) e-2pimn/N
(24)
while the inverse discrete Hankel mode v-n(m,q) is the discrete Fourier transform of the reciprocal 1/u-m(m,q) ,
v-n(m,q) = 1
N		 [N/2]
S´
m = -[N/2] 
 1
u-m(m,q)
  e-2pimn/N .
(25)
The Hankel transform matrices v+m+n(m,q) and v-m+n(m,q) are mutually inverse
[N/2]
S´
l = -[N/2] 
 v+m+l(m,q) v-l+n(m,q) = dmn
(26)
where dmn denotes the Kronecker delta. The forward and inverse Hankel modes have the interesting property of being self-similar; that is, Hankel modes v+m+n(m,q) [or v-m+n(m,q)] with different indices m consist of the same periodic sequence v+n(m,q) [or v-n(m,q)] cyclically shifted by m notches.

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 um(m,-q), according to equations (16) and  (18),

1
u-m(m,q)
 = um(m,-q)        (m ¹ N/2)
(27)
except in the case m = ±N/2 for even N, when the replacement (19) generally invalidates equation (27). However, in the special case where u±N/2 are already real, then equation (19) leaves u±N/2 unchanged, and equation (27) remains valid also at m = ±N/2. This special case is of particular interest, and is discussed further in §8 below.

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

7. FFTLog algorithm

The FFTLog algorithm for taking the discrete Hankel transform, equation (22), of a sequence an of N logarithmically spaced points is:

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:

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.

8. Low-ringing condition on k0 r0

The central values lnr0 and lnk0 of the periodic intervals in lnr and lnk may be chosen arbitrarily. However, ringing of the discrete Hankel transform may be reduced, for either even or odd N, if the product k0 r0 is chosen in such a way that the boundary points of the sequence um, equation (18), are equal
u-N/2 = uN/2 .
(28)
Recall that the general procedure, for even N, was to replace u±N/2 by their real part, equation (19). The condition (28) requires that u±N/2 are already real. The condition (28) reduces ringing because it makes the periodic sequence um fold smoothly across the period boundary at m = ±N/2.

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

v-n(m,q) = v+n(m,-q) = 1
N		 [N/2]
S´
m = -[N/2] 
 um(m,-q) e-2pimn/N .
(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 k0 r0

ln(k0 r0 ) = L
N
 ì
í
î		 1
p
 Argé
ê
ë		 Um æ
ç
è		 q + piN
L
 ö
÷
ø		 ù
ú
û		 + integer ü
ý
þ
(30)
where Argz º Im lnz denotes the argument of a complex number, and integer is any integer. In other words, to reduce ringing, it may help to choose k0 r0 so as to satisfy the condition (30). This is not too much of a restriction, since L/N is the logarithmic spacing between points (= one notch), so the low-ringing condition (30) allows k0 r0 to be chosen to lie within half a notch [= L/(2 N)] of whatever number one chooses, for example within half a notch of k0 r0 = 1.

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

9. Unitary Hankel transform

The discrete Hankel transform with both low-ringing k0 r0 and no power law bias, q = 0, is of particular interest because it is unitary, like the Fourier transform. Indeed, being also real, the low-ringing unbiased Hankel transform is orthogonal, i.e. self-inverse, like the Fourier sine and cosine transforms. This is like the continuous unbiased (q = 0) Hankel transform, equations (2), which is self-inverse.

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, equation (26),

[N/2]
S´
l = -[N/2] 
 vm+l(m,0) vl+n(m,0) = dmn .
(31)

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.

10. Ringing and aliasing

FFTLog suffers from the same problems of ringing (response to sudden steps) and aliasing (periodic folding of frequencies) as the normal FFT.

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.

Ringing and 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

a(r) = exp[-(lnr)2/2]  .
(32)
Lines are red where values are negative.

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.

11. Troubleshooting

FFTLog does not work well with my function. What should I do?

12. References, Links

File translated from TEX by TTH, version 2.01.
On 13 Mar 1999, 21:17.

Last updated 9 Dec 2005.