/*
   Random Phase test

   find the significance of a correlation with
   serial correlation (i.e. autocorrelated)

   v1.0 n = power of 2  (yes, I used a power of 2 fft)
              July 2000, w. ebisuzaki, CPC/NCEP/NOAA

*/

/* number of random phase samples to generate the statistics */
#define NSAMPLE 50000

/* fftutil.c */
void do_fft(float *series, float *spectral, int n);
void undo_fft(float *spectral, float *series, int n);
double s_variance(float *spectral, int n);
double s_corr(float *spectral1, float *spectral2, int n);
double s_covar(float *spectral1, float *spectral2, int n);
double variance(float *x, int n);
double corr(float *x, float *y, int n);
void s_random_phase(float *a, float *ran, int n);

#include <stdio.h>
#include <stdlib.h>
#include <stddef.h>
#include <math.h>

int flt_cmp(const void *, const void *);


long int iseed = -22459;

void main(int argc, char **argv) {

    int n;
    float *x, *y, *s_x, *s_y, *s_ran;

    float ac_n[NSAMPLE];
    double sig, ac_0;
    int count, i;

    FILE *in1, *in2;

    if (argc <= 4) {
	fprintf(stderr,"usage %s sig n [time series 1] [time series 2]\n", 
 	    argv[0]);
	exit(8);
    }
    sig = 1.0 - atof(argv[1]);
    n = atoi(argv[2]);

    /* check if n is a power of 2 */
    i = 2;
    while (i < n) i = i + i;
    if (i != n) {
	fprintf(stderr, "n (%d) is not a power of 2\n", n);
    }

    in1 = fopen(argv[3],"r");    
    in2 = fopen(argv[4],"r");    

    x = (float *) malloc( n * sizeof(float));
    y = (float *) malloc( n * sizeof(float));
    s_x = (float *) malloc( n * sizeof(float));
    s_y = (float *) malloc( n * sizeof(float));
    s_ran = (float *) malloc( n * sizeof(float));

    if (x == NULL || y == NULL || s_x == NULL || s_x == NULL || s_ran == NULL) {
        fprintf(stderr,"lack of memory\n");
	exit(9);
    }

    printf("reading %s\n", argv[3]);
    for (i = 0; i < n; i++) {
	if (fscanf(in1, "%f", x+i) != 1) {
	    printf("error reading %s, element %d\n", argv[3], i);
	    exit(9);
	}
    }
    printf("reading %s\n", argv[4]);
    for (i = 0; i < n; i++) {
	if (fscanf(in2, "%f", y+i) != 1) {
	    printf("error reading %s, element %d\n", argv[4], i);
	    exit(9);
	}
    }
    printf("calculating\n");

    ac_0 = corr(x, y, n);
    ac_0 = fabs((double) ac_0);

    /* fft of x and y time series */

    do_fft(x, s_x, n);
    do_fft(y, s_y, n);

    count = 0;
    for (i = 0; i < NSAMPLE; i++) {
	s_random_phase(s_y, s_ran, n);
	ac_n[i] = fabs(s_corr(s_x, s_ran, n));
	if (ac_n[i] > ac_0) count++;
    }

    qsort(ac_n, (size_t) NSAMPLE, sizeof (float), flt_cmp);

    printf("testing %s and %s (%d points) for %f significance\n",
	argv[3], argv[4], n, 1.0-sig);

    printf("sample |corr| %f,  fraction of samples with larger |corr| %f\n",
	ac_0, (double) count / NSAMPLE);
    printf("critical |corr| %f at %f sig level, %d random samples used\n", 
	ac_n[(int) ((NSAMPLE-1)*sig)], 1.0-sig, NSAMPLE);
}

int flt_cmp(const void *a, const void *b) {
	float dif;

	dif = *((float *) a) -  *((float *) b);

	if (dif < 0.0) return -1;
	if (dif > 0.0) return +1;
	return 0;
}