Discrete Fourier Transform

  

 

Input description: A sequence of n real or complex values , , sampled at uniform intervals from a function h.

Problem description: The discrete Fourier transform H of h, , for .

Discussion: Although computer scientists tend to be relatively unfamiliar with Fourier transforms, electrical engineers and signal processors eat them for breakfast. Functionally, Fourier transforms provide a way to convert samples of a standard time-series into the ``frequency domain''. This provides a ``dual'' representation of the function, in which certain operations become easier than in the time domain. Applications of Fourier transforms include:         

• Filtering - Taking the Fourier transform of a function is equivalent to representing it as the sum of sine functions.       By eliminating undesirable high- and/or low-frequency components (i.e. dropping some of the sine functions) and taking an inverse Fourier transform to get us back into the time domain, we can filter an image to remove noise and other artifacts.   For example, the sharp spike in the figure above describes the period of a single sine function that closely models the input data.
• Image Compression - A smoothed, filtered image contains less information than a noisy image, while retaining a similar appearance. Thus encoding the smoothed image will require fewer bits to represent than the original image.     By eliminating the coefficients of sine functions that contribute relatively little to the image, we can further reduce the size of the image, at little cost in image fidelity.
• Convolution and Deconvolution - Fourier transforms can be used to efficiently compute convolutions of two sequences.    A convolution is the pairwise product of elements from two different sequences, such as in multiplying two n-variable polynomials f and g or multiplying two long integers.     Implementing the product directly takes , while suffices using the fast Fourier transform. Another example comes from image processing.   Because a scanner measures the darkness of an image patch instead of a single point, the scanned input is always blurred. A reconstruction of the original signal can be obtained by deconvoluting the input signal with a Gaussian point-spread function.   
• Computing the correlation of functions -   The correlation function of two functions f(t) and g(t) is defined by

  

  and can be easily computed using Fourier transforms. Note that if the two functions are similar in shape but one is shifted relative to the other (such as and ), the value of will be large at this shift offset .   As an application, suppose that we want to detect whether there are any funny periodicities in our random number generator.   We can generate a large series of random numbers, turn it into a time series (the ith number at time i), and take the Fourier transform of this series. Any funny spikes will correspond to potential periodicities.  

The discrete Fourier transform takes as input n   complex numbers , , corresponding to equally spaced points in a time series, and outputs n complex numbers , , each describing a sine function of given frequency. The discrete Fourier transform is defined by

and the inverse Fourier transform is defined by

which enables us move easily between h and H.

Since the output of the discrete Fourier transform consists of n numbers, each of which can be computed using a formula on n numbers, they can be computed in time. The fast Fourier transform (FFT) is an algorithm that computes the discrete Fourier transform in . This is arguably the most important algorithm known, as measured by practical impact, for it opened the door to modern image processing. There are several different algorithms that call themselves FFTs, all of which are based on a divide-and-conquer approach. Essentially, the problem of computing the discrete Fourier transform on n points is reduced to computing two transforms on n/2 points each and is then applied recursively     

The FFT usually assumes that n is a power of two. If this is not the case for your data, you are usually better off padding your data with zeros to create elements rather than hunting for a more general code.

Since many image processing systems have strong real-time constraints, FFTs are often implemented in hardware, or at least in assembly language tuned to the particular machine. Be aware of this possibility if the codes below prove too slow.    

Implementations: FFTPACK is a package of Fortran subprograms for the fast Fourier transform of periodic and other     symmetric sequences, written by P. Swartzrauber. It includes complex, real, sine, cosine, and quarter-wave transforms. A C language translation of the main routines is also provided.   FFTPACK resides on Netlib (see Section ) at http://www.netlib.org/fftpack.  

Algorithm 545 [Fra79] of the Collected Algorithms of the ACM is an implementation of the fast Fourier transform optimizing virtual memory performance and written in Fortran.     See Section for further information.

XTango (see Section ) is an algorithm animation system for UNIX and X-windows, which includes an interesting animation of the fast Fourier transform.  

A Pascal implementation of the fast Fourier transform for points appears in [MS91].   For more details, see Section . Sedgewick [Sed92] provides a bare bones implementation of the fast Fourier transform in C++.   See Section for details.

Notes: Brigham [Bri74] is an excellent introduction to Fourier transforms and the FFT and is strongly recommended, as is the exposition in [PFTV86]. Expositions in algorithms texts on the fast Fourier transform include [AHU74, Baa88, CLR90, Man89].

Credit for inventing the fast Fourier transform is usually given to Cooley and Tukey [CT65], although it is not completely deserved. See [Bri74] for a complete history.

An interesting divide-and-conquer algorithm for polynomial multiplication [KO63] does the job in time and is discussed in [AHU74, Man89]. An FFT-based algorithm that multiplies two n-bit numbers in time is due to Schönhage and Strassen [SS71] and is presented in [AHU74].

In recent years, wavelets have been proposed to replace Fourier transforms in filtering.   See [Dau92] for an introduction to wavelets.

Related Problems: Data compression (see page ), high-precision arithmetic (see page ).