Recursive Fast Fourier Transform


For this I essentially use the same equation as I used for the DFT. But the equation is split into two separate and smaller DFTs. One for the Even indices and the other for the Odd indices. 
\[\color{silver} { \huge X(k) = \sum_{t=0}^{N-1}x(t)e^{\frac{-{i2\pi tk}}{N}}}\]

Even 

Odd
\[\color{silver} { \huge = \sum_{t=0}^{N/2-1}x(2t)e^{\frac{-i2\pi (2t)k}{N}} +\sum_{t=0}^{N/2-1}x(2t+1)e^{\frac{-i2\pi (2t+1)k}{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 utilizes the properties of recursion to calculate both the Even and Odd DFTs. 

Input Samples (N = 8)

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

FFT 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 2 millseconds to calculate the FFT for 1024 samples. This is much more efficient than the DFT algorithm mainly because of splitting the equation into two separate DFTs.

C++ Source Code

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

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

using namespace std;

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

/***************************************************************************/
/*!
\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  Function to calculate the FFT for given N samples
*/
/***************************************************************************/
void fft(complex<double> *input, complex<double> *output, int N, int step)
{
  if (step < N) 
  {
    // Recursive Functionality 
    fft(output, input, N, step * 2);
    fft(output + step, input + step, N, step * 2);
    
    // For loop for k cycles
    for (int k = 0; k < N; k += 2 * step)
    {
      double angle = -PI * k / N;
      double realPart, imagPart;      

      // Calculating the Complex Sum
      realPart = output[k + step].real() * cos(angle) - output[k + step].imag() * sin(angle);
      imagPart = output[k + step].real() * sin(angle) + output[k + step].imag() * cos(angle);

      complex<double> out(realPart, imagPart);

      input[k / 2] = output[k] + out;
      input[(k + N) / 2] = output[k] - out;
    }
  }
}

/***************************************************************************/
/*!
\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);

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

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

  ifstream f1(argv[2]);

  // For loop to read in Complex Number Samples
  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);
      buffer[i] = c1;
    }
  }

  cout << setprecision(4);

  for (int i = 0; i < N; i++)
  {
    outputs[i] = buffer[i];
  }

  // Calculate FFT
  fft(buffer, outputs, N, 1);

  // Print out the output
  for (int i = 0; i < N; i++)
  {
    std::cout << buffer[i] << std::endl;
  }

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

  return 0;
}