/* fourier.c  for unequally spaced points, FFT may not work,   */
/*            here we do a DFT Discrete Fourier Transform to   */
/*            get the Fourier Series to approximate y = f(x)   */
/*            from a data file with x,y values not equally     */
/*            spaced. x in increasing order.                   */
/*                                                             */
/* y_approx(xj) = a0/2 + sum_i=j,n( ai*cos(i*2*Pi*xj/T)+       */
/*                                  bi*sin(i*2*Pi*xj/T) )      */
/*            ai = 2/T * integral_0,T( f(x)*cos(i*2*Pi*xj/T)dx */
/*            bi = 2/T * integral_0,T( f(x)*sin(i*2*Pi*xj/T)dx */
/*                                                             */
/* for numerical computation, x is scaled to be in range [0,1] */
/* xs = (x-xmin)/(xmax-xmin)  T=1                              */
/* then y_approx is computed  y_approx((x-xmin)/(xmax-xmin))   */
/*                                                             */
/* Real data has coefficients aliased at and above Nyquist     */
/* thus n max is (npoints-2)/2                                 */ 
/* a "smoother" fit can be obtained by a Fejer approximation.  */
/* many coefficients may be near zero.                         */

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

#undef max
#define max(x,xmax) ((x)>(xmax)?(x):(xmax))
#undef min
#define min(x,xmin) ((x)<(xmin)?(x):(xmin))
#undef abs
#define abs(x) ((x)>0.0?(x):(-(x)))

#define Pi 3.1415926535897932384626433832795028841971

static int debug = 1;

int main(int argc, char* argv[])
{
  FILE *inp;
  double x[1000], y[1000], a[1000], b[1000];
  double xmin, xmax, ymin, ymax, xs, xs1, xs2;
  double err, avgerr=0.0, rmserr=0.0, maxerr=0.0;
  double sum, sumsq=0.0;
  int i, j, npts, ntrm;
  double tpi = 2.0*Pi;
  double dx, ti;

  printf("fourier.c running \n");
  if(argc>2) /* file and ntrm */
  {
    inp = fopen(argv[1],"r");
    if(inp==NULL)
    {
      printf("can not open %s, bye. \n", argv[1]);
      return 1;
    }
    ntrm = atoi(argv[2]);
    printf("fourier series of %s for %d terms \n", argv[1], ntrm);
  }
  else if(argc>1) /* just file, all terms */
  {
    inp = fopen(argv[1],"r");
    if(inp==NULL)
    {
      printf("can not open %s, bye. \n", argv[1]);
      return 1;
    }
    ntrm = -1;
    printf("fourier series of %s for all terms \n", argv[1]);
  }
  else
  {
    inp = fopen("curve.dat","r");
    if(inp==NULL)
    {
      printf("can not open %s, bye. \n", "curve.dat");
      return 1;
    }
    ntrm = 6;
    printf("fourier series of %s for %d terms \n", "curve.dat", ntrm);
  }


  /* read data */
  npts = 0;
  while(!feof(inp))
  {
    fscanf(inp, "%lf %lf\n", &x[npts], &y[npts]);
    if(debug) printf("npts=%d, x=%f, y=%f \n", npts, x[npts], y[npts]);
    npts++;
  }
  if(ntrm==-1) ntrm = npts;
  ntrm = min(ntrm,(npts-2)/2);
  printf("%d points read, using %d terms \n", npts, ntrm);

  xmin = x[0];
  xmax = x[0];
  ymin = y[0];
  ymax = y[0];
  for(j=1; j<npts; j++)
  {
    xmin = min(x[j],xmin);
    xmax = max(x[j],xmax);
    ymin = min(y[j],ymin);
    ymax = max(y[j],ymax);
  }
  printf("xmin=%f, xmax=%f, ymin=%f, ymax=%f \n", xmin, xmax, ymin, ymax);

  /* integration, simple trapezoidal rule */
  printf("coefficients \n");
  for(i=0; i<ntrm; i++)
  {
    a[i] = 0.0;
    b[i] = 0.0;
    ti = (double)i;
    for(j=0; j<npts-1; j++)
    {
      xs1 = (x[j]-xmin)/(xmax-xmin);
      xs2 = (x[j+1]-xmin)/(xmax-xmin);
      dx = xs2-xs1;
      a[i] = a[i] + dx*(y[j]*cos(ti*tpi*xs1)+y[j+1]*cos(ti*tpi*xs2));
      b[i] = b[i] + dx*(y[j]*sin(ti*tpi*xs1)+y[j+1]*sin(ti*tpi*xs2));
    }
    printf("a[%d]=%f, b[%d]=%f \n", i, a[i], i, b[i]);
  }
  printf("\n");
  /* check fit */

  for(j=0; j<npts; j++)
  {
    sum = 0.0;
    sumsq = 0.0;
    xs = (x[j]-xmin)/(xmax-xmin);
    for(i=0; i<ntrm; i++)
    {
      ti = (double)i;
      sum = sum + a[i]*cos(ti*tpi*xs) + b[i]*sin(ti*tpi*xs);
    }
    err = y[j]-sum;
    sumsq = sumsq + err*err;
    avgerr = avgerr + abs(err);
    maxerr = max(abs(err),maxerr);
    if(debug) printf("y[%d]=%f, approx=%f, err=%g \n", j, y[j], sum, err);
  }
  rmserr = sqrt(sumsq/(double)npts);
  avgerr = avgerr/(double)npts;
  printf("avgerr=%g, rmserr=%g, maxerr=%g \n", avgerr, rmserr, maxerr);

  printf("fourier.c finished \n");
  return 0;
}