c13-2 (Numerical Recipes in C)

54513.2 Correlation and Autocorrelation Using the FFTElliott, D.F., and Rao, K.R. 1982, Fast Transforms: Algorithms, Analyses, Applications (NewYork: Academic Press).Brigham, E.O. 1974, The Fast Fourier Transform (Englewood Cliffs, NJ: Prentice-Hall), Chapter 13.Correlation is the close mathematical cousin of convolution. It is in someways simpler, however, because the two functions that go into a correlation are notas conceptually distinct as were the data and response functions that entered intoconvolution.

Rather, in correlation, the functions are represented by different, butgenerally similar, data sets. We investigate their “correlation,” by comparing themboth directly superposed, and with one of them shifted left or right.We have already defined in equation (12.0.10) the correlation between twocontinuous functions g(t) and h(t), which is denoted Corr(g, h), and is a functionof lag t. We will occasionally show this time dependence explicitly, with the ratherawkward notation Corr(g, h)(t). The correlation will be large at some value oft if the first function (g) is a close copy of the second (h) but lags it in time byt, i.e., if the first function is shifted to the right of the second. Likewise, thecorrelation will be large for some negative value of t if the first function leads thesecond, i.e., is shifted to the left of the second.

The relation that holds when thetwo functions are interchanged isCorr(g, h)(t) = Corr(h, g)(−t)(13.2.1)The discrete correlation of two sampled functions gk and hk , each periodicwith period N , is defined byCorr(g, h)j ≡N−1Xgj+k hk(13.2.2)k=0The discrete correlation theorem says that this discrete correlation of two realfunctions g and h is one member of the discrete Fourier transform pairCorr(g, h)j ⇐⇒ Gk Hk *(13.2.3)where Gk and Hk are the discrete Fourier transforms of gj and hj , and the asteriskdenotes complex conjugation.

This theorem makes the same presumptions about thefunctions 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 transform by the complex conjugate of the other, and inversetransform the product. The result (call it rk ) will formally be a complex vectorof length N .

However, it will turn out to have all its imaginary parts zero sincethe original data sets were both real. The components of rk are the values of 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.

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).13.2 Correlation and Autocorrelation Usingthe FFT546Chapter 13.Fourier and Spectral Applications#include "nrutil.h"void correl(float data1[], float data2[], unsigned long n, float ans[])Computes the correlation of two real data sets data1[1..n] and data2[1..n] (including anyuser-supplied zero padding).

n MUST be an integer power of two. The answer is returned asthe first n points in ans[1..2*n] stored in wrap-around order, i.e., correlations at increasinglynegative lags are in ans[n] on down to ans[n/2+1], while correlations at increasingly positivelags are in ans[1] (zero lag) on up to ans[n/2]. Note that ans must be supplied in the callingprogram with length at least 2*n, since it is also used as working space.

Sign convention ofthis routine: if data1 lags data2, i.e., is shifted to the right of it, then ans will show a peakat positive lags.{void realft(float data[], unsigned long n, int isign);void twofft(float data1[], float data2[], float fft1[], float fft2[],unsigned long n);unsigned long no2,i;float dum,*fft;fft=vector(1,n<<1);twofft(data1,data2,fft,ans,n);Transform both data vectors at once.no2=n>>1;Normalization for inverse FFT.for (i=2;i<=n+2;i+=2) {ans[i-1]=(fft[i-1]*(dum=ans[i-1])+fft[i]*ans[i])/no2;Multiply to findans[i]=(fft[i]*dum-fft[i-1]*ans[i])/no2;FFT of their cor}relation.ans[2]=ans[n+1];Pack first and last into one element.realft(ans,n,-1);Inverse transform gives correlation.free_vector(fft,1,n<<1);}As in convlv, it would be better to substitute two calls to realft for the onecall to twofft, if data1 and data2 have very different magnitudes, to minimizeroundoff error.The discrete autocorrelation of a sampled function gj is just the discretecorrelation of the function with itself.

Obviously this is always symmetric withrespect to positive and negative lags. Feel free to use the above routine correlto obtain autocorrelations, simply calling it with the same data vector in botharguments. If the inefficiency bothers you, routine realft can, of course, be usedto transform the data vector instead.CITED REFERENCES AND FURTHER READING:Brigham, E.O. 1974, The Fast Fourier Transform (Englewood Cliffs, NJ: Prentice-Hall), §13–2.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).correlation at different lags, with positive and negative lags stored in the by nowfamiliar wrap-around order: The correlation at zero lag is in r0 , the first component;the correlation at lag 1 is in r1 , the second component; the correlation at lag −1is in rN−1 , the last component; etc.Just as in the case of convolution we have to consider end effects, since ourdata will not, in general, be periodic as intended by the correlation theorem.

Hereagain, we can use zero padding. If you are interested in the correlation for lags aslarge as ±K, then you must append a buffer zone of K zeros at the end of bothinput data sets. If you want all possible lags from N data points (not a usual thing),then you will need to pad the data with an equal number of zeros; this is the extremecase.

So here is the program:54713.3 Optimal (Wiener) Filtering with the FFT13.3 Optimal (Wiener) Filtering with the FFTZ∞s(t) =−∞r(t − τ )u(τ ) dτor S(f) = R(f)U (f)(13.3.1)where S, R, U are the Fourier transforms of s, r, u, respectively. Second, themeasured signal c(t) may contain an additional component of noise n(t),c(t) = s(t) + n(t)(13.3.2)We already know how to deconvolve the effects of the response function r inthe absence of any noise (§13.1); we just divide C(f) by R(f) to get a deconvolvedsignal.

We now want to treat the analogous problem when noise is present. Ourtask is to find the optimal filter, φ(t) or Φ(f), which, when applied to the measuredsignal c(t) or C(f), and then deconvolved by r(t) or R(f), produces a signal ue(t)eor U(f) that is as close as possible to the uncorrupted signal u(t) or U (f). In otherwords we will estimate the true signal U bye (f) = C(f)Φ(f)UR(f)e to be close to U ?In what sense is Uleast-square senseZ∞−∞Z2|eu(t) − u(t)| dt =∞−∞(13.3.3)We ask that they be close in the2eU(f) − U (f) dfis minimized.(13.3.4)Substituting equations (13.3.3) and (13.3.2), the right-hand side of (13.3.4) becomes2 [S(f) + N (f)]Φ(f)S(f) −dfR(f)R(f) −∞Z ∞no−22222|S(f)| |1 − Φ(f)| + |N (f)| |Φ(f)|df|R(f)|=Z∞(13.3.5)−∞The signal S and the noise N are uncorrelated, so their cross product, whenintegrated over frequency f, gave zero. (This is practically the definition of what wemean by noise!) Obviously (13.3.5) will be a minimum if and only if the integrandSample 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).There are a number of other tasks in numerical processing that are routinelyhandled with Fourier techniques.

One of these is filtering for the removal of noisefrom a “corrupted” signal. The particular situation we consider is this: There is someunderlying, uncorrupted signal u(t) that we want to measure. The measurementprocess is imperfect, however, and what comes out of our measurement device is acorrupted signal c(t). The signal c(t) may be less than perfect in either or both oftwo respects.

First, the apparatus may not have a perfect “delta-function” response,so that the true signal u(t) is convolved with (smeared out by) some known responsefunction r(t) to give a smeared signal s(t),.