f13-2

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 defined 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.,ifthefirst 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, is defined 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 13.3Optimal(Wiener)FilteringwiththeFFT 539Sample 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).input data sets. If youwant all possible lags from Ndata points (nota usual thing), thenyouwillneedtopadthedatawithanequalnumberofzeros;thisis theextremecase. So here is the program: SUBROUTINE correl(data1,data2,n,ans) INTEGER n,NMAX REAL data1(n),data2(n)COMPLEX ans(n) PARAMETER (NMAX=4096) Maximum anticipated FFT size. C USES realft,twofft Computes the correlation of two real data sets data1(1:n) anddata2(1:n) (includ- ing any user-supplied zero padding). nMUST be an integer power of two. The answer is returned as the first npoints in ans stored in wrap-around order, i.e., correlations at increasingly negative lags are in ans(n) on down to ans(n/2+1) , while correlations at increasingly positive lags are in ans(1) ( z e r ol a g )o nu pt o ans(n/2) .N o t et h a t ans must be supplied in the calling program with length at least 2*n , since it is also used as working space. Sign convention of this routine: if data1 lagsdata2 , i.e., is shifted to the right of it, then ans will show a peak at positive lags. INTEGER i,no2 COMPLEX fft(NMAX) call twofft(data1,data2,fft,ans,n) Transform b oth data vectors at once. no2=n/2 Normalization for inverse FFT. do11i=1,no2+1 ans(i)=fft(i)*conjg(ans(i))/float(no2) Multiply to find FFT of their corre- lation. enddo 11 ans(1)=cmplx(real(ans(1)),real(ans(no2+1))) Pack first and last into one element. call realft(ans,n,-1) Inverse transform gives correlation. returnEND As in convlv, it would be better to substitute two calls to realftfor the one call to twofft,i fdata1anddata2have very different magnitudes, to minimize roundoff error. Thediscrete autocorrelation of a sampled function gjis just the discrete correlation of the function with itself. Obviously this is always symmetric withrespect to positive and negative lags. Feel free to use the above routine correl to obtain autocorrelations, simply calling it with the same datavector in both arguments. If the inefficiency bothers you, routine realftcan, of course, be used to transform the datavector instead. CITED REFERENCES AND FURTHER READING: Brigham, E.O. 1974, TheFast FourierTransform (EnglewoodCliffs, NJ: Prentice-Hall), §13–2. 13.3 Optimal (Wiener) Filtering with the FFT There are a number of other tasks in numerical processing that are routinely handled with Fourier techniques. One of these is filtering for the removal of noise froma“corrupted”signal. Theparticularsituationweconsideristhis: Thereissomeunderlying, uncorrupted signal u(t)that we want to measure. The measurement process is imperfect, however, and what comes out of our measurement device is a corrupted signal c(t). The signal c(t)may be less than perfect in either or both of two respects. First, the apparatus may not have a perfect “delta-function”response,