c13-1 (Numerical Recipes in C)

538Chapter 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 :Figure 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.N−10N−1Figure 13.1.2. Convolution of discretely sampled functions. Note how the response function for negativetimes is wrapped around and stored at the extreme right end of the array rk .Sample 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).0(r* s)jN−1rk0sjtr* s(t)tr (t)ts(t)53913.1 Convolution and Deconvolution Using the FFT540Chapter 13.Fourier and Spectral Applicationsdiscrete convolution with a response function of finite duration N is a member ofthe discrete Fourier transform pair,XN/2sj−k rk⇐⇒Sn Rn(13.1.2)k=−N/2+1Treatment of End Effects by Zero PaddingThe discrete convolution theorem presumes a set of two circumstances thatare not universal.

First, it assumes that the input signal is periodic, whereas realdata often either go forever without repetition or else consist of one nonperiodicstretch of finite length. Second, the convolution theorem takes the duration of theresponse to be the same as the period of the data; they are both N . We need towork around these two constraints.The second is very straightforward. Almost always, one is interested in aresponse function whose duration M is much shorter than the length of the dataset N .

In this case, you simply extend the response function to length N bypadding it with zeros, i.e., define rk = 0 for M/2 ≤ k ≤ N/2 and also for−N/2 + 1 ≤ k ≤ −M/2 + 1. Dealing with the first constraint is more challenging.Since the convolution theorem rashly assumes that the data are periodic, it willfalsely “pollute” the first output channel (r ∗ s)0 with some wrapped-around datafrom 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 sj vector,in order to make this pollution zero. How many zero values do we need in thisbuffer? Exactly as many as the most negative index for which the response functionis nonzero. For example, if r−3 is nonzero, while r−4 , r−5 , . . . are all zero, then weneed three zero pads at the end of the data: sN−3 = sN−2 = sN−1 = 0. Thesezeros will protect the first output channel (r ∗ s)0 from wrap-around pollution. Itshould be obvious that the second output channel (r ∗ s)1 and subsequent ones willalso be protected by these same zeros.

Let K denote the number of padding zeros,so that the last actual input data point is sN−K−1 .What now about pollution of the very last output channel? Since the datanow end with sN−K−1 , the last output channel of interest is (r ∗ s)N−K−1 . Thischannel can be polluted by wrap-around from input channel s0 unless the numberK is also large enough to take care of the most positive index k for which theresponse function rk is nonzero.

For example, if r0 through r6 are nonzero, whiler7 , r8 . . . are all zero, then we need at least K = 6 padding zeros at the end of thedata: sN−6 = . . . = sN−1 = 0.To summarize — we need to pad the data with a number of zeros on oneend equal to the maximum positive duration or maximum negative duration ofthe response function, whichever is larger. (For a symmetric response function ofduration M , you will need only M/2 zero pads.) Combining this operation with theSample 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.