538Chapter 13. Fourier and Spectral ApplicationsSample 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 machine-readable files (including this one) to any servercomputer, 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 [email protected] (outside North America).13.1 Convolution and Deconvolution Usingthe FFTWe 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 convolutionof two functionsr(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-functionof unit area ins which occursat some time t0is supposed to be smeared into the shape of the response functionitself, but translated from time 0 to time t0as r(t − t0).In the discrete case, the signal s(t) is represented by its sampled values at equaltime intervals sj. The response functionis also a discrete set of numbers rk, with thefollowinginterpretation: r0tellswhat multipleof theinputsignalinonechannel (oneparticular value of j) is copied into the identical output channel (same value of j);r1tells what multiple of input signal in channel j is additionally copied into outputchannel j +1;r−1tells the multiplethat 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=1and 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.5andall 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:(r ∗ s)j≡M/2Xk=−M/2+1sj−krk(13.1.1)If 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 responsereally has a finiteduration, or because we chooseto truncate it at some point and approximate it by a finite-durationresponse function.The discrete convolution theorem is this: If a signal sjis periodic with periodN, so that it is completely determined by the N values s0,...,sN−1, then its13.1 Convolution and Deconvolution Using the FFT539Sample 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 machine-readable files (including this one) to any servercomputer, 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 [email protected] (outside North America).s(t)r(t)r*s(t)tttFigure 13.1.1. Example of the convolution of two functions. A signal s(t) is convolved with aresponse function r(t). Since the response function is broader than some features in the original signal,these are “washed out” in the convolution. In the absence of any additional noise, the process can bereversed by deconvolution.sj000N − 1rk(r*s)jN − 1N − 1Figure 13.1.2. Convolution of discretely sampledfunctions. Notehow the responsefunctionfornegativetimes is wrapped around and stored at the extreme right end of the array rk.540Chapter 13. Fourier and Spectral ApplicationsSample 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 machine-readable files (including this one) to any servercomputer, 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 [email protected] (outside North America).discrete convolution with a response function of finite duration N is a member ofthe discrete Fourier transform pair,N/2Xk=−N/2+1sj−krk⇐⇒ SnRn( 13.1.2)Here Sn, (n =0,...,N − 1) is the discrete Fourier transform of the valuessj, (j =0,...,N −1), while Rn, (n =0,...,N − 1) is the discrete Fouriertransform of the values rk, (k =0,...,N −1). These values of rkare the sameones as for the range k = −N/2+1,...,N/2, but in wrap-around order, exactlyas was described at the end of §12.2.Treatment of End Effects by Zero PaddingThe discrete convolution theorem presumes