/*
   misc routines
*/

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include "nr.h"

#ifndef PI
#define PI 	3.14159265358979323846
#endif

#ifndef abs
#define abs(x)	((x) < 0 ? (-(x)) : (x))
#endif

#define Window(j) (1.0 - fabs((2.0*j - (n-1.0))/(n-1.0)))

#ifdef AMIGA
#include <m68881.h>
#endif



/*
 * do_fft:
 *          series -> spectral
 */
void do_fft(float *series, float *spectral, int n) {

	int i;
	double factor;

	for (i = 0; i < n; i++) {
		spectral[i] = series[i];
	}

	realft(spectral - 1, (unsigned long) n, 1);

	factor = 2.0 / n;

	for (i = 0; i < n; i++) {
		spectral[i] *= factor;
	}
}

/*
 * undo_fft:
 *           spectral -> series
 */
void undo_fft(float *spectral, float *series, int n) {

	int i;
	for (i = 0; i < n; i++) {
		series[i] = spectral[i];
	}
        realft(series - 1, (unsigned long) n, -1);
}

/*
 * window the data
 */
void do_window(float *series, int n) {
	int i;
	double wt, window;

	wt = 0.0;
	for (i = 0; i < n; i++) {
		window = Window(i);
		wt += window*window;
		series[i] *= window;
	}
	wt = n / wt;
	for (i = 0; i < n; i++) {
		series[i] *= wt;
	}
}

/*
 * smooth the spectral data
 */
void do_smooth(float *spectral, int smooth, int n) {

	int i, j, low, high, m;
	double *spec, sum, spectral0;

	if (smooth < 1) return;

	spectral0 = spectral[0];
	if ((spec = (double *) malloc(n * sizeof(double))) == NULL) {
		fprintf(stderr,"ran out of memory in do_smooth\n");
		exit(9);
	}
	m = 0;
	spec[m++] = 0.0;
	for (i = 2; i < n; i += 2) {
		spec[m++] = spectral[i]*spectral[i] + 
			spectral[i+1]*spectral[i+1];
	}

	/* nyquist frequency */
	spec[m] = 2.0 * spectral[1] * spectral[1];

	/* zero frequency  see p353 Time Series: Theory and Methods */
	sum = spec[1];
	for (i = 1; i <= smooth/2; i++) {
		sum += 2.0 * spec[i+1];
	}
	spec[0] = sum / smooth;

	/* spec[] now contains the raw spectrum */

	/* spectral <- spec */
	for (i = 0; i <= m; i++) {
		spectral[i] = spec[i];
	}

	/* spec = smooth(spectral) */
	for (i = 1; i <= m; i++) {
		low = i - (smooth/2);
		high = i + (smooth/2);
		if (high > m) high = m;

		sum = 0.0;
		for (j = low; j <= high; j++) {
			sum += spectral[abs(j)];
		}
		sum /= (high - low + 1);
		spec[i] = sqrt(sum);
	}

	spectral[0] = spectral0;
	spectral[1] = 0.5 * spec[m];
	j = 2;
	for (i = 1; i < m; i++, j += 2) {
		spectral[j] = spec[i];
		spectral[j+1] = 0.0;
	}

	free(spec);
	return;
}

/*
 * s_variance:
 *            returns variance of spectral data
 */
double s_variance(float *spectral, int n) {

	int i;
	double var;

	var = 0.5 * spectral[1]*spectral[1];

	for (i = 2; i < n; i++) {
		var += spectral[i] * spectral[i];
	}

	return 0.5 * var;
}

/*
 * s_corr:
 *        returns the correlation of two spectral data sets
 */
double s_corr(float *spectral1, float *spectral2, int n) {
 
        int i;
        double var1, var2, cross;
 
        var1 = 0.5 * spectral1[1]*spectral1[1];
        var2 = 0.5 * spectral2[1]*spectral2[1];
        cross = 0.5 * spectral1[1]*spectral2[1];
 
        for (i = 2; i < n; i++) {
                var1 += spectral1[i] * spectral1[i];
                var2 += spectral2[i] * spectral2[i];
        	cross += spectral1[i] * spectral2[i];
        }
 
        return cross/sqrt(var1*var2);
}

/*
 * s_covar:
 *        returns the covariance of two spectral data sets
 */
double s_covar(float *spectral1, float *spectral2, int n) {
 
        int i;
        double covar;
 
        covar = 0.5 * spectral1[1]*spectral2[1];
 
        for (i = 2; i < n; i++) {
                covar += spectral1[i] * spectral2[i];
        }
 
        return 0.5 * covar;
}

/*
 * variance:
 *          returns the variance of a series
 */
double variance(float *x, int n) {

	int i;
	double var, mean;

	var = mean = 0.0;

	for (i = 0; i < n; i++) {
		mean += x[i];
		var += x[i]*x[i];
	}
	mean /= (double) n;
	return var/n - mean*mean;
}

/*
 * corr:
 *      returns the correlation of two series
 */
double corr(float *x, float *y, int n) {
	int i;
	double mx, xx, my, yy, xy;
	mx = my = xx = yy = xy = 0.0;

	for (i = 0; i < n; i++) {
		mx += x[i];
		my += y[i];
		xx += x[i]*x[i];
		yy += y[i]*y[i];
		xy += x[i]*y[i];
	}

	mx /= n;
	my /= n;
	xx = xx/n - mx*mx;
	yy = yy/n - my*my;
	xy = xy/n - mx*my;
	return (xy)/sqrt(xx*yy);
}

void s_random_phase(float *a, float *ran, int n) {

	int i;
	extern long int iseed;
	double angle, cos_ang, sin_ang;

	/* retain means */
	ran[0] = a[0];

	/* Nyquist frequency */
        ran[1] = 1.414213562 * a[1] * sin(2.0*PI*ran1(&iseed));

	/* rest of the frequencies */
	for (i = 2; i < n; i += 2) {
		angle = 2.0 * PI * ran1(&iseed);
		cos_ang = cos(angle);
		sin_ang = sin(angle);
		ran[i]  =   cos_ang*a[i] + sin_ang*a[i+1];
		ran[i+1] = -sin_ang*a[i] + cos_ang*a[i+1];
	}
}

#ifdef MAIN
long int iseed=-123225;
#define N 32
main() {

int i;
float s[N], spec[N], v;
float t[N], tpec[N];

v=0;
for (i=0; i < N; i++) {
	s[i] = 2.0*sin(PI/9.0*i);
	t[i] = i - (N-1)/2.0;
	s[i] += 10;
}

do_fft(s, spec, N);



do_fft(t, tpec, N);
printf("corr %f %f\n", corr(s,t,N), s_corr(spec,tpec,N));
printf("var %f %f   %f %f\n",variance(s,N),variance(t,N),
	s_variance(spec,N), s_variance(tpec,N) );

}

#endif