f13-1

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)and s(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 532 Chapter13. FourierandSpectralApplicationsSample 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).s(t) r(t) r*s(t)t t t Figure 13.1.1. Example of the convolution of two functions. A signal s(t)is convolved with a response 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. sj0 0 0N − 1rk (r*s)jN − 1 N − 1 Figure13.1.2. Convolution ofdiscretely sampledfunctions. Notehowtheresponsefunction fornegative times is wrapped around and stored at the extreme right end of the array rk. 13.1ConvolutionandDeconvolutionUsingtheFFT 533Sample 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).discrete convolution with a response function offinite duration Nis a member of the discrete Fourier transform pair, N/2/summationdisplay k=−N/2+1sj−krk⇐⇒ SnRn (13.1.2 ) Here Sn,(n=0 ,...,N −1)is the discrete Fourier transform of the values sj,(j=0 ,...,N −1), while Rn,(n=0 ,...,N −1)is the discrete Fourier transform of the values rk,(k=0 ,...,N −1). These values of rkare the same ones as for the range k=−N/2+1 ,...,N / 2, but in wrap-around order, exactly as was described at the end of §12.2. Treatment of End Effects by ZeroPadding The discrete convolution theorem presumes a set of two circumstances that are not universal. First, it assumes that the input signal is periodic, whereas real data often either go forever without repetition or else consist of one nonperiodicstretch of finite length. Second, the convolution theorem takes the duration of the response to be the same as the period of the data; they are both N. We need to work around these two constraints. The second is very straightforward. Almost always, one is interested in a response function whose duration Mis much shorter than the length of the data set N. In this case, you simply extend the response function to length Nby padding it with zeros, i.e., de finer k=0for M/ 2≤k≤N/2and also for −N/2+1≤k≤− M/ 2+1. Dealingwiththe firstconstraintismorechallenging. Since the convolution theorem rashly assumes that the data are periodic, it will falsely“pollute”thefirst output channel (r∗s)0with some wrapped-around data from the far end of the data stream sN−1,sN−2, etc. (See Figure 13.1.3.) So, we need to set up a buffer zone of zero-padded values at the end of the sjvector, in order to make this pollution zero. How many zero values do we need in this buffer? Exactlyas manyas the most negativeindexforwhichthe responsefunction is nonzero. Forexample,if r−3is nonzero,while r−4,r−5,...are all zero,then we need three zero pads at the end of the data: sN−3=sN−2=sN−1=0. These zeros will protect the first output channel (r∗s)0from wrap-around pollution. It should be obviousthat the second output channel (r∗s)1and subsequent ones will also be protected by these same zeros. Let Kdenote the number of padding zeros, so that the last actual input data point is sN−K−1. What now about pollution of the very lastoutput channel? Since the data now end with sN−K−1, the last output channel of interest is (r∗s)N−K−1. This channel can be polluted by wrap-around from input channel s0unless the number Kis also large enough to take care of the most positive index kfor which the response function rkis nonzero. For example, if r0through r6are nonzero, while r7,r8...are all zero, then we need at least K=6padding zeros at the end of the data: sN−6=... =sN−1=0. To summarize —we need to pad the data with a number of zeros on one endequal to the maximum positive duration ormaximum negative duration of the response function, whichever is larger . (For a symmetric response function of duration M,youwill needonly M/ 2zeropads.) Combiningthis operationwiththe 534 Chapter13. FourierandSpectralApplicationsSample 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).m+ spoiled spoiled unspoiledm−response function sample of original function convolutionm+m− Figure 13.1.3. The wrap-around problem in convolving finite segments of a function. Not only must the response function wrap be viewed as cyclic, but so must the sampled original function. Thereforea portion at each end of the original function is erroneously wrapped around by convolution with the response function. response function m+ m− m− m+ m−m+zero padding original function spoiled but irrelevantunspoilednot spoiled because zero Figure 13.1.4. Zero padding as solution to the wrap-around problem. The original function is extended by zeros, serving a dual purpose: When the zeros wrap around, they do not disturb the true convolution; and while the original function wraps around onto the zero region, that region can be discarded. 13.1ConvolutionandDeconvolutionUsingtheFFT 535Sample 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).padding of the response rkdescribed above, we effectively insulate the data from artifacts of undesired periodicity. Figure 13.1.4 illustrates matters. Use of FFTfor Convolution The data, complete with zero padding, are now a set of real numbers sj,j = 0,...,N −1, and the response function is zero padded out to duration Nand arrangedinwrap-aroundorder. (Generallythismeansthatalargecontiguoussection of the rk’s, in the middle of that array, is zero, with nonzero values clustered at the two extreme ends of the array.) You now compute the discrete convolution as follows: UsetheFFTalgorithmtocomputethediscreteFouriertransformof sandof r. Multiplythetwotransformstogethercomponentbycomponent,rememberingthat the transformsconsist ofcomplexnumbers. Thenuse the FFT algorithmto take the inversediscreteFouriertransformoftheproducts. Theansweristheconvolution r∗s. What about deconvolution ? Deconvolution is the process of undoingthe smearing in a data set that has occurred under the in fluence of a known response function,for example,because of the known effect of a less-than-perfectmeasuring apparatus. Thede finingequationofdeconvolutionisthesameasthatforconvolution, namely (13.1.1),exceptnow the left-handside is taken to be known,and (13.1.1)is tobeconsideredasasetof Nlinearequationsfortheunknownquantities sj. Solving these simultaneous linear equations in the time domain of (13.1.1) is unrealistic in most cases, but the FFT renders the problem almost trivial. Instead of multiplying the transformof the signal and response to get the transform of the convolution,wejustdividethetransformofthe(known)convolutionbythetransformoftheresponse to get the transform of the deconvolved signal. This procedure can go wrong mathematically if the transform of the response function is exactly zero for some value R n, so that we can ’t divide by it. This indicates that the original convolution has truly lost all information at that onefrequency, so that a reconstruction of that frequency component is not possible. You should be aware, however, that apart from mathematical problems, the process of deconvolution has other practical shortcomings. The process is generally quitesensitivetonoiseintheinputdata,andtotheaccuracytowhichtheresponsefunction r kisknown. Perfectlyreasonableattemptsatdeconvolutioncansometimesproduce nonsenseforthesereasons. Insuchcasesyoumaywanttomakeuseoftheadditional process of optimalfiltering, which is discussed in §13.3. Here is our routine for convolution and deconvolution, using the FFT as implemented in four1of§12.2. Since the data and response functions are real, not complex, both of their transforms can be taken simultaneously using twofft. Note, however, that two calls to realftshould be substituted if dataandrespns have very different magnitudes, to minimize roundoff. The data are assumed to be stored in a real array dataof length n, which must be an integer power of two. The response function is assumed to be stored in wrap-around order in a real array respnsof length m. The value of mcan be any oddinteger less than or equal to n, since the first thing the program does is to recopy the response function into the appropriate wrap-around order in an array of length n. The answer is returned in ans, which is also used as working space. 536 Chapter13. FourierandSpectralApplicationsSample 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).SUBROUTINE convlv(data,n,respns,m,isign,ans) INTEGER isign,m,n,NMAX REAL data(n),respns(n) COMPLEX ans(n)PARAMETER (NMAX=4096) Maximum anticipated size of FFT. C USES realft,twofft Convolves or deconvolves a real data set data(1:n) (including any user-supplied zero padding) with a response function respns , stored in wrap-around order in a real array of length m≤n.(mshould be an odd integer.) Wrap-around order means that the first half of the array respns contains the impulse response function at positive times, while the second half of the array contains the impulse response function at negative times, countingdown from the highest element respns(m) . On input isign is+1for convolution, −1 for deconvolution. The answer is returned in the first ncomponents of ans. However, ans must be supplied in the calling program with length at least 2*n, for consistency with twofft .nMUST be an integer power of two. INTEGER i,no2 COMPLEX fft(NMAX) do11i=1,(m-1)/2 Putrespns in array of length n. respns(n+1-i)=respns(m+1-i) enddo 11 do12i=(m+3)/2,n-(m-1)/2 P a dw i t hz e r o s . respns(i)=0.0 enddo 12 call twofft(data,respns,fft,ans,n) FFT both at once. no2=n/2do 13i=1,no2+1 if (isign.eq.1) then ans(i)=fft(i)*ans(i)/no2 Multiply FFTs to convolve. else if (isign.eq.-1) then if (abs(ans(i)).eq.0.0) pause ’deconvolving at response zero in convlv’ ans(i)=fft(i)/ans(i)/no2 Divide FFTs to deconvolve. else pause ’no meaning for isign in convlv’ endif enddo 13 ans(1)=cmplx(real(ans(1)),real(ans(no2+1))) Pack last element with first for realft . call realft(ans,n,-1) Inverse transform back to time domain. return END Convolvingor DeconvolvingVery Large Data Sets If your data set is so long that you do not want to fit it into memory all at once, then you must break it up into sections and convolve each section separately. Now, however, the treatment of end effects is a bit different. You have to worry not only about spuriouswrap-aroundeffects, but also aboutthe fact that the ends ofeach section of data shouldhave been in fluenced by data at the nearby ends of the immediately preceding and following sections of data, but were not so in fluenced since only one section of data is in the machine at a time. There are two, related, standard solutions to this problem. Both are fairly obvious,sowitha fewwordsofdescriptionhere,yououghtto beabletoimplementthem for yourself. The first solution is called the overlap-save method . In this technique you pad only the very beginning of the data with enough zeros to avoid wrap-around pollution. After this initial padding, you forget about zero paddingaltogether. Bringina sectionofdataandconvolveordeconvolveit. Thenthrowout the pointsat each endthat arepollutedbywrap-aroundendeffects. Outputonlythe remaininggood points in the middle. Now bring in the next section of data, but not all new data. The first points in each new section overlap the last points from the 13.1ConvolutionandDeconvolutionUsingtheFFT 537Sample 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).convolution (out)AA + BB B + CCCcc b a A b B00 a 0 0 00data (in) Figure 13.1.5. The overlap-add method for convolving a response with a very long signal. The signal data is broken up into smaller pieces. Each is zero padded at both ends and convolved (denoted by bold arrows in the figure). Finally the pieces are added back together, including the overlapping regions formed by the zero pads. preceding section of data. The sections must be overlapped suf ficiently so that the polluted output points at the end of one section are recomputed as the first of the unpollutedoutputpointsfromthe subsequentsection. With a bitofthoughtyoucan easily determine how many points to overlap and save. The second solution, called the overlap-add method , is illustrated in Figure 13.1.5. Here you don’toverlapthe input data. Each section of data is disjoint from the othersandis used exactlyonce. However,youcarefullyzero-padit at bothends sothatthereisnowrap-aroundambiguityintheoutputconvolutionordeconvolution. Now you overlap and addthese sections of output. Thus, an output point near the end of one section will have the response due to the input points at the beginning of the next section of data properly added in to it, and likewise for an output point near the beginning of a section, mutatis mutandis . Evenwhencomputermemoryisavailable,thereissomeslightgainincomputing speedinsegmentingalongdataset,sincetheFFTs ’Nlog2Nisslightlyslowerthan linearin N. However,the logterm is so slowlyvaryingthat youwill oftenbe much happier to avoid the bookkeeping complexities of the overlap-add or overlap-save methods: If it is practical to do so, just cram the whole data set into memory andFFT away. Then you will have moretime for the finer things in life, some of which are described in succeeding sections of this chapter. CITED REFERENCES AND FURTHER READING: Nussbaumer,H.J.1982, FastFourierTransformandConvolutionAlgorithms (NewYork:Springer- Verlag). Elliott,D.F.,andRao,K.R.1982, FastTransforms:Algorithms,Analyses,Applications (NewYork: Academic Press). Brigham, E.O. 1974, The Fast Fourier Transform (Englewood Cliffs, NJ: Prentice-Hall), Chap- ter 13. 538 Chapter13. FourierandSpectralApplicationsSample 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.2 Correlation and Autocorrelation Using the FFT Correlation is the close mathematical cousin of convolution. It is in some ways simpler, however, because the two functions that go into a correlation are notas conceptually distinct as were the data and response functions that entered into convolution. Rather, in correlation, the functions are represented by different, but generally similar, data sets. We investigate their “correlation, ”by comparing them both directly superposed, and with one of them shifted left or right. We have already de fined in equation (12.0.10) the correlation between two continuous functions g(t)and h(t), which is denoted Corr (g, h ), and is a function oflag t. We will occasionally show this time dependenceexplicitly,with the rather awkwardnotationCorr (g, h )(t). Thecorrelationwillbelargeatsomevalueof tifthe firstfunction( g)isaclosecopyofthesecond( h)butlagsitintimeby t,i.e.,ifthe first functionis shifted to the right of the second. Likewise, the correlationwill be large forsomenegativevalueof tifthefirstfunction leadsthesecond,i.e.,isshiftedtothe leftofthesecond. Therelationthatholdswhenthetwofunctionsareinterchangedis Corr (g, h )(t)=Corr (h, g )(−t)( 13.2.1 ) The discrete correlation of two sampled functions g kand hk, each periodic with period N,i sd efined by Corr (g, h )j≡N−1/summationdisplay k=0gj+khk (13.2.2 ) Thediscrete correlation theorem says that this discrete correlation of two real functions gand his one member of the discrete Fourier transform pair Corr (g, h )j⇐⇒ GkHk*( 13.2.3 ) where Gkand Hkare the discrete Fourier transforms of gjand hj, and the asterisk denotescomplexconjugation. Thistheoremmakesthesamepresumptionsaboutthe functions as those encountered for the discrete convolution theorem. We can compute correlations using the FFT as follows: FFT the two data sets, multiply one resulting transformby the complexconjugateof the other, and inverse transform the product. The result (call it rk) will formally be a complex vector of length N. However, it will turn out to have all its imaginary parts zero since the original data sets were both real. The components of rkare the values of the correlation at different lags, with positive and negative lags stored in the by nowfamiliarwrap-aroundorder: Thecorrelationat zerolagisin r 0,thefirst component; the correlation at lag 1 is in r1, the second component; the correlation at lag −1 is in rN−1, the last component; etc. Just as in the case of convolution we have to consider end effects, since our data will not, in general, be periodic as intended by the correlation theorem. Here again, we can use zero padding. If you are interested in the correlation for lags as large as ±K, then you must append a buffer zone of Kzeros at the end of both