/* To get some feel for Fourier transforms, try changing the parameters
NMAX (total number of points) and dt (interval between points), or
the function being transformed (defined in the first loop).  Here's a
refresher on some possibly useful Unix commands:
g++ discft.c -o discft      (if you save this file as discft.c, this
			     compiles and creates an executable, discft)
./discft                    (runs the program; the first two characters
			     may not be necessary)
./discft > outputfile       (runs the program, putting results in a file
			     instead of on the screen)
gnuplot                     (runs gnuplot, a simple graphing package)
  plot [xrange=0:10] "outputfile" with lines      
			    (within gnuplot, plots the first
      two columns of outputfile, using the range from 0 to 10 for the
      horizontal axis.  These columns are frequency and the real part
      of the Fourier transform.)
  plot [0:5][0:0.5] "outputfile" with linespoints
			     (within gnuplot, plots first two columns
      of outputfile, using range 0 to 10 for horizontal axis and 0 to 0.5
      for vertical axis, with both symbols and lines connecting the data 
      points)	
*/
 
#include <cmath>
#include <iostream>
#include <fstream>
#include <cstdlib>
using namespace std;
#define NMAX 80000        // You may need to change the number of points!
#define NFILE 80000      // Number of data points in file

main ()
{
   int i, k, n;
   double creal[NMAX], cimag[NMAX], freal[NFILE], fimag[NFILE];
   double expreal[NMAX], expimag[NMAX];
   double dt, tval, const1, const2, freq;
      
   /* Choose an appropriate dt. */
   dt = 1;  //one of the questions on Set 3 is to find the "real" dt value
   
   /* Parameters for the truncated cosine below. */
   double cosfreq = 6.6, numperiods = 2.5, tmax;
   
   /* This loop, presently commented out, calculates values of a cosine 
   function at evenly spaced points, in preparation for a discrete Fourier
   transform. */ 
/*   for (k=0; k<NMAX; k++) {
       freal[k] = cos(2*M_PI*dt*k);
       fimag[k] = 0;
   }
*/

   /* This loop, also commented out, calculates values of a single triangular
   peak, as on problem set 1.  With my numbers for a, dt, and NMAX, the flat 
   portion of the function extends from -5a to 5a; then the function
   repeats. (So there are additional triangle peaks from 9a to 11a, 19a to 21a,
   -11a to -9a, etc.) */
/*
   double a=10;
   for (k=0; k<=NMAX/2; k++) {
       tval = k*dt;
       freal[k] = 0;
       fimag[k]=0;
       if (tval<=a) {
           freal[k] = 2*(a-tval); 
       }
   }
   for (k=1; k<NMAX/2; k++) {
       freal[NMAX-k] = freal[k];
       fimag[NMAX-k] = 0;
   }
*/

   /* This is a new piece that reads data points from a file. */
   ifstream from("hw3signal.dat");
   for (k=0; k<NFILE; k++) {
       from >> freal[k];
       fimag[k]=0;
   }

     
   /* The next two loops, also commented out, do the same job as the
   previous two function assignments, but for a truncated cosine wave.
   (``Truncated" because the function is a cosine for small enough |k|
   but zero outside.) */
/*   const1 = NMAX*dt; 
   tmax = M_PI * numperiods / cosfreq;
   for (k=0; k<=NMAX/2; k++) {    //the <= will work whether NMAX is even or odd
        tval = k * dt;
        freal[k] = 0;
        if (fabs(tval)<=tmax) { 
            freal[k] = cos(cosfreq*tval);
        }
        fimag[k] = 0;
    }
*/    
    /* assign the second half of array by using evenness of cosine; note that
    periodicity means that freal[NMAX-k] is completely equivalent to
    freal[-k] */
/*    for (k=1; k<NMAX/2; k++) {
        freal[NMAX-k] = freal[k];
        fimag[NMAX-k] = 0;
    }      
*/                         
   
   /* This loop calculates values of the exponential factors.  It is more
   efficient to calculate the exponentials ONCE, and to save them, than
   to recalculate them for each new value of n.  This part should not need
   changing.  */
   const1 = 2. * M_PI/NMAX;
   for (i=0; i<NMAX; i++) {   
       const2 = const1 * i;    
       expreal[i] = cos(const2);
       expimag[i] = (-1)*sin(const2);
   }
   
   /* The following double loop calculates the Fourier coefficients, and
   does not need to be changed when choosing a new function. */
   for (n=0; n<NMAX; n++) {
       creal[n] = 0; cimag[n] = 0;
       for (k=0; k<NMAX; k++) {
           i = (k*n) % NMAX;
	   if (i<0) {i += NMAX;}
           creal[n] += freal[k] * expreal[i] - fimag[k] * expimag[i];
           cimag[n] += fimag[k] * expreal[i] + freal[k] * expimag[i];
       }
       cimag[n] = cimag[n]/NMAX;
       creal[n] = creal[n]/NMAX;
   }
       
   /* Finally, print the results.  */
   const1 = 2*M_PI/NMAX/dt;
   for (n=0; n<NMAX; n++) {
       freq = n*const1;
       cout << freq << "   "<<creal[n]<<"  "<<cimag[n]<<"\n";
   }
    
   return (0);
}