546 Chapter 13.Fourier and Spectral Applications the correlation at lag I 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 K zeros at the end of both input 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 extreme case.So here is the program: #include "nrutil.h" void correl(float datai[],float data2[],unsigned long n,float ans]) Computes the correlation of two real data sets datal[1..n]and data2[1..n](including any user-supplied zero padding).n MUST be an integer power of two.The answer is returned as the first n points in ans [1..2*n]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](zero lag)on up to ans [n/2].Note that 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 datal lags data2,i.e.,is shifted to the right of it,then ans will show a peak at positive lags void realft(float data],unsigned long n,int isign); Press. void twofft(float datal[],float data2[],float ffti[], float fft2], unsig即ed long n); unsigned long no2,i; 9 float dum,*fft; fft=vector(1.n<<1); SCIENTIFIC twofft(datal,data2,fft,ans,n); Transform both data vectors at once. no2=n>>1; Normalization for inverse FFT. for(1=2;1<=n+2;1+=2){ ans [i-1]=(fft[i-1]*(dum=ans[i-1])+fft[i]*ans[i])/no2; Multiply to find ans [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 1920 realft(ans,n,-1); Inverse transform gives correlation. (ISBN free_vector(fft,1,n<<1); Numerical 色 43108 As in convlv,it would be better to substitute two calls to realft for the one uction call to twofft,if datal and data2 have very different magnitudes,to minimize Recipes roundoff error. (outside The discrete autocorrelation of a sampled function gi is just the discrete North Software. correlation of the function with itself.Obviously this is always symmetric with respect to positive and negative lags.Feel free to use the above routine correl to obtain autocorrelations,simply calling it with the same data vector in both arguments.If the inefficiency bothers you,routine realft can,of course,be used to 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.546 Chapter 13. Fourier and Spectral Applications Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machine￾readable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). 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 K zeros at the end of both input 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 extreme case. So here is the program: #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 any user-supplied zero padding). n MUST be an integer power of two. The answer is returned as the first n points in ans[1..2*n] 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] (zero lag) on up to ans[n/2]. Note that 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 lags data2, i.e., is shifted to the right of it, then ans will show a peak at 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 find FFT of their cor￾relation. ans[i]=(fft[i]*dum-fft[i-1]*ans[i])/no2; } 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 one call to twofft, if data1 and data2 have very different magnitudes, to minimize roundoff error. The discrete autocorrelation of a sampled function gj is just the discrete correlation of the function with itself. Obviously this is always symmetric with respect to positive and negative lags. Feel free to use the above routine correl to obtain autocorrelations, simply calling it with the same data vector in both arguments. If the inefficiency bothers you, routine realft can, of course, be used to 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