#include <math.h>
#include "nrutil.h"
#include "four1.c"

/* correl.c								*/
/* Cross-correlation							*/
/* From W. H.Press, S. A. Teukolsky, W. T. Vettering and B. P. Flannery */
/* (1992) Numerical Recipes in C: The Art of Scientific Computing	*/
/* (2nd edition). Cambridge University Press.				*/

/* Copyright (c) Numerical Recipes Software 1987, 1988, 1992,1997	*/

/* p. 546								*/

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-specified 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);
	no2=n>>1;
	for (i=2;i<=n+2;i+=2) {
		ans[i-1]=(fft[i-1]*(dum=ans[i-1])+fft[i]*ans[i])/no2;
		ans[i]=(fft[i]*dum-fft[i-1]*ans[i])/no2;
	}
	ans[2]=ans[n+1];
	realft(ans,n,-1);
	free_vector(fft,1,n<<1);
}

/* FFT of two real functions simultaneously				*/
/* Press, Teukolsky, Vettering and Flannery (1992: 511)			*/
/* Copyright (c) Numerical Recipes Software 1987, 1988, 1992,1997	*/

void twofft(float data1[], float data2[], float fft1[], float fft2[], unsigned long n)
/* Given two real input arrays data1[1..n] and data2[1..n], this routine calls four1
and returns two complex output arrays, fft1[1..2n] and fft2[1..2n], each of complex
length n (i.e. real length 2*n), which contains the discrete Fourier transforms of
the respective data arrays. n MUST be an integer power of 2. */
{
	void four1(float data[], unsigned long nn, int isign);
	unsigned long nn3, nn2, jj, j;
	float rep, rem, aip, aim;

	nn3=1+(nn2=2+n+n);
	for (j=1,jj=2;j<=n;j++,jj+=2) {
		fft1[jj-1]=data1[j];
		fft1[jj]=data2[j];
	}
	four1(fft1,n,1);	
	fft2[1]=fft1[2];
	fft1[2]=fft2[2]=0.0;
	for (j=3;j<=n+1;j+=2) {
		rep=0.5*(fft1[j]+fft1[nn2-j]);
		rem=0.5*(fft1[j]-fft1[nn2-j]);
		aip=0.5*(fft1[j+1]+fft1[nn3-j]);
		aim=0.5*(fft1[j+1]-fft1[nn3-j]);
		fft1[j]=rep;
		fft1[j+1]=aim;
		fft1[nn2-j]=rep;
		fft1[nn3-j] = -aim;
		fft2[j]=aip;
		fft2[j+1] = -rem;
		fft2[nn2-j]=aip;
		fft2[nn3-j]=rem;
	}
}

/* FFT of single real function						*/
/* Press, Teukolsky, Vettering and Flannery (1992: 513-4)		*/
/* Copyright (c) Numerical Recipes Software 1987, 1988, 1992,1997	*/

void realft(float data[], unsigned long n, int isign)
/* Calculates the Fourier transform of a set of n real-valued data points.
Replaces this data (which is stored in array data[1..n]) by the positive
frequency half of its complex Fourier transform. The real-valued first and
last components of the complex transform are returned as elements data[1] and
data[2], respectively. n must be a power of 2. This routine also calculates
the inverse transform of a complex data array if it is the transform of real
data. (Result in this case must be multiplied by 2/n.) */
{
	void four1(float data[], unsigned long nn, int isign);
	unsigned long i, i1, i2, i3, i4, np3;
	float c1=0.5, c2, h1r, h1i, h2r, h2i;
	double wr, wi, wpr, wpi, wtemp, theta;

	theta=3.141592653589793/(double) (n>>1);
	if (isign == 1) {
		c2 = -0.5;
		four1(data,n>>1,1);
	} else {
		c2 = 0.5;
		theta = -theta;
	}
	wtemp=sin(0.5*theta);
	wpr = -2.0*wtemp*wtemp;
	wpi=sin(theta);
	wr=1.0+wpr;
	wi=wpi;
	np3=n+3;
	for (i=2;i<=(n>>2);i++) {
		i4=1+(i3=np3-(i2=1+(i1=i+i-1)));
		h1r=c1*(data[i1]+data[i3]);
		h1i=c1*(data[i2]-data[i4]);
		h2r = -c2*(data[i2]+data[i4]);
		h2i=c2*(data[i1]-data[i3]);
		data[i1]=h1r+wr*h2r-wi*h2i;
		data[i2]=h1i+wr*h2i+wi*h2r;
		data[i3]=h1r-wr*h2r+wi*h2i;
		data[i4] = -h1i+wr*h2i+wi*h2r;
		wr=(wtemp=wr)*wpr-wi*wpi+wr;
		wi=wi*wpr+wtemp*wpi+wi;
	}
	if (isign == 1) {
		data[1] = (h1r=data[1])+data[2];
		data[2] = h1r-data[2];
	} else {
		data[1]=c1*((h1r=data[1])+data[2]);
		data[2]=c1*(h1r-data[2]);
		four1(data,n>>1,-1);
	}
}