f13-0

Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).Chapter 13. Fourier and Spectral Applications 13.0 Introduction Fourier methods have revolutionized fields of science and engineering, from radio astronomy to medical imaging, from seismology to spectroscopy. In this chapter, we present some of the basic applications of Fourier and spectral methods that have made these revolutions possible. Say the word “Fourier” to a numericist, and the response, as if by Pavlovian conditioning,will likely be “FFT.” Indeed, the wide application of Fourier methodsmust be credited principally to the existence of the fast Fourier transform. Better mousetraps stand aside: If you speed up anynontrivial algorithm by a factor of a million 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 of data(§13.1),correlationandautocorrelation( §13.2),optimalfiltering( §13.3),power spectrum estimation ( §13.4),and the computationof Fourierintegrals ( §13.9). As important as they are, however, FFT methods are not the be-all and end-all of spectral analysis. Section 13.5 is a brief introductionto the field of time-domain digital filters. In the spectral domain, one limitation of the FFT is that it always represents a function’s Fourier transform as a polynomial in z=e x p ( 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 to have poles in z, is used in the techniquesof linear prediction( §13.6)and maximum entropy spectral estimation ( §13.7). Another significant limitation of all FFT methods is that they require the input data to be sampled at evenly spaced intervals. For irregularly or incompletely sampled data, other(albeit slower) methodsare available, as discussed in §13.8. So-called wavelet methods inhabit a representation of function space that is neitherinthetemporal,norinthespectral,domain,butrathersomethingin-between.Section 13.10 is an introduction to this subject. Finally §13.11 is an excursion into numerical use of the Fourier sampling theorem. 530 13.1ConvolutionandDeconvolutionUsingtheFFT 531Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).13.1 Convolution and Deconvolution Using the FFT We have defined the convolution of two functions for the continuous case in equation(12.0.8),andhavegiventhe convolutiontheorem as equation(12.0.9). The theoremsays that the Fouriertransformof the convolutionof two functionsis equal to the product of their individual Fourier transforms. Now, we want to deal with the discrete case. We will mention first the context in which convolutionis a usefulprocedure, and then discuss how to compute it efficiently using the FFT. Theconvolutionoftwofunctions r(t)ands(t),denoted r∗s,ismathematically equal to their convolution in the opposite order, s∗r. Nevertheless, in most applicationsthe two functionshavequite differentmeaningsandcharacters. Oneof the functions, say s, is typically a signal or data stream, which goes on indefinitely in time (or in whatever the appropriate independent variable may be). The other function ris a “response function,” typically a peaked function that falls to zero in both directions from its maximum. The effect of convolutionis to smear the signals(t)intimeaccordingtotherecipeprovidedbytheresponsefunction r(t),asshown inFigure13.1.1. Inparticular,aspikeordelta-functionofunitareain swhichoccurs at some time t 0is supposed to be smeared into the shape of the response function itself, but translated from time 0 to time t0asr(t−t0). Inthediscretecase,thesignal s(t)isrepresentedbyitssampledvaluesatequal timeintervals sj. Theresponsefunctionisalsoadiscretesetofnumbers rk,withthe followinginterpretation: r0tellswhatmultipleoftheinputsignalinonechannel(one particular value of j) is copied into the identical output channel (same value of j); r1tells what multiple of input signal in channel jis additionally copied into output channel j+1;r−1tells themultiplethatis copiedintochannel j−1; andso onfor bothpositiveandnegativevalues of kinrk. Figure13.1.2illustrates thesituation. Example: a response function with r0=1and all other rk’s equal to zero is just the identity filter: convolution of a signal with this response function givesidenticallythe signal. Anotherexampleis the responsefunctionwith r 14=1.5and all other rk’s equal to zero. This produces convolvedoutput that is the input signal multiplied by 1.5and delayed by 14sample intervals. Evidently, we have just described in words the following definition of discrete convolution with a response function of finite duration M: (r∗s)j≡M/ 2/summationdisplay k=−M/ 2+1sj−krk (13.1.1 ) If a discrete response function is nonzero only in some range −M/2<k≤M/2, where Mis a sufficiently large even integer, then the response function is called a finiteimpulseresponse(FIR) ,anditsdurationisM. (Noticethatwearedefining M as the number of nonzero valuesofrk; these values span a time interval of M−1 sampling times.) In most practical circumstances the case of finite Mis the case of interest,eitherbecausetheresponsereallyhasafiniteduration,orbecausewechoose totruncateitatsomepointandapproximateitbyafinite-durationresponsefunction. Thediscrete convolutiontheorem is this: If a signal sjisperiodicwith period N, so that it is completely determined by the Nvalues s0,...,s N−1, then its