c06 Chapter Introduction : NAG Library, Mark 23

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)	Direct summation of orthogonal series.

Most of the functions 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}}) .

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 one scheme being used in this chapter. In this 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 [n]

in your calling function, 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}

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.

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

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

The choice of function 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 function.

Two groups, each of three functions, are provided in real storage format.

Group 1 functions each compute a single transform of length

n

, without requiring any extra working storage. The Group 1 functions 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 functions are all designed to perform several transforms in a single call, all with the same value of

n

. They are designed to be much faster than the Group 1 functions on vector-processing machines. They do however require more working storage. Even on scalar processors, they may be somewhat faster than repeated calls to Group 1 functions because of reduced overheads and because they pre-compute and store the required values of trigonometric functions. Group 2 functions 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 functions are particularly efficient if the only prime factors of

n

are

2

,

3

or

5

.

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

To compute inverse (backward) discrete Fourier transforms the functions should be used in conjunction with the utility functions nag_conjugate_hermitian (c06gbc), nag_conjugate_complex (c06gcc) and nag_multiple_conjugate_hermitian (c06gqc) which form the complex conjugate of a Hermitian or general sequence of complex data values.

Four functions are provided for computing fast Fourier transforms (FFTs) of real symmetric sequences. nag_fft_multiple_sine (c06hac) computes compute multiple Fourier sine transforms, nag_fft_multiple_cosine (c06hbc) computes multiple Fourier cosine transforms, nag_fft_multiple_qtr_sine (c06hcc) computes multiple quarter-wave Fourier sine transforms, and nag_fft_multiple_qtr_cosine (c06hdc) computes multiple quarter-wave Fourier cosine transforms.

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

nag_fft_multiple_sine (c06hac) may be used to solve equations where the solution is specified along the boundary.

nag_fft_multiple_cosine (c06hbc) may be used to solve equations where the derivative of the solution is specified along the boundary.

nag_fft_multiple_qtr_sine (c06hcc) 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.

nag_fft_multiple_qtr_cosine (c06hdc) 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 nag_fft_multiple_real (c06fpc) and nag_fft_multiple_hermitian (c06fqc) are appropriate.

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

	Real storage	Complex storage
2 dimensions	nag_fft_2d_complex (c06fuc)
3 dimensions		nag_fft_3d (c06pxc)
any number of dimensions		nag_fft_multid_full (c06pjc)

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

nag_fft_2d_complex (c06fuc) and nag_fft_3d (c06pxc) should be used in preference to nag_fft_multid_full (c06pjc) for two- and three-dimensional transforms, as they are easier to use and are likely to be more efficient, especially on vector processors.

nag_convolution_real (c06ekc) computes either the discrete convolution or the discrete correlation of two real vectors.

The only function available is nag_sum_cheby_series (c06dcc), 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 an argument.

None.

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
Real sequences	nag_fft_real (c06eac)	nag_fft_multiple_real (c06fpc)
Hermitian sequences	nag_fft_hermitian (c06ebc)	nag_fft_multiple_hermitian (c06fqc)
General complex sequences	nag_fft_complex (c06ecc)	nag_fft_multiple_complex (c06frc)

NAG Library Chapter Introduction

c06 – Fourier Transforms

+− 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 Direct Summation of Orthogonal Series

3 Recommendations on Choice and Use of Available Functions

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.5 Convolution and Correlation

3.6 Direct Summation of Orthogonal Series

4 Functionality Index

5 Functions Withdrawn or Scheduled for Withdrawal

6 References

NAG Library Chapter Introductionc06 – Fourier Transforms