NAG C Library Chapter Introduction
c06 – Fourier Transforms
1
Scope of the Chapter
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 fast Gauss transform approximation to the discrete Gauss transform. |
(d) |
Direct summation of orthogonal series. |
2
Background to the Problems
2.1
Discrete Fourier Transforms
2.1.1
Complex transforms
Most of the functions in this chapter calculate the finite
discrete Fourier transform (DFT) of a sequence of
complex numbers
, for
. The direct transform is defined by
for
. Note that equation
(1) makes sense for all integral
and with this extension
is periodic with period
, i.e.,
, and in particular
. Note also that the scale-factor of
may be omitted in the definition of the DFT, and replaced by
in the definition of the inverse.
If we write
and
, then the definition of
may be written in terms of sines and cosines as
The original data values
may conversely be recovered from the transform
by an
inverse discrete Fourier transform:
for
. If we take the complex conjugate of
(2), we find that the sequence
is the DFT of the sequence
. Hence the inverse DFT of the sequence
may be obtained by taking the complex conjugates of the
; 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
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.
2.1.2
Real transforms
If the original sequence is purely real valued, i.e.,
, then
and
is the complex conjugate of
. 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
and, if
is even,
.
Thus a Hermitian sequence of
complex data values can be represented by only
, rather than
, independent real values. This can obviously lead to economies in storage, with two schemes being used in this chapter. In
the first (deprecated) scheme, which will be referred to as the
real storage format for Hermitian sequences, the real parts
for
are stored in normal order in the first
locations of an array
x of length
; the corresponding nonzero imaginary parts are stored in reverse order in the remaining locations of
x. To clarify,
the following two tables illustrate the storage of the real and imaginary parts of
for the two cases:
even and
odd.
If is even then the sequence has two purely real elements and is stored as follows:
Index of x |
0 |
1 |
2 |
|
|
|
|
|
Sequence |
|
|
|
|
|
|
|
|
Stored values |
|
|
|
|
|
|
|
|
If is odd then the sequence has one purely real element and, letting , is stored as follows:
Index of x |
0 |
1 |
2 |
|
|
|
|
|
|
Sequence |
|
|
|
|
|
|
|
|
|
Stored values |
|
|
|
|
|
|
|
|
|
The second (recommended) storage scheme, referred to in this chapter as the
complex storage format for Hermitian sequences, stores the real and imaginary parts
, for
, in consecutive locations of an array
x of length
.
The
following two tables illustrate the storage of the real and imaginary parts of
for the two cases:
even and
odd.
If is even then the sequence has two purely real elements and is stored as follows:
Index of x |
0 |
1 |
2 |
3 |
|
|
|
|
|
Stored values |
|
|
|
|
|
|
|
|
|
If is odd then the sequence has one purely real element and, letting , is stored as follows:
Index of x |
0 |
1 |
2 |
3 |
|
|
|
|
|
Stored values |
|
|
|
|
|
|
|
|
|
Also, given a Hermitian sequence, the inverse (or backward) discrete transform produces a real sequence. That is,
where
if
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
is the complex conjugate of
. It is more convenient for transformed data of two or more dimensions to be stored as a complex sequence of length
where
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 (
).
2.1.3
Real symmetric transforms
In many applications the sequence
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
is
odd,
, then the discrete Fourier transform of
contains only sine terms. Rather than compute the transform of an odd sequence, we define the
sine transform of a real sequence by
which could have been computed using the Fourier transform of a real odd sequence of length
. In this case the
are arbitrary, and the symmetry only becomes apparent when the sequence is extended. Similarly we define the
cosine transform of a real sequence by
which could have been computed using the Fourier transform of a real
even sequence of length
.
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
which could have been computed using the Fourier transform of a real sequence of length
of the form
Similarly we may define the
quarter-wave cosine transform by
which could have been computed using the Fourier transform of a real sequence of length
of the form
2.1.4
Fourier integral transforms
The usual application of the discrete Fourier transform is that of obtaining an approximation of the
Fourier integral transform
when
is negligible outside some region
. Dividing the region into
equal intervals we have
and so
for
, where
and
.
Hence the discrete Fourier transform gives an approximation to the Fourier integral transform in the region to .
If the function is defined over some more general interval , then the integral transform can still be approximated by the discrete transform provided a shift is applied to move the point to the origin.
2.1.5
Convolutions and correlations
One of the most important applications of the discrete Fourier transform is to the computation of the discrete
convolution or
correlation of two vectors
and
defined (as in
Brigham (1974)) by
- convolution:
- correlation:
(Here
and
are assumed to be periodic with period
.)
Under certain circumstances (see
Brigham (1974)) these can be used as approximations to the convolution or correlation integrals defined by
and
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).
2.1.6
Applications to solving partial differential equations (PDEs)
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).
2.2
Fast Gauss Transform
Gauss transforms have applications in areas including statistics, machine learning, and numerical solution of the heat equation. The discrete Gauss transform (DGT),
, evaluated at a set of target points
, for
, is defined as:
where
, for
, are the Gaussian source points,
, for
, are the source weights and
, for
, are the source standard deviations (alternatively source scales or source bandwidths).
The fast Gauss transform (FGT) algorithm presented in
Raykar and Duraiswami (2005) approximates the DGT by using two Taylor series and clustering of the source points.
2.3
Direct Summation of Orthogonal Series
For any series of functions
which satisfy a recurrence
the sum
is given by
where
This may be used to compute the sum of the series. For further reading, see
Hamming (1962).
3
Recommendations on Choice and Use of Available Functions
The fast Fourier transform algorithm ceases to be ‘fast’ if applied to values of which cannot be expressed as a product of small prime factors. All the FFT functions in this chapter are particularly efficient if the only prime factors of are , or .
3.1
One-dimensional Fourier Transforms
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.
3.1.1
Real and Hermitian data
nag_sum_fft_realherm_1d (c06pac) transforms a single sequence of real data onto (and in-place) a representation of the transformed Hermitian sequence using the
complex storage scheme described in
Section 2.1.2.
nag_sum_fft_realherm_1d (c06pac) also performs the inverse transform using the representation of Hermitian data and transforming back to a real data sequence.
Alternatively, the two-dimensional function
nag_sum_fft_real_2d (c06pvc) can be used (on setting the second dimension to 1) to transform a sequence of real data onto an Hermitian sequence whose first half is stored in a separate Complex array. The second half need not be stored since these are the complex conjugate of the first half in reverse order.
nag_sum_fft_hermitian_2d (c06pwc) performs the inverse operation, transforming the the Hermitian sequence (half-)stored in a Complex array onto a separate real array.
3.1.2
Complex data
nag_sum_fft_complex_1d (c06pcc) transforms a single complex sequence in-place; it also performs the inverse transform.
nag_sum_fft_complex_1d_multi (c06psc) transforms multiple complex sequences, each stored sequentially; it also performs the inverse transform on multiple complex sequences. This function is designed to perform several transforms in a single call, all with the same value of
.
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.
3.2
Half- and Quarter-wave Transforms
Four functions are provided for computing fast Fourier transforms (FFTs) of real symmetric sequences.
nag_sum_fft_sine (c06rec) computes multiple Fourier sine transforms,
nag_sum_fft_cosine (c06rfc) computes multiple Fourier cosine transforms,
nag_sum_fft_qtrsine (c06rgc) computes multiple quarter-wave Fourier sine transforms, and
nag_sum_fft_qtrcosine (c06rhc) computes multiple quarter-wave Fourier cosine transforms.
3.3
Application to Elliptic Partial Differential Equations
As described in
Section 2.1.6, Fourier transforms may be used in the solution of elliptic PDEs.
nag_sum_fft_sine (c06rec) may be used to solve equations where the solution is specified along the boundary.
nag_sum_fft_cosine (c06rfc) may be used to solve equations where the derivative of the solution is specified along the boundary.
nag_sum_fft_qtrsine (c06rgc) 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_sum_fft_qtrcosine (c06rhc) 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_sum_fft_realherm_1d (c06pac) are appropriate.
3.4
Multidimensional Fourier Transforms
The following functions compute multidimensional discrete Fourier transforms of real, Hermitian and complex data stored in Complex arrays:
The Hermitian data, either transformed from or being transformed to real data, is compacted (due to symmetry) along its first dimension when stored in Complex arrays; thus approximately half the full Hermitian data is stored.
nag_sum_fft_complex_2d (c06puc) 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.
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 functions are provided for performing forward and inverse transforms.
nag_sum_fft_real_2d (c06pvc) performs the forward two-dimensionsal transform while
nag_sum_fft_hermitian_2d (c06pwc) performs the inverse of this transform.
nag_sum_fft_real_3d (c06pyc) performs the forward three-dimensional transform while
nag_sum_fft_hermitian_3d (c06pzc) performs the inverse of this transform.
The complex sequences computed by
nag_sum_fft_real_2d (c06pvc) and
nag_sum_fft_real_3d (c06pyc) 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
are the complex conjugate of
for
, and
.
3.5
Convolution and Correlation
nag_sum_convcorr_real (c06fkc) computes either the discrete convolution or the discrete correlation of two real vectors.
3.6
Fast Gauss Transform
The only function available is
nag_sum_fast_gauss (c06sac). If the dimensionality of the data is low or the number of source and target points is small, however, it may be more efficient to evaluate the discrete Gauss transform directly.
3.7
Direct Summation of Orthogonal Series
The only function available is
nag_sum_cheby_series (c06dcc), which sums a finite Chebyshev series
depending on the choice of argument.
4
Decision Trees
Tree 1: Fourier Transform of Discrete Complex Data
Is the data one-dimensional? |
|
Multiple vectors? |
|
c06psc |
yes | yes |
| no | | | no | |
|
c06pcc |
|
|
Is the data two-dimensional? |
|
c06puc |
yes |
| no | |
Is the data three-dimensional? |
|
c06pxc |
yes |
| no | |
Transform on one dimension only? |
|
c06pfc |
yes |
| no | |
Transform on all dimensions? |
|
c06pjc |
yes |
Tree 2: Fourier Transform of Real Data or Data in Complex Hermitian Form Resulting from the Transform of Real Data
Quarter-wave sine (inverse) transform? |
|
c06rgc |
yes |
| no | |
Quarter-wave cosine (inverse) transform? |
|
c06rhc |
yes |
| no | |
Sine (inverse) transform? |
|
c06rec |
yes |
| no | |
Cosine (inverse) transform? |
|
c06rfc |
yes |
| no | |
Is the data three-dimensional? |
|
Forward transform on real data? |
|
c06pyc |
yes | yes |
| no | | | no | |
|
Inverse transform on Hermitian data? |
|
c06pzc |
| yes |
|
Is the data two-dimensional? |
|
Forward transform on real data? |
|
c06pvc |
yes | yes |
| no | | | no | |
|
Inverse transform on Hermitian data? |
|
c06pwc |
| yes |
|
Is the data multi one-dimensional? |
|
Sequences stored by row? |
|
c06fpc, c06fqc |
yes | yes |
| no | | | no | |
|
Sequences stored by column? |
|
c06pac (repeated calls) |
| yes |
|
c06pac |
|
5
Functionality Index
Convolution or Correlation, | | |
Discrete Fourier Transform, | | |
multiple half- and quarter-wave transforms, | | |
6
Auxiliary Functions Associated with Library Function Arguments
None.
7
Functions Withdrawn or Scheduled for Withdrawal
The following lists all those functions that have been withdrawn since Mark 23 of the Library or are scheduled for withdrawal at one of the next two marks.
8
References
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
Raykar V C and Duraiswami R (2005) Improved Fast Gauss Transform With Variable Source Scales University of Maryland Technical Report CS-TR-4727/UMIACS-TR-2005-34
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 transformation Math. Tables Aids Comput. 10 91–96