c13-0 (Numerical Recipes in C)

13.0 IntroductionFourier methods have revolutionized fields of science and engineering, fromradio astronomy to medical imaging, from seismology to spectroscopy. In thischapter, we present some of the basic applications of Fourier and spectral methodsthat have made these revolutions possible.Say the word “Fourier” to a numericist, and the response, as if by Pavlovianconditioning, will likely be “FFT.” Indeed, the wide application of Fourier methodsmust be credited principally to the existence of the fast Fourier transform. Bettermousetraps stand aside: If you speed up any nontrivial algorithm by a factor of amillion or so, the world will beat a path towards finding useful applications for it.The most direct applications of the FFT are to the convolution or deconvolution ofdata (§13.1), correlation and autocorrelation (§13.2), optimal filtering (§13.3), powerspectrum estimation (§13.4), and the computation of Fourier integrals (§13.9).As important as they are, however, FFT methods are not the be-all and end-allof spectral analysis.

Section 13.5 is a brief introduction to the field of time-domaindigital filters. In the spectral domain, one limitation of the FFT is that it alwaysrepresents a function’s Fourier transform as a polynomial in z = exp(2πif∆)(cf. equation 12.1.7).

Sometimes, processes have spectra whose shapes are notwell represented by this form. An alternative form, which allows the spectrum tohave poles in z, is used in the techniques of linear prediction (§13.6) and maximumentropy spectral estimation (§13.7).Another significant limitation of all FFT methods is that they require the inputdata to be sampled at evenly spaced intervals. For irregularly or incompletelysampled data, other (albeit slower) methods are available, as discussed in §13.8.So-called wavelet methods inhabit a representation of function space that isneither in the temporal, nor in the spectral, domain, but rather something in-between.Section 13.10 is an introduction to this subject. Finally §13.11 is an excursion intonumerical use of the Fourier sampling theorem.537Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use.

Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).Chapter 13.

Fourier and SpectralApplications538Chapter 13.Fourier and Spectral Applications13.1 Convolution and Deconvolution Usingthe FFTXM/2(r ∗ s)j ≡sj−k rk(13.1.1)k=−M/2+1If a discrete response function is nonzero only in some range −M/2 < k ≤ M/2,where M is a sufficiently large even integer, then the response function is called afinite impulse response (FIR), and its duration is M . (Notice that we are defining Mas the number of nonzero values of rk ; these values span a time interval of M − 1sampling times.) In most practical circumstances the case of finite M is the case ofinterest, either because the response really has a finite duration, or because we chooseto truncate it at some point and approximate it by a finite-duration response function.The discrete convolution theorem is this: If a signal sj is periodic with periodN , so that it is completely determined by the N values s0 , .

. . , sN−1 , then itsSample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-1992 by Numerical Recipes Software.Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).We have defined the convolution of two functions for the continuous case inequation (12.0.8), and have given the convolution theorem as equation (12.0.9).

Thetheorem says that the Fourier transform of the convolution of two functions is equalto the product of their individual Fourier transforms. Now, we want to deal withthe discrete case. We will mention first the context in which convolution is a usefulprocedure, and then discuss how to compute it efficiently using the FFT.The convolution of two functions r(t) and s(t), denoted r ∗ s, is mathematicallyequal to their convolution in the opposite order, s ∗ r. Nevertheless, in mostapplications the two functions have quite different meanings and characters.

One ofthe functions, say s, is typically a signal or data stream, which goes on indefinitelyin time (or in whatever the appropriate independent variable may be). The otherfunction r is a “response function,” typically a peaked function that falls to zero inboth directions from its maximum. The effect of convolution is to smear the signals(t) in time according to the recipe provided by the response function r(t), as shownin Figure 13.1.1. In particular, a spike or delta-function of unit area in s which occursat some time t0 is supposed to be smeared into the shape of the response functionitself, but translated from time 0 to time t0 as r(t − t0 ).In the discrete case, the signal s(t) is represented by its sampled values at equaltime intervals sj .

The response function is also a discrete set of numbers rk , with thefollowing interpretation: r0 tells what multiple of the input signal in one channel (oneparticular value of j) is copied into the identical output channel (same value of j);r1 tells what multiple of input signal in channel j is additionally copied into outputchannel j + 1; r−1 tells the multiple that is copied into channel j − 1; and so on forboth positive and negative values of k in rk .

Figure 13.1.2 illustrates the situation.Example: a response function with r0 = 1 and all other rk ’s equal to zerois just the identity filter: convolution of a signal with this response function givesidentically the signal. Another example is the response function with r14 = 1.5 andall other rk ’s equal to zero. This produces convolved output that is the input signalmultiplied by 1.5 and delayed by 14 sample intervals.Evidently, we have just described in words the following definition of discreteconvolution with a response function of finite duration M :.