// fourier.java  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.                         

import java.io.*;

class fourier
{
  fourier(String filename, int ntrm)
  {
    boolean debug = true;
    double Pi = 3.1415926535897932384626433832795028841971;
    double x[] = new double[1000];
    double y[] = new double[1000];
    double a[] = new double[1000];
    double b[] = new double[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 npts, index, last;
    double tpi = 2.0*Pi;
    double dx, ti;
    FileReader file_in;
    BufferedReader in;
    String input_line, intStr;

    // read data 
    try
    {
      file_in = new FileReader(filename);
      in = new BufferedReader(file_in);
    }
    catch(Exception e)
    {
      System.out.println("read unable to open file "+filename);
      return;
    }
    npts = 0;
    try
    {
      input_line = in.readLine();
      // System.out.println("line="+input_line);
      while(input_line!=null)
      {
	index = 2; // x allow for two blanks
	last = input_line.indexOf(' ',index);
        // System.out.println("xindex="+0+", xlast="+last);
        intStr = input_line.substring(0,last);
        // System.out.println("xstring="+intStr);
	x[npts] = Double.parseDouble(intStr);
	index = last+1; // y
	last = input_line.length();
        // System.out.println("yindex="+index+", ylast="+last);
        intStr = input_line.substring(index,last);
        // System.out.println("ystring="+intStr);
	y[npts] = Double.parseDouble(intStr);
        if(debug)
           System.out.println("npts="+npts+", x="+x[npts]+", y="+y[npts]);
        npts++;
        input_line = in.readLine();
        // System.out.println("line="+input_line);
      }
      file_in.close();
    }
    catch(Exception e)
    {
      System.out.println("read exception"+e);
      return;
    }

    if(ntrm==-1) ntrm = npts;
    ntrm = Math.min(ntrm,(npts-2)/2);
    System.out.println(npts+" points read, using "+ntrm+" terms");

    xmin = x[0];
    xmax = x[0];
    ymin = y[0];
    ymax = y[0];
    for(int j=1; j<npts; j++)
    {
      xmin = Math.min(x[j],xmin);
      xmax = Math.max(x[j],xmax);
      ymin = Math.min(y[j],ymin);
      ymax = Math.max(y[j],ymax);
    }
    System.out.println("xmin="+xmin+", xmax="+xmax+
                       ", ymin="+ymin+", ymax="+ymax);

    // integration, simple trapezoidal rule 
    System.out.println("coefficients");
    for(int i=0; i<ntrm; i++)
    {
      a[i] = 0.0;
      b[i] = 0.0;
      ti = (double)i;
      for(int 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]*Math.cos(ti*tpi*xs1)+y[j+1]*
                               Math.cos(ti*tpi*xs2));
        b[i] = b[i] + dx*(y[j]*Math.sin(ti*tpi*xs1)+y[j+1]*
                               Math.sin(ti*tpi*xs2));
      }
      System.out.println("a["+i+"]="+a[i]+", b["+i+"]="+b[i]);
    }
    System.out.println(" ");
    // check fit 

    for(int j=0; j<npts; j++)
    {
      sum = 0.0;
      sumsq = 0.0;
      xs = (x[j]-xmin)/(xmax-xmin);
      for(int i=0; i<ntrm; i++)
      {
        ti = (double)i;
        sum = sum + a[i]*Math.cos(ti*tpi*xs) + b[i]*Math.sin(ti*tpi*xs);
      }
      err = y[j]-sum;
      sumsq = sumsq + err*err;
      avgerr = avgerr + Math.abs(err);
      maxerr = Math.max(Math.abs(err),maxerr);
      if(debug) System.out.println("y["+j+"]="+y[j]+", approx="+sum+
                                   ", err"+err);
    }
    rmserr = Math.sqrt(sumsq/(double)npts);
    avgerr = avgerr/(double)npts;
    System.out.println("avgerr="+avgerr+", rmserr="+rmserr+", maxerr="+maxerr);
    System.out.println("fourier.java finished ");
  } // end constructor fourier

  public static void main (String[] args)
  {
    String filename;
    int ntrm = -1;
    int argc = args.length;

    System.out.println("fourier.java running  ");
    if(argc>1) // file and ntrm 
    {
      filename = new String(args[0]);
      ntrm = Integer.parseInt(args[1]);
      System.out.println("fourier series of "+filename+" for "+ntrm+" terms");
      new fourier(filename, ntrm);
    }
    else if(argc>0) // just file, all terms 
    {
      filename = new String(args[0]);
      System.out.println("fourier series of "+filename+" for all terms");
      new fourier(filename, ntrm);
    }
    else
    {
      ntrm = 6;
      filename = new String("curve.dat");
      System.out.println("fourier series of "+filename+" for "+ntrm);
      new fourier(filename, ntrm);
    }
  } // end main  
} // end class fourier