trejo@nprdc.arpa (Leonard J. Trejo) (08/22/89)
We are experiencing problems with the subroutine "correl()" on page
433 of "Numerical Recipes in C" by Press et al. Specifically, we find
that the correlation function returned by a program using this
subroutine and the supporting subroutines (vector(), twofft(),
realft(), free_vector()) is incorrect. For example, the correlation of
two identical cosine bursts is asymmetrical.
If other have detected problems with this code, please let us know
about it by e-mail.
I'm attaching a copy of the program we wrote to this message.
It attempts to correlate two functions in ascii files "mask" and
"wave" and puts the result in "ans".
L. J. T.
===============================================================================
/* Maria 8/10/89 */
/* This program takes bogus files <data1> and <data2>, */
/* cross-correlates them, and stores the corr. values */
/* in bogus file <ans>. */
#include <stdio.h>
#include <fcntl.h>
#include <math.h>
/* for function four1.... */
#define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr
FILE *ptrdata1, *ptrdata2, *ptrout;
/* file pointers to <data1>, <data2>, and <ans> */
/* ##### BEGIN function free_vector ##### */
/* Frees a float vector allocated by vector(). */
void free_vector(v,nl,nh)
float *v;
int nl,nh;
{
free((char*) (v+nl));
}
/* ##### END function free_vector ##### */
/* ##### BEGIN function four1 ##### */
/* Replaces data by its discrete Fourier transform, */
/* if isign is input as 1; or replaces data by nn */
/* times its inverse discrete Fourier transform, if */
/* isign is input as -1. data is a complex array of */
/* length nn, input as a real array data[1..2*nn]. */
/* nn MUST be an integer power of 2 (this is not */
/* checked for!!!!!). */
void four1(data,nn,isign)
float data[];
int nn,isign;
{
int n,mmax,i,m,j,istep;
float tempr,tempi;
double wtemp,wr,wpr,wpi,wi,theta;
/* double precision for the trigonometric recurrences. */
n = nn << 1;
j = 1;
/* This is the bit-reversal section of the routine. */
for (i = 1; i < n; i += 2) {
if (j>i) {
/* Exchange the two complex numbers. */
SWAP(data[j],data[i]);
SWAP(data[j+1],data[i+1]);
}
m = n >> 1;
while (m >= 2 && j > m) {
j -= m;
m >>= 1;
}
j += m;
}
/* Here begins the Danielson-Lanczos section of the routine. */
mmax = 2;
while (n > mmax) {
/* Outer loop executed log(base 2) nn times. */
istep = 2*mmax;
/* Initialize for the trigonometric recurrence. */
theta = 6.28318530717959/(isign*mmax);
wtemp = sin(0.5*theta);
wpr = -2.0*wtemp*wtemp;
wpi = sin(theta);
wr = 1.0;
wi = 0.0;
/* Here are the two nested inner loops. */
for (m = 1; m < mmax; m += 2) {
for (i = m; i <= n; i += istep) {
/* This is the Danielson_Lanczos: */
j = i+mmax;
tempr = wr*data[j]-wi*data[j+1];
tempi = wr*data[j+1]+wi*data[j];
data[j] = data[i]-tempr;
data[j+1] = data[i+1]-tempi;
data[i] += tempr;
data[i+1] += tempi;
}
/* Trigonometric recurrence. */
wr = (wtemp = wr)*wpr-wi*wpi+wr;
wi = wi*wpr+wtemp*wpi+wi;
}
mmax = istep;
}
}
/* ##### END function four1 ##### */
/* ##### BEGIN function realft ##### */
/* This procedure calculates the Fourier Transform */
/* of a set of 2n real-valued data points. */
void realft(data,n,isign)
float data[];
int n, isign;
{
int i,i1,i2,i3,i4,n2p3;
float c1=0.5,c2,h1r,h1i,h2r,h2i;
double wr,wi,wpr,wpi,wtemp,theta;
/* double precision for the trigonometric recurrences */
void four1();
/* initialize the recurrence */
theta=3.141592653589793/(double) n;
if (isign == 1) {
/* the forward transform is here */
c2 = -0.5;
four1(data,n,1);
} else {
/* otherwise set up for an inverse transform */
c2 = 0.5;
theta = -theta;
}
wtemp = sin(0.5 * theta);
wpr = -2.0 * wtemp * wtemp;
wpi = sin(theta);
wr = 1.0 + wpr;
wi = wpi;
n2p3 = 2 * n + 3;
for ( i = 2; i <= n/2; i++) {
/* case i=1 done separately below */
i4 = 1+(i3 = n2p3-(i2 = 1+(i1 = i+i-1)));
/* The two separate transforms are */
/* separated out of data. */
h1r = c1*(data[i1]+data[i3]);
h1i = c1*(data[i2]-data[i4]);
h2r = -c2*(data[i2]+data[i4]);
h2i = c2*(data[i1]-data[i3]);
/* Here they are recombined to form the true */
/* transform of the original real data. */
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;
/* The recurrence. */
wr = (wtemp = wr)*wpr-wi*wpi+wr;
wi = wi*wpr+wtemp*wpi+wi;
}
if (isign == 1) {
/* squeeze the first and last data together */
/* to get them all in the original array. */
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]);
/* This is the inverse transform for the case isign=-1 */
four1(data,n,-1);
}
}
/* ##### END function realft ##### */
/* ##### BEGIN function twofft ##### */
/* Given two real input arrays data[1..n] and */
/* data2[1..n], this routine calls four1 and returns */
/* two complex output arrays, fft1 and fft2, each of */
/* complex length n (ie. real dimensions [1..2n], which*/
/* contain the discrete Fourier transforms of the */
/* respective datas. n MUST be an integer power of 2. */
void twofft(data1,data2,fft1,fft2,n)
float data1[],data2[],fft1[],fft2[];
int n;
{
int nn3,nn2,jj,j;
float rep,rem,aip,aim;
void four1();
nn3 = 1+(nn2 = 2+n+n);
for ( j = 1, jj = 2; j <= n; j++, jj += 2) {
/* Pack the two real arrays into one complex array. */
fft1[jj-1] = data1[j];
fft1[jj] = data2[j];
}
/* Transform the complex array. */
four1(fft1,n,1);
fft2[1] = fft1[2];
fft1[2] = fft2[2] = 0.0;
for (j = 3; j <= n+1; j += 2) {
/* Use symmetries to separate the two transforms. */
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]);
/* Ship them out in the complex arrays. */
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;
}
}
/* ##### END function twofft ##### */
/* ##### BEGIN function nrerror ##### */
/* The procedure is the numerical recipes standard error handler */
void nrerror(error_text)
char error_text[];
{
void exit();
fprintf(stderr,"Numerical Recipes run-time error...\n");
fprintf(stderr,"%s\n",error_text);
fprintf(stderr,"...now exiting to system...\n");
exit(1);
}
/* ##### END function nrerror ##### */
/* ##### BEGIN function *vector ##### */
/* This procedure allocates a float vector with range */
/* [nl..nh] */
float *vector(nl,nh)
int nl,nh;
{
float *v;
v = (float *)malloc((unsigned) (nh-nl+1)*sizeof(float));
if (!v) nrerror("allocation failure in vector()");
return v-nl;
}
/* ##### END function *vector ##### */
/* ##### BEGIN function correl ##### */
/* This procedure computes the correlation of the two */
/* real data sets data1[1..n] and data2[1..n] each of */
/* length n (including any user supplied padding). */
/* ******* n MUST be a power of 2!! ******* */
/* The answer is returned as the first n points in */
/* ans[1..2*n] stored in wraparound order, ie. */
/* correlations at increasingly negative lags are in */
/* ans[n] on down to ans[n/2+1], while correlations */
/* increasingly positive lags are in and[1] on up to */
/* ans[n/2]. */
void correl(data1,data2,n,ans)
float data1[],data2[],ans[];
int n;
{
int no2, i;
float dum, *fft, *vector();
void twofft(), realft(), free_vector();
fft = vector(1,2*n);
/* transforms both data vectors at one. */
twofft(data1,data2,fft,ans,n);
/* normalization for inverse FFT */
no2 = n/2;
for (i = 2; i <= n + 2; i += 2) {
/* multiply to find FFT of their correlation */
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;
}
/* pack first and last into one element */
ans[2] = ans[n+1];
/* inverse transform gives correlation */
realft(ans,no2,-1);
free_vector(fft,1,2*n);
}
/* ##### END function correl ##### */
/* ##### BEGIN main body ##### */
main (argc, argv)
int argc;
char *argv[];
{
float data1[1024],data2[1024],datat[1024],ans[2048],hold=0.0,avg;
int n,i=0,j,k,bot,top,c1=0,c2=0,bigger,twopwr=2,upper,lower;
char mask[10], *wave, *out, command[100], c, npts;
/* check command line for correct # of args, then */
/* identify the three files: mask, wave and out. */
if (argc != 7) {
printf("Command line has incorrect number of entries!\n");
exit();
}
while (--argc > 0 && (*++argv)[0]== '-') {
while (c = *++argv[0]) {
switch(c) {
case 't':
strcpy(mask, argv[1]);
break;
case 'n':
npts = atoi(argv[1]);
sprintf(command,"series 0 %d | dm \"sin(x1*6.283/%d)\" > mask",npts,2*npts);
system(command);
strcpy(mask,"mask");
break;
case 'd':
wave = argv[1];
break;
case 'r':
out = argv[1];
break;
default:
printf("This is an illegal command: %c\n",c);
exit();
}
}
++argv;
--argc;
}
/* open files wave and mask and check to make sure */
/* that they are not NULL files. If they are, exit */
/* this program and report NULL file to screen. */
if ((ptrdata1 = fopen(wave,"r")) == NULL) {
printf("results file can't be opened!\n");
exit();
}
if ((ptrdata2 = fopen(mask,"r")) == NULL) {
printf("mask file can't be opened!\n");
exit();
}
if ((ptrout = fopen(out,"w")) == NULL) {
printf("results file can't be opened!\n");
exit();
}
/* Fill array data1 from file wave. */
while ( fscanf(ptrdata1,"%f", &data1[c1]) != 0 && !feof(ptrdata1)) {
c1 ++;
}
c1 --;
/* Fill array data2 from file mask. */
while ( fscanf(ptrdata2,"%f", &data2[c2]) != 0 && !feof(ptrdata2)) {
c2 ++;
}
c2 --;
/* Normalize mask. */
for ( i = 0; i <= c2; i++) {
hold = hold + data2[i];
}
avg = hold/(c2+1);
hold=0;
for ( i = 0; i <= c2; i++) {
datat[i] = data2[i] - avg;
hold += datat[i] * data2[i];
}
for (i=0; i <= c2; i++) {
data2[i]=datat[i]/hold;
}
/* Find out which array is larger. */
if ( c1 >= c2 )
bigger = c1+1;
else {
printf("Mask file length cannot exceed ");
printf("data file length!!!\n");
exit();
}
/* Compare the larger array with a power of two */
/* until the smallest power of two that fits the */
/* 'biggest' array size is found. */
/* Assign that power of two (- 1) to var n. */
while ( bigger > twopwr)
twopwr = twopwr * 2;
/* Double array sizes, so convolution is not */
/* circular...but linear. */
n = twopwr * 2;
correl(data1,data2,n,ans);
/* Move first half of mask to end portion of new */
/* mask to remove delay from ans. */
/* Check to see if mask size is even or odd... */
/* if (c2 % 2 == 0) { */
/* odd */
/* upper = c2/2;
lower = upper - 1;
for (i=0; i<n; i++)
datat[i]=0;
k = n - upper;
for (i=0; i<=upper; i++) {
j=i+(c2/2);
datat[i]=data2[j];
}
for (i=0; i<=lower; i++) {
datat[k]=data2[i];
k++;
}
}
else { */
/* even */
/* upper = c2/2 + 1;
lower = upper - 1;
for (i=0; i<n; i++)
datat[i]=0;
k = n - upper;
j = 0;
for (i=upper; i<=c2; i++) {
datat[j]=data2[i];
j++;
}
for (i=0; i<=lower; i++) {
datat[k]=data2[i];
k++;
}
}
for (i=0; i<n; i++)
data2[i]=datat[i];*/
for (i = 1; i < n; i++) {
/* write array ans to file <ans> */
fprintf(ptrout,"%f\n", ans[i]);
}
close(out);
}
============================================================================
ARPANET : trejo@nprdc.navy.mil UUCP: ucsd!nprdc!trejo
U.S. Mail: Leonard J. Trejo, Ph. D. Phone: (619) 553-7711
Neuroscience Laboratory (AV) 553-7711
NPRDC, Code 141
San Diego, CA 92152-6800