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

(d) | Direct summation of orthogonal series. |

(e) | Acceleration of convergence of a seuqnce of real values. |

# Syntax

C# |
---|

public static class C06 |

Visual Basic |
---|

Public NotInheritable Class C06 |

Visual C++ |
---|

public ref class C06 abstract sealed |

F# |
---|

[<AbstractClassAttribute>] [<SealedAttribute>] type C06 = class end |

# Background to the Problems

# Discrete Fourier Transforms

# Complex transforms

Most of the methods in this chapter calculate the finite

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.

**discrete Fourier transform**(DFT) of a sequence of $n$ complex numbers ${z}_{\mathit{j}}$, for $\mathit{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}\mathrm{exp}\left(-i\frac{2\pi jk}{n}\right)$$ | (1) |

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

The original data values ${z}_{j}$ may conversely be recovered from the transform ${\hat{z}}_{k}$ by an

for $j=0,1,\dots ,n-1$. If we take the complex conjugate of (2), we find that the sequence ${\stackrel{-}{z}}_{j}$ is the DFT of the sequence ${\stackrel{-}{\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

$${a}_{k}=\frac{1}{\sqrt{n}}\sum _{j=0}^{n-1}\left({x}_{j}\mathrm{cos}\left(\frac{2\pi jk}{n}\right)+{y}_{j}\mathrm{sin}\left(\frac{2\pi jk}{n}\right)\right)$$ |

$${b}_{k}=\frac{1}{\sqrt{n}}\sum _{j=0}^{n-1}\left({y}_{j}\mathrm{cos}\left(\frac{2\pi jk}{n}\right)-{x}_{j}\mathrm{sin}\left(\frac{2\pi jk}{n}\right)\right)\text{.}$$ |

**inverse discrete Fourier transform**:$${z}_{j}=\frac{1}{\sqrt{n}}\sum _{k=0}^{n-1}{\hat{z}}_{k}\mathrm{exp}\left(+i\frac{2\pi jk}{n}\right)$$ | (2) |

**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}}\mathrm{exp}\left(-i\frac{2\pi {j}_{1}{k}_{1}}{{n}_{1}}\right)\mathrm{exp}\left(-i\frac{2\pi {j}_{2}{k}_{2}}{{n}_{2}}\right)\text{.}$$ | (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.# Real transforms

If the original sequence is purely real valued, i.e., ${z}_{j}={x}_{j}$, then

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

and, if $n$ is even, ${b}_{n/2}=0$.

$${\hat{z}}_{k}={a}_{k}+i{b}_{k}=\frac{1}{\sqrt{n}}\sum _{j=0}^{n-1}{x}_{j}\mathrm{exp}\left(-i\frac{2\pi jk}{n}\right)$$ |

**Hermitian**sequence, or**half-complex**or**conjugate symmetric**, with the properties$${a}_{n-k}={a}_{k}\text{\hspace{1em}}{b}_{n-k}={-b}_{k}\text{\hspace{1em}}{b}_{0}=0$$ |

Thus a Hermitian sequence of $n$ complex data values can be represented by only $n$, rather than $2n$, 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 ${a}_{k}$ for $0\le k\le 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 $\left(0:n-1\right)$ in your calling method, 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}+{ib}_{1}$ | ${a}_{2}+{ib}_{2}$ | $\dots $ | ${a}_{n/2}$ | $\dots $ | ${a}_{2}-{ib}_{2}$ | ${a}_{1}-{ib}_{1}$ |

Stored values | ${a}_{0}$ | ${a}_{1}$ | ${a}_{2}$ | $\dots $ | ${a}_{n/2}$ | $\dots $ | ${b}_{2}$ | ${b}_{1}$ |

$$\begin{array}{cc}\mathbf{x}\left[k\right]={a}_{k}\text{,}\hfill & \text{for}k=0,1,\dots ,n/2\text{, and}\hfill \\ \mathbf{x}\left[n-k\right]={b}_{k}\text{,}\hfill & \text{for}k=1,2,\dots ,n/2-1\text{.}\hfill \end{array}$$ |

If $n$ is odd then the sequence has one purely real element and, letting $n=2s+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}+{ib}_{1}$ | ${a}_{2}+{ib}_{2}$ | $\dots $ | ${a}_{s}+i{b}_{s}$ | ${a}_{s}-i{b}_{s}$ | $\dots $ | ${a}_{2}-{ib}_{2}$ | ${a}_{1}-{ib}_{1}$ |

Stored values | ${a}_{0}$ | ${a}_{1}$ | ${a}_{2}$ | $\dots $ | ${a}_{s}$ | ${b}_{s}$ | $\dots $ | ${b}_{2}$ | ${b}_{1}$ |

$$\begin{array}{cc}\mathbf{x}\left[k\right]={a}_{k}\text{,}\hfill & \text{for}k=0,1,\dots ,s\text{, and}\hfill \\ \mathbf{x}\left[n-k\right]={b}_{k}\text{,}\hfill & \text{for}k=1,2,\dots ,s\text{.}\hfill \end{array}$$ |

The second (recommended) 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\le k\le n/2$, in consecutive locations of an array x of length $n+2$. If x is declared with bounds $\left(0:n+1\right)$ in your calling method, 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{array}{cc}\mathbf{x}\left[2\times k-1\right]={a}_{k}\text{,}\hfill & \text{for}k=0,1,\dots ,n/2\text{, and}\hfill \\ \mathbf{x}\left[2\times k+1-1\right]={b}_{k}\text{,}\hfill & \text{for}k=0,1,\dots ,n/2\text{.}\hfill \end{array}$$ |

If $n$ is odd then the sequence has one purely real element and, letting $n=2s+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{array}{cc}\mathbf{x}\left[2\times k-1\right]={a}_{k}\text{,}\hfill & \text{for}k=0,1,\dots ,s\text{, and}\hfill \\ \mathbf{x}\left[2\times k+1-1\right]={b}_{k}\text{,}\hfill & \text{for}k=0,1,\dots ,s\text{.}\hfill \end{array}$$ |

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

where ${a}_{n/2}=0$ if $n$ is odd.

$${x}_{j}=\frac{1}{\sqrt{n}}\left({a}_{0}+2\sum _{k=1}^{n/2-1}\left({a}_{k}\mathrm{cos}\left(\frac{2\pi jk}{n}\right)-{b}_{k}\mathrm{sin}\left(\frac{2\pi jk}{n}\right)\right)+{a}_{n/2}\right)$$ |

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}}_{\left({n}_{1}-{k}_{1}\right){k}_{2}}$. It is more convenient for transformed data of two or more dimensions to be stored as a complex sequence of length $\left({n}_{1}/2+1\right)\times {n}_{2}\times \cdots \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 \cdots \times {n}_{d}$).

# Real symmetric transforms

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

which could have been computed using the Fourier transform of a real odd sequence of length $2n$. In this case the ${x}_{j}$ are arbitrary, and the symmetry only becomes apparent when the sequence is extended. Similarly we define the

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

**odd**, $\left({x}_{j}={-x}_{n-j}\right)$, 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}\mathrm{sin}\left(\frac{\pi jk}{n}\right)\text{,}$$ |

**cosine transform**of a real sequence by$${\hat{x}}_{k}=\sqrt{\frac{2}{n}}\left(\frac{1}{2}{x}_{0}+\sum _{j=1}^{n-1}{x}_{j}\mathrm{cos}\left(\frac{\pi jk}{n}\right)+\frac{1}{2}{\left(-1\right)}^{k}{x}_{n}\right)$$ |

**even**sequence of length $2n$.In addition to these ‘half-wave’ symmetries described above, sequences arise in practice with ‘quarter-wave’ symmetries. We define the

which could have been computed using the Fourier transform of a real sequence of length $4n$ of the form

Similarly we may define the

which could have been computed using the Fourier transform of a real sequence of length $4n$ of the form

**quarter-wave sine transform**by$${\hat{x}}_{k}=\frac{1}{\sqrt{n}}\left(\sum _{j=1}^{n-1}{x}_{j}\mathrm{sin}\left(\frac{\pi j\left(2k-1\right)}{2n}\right)+\frac{1}{2}{\left(-1\right)}^{k-1}{x}_{n}\right)$$ |

$$\left(0,{x}_{1},\dots ,{x}_{n},{x}_{n-1},\dots ,{x}_{1},0,{-x}_{1},\dots ,{-x}_{n},{-x}_{n-1},\dots ,{-x}_{1}\right)\text{.}$$ |

**quarter-wave cosine transform**by$${\hat{x}}_{k}=\frac{1}{\sqrt{n}}\left(\frac{1}{2}{x}_{0}+\sum _{j=1}^{n-1}{x}_{j}\mathrm{cos}\left(\frac{\pi j\left(2k-1\right)}{2n}\right)\right)$$ |

$$\left({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}\right)\text{.}$$ |

# Fourier integral transforms

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

when $f\left(t\right)$ is negligible outside some region $\left(0,c\right)$. Dividing the region into $n$ equal intervals we have

and so

for $k=0,1,\dots ,n-1$, where ${f}_{j}=f\left(jc/n\right)$ and ${F}_{k}=F\left(k/c\right)$.

**Fourier integral transform**$$F\left(s\right)=\underset{-\infty}{\overset{\infty}{\int}}f\left(t\right)\mathrm{exp}\left(-i2\pi st\right)dt$$ |

$$F\left(s\right)\cong \frac{c}{n}\sum _{j=0}^{n-1}{f}_{j}\mathrm{exp}\left(\frac{-i2\pi sjc}{n}\right)$$ |

$${F}_{k}\cong \frac{c}{n}\sum _{j=0}^{n-1}{f}_{j}\mathrm{exp}\left(\frac{-i2\pi jk}{n}\right)$$ |

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\left(t\right)$ is defined over some more general interval $\left(a,b\right)$, 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.

# 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 $x$ and $y$ defined (as in Brigham (1974)) by- convolution: ${z}_{k}={\displaystyle \sum _{j=0}^{n-1}}{x}_{j}{y}_{k-j}$
- correlation: ${w}_{k}={\displaystyle \sum _{j=0}^{n-1}}{\stackrel{-}{x}}_{j}{y}_{k+j}$

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

$$z\left(s\right)=\underset{-\infty}{\overset{\infty}{\int}}x\left(t\right)y\left(s-t\right)dt$$ |

$$w\left(s\right)=\underset{-\infty}{\overset{\infty}{\int}}\stackrel{-}{x}\left(t\right)y\left(s+t\right)dt\text{, \hspace{1em}}-\infty <s<\infty \text{.}$$ |

# 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 [Application to Elliptic Partial Differential Equations]. Full details of the Fourier method for the solution of PDEs may be found in Swarztrauber (1977) and Swarztrauber (1984).

# Inverse Laplace Transforms

Let $f\left(t\right)$ be a real function of $t$, with $f\left(t\right)=0$ for $t<0$, and be piecewise continuous and of exponential order $\alpha $, i.e.,

for large $t$, where $\alpha $ is the minimal such exponent.

$$\left|f\left(t\right)\right|\le M{e}^{\alpha t}$$ |

The Laplace transform of $f\left(t\right)$ is given by

where $F\left(s\right)$ is defined for $\mathrm{Re}\left(s\right)>\alpha $.

$$F\left(s\right)=\underset{0}{\overset{\infty}{\int}}{e}^{-st}f\left(t\right)dt\text{, \hspace{1em}}t>0$$ |

The inverse transform is defined by the Bromwich integral

The integration is performed along the line $s=a$ in the complex plane, where $a>\alpha $. This is equivalent to saying that the line $s=a$ lies to the right of all singularities of $F\left(s\right)$. For this reason, the value of $\alpha $ is crucial to the correct evaluation of the inverse. It is not essential to know $\alpha $ exactly, but an upper bound must be known.

$$f\left(t\right)=\frac{1}{2\pi i}\underset{a-i\infty}{\overset{a+i\infty}{\int}}{e}^{st}F\left(s\right)ds\text{, \hspace{1em}}t>0\text{.}$$ |

The problem of determining an inverse Laplace transform may be classified according to whether (a) $F\left(s\right)$ is known for real values only, or (b) $F\left(s\right)$ is known in functional form and can therefore be calculated for complex values of $s$. Problem (a) is very ill-defined and no methods are provided. Two methods are provided for problem (b).

# Direct Summation of Orthogonal Series

For any series of functions ${\varphi}_{i}$ 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).

$${\varphi}_{r+1}\left(x\right)+{\alpha}_{r}\left(x\right){\varphi}_{r}\left(x\right)+{\beta}_{r}\left(x\right){\varphi}_{r-1}\left(x\right)=0$$ |

$$\sum _{r=0}^{n}{a}_{r}{\varphi}_{r}\left(x\right)$$ |

$$\sum _{r=0}^{n}{a}_{r}{\varphi}_{r}\left(x\right)={b}_{0}\left(x\right){\varphi}_{0}\left(x\right)+{b}_{1}\left(x\right)\left({\varphi}_{1}\left(x\right)+{\alpha}_{0}\left(x\right){\varphi}_{0}\left(x\right)\right)$$ |

$${b}_{r}\left(x\right)+{\alpha}_{r}\left(x\right){b}_{r+1}\left(x\right)+{\beta}_{r+1}\left(x\right){b}_{r+2}\left(x\right)={a}_{r}{b}_{n+1}\left(x\right)={b}_{n+2}\left(x\right)=0\text{.}$$ |

# Acceleration of Convergence

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 $\left\{{s}_{n}\right\}$, for $\mathit{n}=m,\dots ,m+2l$, then, except in certain singular cases, parameters, $a$, ${b}_{i}$, ${c}_{i}$ may be determined such that

If the sequence $\left\{{s}_{n}\right\}$ 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.

$${s}_{n}=a+\sum _{i=1}^{l}{b}_{i}{c}_{i}^{n}\text{.}$$ |

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}={\displaystyle \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

and to consider the sequence of values

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 not in this release) 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\left(x\right)$, so that successive contributions to the integral are of opposite sign. As an example, consider the case where $f\left(x\right)=M\left(x\right)\mathrm{sin}\u200ax$ and $M\left(x\right)>0$: convergence will be much improved if ${a}_{k}=k\pi $ rather than ${a}_{k}=2k\pi $.

$$\underset{0}{\overset{\infty}{\int}}f\left(x\right)dx=\underset{0}{\overset{{a}_{1}}{\int}}f\left(x\right)dx+\underset{{a}_{1}}{\overset{{a}_{2}}{\int}}f\left(x\right)dx+\underset{{a}_{2}}{\overset{{a}_{3}}{\int}}f\left(x\right)dx+\dots $$ |

$${s}_{1}=\underset{0}{\overset{{a}_{1}}{\int}}f\left(x\right)dx\text{, \hspace{1em}}{s}_{2}=\underset{0}{\overset{{a}_{2}}{\int}}f\left(x\right)dx={s}_{1}+\underset{{a}_{1}}{\overset{{a}_{2}}{\int}}f\left(x\right)dx\text{, etc.,}$$ |

# References

Brigham E O (1974)

*The Fast Fourier Transform*Prentice–HallDavies S B and Martin B (1979) Numerical inversion of the Laplace transform: A survey and comparison of methods

*J. Comput. Phys.***33**1–32Fox L and Parker I B (1968)

*Chebyshev Polynomials in Numerical Analysis*Oxford University PressGentleman W S and Sande G (1966) Fast Fourier transforms for fun and profit

*Proc. Joint Computer Conference, AFIPS***29**563–578Hamming R W (1962)

*Numerical Methods for Scientists and Engineers*McGraw–HillShanks D (1955) Nonlinear transformations of divergent and slowly convergent sequences

*J. Math. Phys.***34**1–42Swarztrauber 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–501Swarztrauber P N (1984) Fast Poisson solvers

*Studies in Numerical Analysis*(ed G H Golub) Mathematical Association of AmericaSwarztrauber P N (1986) Symmetric FFT's

*Math. Comput.***47(175)**323–346Wynn P (1956) On a device for computing the ${e}_{m}\left({S}_{n}\right)$ transformation

*Math. Tables Aids Comput.***10**91–96