/* =======================================================================
   ************************* RectangularToPolar **************************
   ======================================================================= */

/* This routine converts the rectangular real and imaginary coefficients generated
   by a DFT into polar mag and phase. */ 

#define M_PI 3.14159265358979323846

void RectangularToPolar( imaginary, real, mag, phase, dft_num_freqs, dft_num_points)
   double *imaginary, *real, *mag, *phase;
   long dft_num_freqs, dft_num_points; 
   {
   int i;

/* Compute the DC component. */

   mag[0] = real[0]/dft_num_points;
   phase[0] = 0;

/* Convert to polar form. */

   for ( i = 1; i < dft_num_freqs; i++ ) 
      {
      imaginary[i] *= 2.0/dft_num_points;
      real[i] *= 2.0/dft_num_points;

      mag[i] = sqrt(imaginary[i]*imaginary[i]+real[i]*real[i]);

/* Standard way -> atan(-b/a). */

      phase[i] = atan2(-imaginary[i], real[i])*180.0/M_PI;
      }
   }


/* =======================================================================
   ******************************* CalcDFT *******************************
   ======================================================================= */

/* CalcDFT performs discrete fourier analysis. y_time_data is the raw time 
   domain data to be used as input. Output consists of 5 arrays (imaginary,
   real, mag, phase and freq). dft_num_points is the number of elements in 
   the y_time_data, dft_fund_freq is the fundamental frequency of the 
   analysis.  dft_num_freqs is number of harmonics to calculate. The arrays 
   (y_time_data, imaginary, real, mag, phase and freq) must all be allocated 
   by the caller. y_time_data is a one-dimensional array of size dft_num_points
   while imaginary, real, mag, phase and freq are one dimensional arrays of 
   size dft_num_freqs. The y_time_data array must be reasonably distributed 
   over at least one full period of the fundamental frequency for the fourier 
   transform to be useful. Also, dft_num_freqs is always <= 1/2*dft_num_points. */

void CalcDFT( y_time_data, dft_num_points, dft_fund_freq, dft_num_freqs, 
   imaginary, real, mag, phase, freq)
   double *y_time_data; 
   long dft_num_points, dft_num_freqs; 
   double dft_fund_freq; 
   double *imaginary, *real, *mag, *phase, *freq; 
   {
   long i, j;

/* We are assuming that the caller has provided exactly one period of the 
   fundamental frequency. Clear output arrays */

   for ( i = 0; i < dft_num_freqs; i++ ) 
      {
      real[i] = 0;
      imaginary[i] = 0;
      }

   for (i = 0; i < dft_num_points; i++ ) 
      {
      for ( j = 0; j < dft_num_freqs; j++ ) 
         {
         imaginary[j] += (y_time_data[i]*sin(j*2.0*M_PI*i/((double) dft_num_points)));
         real[j] += (y_time_data[i]*cos(j*2.0*M_PI*i/((double) dft_num_points)));
         }
      }

/* Construct the frequency values, which are just multiples of the fundamental. */

   for ( i = 0; i < dft_num_freqs; i++ ) 
      freq[i] = i * dft_fund_freq;

   RectangularToPolar(imaginary, real, mag, phase, dft_num_freqs, dft_num_points);
   }