/*
 *
 *  COMPUTE VARIANCE OF THE CURVE RESULTING FROM FFT FILTER
 *
 *  CUTOFF  - FILTER CUTOFF VALUE DEFINED IN SAME UNITS AS X(*)
 *  DINTERV - SAMPLING INTERVAL SPECIFIED IN SAME UNITS AS X(*)
 *  XFILT   - TIME VALUES OF FILTERED CURVE
 *  FILT    - ARRAY OF FILTERED VALUES
 *  NFILT   - NUMBER OF POINTS IN XFILT(*), FILT(*)
 *  X       - ARRAY SPECIFYING TIME VALUES OF ORIGINAL DATA POINTS
 *  Y       - ARRAY OF CONCENTRATIONS OF ORIGINAL DATA POINTS
 *
 */
#include <stdio.h>
#include <stdlib.h>

#include <math.h>

#include "filter.h"

#define NSIZE 2048

/***************************************************************/
float filtvar (
float cutoff,
float dinterv,
double xfilt[],
double filt[],
int nfilt,
double x[],
double y[],
int n
)
{
	int n0, z, i, j;
	int nw;
	double *xtemp, *ytemp, weights[2*NSIZE+2], xtemp2[2*NSIZE+2];
	double ssw, yv, yd, sm, cor, r;
	float rmean, rsd, var;

/*
 * First step: Compute weights of filter by filtering a single
 * point in the middle of zero values (impulse response)
 */

#ifdef DEBUG
printf ("Start filtvar, cutoff = %f, dinterv = %e\n", cutoff, dinterv);
printf ("n = %d, nfilt = %d\n", n, nfilt);
#endif

	n0 = 4* (int) (cutoff/dinterv);

	z=1;
	while (n0 > z) z = z<<1;
	n0 = z;

	if ( NSIZE < z) n0 = NSIZE;
	if (n0 <= 1) {
		return (0.0);
	}

#ifdef DEBUG
printf ("In filtvar, n0 = %d\n", n0);
#endif

	ytemp = (double *) calloc (n0, sizeof(double));
	xtemp = (double *) calloc (n0, sizeof(double));
if (!xtemp) {
fprintf (stderr, "In filtvar: Can't allocate double array of %d elements.\n", n0);
exit (1);
}
	for (i=0; i<n0; i++) {
		xtemp[i] = i*dinterv;
		ytemp[i] = 0.0;
	}
	ytemp[n0/2] = 1;

	filter (xtemp, ytemp, n0, dinterv, cutoff, xtemp2, weights, &nw);
#ifdef DEBUG
printf ("finished filter, nw = %d\n", nw);
#endif
	free(ytemp);
	free(xtemp);

/*
 *  SECOND STEP: COMPUTE SUM OF SQUARES OF WEIGHTS
 */

	ssw = 0.0;
	for (i=0; i< nw; i++) {
		ssw += weights[i]*weights[i];
	}

#ifdef DEBUG
printf ("ssw = %e\n", ssw);
#endif

/*
 *  THIRD STEP: COMPUTE RESIDUALS OF DATA ABOUT FILTERED CURVE.
 *              INTERPOLATE FILTERED POINTS SINCE ORIGINAL RAW
 *              DATA IS NOT NECESSARILY AT SAME LOCATION IN TIME
 */

	ytemp = (double *) calloc (nfilt, sizeof(double));
	xtemp = (double *) calloc (n, sizeof(double));

	spline (xfilt, filt, nfilt, ytemp);
/*
	spline_f (xfilt, filt, nfilt, x, ytemp, dummy, n);
*/
	for (i=0; i<n; i++) {
		splint (xfilt, filt, ytemp, nfilt, x[i], &yv, &yd);
		xtemp[i] = y[i] - yv;
/*
		xtemp[i] = y[i] - ytemp[i];
*/
	}

/*
 *  FOURTH STEP: COMPUTE RESIDUAL STANDARD DEVIATION
 */

	means (xtemp, &rmean, &rsd, n);

#ifdef DEBUG
printf ("rsd = %e\n", rsd);
#endif

/*
 *  FIFTH STEP: COMPUTE LAG 1 AUTO COVARIANCE
 */

	sm = 0.0;
	for (i=0; i<n-1; i++) {
		sm += (xtemp[i]-rmean)*(xtemp[i+1]-rmean);
	}
	cor = sm/(n-1)/(rsd*rsd);

#ifdef DEBUG
printf ("cor = %e\n", cor);
#endif

/*
 *  SIXTH STEP: COMPUTE AUTO COVARIANCES
 */

	sm = 0.0;
	for (i=0; i<nw-1; i++) {
		for (j=i+1; j<nw; j++) {
			r = pow(cor, (double) (j-i));
			if (r < 1e-3) goto out;
			r = r*weights[i]*weights[j];
			sm += r;
		}
	}


out:
	var = rsd*rsd*ssw+sm;

#ifdef DEBUG
printf ("sm = %e, var = %e\n", sm, var);
#endif

	free(ytemp);
	free(xtemp);

	return (var);
}