c13-1 (779564), страница 2
Текст из файла (страница 2)
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).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 rk are 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.54113.1 Convolution and Deconvolution Using the FFTresponse functionm−m+m−convolutionspoiledunspoiledspoiledFigure 13.1.3.
The wrap-around problem in convolving finite segments of a function. Not only mustthe 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 theresponse function.response functionm+m−original functionzero paddingm−m+not spoiled because zerom+m−unspoiledspoiledbut irrelevantFigure 13.1.4. Zero padding as solution to the wrap-around problem. The original function is extendedby 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.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).m+sample of original function542Chapter 13.Fourier and Spectral Applicationspadding of the response rk described above, we effectively insulate the data fromartifacts of undesired periodicity.
Figure 13.1.4 illustrates matters.Use of FFT for ConvolutionSample 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).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 N andarranged in wrap-around order. (Generally this means that a large contiguous sectionof the rk ’s, in the middle of that array, is zero, with nonzero values clustered atthe two extreme ends of the array.) You now compute the discrete convolution asfollows: Use the FFT algorithm to compute the discrete Fourier transform of s andof r.
Multiply the two transforms together component by component, rememberingthat the transforms consist of complex numbers. Then use the FFT algorithm totake the inverse discrete Fourier transform of the products. The answer is theconvolution r ∗ s.What about deconvolution? Deconvolution is the process of undoing thesmearing in a data set that has occurred under the influence of a known responsefunction, for example, because of the known effect of a less-than-perfect measuringapparatus. The defining equation of deconvolution is the same as that for convolution,namely (13.1.1), except now the left-hand side is taken to be known, and (13.1.1) isto be considered as a set of N linear equations for the unknown quantities sj .
Solvingthese simultaneous linear equations in the time domain of (13.1.1) is unrealistic inmost cases, but the FFT renders the problem almost trivial. Instead of multiplyingthe transform of the signal and response to get the transform of the convolution, wejust divide the transform of the (known) convolution by the transform of the responseto get the transform of the deconvolved signal.This procedure can go wrong mathematically if the transform of the responsefunction is exactly zero for some value Rn , so that we can’t divide by it. Thisindicates 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 processof deconvolution has other practical shortcomings.
The process is generally quitesensitive to noise in the input data, and to the accuracy to which the response functionrk is known. Perfectly reasonable attempts at deconvolution can sometimes producenonsense for these reasons. In such cases you may want to make use of the additionalprocess of optimal filtering, which is discussed in §13.3.Here is our routine for convolution and deconvolution, using the FFT asimplemented in four1 of §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 realft should be substituted if data and respnshave very different magnitudes, to minimize roundoff. The data are assumed to bestored in a float array data[1..n], with n an integer power of two. The responsefunction is assumed to be stored in wrap-around order in a sub-array respns[1..m]of the array respns[1..n]. The value of m can be any odd integer less than or equalto n, since the first thing the program does is to recopy the response function into theappropriate wrap-around order in respns[1..n].
The answer is provided in ans.13.1 Convolution and Deconvolution Using the FFT543#include "nrutil.h"fft=vector(1,n<<1);for (i=1;i<=(m-1)/2;i++)Put respns in array of length n.respns[n+1-i]=respns[m+1-i];for (i=(m+3)/2;i<=n-(m-1)/2;i++)Pad with zeros.respns[i]=0.0;twofft(data,respns,fft,ans,n);FFT both at once.no2=n>>1;for (i=2;i<=n+2;i+=2) {if (isign == 1) {ans[i-1]=(fft[i-1]*(dum=ans[i-1])-fft[i]*ans[i])/no2;Multiply FFTsans[i]=(fft[i]*dum+fft[i-1]*ans[i])/no2;to convolve.} else if (isign == -1) {if ((mag2=SQR(ans[i-1])+SQR(ans[i])) == 0.0)nrerror("Deconvolving at response zero in convlv");ans[i-1]=(fft[i-1]*(dum=ans[i-1])+fft[i]*ans[i])/mag2/no2;Divide FFTsans[i]=(fft[i]*dum-fft[i-1]*ans[i])/mag2/no2;to deconvolve.} else nrerror("No meaning for isign in convlv");}ans[2]=ans[n+1];Pack last element with first for realft.realft(ans,n,-1);Inverse transform back to time domain.free_vector(fft,1,n<<1);}Convolving or Deconvolving Very Large Data SetsIf your data set is so long that you do not want to fit it into memory all atonce, 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 worrynot only about spurious wrap-around effects, but also about the fact that the ends ofeach section of data should have been influenced by data at the nearby ends of theimmediately preceding and following sections of data, but were not so influencedsince only one section of data is in the machine at a time.There are two, related, standard solutions to this problem.
Both are fairlyobvious, so with a few words of description here, you ought to be able to implementthem for yourself. The first solution is called the overlap-save method. In thistechnique you pad only the very beginning of the data with enough zeros to avoidwrap-around pollution. After this initial padding, you forget about zero paddingSample 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).void convlv(float data[], unsigned long n, float respns[], unsigned long m,int isign, float ans[])Convolves or deconvolves a real data set data[1..n] (including any user-supplied zero padding)with a response function respns[1..n]. The response function must be stored in wrap-aroundorder in the first m elements of respns, where m is an odd integer ≤ n. Wrap-around ordermeans that the first half of the array respns contains the impulse response function at positivetimes, while the second half of the array contains the impulse response function at negative times,counting down from the highest element respns[m].
On input isign is +1 for convolution,−1 for deconvolution. The answer is returned in the first n components of ans. However,ans must be supplied in the calling program with dimensions [1..2*n], for consistency withtwofft. n MUST be an integer power of two.{void realft(float data[], unsigned long n, int isign);void twofft(float data1[], float data2[], float fft1[], float fft2[],unsigned long n);unsigned long i,no2;float dum,mag2,*fft;544Chapter 13.a0baAFourier and Spectral Applicationscdata (in)00b000cCAA+BBB+CCconvolution (out)Figure 13.1.5.The overlap-add method for convolving a response with a very long signal. Thesignal data is broken up into smaller pieces. Each is zero padded at both ends and convolved (denotedby bold arrows in the figure).
Finally the pieces are added back together, including the overlappingregions formed by the zero pads.altogether. Bring in a section of data and convolve or deconvolve it. Then throwout the points at each end that are polluted by wrap-around end effects. Output onlythe remaining good points in the middle. Now bring in the next section of data, butnot all new data. The first points in each new section overlap the last points fromthe preceding section of data.
The sections must be overlapped sufficiently so thatthe polluted output points at the end of one section are recomputed as the first of theunpolluted output points from the subsequent section. With a bit of thought you caneasily determine how many points to overlap and save.The second solution, called the overlap-add method, is illustrated in Figure13.1.5. Here you don’t overlap the input data. Each section of data is disjointfrom the others and is used exactly once. However, you carefully zero-pad it atboth ends so that there is no wrap-around ambiguity in the output convolution ordeconvolution.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
