/**********************************************************

     DETREND is intended for detrending of simple
     vectors consisting of equally spaced data.  The
     index serves as the x-coordinate.

*/
#include "common.h"
#include "vector.h"

#if PROTOTYPE_ALLOWED
void detrend(float array[], int n)
#else
void detrend(array, n)
float array[];
int n;
#endif
{
     double x, y, a, b;
     double sy = 0.0,
            sxy = 0.0,
            sxx = 0.0;
     int i;

     for (i=0, x=(-n/2.0+0.5); i<n; i++, x+=1.0)
     {
          y = array[i];
          sy += y;
          sxy += x * y;
          sxx += x * x;
     }
     b = sxy / sxx;
     a = sy / n;

     for (i=0, x=(-n/2.0+0.5); i<n; i++, x+=1.0)   
          array[i] -= (a+b*x);
}
/*********** end: detrend() **********************************/

/**************************************************************

          Simple testing function for detrend:

****************************************************************/

#if 0

#include <stdio.h>
#include <stdlib.h>
#include "prtarray.h"  /* print_array_*() */

#define N 512
#if PROTOTYPE_ALLOWED
void main(int argc, char *argv[])
#else
void main(argc, argv)
int argc;
char *argv[];
#endif
{
     float data[N];
     int i;
     int n;
     double a, b;

     n = atoi(argv[1]);
     a = atof(argv[2]);
     b = atof(argv[3]);

     for (i=0; i<n; i++)  data[i] = a + i * b;
     print_array_float(stdout, data, n, 5, "%12.3f ");
     detrend(data, n);
     print_array_float(stdout, data, n, 5, "%12.3f ");
     return;
}
#endif