C06 Chapter Introduction : NAG Library, Mark 24

This chapter is concerned with the following tasks.

(a)	Calculating the discrete Fourier transform of a sequence of real or complex data values.
(b)	Calculating the discrete convolution or the discrete correlation of two sequences of real or complex data values using discrete Fourier transforms.
(c)	Calculating the inverse Laplace transform of a user-supplied function.
(d)	Direct summation of orthogonal series.
(e)	Acceleration of convergence of a sequence of real values.

Most of the routines in this chapter calculate the finite discrete Fourier transform (DFT) of a sequence of

n

complex numbers

z_{j}

, for

j = 0, 1, \dots, n - 1

. The direct transform is defined by

{\hat{z}}_{k} = \frac{1}{\sqrt{n}} \sum_{j = 0}^{n - 1} z_{j} \exp (- i \frac{2 π j k}{n})

(1)

for

k = 0, 1, \dots, n - 1

. Note that equation (1) makes sense for all integral

k

and with this extension

{\hat{z}}_{k}

is periodic with period

n

, i.e.,

{\hat{z}}_{k} = {\hat{z}}_{k \pm n}

, and in particular

{\hat{z}}_{- k} = {\hat{z}}_{n - k}

. Note also that the scale-factor of

\frac{1}{\sqrt{n}}

may be omitted in the definition of the DFT, and replaced by

\frac{1}{n}

in the definition of the inverse.

If we write

z_{j} = x_{j} + i y_{j}

and

{\hat{z}}_{k} = a_{k} + i b_{k}

, then the definition of

{\hat{z}}_{k}

may be written in terms of sines and cosines as

a_{k} = \frac{1}{\sqrt{n}} \sum_{j = 0}^{n - 1} (x_{j} \cos (\frac{2 π j k}{n}) + y_{j} \sin (\frac{2 π j k}{n}))

b_{k} = \frac{1}{\sqrt{n}} \sum_{j = 0}^{n - 1} (y_{j} \cos (\frac{2 π j k}{n}) - x_{j} \sin (\frac{2 π j k}{n})) .

The original data values

z_{j}

may conversely be recovered from the transform

{\hat{z}}_{k}

by an inverse discrete Fourier transform:

z_{j} = \frac{1}{\sqrt{n}} \sum_{k = 0}^{n - 1} {\hat{z}}_{k} \exp (+ i \frac{2 π j k}{n})

(2)

for

j = 0, 1, \dots, n - 1

. If we take the complex conjugate of (2), we find that the sequence

{\bar{z}}_{j}

is the DFT of the sequence

{\bar{\hat{z}}}_{k}

. Hence the inverse DFT of the sequence

{\hat{z}}_{k}

may be obtained by taking the complex conjugates of the

{\hat{z}}_{k}

; performing a DFT, and taking the complex conjugates of the result. (Note that the terms forward transform and backward transform are also used to mean the direct and inverse transforms respectively.)

The definition (1) of a one-dimensional transform can easily be extended to multidimensional transforms. For example, in two dimensions we have

{\hat{z}}_{k_{1} k_{2}} = \frac{1}{\sqrt{n_{1} n_{2}}} \sum_{j_{1} = 0}^{n_{1} - 1} \sum_{j_{2} = 0}^{n_{2} - 1} z_{j_{1} j_{2}} \exp (- i \frac{2 π j_{1} k_{1}}{n_{1}}) \exp (- i \frac{2 π j_{2} k_{2}}{n_{2}}) .

(3)

Note: definitions of the discrete Fourier transform vary. Sometimes (2) is used as the definition of the DFT, and (1) as the definition of the inverse.

If the original sequence is purely real valued, i.e.,

z_{j} = x_{j}

, then

{\hat{z}}_{k} = a_{k} + i b_{k} = \frac{1}{\sqrt{n}} \sum_{j = 0}^{n - 1} x_{j} \exp (- i \frac{2 π j k}{n})

and

{\hat{z}}_{n - k}

is the complex conjugate of

{\hat{z}}_{k}

. Thus the DFT of a real sequence is a particular type of complex sequence, called a Hermitian sequence, or half-complex or conjugate symmetric, with the properties

a_{n - k} = a_{k} b_{n - k} = {- b}_{k} b_{0} = 0

and, if

n

is even,

b_{n / 2} = 0

.

Thus a Hermitian sequence of

n

complex data values can be represented by only

n

, rather than

2 n

, independent real values. This can obviously lead to economies in storage, with two schemes being used in this chapter. In the first scheme, which will be referred to as the real storage format for Hermitian sequences, the real parts

a_{k}

for

0 \leq k \leq n / 2

are stored in normal order in the first

n / 2 + 1

locations of an array X of length

n

; the corresponding nonzero imaginary parts are stored in reverse order in the remaining locations of X. To clarify, if X is declared with bounds

(0 : n - 1)

in your calling subroutine, the following two tables illustrate the storage of the real and imaginary parts of

{\hat{z}}_{k}

for the two cases:

n

even and

n

odd.

If

n

is even then the sequence has two purely real elements and is stored as follows:

Index of X	0	1	2	$\dots$	$n / 2$	$\dots$	$n - 2$	$n - 1$
Sequence	$a_{0}$	$a_{1} + {i b}_{1}$	$a_{2} + {i b}_{2}$	$\dots$	$a_{n / 2}$	$\dots$	$a_{2} - {i b}_{2}$	$a_{1} - {i b}_{1}$
Stored values	$a_{0}$	$a_{1}$	$a_{2}$	$\dots$	$a_{n / 2}$	$\dots$	$b_{2}$	$b_{1}$

\begin{matrix} X (k) = a_{k}, & for ​ k = 0, 1, \dots, n / 2, and \\ X (n - k) = b_{k}, & for ​ k = 1, 2, \dots, n / 2 - 1 . \end{matrix}

If

n

is odd then the sequence has one purely real element and, letting

n = 2 s + 1

, is stored as follows:

Index of X	0	1	2	$\dots$	$s$	$s + 1$	$\dots$	$n - 2$	$n - 1$
Sequence	$a_{0}$	$a_{1} + {i b}_{1}$	$a_{2} + {i b}_{2}$	$\dots$	$a_{s} + i b_{s}$	$a_{s} - i b_{s}$	$\dots$	$a_{2} - {i b}_{2}$	$a_{1} - {i b}_{1}$
Stored values	$a_{0}$	$a_{1}$	$a_{2}$	$\dots$	$a_{s}$	$b_{s}$	$\dots$	$b_{2}$	$b_{1}$

\begin{matrix} X (k) = a_{k}, & for ​ k = 0, 1, \dots, s, and \\ X (n - k) = b_{k}, & for ​ k = 1, 2, \dots, s . \end{matrix}

The second storage scheme, referred to in this chapter as the complex storage format for Hermitian sequences, stores the real and imaginary parts

a_{k}, b_{k}

, for

0 \leq k \leq n / 2

, in consecutive locations of an array X of length

n + 2

. If X is declared with bounds

(0 : n + 1)

in your calling subroutine, the following two tables illustrate the storage of the real and imaginary parts of

{\hat{z}}_{k}

for the two cases:

n

even and

n

odd.

If

n

is even then the sequence has two purely real elements and is stored as follows:

Index of X	0	1	2	3	$\dots$	$n - 2$	$n - 1$	$n$	$n + 1$
Stored values	$a_{0}$	$b_{0} = 0$	$a_{1}$	$b_{1}$	$\dots$	$a_{n / 2 - 1}$	$b_{n / 2 - 1}$	$a_{n / 2}$	$b_{n / 2} = 0$

\begin{matrix} X (2 \times k) = a_{k}, & for ​ k = 0, 1, \dots, n / 2, and \\ X (2 \times k + 1) = b_{k}, & for ​ k = 0, 1, \dots, n / 2 . \end{matrix}

If

n

is odd then the sequence has one purely real element and, letting

n = 2 s + 1

, is stored as follows:

Index of X	0	1	2	3	$\dots$	$n - 2$	$n - 1$	$n$	$n + 1$
Stored values	$a_{0}$	$b_{0} = 0$	$a_{1}$	$b_{1}$	$\dots$	$b_{s - 1}$	$a_{s}$	$b_{s}$	$0$

\begin{matrix} X (2 \times k) = a_{k}, & for ​ k = 0, 1, \dots, s, and \\ X (2 \times k + 1) = b_{k}, & for ​ k = 0, 1, \dots, s . \end{matrix}

Also, given a Hermitian sequence, the inverse (or backward) discrete transform produces a real sequence. That is,

x_{j} = \frac{1}{\sqrt{n}} (a_{0} + 2 \sum_{k = 1}^{n / 2 - 1} (a_{k} \cos (\frac{2 π j k}{n}) - b_{k} \sin (\frac{2 π j k}{n})) + a_{n / 2})

where

a_{n / 2} = 0

if

n

is odd.

For real data that is two-dimensional or higher, the symmetry in the transform persists for the leading dimension only. So, using the notation of equation (3) for the complex two-dimensional discrete transform, we have that

{\hat{z}}_{k_{1} k_{2}}

is the complex conjugate of

{\hat{z}}_{(n_{1} - k_{1}) k_{2}}

. It is more convenient for transformed data of two or more dimensions to be stored as a complex sequence of length

(n_{1} / 2 + 1) \times n_{2} \times \dots \times n_{d}

where

d

is the number of dimensions. The inverse discrete Fourier transform operating on such a complex sequence (Hermitian in the leading dimension) returns a real array of full dimension (

n_{1} \times n_{2} \times \dots \times n_{d}

).

In many applications the sequence

x_{j}

will not only be real, but may also possess additional symmetries which we may exploit to reduce further the computing time and storage requirements. For example, if the sequence

x_{j}

is odd,

(x_{j} = {- x}_{n - j})

, then the discrete Fourier transform of

x_{j}

contains only sine terms. Rather than compute the transform of an odd sequence, we define the sine transform of a real sequence by

{\hat{x}}_{k} = \sqrt{\frac{2}{n}} \sum_{j = 1}^{n - 1} x_{j} \sin (\frac{π j k}{n}),

which could have been computed using the Fourier transform of a real odd sequence of length

2 n

. In this case the

x_{j}

are arbitrary, and the symmetry only becomes apparent when the sequence is extended. Similarly we define the cosine transform of a real sequence by

{\hat{x}}_{k} = \sqrt{\frac{2}{n}} (\frac{1}{2} x_{0} + \sum_{j = 1}^{n - 1} x_{j} \cos (\frac{π j k}{n}) + \frac{1}{2} {(- 1)}^{k} x_{n})

which could have been computed using the Fourier transform of a real even sequence of length

2 n

.

In addition to these ‘half-wave’ symmetries described above, sequences arise in practice with ‘quarter-wave’ symmetries. We define the quarter-wave sine transform by

{\hat{x}}_{k} = \frac{1}{\sqrt{n}} (\sum_{j = 1}^{n - 1} x_{j} \sin (\frac{π j (2 k - 1)}{2 n}) + \frac{1}{2} {(- 1)}^{k - 1} x_{n})

which could have been computed using the Fourier transform of a real sequence of length

4 n

of the form

(0, x_{1}, \dots, x_{n}, x_{n - 1}, \dots, x_{1}, 0, {- x}_{1}, \dots, {- x}_{n}, {- x}_{n - 1}, \dots, {- x}_{1}) .

Similarly we may define the quarter-wave cosine transform by

{\hat{x}}_{k} = \frac{1}{\sqrt{n}} (\frac{1}{2} x_{0} + \sum_{j = 1}^{n - 1} x_{j} \cos (\frac{π j (2 k - 1)}{2 n}))

which could have been computed using the Fourier transform of a real sequence of length

4 n

of the form

(x_{0}, x_{1}, \dots, x_{n - 1}, 0, {- x}_{n - 1}, \dots, {- x}_{0}, {- x}_{1}, \dots, {- x}_{n - 1}, 0, x_{n - 1}, \dots, x_{1}) .

The usual application of the discrete Fourier transform is that of obtaining an approximation of the Fourier integral transform

F (s) = \int_{- \infty}^{\infty} f (t) \exp (- i 2 π s t) d t

when

f (t)

is negligible outside some region

(0, c)

. Dividing the region into

n

equal intervals we have

F (s) ≅ \frac{c}{n} \sum_{j = 0}^{n - 1} f_{j} \exp (\frac{- i 2 π s j c}{n})

and so

F_{k} ≅ \frac{c}{n} \sum_{j = 0}^{n - 1} f_{j} \exp (\frac{- i 2 π j k}{n})

for

k = 0, 1, \dots, n - 1

, where

f_{j} = f (j c / n)

and

F_{k} = F (k / c)

.

Hence the discrete Fourier transform gives an approximation to the Fourier integral transform in the region

s = 0

to

s = n / c

.

If the function

f (t)

is defined over some more general interval

(a, b)

, then the integral transform can still be approximated by the discrete transform provided a shift is applied to move the point

a

to the origin.

One of the most important applications of the discrete Fourier transform is to the computation of the discrete convolution or correlation of two vectors

x

and

y

defined (as in Brigham (1974)) by

convolution: $z_{k} = \sum_{j = 0}^{n - 1} x_{j} y_{k - j}$
correlation: $w_{k} = \sum_{j = 0}^{n - 1} {\bar{x}}_{j} y_{k + j}$

(Here

x

and

y

are assumed to be periodic with period

n

.)

Under certain circumstances (see Brigham (1974)) these can be used as approximations to the convolution or correlation integrals defined by

z (s) = \int_{- \infty}^{\infty} x (t) y (s - t) d t

and

w (s) = \int_{- \infty}^{\infty} \bar{x} (t) y (s + t) d t, - \infty < s < \infty .

For more general advice on the use of Fourier transforms, see Hamming (1962); more detailed information on the fast Fourier transform algorithm can be found in Gentleman and Sande (1966) and Brigham (1974).

A further application of the fast Fourier transform, and in particular of the Fourier transforms of symmetric sequences, is in the solution of elliptic PDEs. If an equation is discretized using finite differences, then it is possible to reduce the problem of solving the resulting large system of linear equations to that of solving a number of tridiagonal systems of linear equations. This is accomplished by uncoupling the equations using Fourier transforms, where the nature of the boundary conditions determines the choice of transforms – see Section 3.3. Full details of the Fourier method for the solution of PDEs may be found in Swarztrauber (1977) and Swarztrauber (1984).

Let

f (t)

be a real function of

t

, with

f (t) = 0

for

t < 0

, and be piecewise continuous and of exponential order

α

, i.e.,

|f (t)| \leq M e^{α t}

for large

t

, where

α

is the minimal such exponent.

The Laplace transform of

f (t)

is given by

F (s) = \int_{0}^{\infty} e^{- s t} f (t) d t, t > 0

where

F (s)

is defined for

Re (s) > α

.

The inverse transform is defined by the Bromwich integral

f (t) = \frac{1}{2 π i} \int_{a - i \infty}^{a + i \infty} e^{s t} F (s) d s, t > 0 .

The integration is performed along the line

s = a

in the complex plane, where

a > α

. This is equivalent to saying that the line

s = a

lies to the right of all singularities of

F (s)

. For this reason, the value of

α

is crucial to the correct evaluation of the inverse. It is not essential to know

α

exactly, but an upper bound must be known.

The problem of determining an inverse Laplace transform may be classified according to whether (a)

F (s)

is known for real values only, or (b)

F (s)

is known in functional form and can therefore be calculated for complex values of

s

. Problem (a) is very ill-defined and no routines are provided. Two methods are provided for problem (b).

For any series of functions

ϕ_{i}

which satisfy a recurrence

ϕ_{r + 1} (x) + α_{r} (x) ϕ_{r} (x) + β_{r} (x) ϕ_{r - 1} (x) = 0

the sum

\sum_{r = 0}^{n} a_{r} ϕ_{r} (x)

is given by

\sum_{r = 0}^{n} a_{r} ϕ_{r} (x) = b_{0} (x) ϕ_{0} (x) + b_{1} (x) (ϕ_{1} (x) + α_{0} (x) ϕ_{0} (x))

where

b_{r} (x) + α_{r} (x) b_{r + 1} (x) + β_{r + 1} (x) b_{r + 2} (x) = a_{r} b_{n + 1} (x) = b_{n + 2} (x) = 0 .

This may be used to compute the sum of the series. For further reading, see Hamming (1962).

This device has applications in a large number of fields, such as summation of series, calculation of integrals with oscillatory integrands (including, for example, Hankel transforms), and root-finding. The mathematical description is as follows. Given a sequence of values

\{s_{n}\}

, for

n = m, \dots, m + 2 l

, then, except in certain singular cases, parameters,

a

,

b_{i}

,

c_{i}

may be determined such that

s_{n} = a + \sum_{i = 1}^{l} b_{i} c_{i}^{n} .

If the sequence

\{s_{n}\}

converges, then

a

may be taken as an estimate of the limit. The method will also find a pseudo-limit of certain divergent sequences – see Shanks (1955) for details.

To use the method to sum a series, the terms

s_{n}

of the sequence should be the partial sums of the series, e.g.,

s_{n} = \sum_{k = 1}^{n} t_{k}

, where

t_{k}

is the

k

th term of the series. The algorithm can also be used to some advantage to evaluate integrals with oscillatory integrands; one approach is to write the integral (in this case over a semi-infinite interval) as

\int_{0}^{\infty} f (x) d x = \int_{0}^{a_{1}} f (x) d x + \int_{a_{1}}^{a_{2}} f (x) d x + \int_{a_{2}}^{a_{3}} f (x) d x + \dots

and to consider the sequence of values

s_{1} = \int_{0}^{a_{1}} f (x) d x, s_{2} = \int_{0}^{a_{2}} f (x) d x = s_{1} + \int_{a_{1}}^{a_{2}} f (x) d x, etc.,

where the integrals are evaluated using standard quadrature methods. In choosing the values of the

a_{k}

, it is worth bearing in mind that C06BAF converges much more rapidly for sequences whose values oscillate about a limit. The

a_{k}

should thus be chosen to be (close to) the zeros of

f (x)

, so that successive contributions to the integral are of opposite sign. As an example, consider the case where

f (x) = M (x) \sin x

and

M (x) > 0

: convergence will be much improved if

a_{k} = k π

rather than

a_{k} = 2 k π

.

The choice of routine is determined first of all by whether the data values constitute a real, Hermitian or general complex sequence. It is wasteful of time and storage to use an inappropriate routine.

The choice is next determined by how you prefer to store complex data. Real storage format routines store general complex sequences and Hermitian sequences in real arrays. In the case of general complex sequences, the real and imaginary parts are stored separately in two real arrays. Since Hermitian sequences contain some symmetries, these can be stored in a compact form in a single real array. Alternatively, complex storage format routines store the corresponding sequence as a complex array for general sequences, and with real and imaginary parts in contiguous locations of a real array for Hermitian sequences.

Two groups, each of three routines, are provided in real storage format and three groups of two routines are provided in complex storage format.

Group 1 routines compute a single transform of length

n

, requiring one additional workspace array. Some of these (C06FAF, C06FBF and C06FCF) impose some restrictions on the value of

n

, namely that no prime factor of

n

may exceed

19

and the total number of prime factors (including repetitions) may not exceed

20

(though the latter restriction only becomes relevant when

n > 10^{6}

).

Group 2 and Group 3 routines are all designed to perform several transforms in a single call, all with the same value of

n

. They are likely to be much faster than the Group 1 routines on modern acrchitectures. They do however require more working storage. It is therefore recommended that, for real storage format, Group 2 routines be used in preference to Group 1 routines, even when only one transform is to be performed. Group 2 and Group 3 routines differ in the way sequences are stored: Group 2 routines store sequences as rows of a two-dimensional array while Group 3 routines store sequences as columns of a two-dimensional array. Group 2 and Group 3 routines impose no practical restrictions on the value of

n

; however, the fast Fourier transform algorithm ceases to be ‘fast’ if applied to values of

n

which cannot be expressed as a product of small prime factors. All the above routines are particularly efficient if the only prime factors of

n

are

2

,

3

or

5

.

If extensive use is to be made of these routines and you are concerned about efficiency, you are advised to conduct your own timing tests.

To compute inverse (backward) discrete Fourier transforms the real storage format routines should be used in conjunction with the complex conjugate of a Hermitian or general sequence of complex data values. In the case of complex storage format routines, there is a direction parameter which determines the direction of the transform; a call to such a routine in the forward direction followed by a call in the backward direction reproduces the original data.

Four routines are provided for computing fast Fourier transforms (FFTs) of real symmetric sequences. C06RAF computes multiple Fourier sine transforms, C06RBF computes multiple Fourier cosine transforms, C06RCF computes multiple quarter-wave Fourier sine transforms, and C06RDF computes multiple quarter-wave Fourier cosine transforms.

As described in Section 2.1, Fourier transforms may be used in the solution of elliptic PDEs.

C06RAF may be used to solve equations where the solution is specified along the boundary.

C06RBF may be used to solve equations where the derivative of the solution is specified along the boundary.

C06RCF may be used to solve equations where the solution is specified on the lower boundary, and the derivative of the solution is specified on the upper boundary.

C06RDF may be used to solve equations where the derivative of the solution is specified on the lower boundary, and the solution is specified on the upper boundary.

For equations with periodic boundary conditions the full-range Fourier transforms computed by C06FPF and C06FQF are appropriate.

The following routines compute multidimensional discrete Fourier transforms of complex data:

	Real storage	Complex storage
2 dimensions		C06PUF
3 dimensions	C06FXF	C06PXF
any number of dimensions	C06FJF	C06PJF

The real storage format routines store sequences of complex data in two real arrays containing the real and imaginary parts of the sequence respectively. The complex storage format routines store the sequences in complex arrays.

C06PUF and C06FXF/C06PXF should be used in preference to C06FJF/C06PJF for two- and three-dimensional transforms, as they are easier to use and are likely to be more efficient.

The transform of multidimensional real data is stored as a complex sequence that is Hermitian in its leading dimension. The inverse transform takes such a complex sequence and computes the real transformed sequence. Consequently, separate routines are provided for performing forward and inverse transforms.

C06PVF performs the forward two-dimensionsal transform while C06PWF performs the inverse of this transform.

C06PYF performs the forward three-dimensional transform while C06PZF performs the inverse of this transform.

The complex sequences computed by C06PVF and C06PYF contain roughly half of the Fourier coefficients; the remainder can be reconstructed by conjugation of those computed. For example, the Fourier coefficients of the two-dimensional transform

{\hat{z}}_{(n_{1} - k_{1}) k_{2}}

are the complex conjugate of

{\hat{z}}_{k_{1} k_{2}}

for

k_{1} = 0, 1, \dots, n_{1} / 2

, and

k_{2} = 0, 1, \dots, n_{2} - 1

.

C06FKF computes either the discrete convolution or the discrete correlation of two real vectors.

C06PKF computes either the discrete convolution or the discrete correlation of two complex vectors.

Two methods are provided: Weeks' method (C06LBF) and Crump's method (C06LAF). Both require the function

F (s)

to be evaluated for complex values of

s

. If in doubt which method to use, try Weeks' method (C06LBF) first; when it is suitable, it is usually much faster.

Typically the inversion of a Laplace transform becomes harder as

t

increases so that all numerical methods tend to have a limit on the range of

t

for which the inverse

f (t)

can be computed. C06LAF is useful for small and moderate values of

t

.

It is often convenient or necessary to scale a problem so that

α

is close to

0

. For this purpose it is useful to remember that the inverse of

F (s + k)

is

\exp (- k t) f (t)

. The method used by C06LAF is not so satisfactory when

f (t)

is close to zero, in which case a term may be added to

F (s)

, e.g.,

k / s + F (s)

has the inverse

k + f (t)

.

Singularities in the inverse function

f (t)

generally cause numerical methods to perform less well. The positions of singularities can often be identified by examination of

F (s)

. If

F (s)

contains a term of the form

\exp (- k s) / s

then a finite discontinuity may be expected in the inverse at

t = k

. C06LAF, for example, is capable of estimating a discontinuous inverse but, as the approximation used is continuous, Gibbs' phenomena (overshoots around the discontinuity) result. If possible, such singularities of

F (s)

should be removed before computing the inverse.

The only routine available is C06DCF, which sums a finite Chebyshev series

\sum_{j = 0}^{n} c_{j} T_{j} (x), \sum_{j = 0}^{n} c_{j} T_{2 j} (x) or \sum_{j = 0}^{n} c_{j} T_{2 j + 1} (x)

depending on the choice of a parameter.

The only routine available is C06BAF.

None.

The following lists all those routines that have been withdrawn since Mark 17 of the Library or are scheduled for withdrawal at one of the next two marks.

Brigham E O (1974) The Fast Fourier Transform Prentice–Hall

Davies S B and Martin B (1979) Numerical inversion of the Laplace transform: A survey and comparison of methods J. Comput. Phys. 33 1–32

Fox L and Parker I B (1968) Chebyshev Polynomials in Numerical Analysis Oxford University Press

Gentleman W S and Sande G (1966) Fast Fourier transforms for fun and profit Proc. Joint Computer Conference, AFIPS 29 563–578

Hamming R W (1962) Numerical Methods for Scientists and Engineers McGraw–Hill

Shanks D (1955) Nonlinear transformations of divergent and slowly convergent sequences J. Math. Phys. 34 1–42

Swarztrauber P N (1977) The methods of cyclic reduction, Fourier analysis and the FACR algorithm for the discrete solution of Poisson's equation on a rectangle SIAM Rev. 19(3) 490–501

Swarztrauber P N (1984) Fast Poisson solvers Studies in Numerical Analysis (ed G H Golub) Mathematical Association of America

Swarztrauber P N (1986) Symmetric FFT's Math. Comput. 47(175) 323–346

Wynn P (1956) On a device for computing the

e_{m} (S_{n})

transformation Math. Tables Aids Comput. 10 91–96

	Group 1	Group 2	Group 3
Real storage format
Real sequences	C06FAF	C06FPF
Hermitian sequences	C06FBF	C06FQF
General complex sequences	C06FCF

Complex storage format/
complex data type
Real/Hermitian sequences	C06PAF	C06PPF	C06PQF
General complex sequences	C06PCF	C06PRF	C06PSF

Is the data one-dimensional?	_ yes	Multiple vectors?	_ yes	Stored as rows?	_ yes	C06PRF
\|		\|		no \|
\|		\|		Stored as columns?	_ yes	C06PSF
\|		no \|
\|		C06PCF
no \|
Is the data two-dimensional?	_ yes	C06PUF
no \|
Is the data three-dimensional?	_ yes	C06PXF
no \|
Transform on one dimension only?	_ yes	C06PFF
no \|
Transform on all dimensions?	_ yes	C06PJF

Quarter-wave sine (inverse) transform	_ yes	C06RCF
no \|
Quarter-wave cosine (inverse) transform	_ yes	C06RDF
no \|
Sine (inverse) transform	_ yes	C06RAF
no \|
Cosine (inverse) transform	_ yes	C06RBF
no \|
Is the data three-dimensional?	_ yes	Forward transform on real data	_ yes	C06PYF
\|		no \|
\|		Inverse transform on Hermitian data	_ yes	C06PZF
no \|
Is the data two-dimensional?	_ yes	Forward transform on real data	_ yes	C06PVF
\|		no \|
\|		Inverse transform on Hermitian data	_ yes	C06PWF
no \|
Is the data multi one-dimensional?	_ yes	Sequences stored by row	_ yes	C06PPF
\|		no \|
\|		Sequences stored by column	_ yes	C06PQF
no \|
Is the data one-dimensional?	_ yes	C06PAF
no \|
Transform in one-dimension only	_ yes	C06FFF
no \|
C06FJF

Withdrawn Routine	Mark of Withdrawal	Replacement Routine(s)
C06DBF	25	C06DCF
C06EAF	26	C06PAF
C06EBF	26	C06PAF
C06ECF	26	C06PCF
C06EKF	26	C06FKF
C06FRF	26	C06PSF
C06FUF	26	C06PUF
C06GBF	26	No replacement required
C06GCF	26	No replacement required
C06GQF	26	No replacement required
C06GSF	26	No replacement required
C06HAF	26	C06RAF
C06HBF	26	C06RAF
C06HCF	26	C06RCF
C06HDF	26	C06RDF

NAG Library Chapter Introduction

C06 – Summation of Series

+− Contents

1 Scope of the Chapter

2 Background to the Problems

2.1 Discrete Fourier Transforms

2.1.1 Complex transforms

2.1.2 Real transforms

2.1.3 Real symmetric transforms

2.1.4 Fourier integral transforms

2.1.5 Convolutions and correlations

2.1.6 Applications to solving partial differential equations (PDEs)

2.2 Inverse Laplace Transforms

2.3 Direct Summation of Orthogonal Series

2.4 Acceleration of Convergence

3 Recommendations on Choice and Use of Available Routines

3.1 One-dimensional Fourier Transforms

3.2 Half- and Quarter-wave Transforms

3.3 Application to Elliptic Partial Differential Equations

3.4 Multidimensional Fourier Transforms

3.4.1 Complex data

3.4.2 Real data

3.5 Convolution and Correlation

3.6 Inverse Laplace Transforms

3.7 Direct Summation of Orthogonal Series

3.8 Acceleration of Convergence

4 Decision Trees

Tree 1: Fourier Transform of Discrete Complex Data

Tree 2: Fourier Transform of Real Data or Data in Complex Hermitian Form Resulting from the Transform of Real Data

5 Functionality Index

6 Auxiliary Routines Associated with Library Routine Parameters

7 Routines Withdrawn or Scheduled for Withdrawal

8 References

NAG Library Chapter IntroductionC06 – Summation of Series

+− Contents

1 Scope of the Chapter

2 Background to the Problems

2.1 Discrete Fourier Transforms

2.1.1 Complex transforms

2.1.2 Real transforms

2.1.3 Real symmetric transforms

2.1.4 Fourier integral transforms

2.1.5 Convolutions and correlations

2.1.6 Applications to solving partial differential equations (PDEs)

2.2 Inverse Laplace Transforms

2.3 Direct Summation of Orthogonal Series

2.4 Acceleration of Convergence

3 Recommendations on Choice and Use of Available Routines

3.1 One-dimensional Fourier Transforms

3.2 Half- and Quarter-wave Transforms

3.3 Application to Elliptic Partial Differential Equations

3.4 Multidimensional Fourier Transforms

3.4.1 Complex data

3.4.2 Real data

3.5 Convolution and Correlation

3.6 Inverse Laplace Transforms

3.7 Direct Summation of Orthogonal Series

3.8 Acceleration of Convergence

4 Decision Trees

Tree 1: Fourier Transform of Discrete Complex Data

Tree 2: Fourier Transform of Real Data or Data in Complex Hermitian Form Resulting from the Transform of Real Data

5 Functionality Index

6 Auxiliary Routines Associated with Library Routine Parameters

7 Routines Withdrawn or Scheduled for Withdrawal

8 References

NAG Library Chapter Introduction

C06 – Summation of Series