Numerical Recipes in C: problems with cross-correlation code
Leonard J. Trejo
trejo at nprdc.arpa
Wed Aug 23 01:56:12 AEST 1989
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 at 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
More information about the Comp.lang.c
mailing list