Usage Notes

Fast Fourier Transforms

A fast Fourier transform (FFT) is simply a discrete Fourier transform that is computed efficiently. The straightforward method for computing the Fourier transform takes approximately $n^2$ operations where n is the number of points in the transform, while the FFT (which computes the same values) takes approximately n log n operations. It uses the system’s high performance library for the computation, if available. The algorithms in this chapter are modeled on the Cooley-Tukey (1965) algorithm. Hence, these functions are most efficient for integers that are highly composite; that is, integers that are a product of small primes.

For the two functions fftReal and fftComplex there is a corresponding initialization function. Use these functions only when repeatedly transforming sequences of the same length. In this situation, the initialization function computes the initial setup once. Subsequently, the user calls the corresponding main function with the appropriate option. This may result in substantial computational savings. For more information on the use of these functions, consult the documentation under the appropriate function name.

In addition to the one-dimensional transformations described above, we also provide a complex two-dimensional FFT and its inverse.

Continuous Versus Discrete Fourier Transform

There is, of course, a close connection between the discrete Fourier transform and the continuous Fourier transform. Recall that the continuous Fourier transform is defined (Brigham 1974) as

$$\hat{f}(\omega) = \int_{-\infty}^{\infty} f(t)e^{-2 \pi i \omega t}dt$$

We begin by making the following approximation:

$$\begin{split}\begin{aligned} \hat{f}(\omega) &\approx \int_{-T/2}^{T/2} f(t)e^{-2{\pi}i{\omega}t}dt \\ &= \int_0^T f(t-T/2)e^{-2{\pi}i{\omega}(t-T/2)}dt \\ &= e^{{\pi}i{\omega}T} \int_0^T f(t-T/2)e^{-2{\pi}i{\omega}t}dt \end{aligned}\end{split}$$

If we approximate the last integral using the rectangle rule with spacing $h = T / n$, we have

$$\hat{f}(\omega) = e^{\pi i \omega T}h \sum_{k=0}^{n-1} e^{-2 \pi i \omega kh} f(kh - T/2)$$

Finally, setting $\omega = j / T$ for $j = 0, ..., n − 1$ yields

$$\hat{f} (j/T) \approx e^{\pi ij} h \sum_{k=0}^{n-1} e^{-2 \pi ijk / n} f(kh - T/2) = (-1)^j \sum_{k=0}^{n-1} e^{-2 \pi ijk/n} f_k^h$$

where the vector $f^h = (f(−T/2), ..., f((n − 1) h − T/2))$. Thus, after scaling the components by $(-1)^j h$, the discrete Fourier transform, as computed in fftComplex (with input $f^h$) is related to an approximation of the continuous Fourier transform by the above formula.

If the function f is expressed as a C function, then the continuous Fourier transform

$$\hat{f}$$

can be approximated using the IMSL function intFcnFourier (Quadrature).