Discrete Fourier Transform

In the field of Mathematics, the Discrete Fourier Transform (DFT) is a way to convert a finite list of samples of a function (Complex Numbers) into a list of coefficients of a finite combination of complex sinusoids, ordered by their frequencies. 

In Digital Signal Processing, the samples can be signals varying over time. For example, pressure of a sound wave, a radio signal etc.

The following is the basic fundamental equation behind the Discrete Fourier Transform:
\[\color{silver} { \huge X(k) = \sum_{t=0}^{N-1}x(t)e^{-{i2\pi tk}/N}}\]
In this equation,

x(t) is the input sample/complex number

i is an imaginary number

N is the total number of input samples

t and k are both indices starting from 0 to N

e^(i * angle) is essentially an exponential that can be expressed as a sum of sinusoids according to Euler's Formula:
\[\color{silver} { \huge e^{i\theta} = cos(\theta) + isin(\theta)}\]
X(k) is the output complex number that encodes both the amplitude and phase of a sinusoidal component of function x(t) for each k cycle.
My program takes in two command line arguments, the number of samples and the name of the text file containing the complex number samples. 

My algorithm is pretty straightforward following the equation directly. One loop to iterate through the k cycles and the second loop to add up all the complex sums.

Sample Input (N = 8)

1+1i
-1+1i
1-1i
-1-1i
-3-1i
4-2i
0-3i
2-6i

DFT Output

3-12i
10.24+5.243i
3+2i
9.071+10.07i
-5+4i
1.757-3.243i
-9+6i
-5.071-4.071i
The program below takes about 235.013 millseconds to calculate the DFT for 1024 samples. This is not as efficient as we would want it to be. This is when the FFT algorithms solves the efficiency problem.

C++ Source Code

/*****************************************************************************/
/*!
\Name   Discrete Fourier Transform
\Author Deepak Chennakkadan
\Input  <N> <FILENAME>
        <N>: Number of input Complex Numbers
        <FILENAME>: Text file containing the Complex Numbers
\Output Prints the DFT Output        
*/      
/*****************************************************************************/

#include <fstream>
#include <complex>
#include <cstdlib>
#include <cmath>
#include <iostream>
#include <iomanip>
#include <string>
#include <complex>
#include <sys/time.h>

using namespace std;

/***************************************************************************/
/*!
\Brief  Overloaded operator to format Complex Number output
*/
/***************************************************************************/
ostream &operator<<(ostream &out, complex<double> c)
{
  if (c.imag() < 0) {
    out << c.real() << c.imag() << "i";
  }
  else {
    out << c.real() << "+" << c.imag() << "i";
  }
  return out;
}

/***************************************************************************/
/*!
\Brief  Function to calculate elapsed program time
*/
/***************************************************************************/
float timedifference_msec(struct timeval t0, struct timeval t1)
{
    return (t1.tv_sec - t0.tv_sec) * 1000.0f + (t1.tv_usec - t0.tv_usec) / 1000.0f;
}

/***************************************************************************/
/*!
\Brief  Main Function
*/
/***************************************************************************/
int main(int argc, char *argv[]) 
{
  struct timeval t0;
  struct timeval t1;
  float elapsed;

  gettimeofday(&t0, 0);

  int N;
  sscanf(argv[1], "%d", &N);

  const double PI = 4.0f*atan(1.0f);
  const double TWOPI = 2.0f*PI;

  complex<double>* inputs;
  inputs = new complex<double>[N];

  ifstream f1(argv[2]);

  // Loop to read in the Complex Numbers from file
  for (int i = 0; i < N; i++)
  {
    if (f1)
    {
      char buf[50];
      f1.getline(buf, 50);
      double d1 = 0, d2 = 0;
      sscanf(buf, "%lf %lf", &d1, &d2);
      complex<double> c1(d1, d2);
      inputs[i] = c1;
    }
  }

  // Loop for each k cycle
  for (int k = 0; k < N; k++)
  {
    double sumreal = 0;
    double sumimag = 0;

    int t;

    // Loop for calculating the sum
    for (t = 0; t < N; t++)
    {
      // e^((2*PI*k*t)/N)
      double angle = - TWOPI * t * k / N;

      // Real Sum after complex dot product
      sumreal += inputs[t].real() * cos(angle) - inputs[t].imag() * sin(angle);

      // Imaginary Sum after complex dot product
      sumimag += inputs[t].real() * sin(angle) + inputs[t].imag() * cos(angle);
    }

    complex<double> out(sumreal, sumimag);

    cout << out << endl;
  }

  delete[] inputs;

  // Time Stuff
  gettimeofday(&t1, 0);  
  elapsed = timedifference_msec(t0, t1);
  printf("\nCode executed in %f milliseconds.\n", elapsed);

  return 0;
}