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:
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.
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 Output3-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; }