Highly useful in digital image processing and signal processing, Fourier Transforms are mathematical models that can transform functions of time (seconds) into functions of frequency (Hz). In computer science, there are 3 methods in which a Fourier transform can be computed. The first is the Discrete Fourier Transform (DFT), the second is a recursive Fast Fourier Transform (FFT1), and the third is an even faster, non-recursive Fast Fourier Transform that uses bit-reversal (FFT2).
Discrete Fourier Transform (DFT)
Discrete Fourier Transform (DFT) is the discrete version of the Fourier Transform (FT) that transforms a signal (or discrete sequence) from the time domain representation to its representation in the frequency domain. This version typically has a time complexity of O(n2).
#include <fstream>
#include <iostream>
#include <iomanip>
#include <string>
#include <vector>
#include <complex>
#define PI 3.14159265359
/**
* @brief Computes the DFT
*
* @param input
* @return std::vector<std::complex<double>>
*/
std::vector<std::complex<double>> DFT(const std::vector<std::complex<double>> &input)
{
std::vector<std::complex<double>> output; // return vector
size_t n = input.size();
for (size_t k = 0; k < n; ++k) // For each output element
{
std::complex<double> sum(0, 0); // exponential sum
for (size_t t = 0; t < n; ++t) // For each input element
{
double angle = 2 * PI * t * k / n;
sum += input[t] * std::exp(std::complex<double>(0, -angle));
}
output.push_back(sum);
}
return output;
}
int main()
{
// Introduction
std::cout << "This program calculates the Discrete Fourier Transform of complex numbers provided by a text file." << std::endl;
std::cout << "======================================================" << std::endl;
// Prompt user for file
std::string filename;
std::cout << "Please enter the name of the text file (with extension .txt) containing the input: " << std::flush;
std::ifstream in;
while (true)
{
in.close();
in.clear();
std::getline(std::cin, filename);
in.open(filename.c_str());
if (in)
{
break;
}
std::cout << "Invalid file. Please enter a valid input file name: " << std::flush;
}
std::cout << std::endl;
// Count the number of complex values (lines) in file
// Assign to N
int N = 0;
std::string line;
std::vector<std::complex<double>> vec;
while (std::getline(in, line))
{
++N;
// Parse the real and imaginary parts of each complex value
std::string re_s, im_s;
double re, im;
std::string::size_type pos = line.find('+');
std::string::size_type neg = line.find('-');
// EDGE CASE: Number begins with '-'
if (neg == 0)
{
std::string newline = line.substr(1);
neg = newline.find('-');
if (neg < 8192)
{
++neg;
}
}
pos = std::min(pos, neg);
char sign = line[pos]; // '+' or '-'
re_s = line.substr(0, pos);
re = stod(re_s);
std::string::size_type i_pos = line.find('i');
im_s = line.substr(pos + 1, i_pos - pos - 1);
im = stod(im_s);
if (sign == '-')
{
im *= -1;
}
std::complex<double> mycomplex(re, im);
// Store them in a vector
vec.push_back(mycomplex);
}
// Close file stream
in.close();
in.clear();
// Do the math
std::vector<std::complex<double>> res = DFT(vec);
// File output
std::string::size_type dot_pos = filename.find('.');
std::string input = filename.substr(0, dot_pos);
std::ofstream output_file("out_dft1_" + std::to_string(N) + "_" + input + ".txt");
for (int i = 0; i < N; ++i)
{
// get sign
std::string sign;
if (res[i].imag() < 0)
{
sign = "";
}
else
{
sign = "+";
}
output_file << std::setprecision(4) << res[i].real() << sign << res[i].imag() << "i" << "\n";
}
std::cout << "==================================================" << std::endl;
std::cout << "Calculation complete. Please check your directory for the output file." << std::endl;
output_file.close();
return 0;
}The following input…
1+1i
-1+1i
1-1i
-1-1i
-3-1i
4-2i
0-3i
2-6i…yields the following output:
3-12i
10.24+5.243i
3+2i
9.071+10.07i
-5+4i
1.757-3.243i
-9+6i
-5.071-4.071iRecursive Fast Fourier Transform (FFT1)
Using recursion, the Discrete Fourier Transform algorithm can be optimized to O(nlog(n)) complexity:
void FFT(std::vector<std::complex<double>>& input)
{
int N = input.size();
// Base case
if (N == 1)
{
return;
}
// Recursive case
std::vector<std::complex<double>> evens(N / 2);
std::vector<std::complex<double>> odds(N / 2);
for (int i = 0; 2 * i < N; i++)
{
evens[i] = input[2 * i];
odds[i] = input[2 * i + 1];
}
FFT(evens, tt);
FFT(odds, tt);
double angle = 2 * PI / N;
std::complex<double> w(1);
std::complex<double> wn(cos(angle), sin(angle));
for (int i = 0; 2 * i < N; i++)
{
input[i] = evens[i] + w * odds[i];
input[i + N/2] = evens[i] - w * odds[i];
w *= wn;
}
}Bit-Reversal Fast Fourier Transform (FFT2)
Finally, using bit-reversal, an even faster FFT algorithm can be written. My first step was to write a function for reversing all bit values from 0 to N:
/**
* @brief Reverses the bits
*
* @param N
* @return std::vector<std::vector<int>>
*/
std::vector<std::vector<std::string>> bitrev(int N)
{
std::vector<std::vector<std::string>> res;
const int k = std::log2(N);
for (int i = 0; i < N; ++i)
{
std::vector<std::string> line;
// 1st: push back i
line.push_back(std::to_string(i));
// 2nd: push back b1 (binary form of i)
std::bitset<16> b1bs(i);
// convert bitset to string
std::string b1String = b1bs.to_string();
// excise the leading zeroes
std::size_t pos = 16 - k;
std::string b1SubStr = b1String.substr(pos);
line.push_back(b1SubStr);
// reverse bits of b1, store as b2
std::bitset<16> b2TMP(b1bs);
std::string b2String = b1SubStr;
std::reverse(b2String.begin(), b2String.end());
std::bitset<16> b2(b2String);
// 3rd: convert b2 to integer ir, push back ir
unsigned long ir = b2.to_ulong();
line.push_back(std::to_string(ir));
// 4th: push back b2
line.push_back(b2String);
res.push_back(line);
}
return res;
}The following table was generated using this value with N = 16. The first column contains the original indices from 0 to N-1, the second column contains the binary representation of those indices, the third column contains the decimal representation of those reversed bits, and the fourth column contains the binary representation:
0 0000 0 0000
1 0001 8 1000
2 0010 4 0100
3 0011 12 1100
4 0100 2 0010
5 0101 10 1010
6 0110 6 0110
7 0111 14 1110
8 1000 1 0001
9 1001 9 1001
10 1010 5 0101
11 1011 13 1101
12 1100 3 0011
13 1101 11 1011
14 1110 7 0111
15 1111 15 1111After the indices are reversed, the algorithm can be called:
void FFT2(std::complex<double>* input
, std::complex<double>* output
, int N)
{
// Shuffle input with reversed bits
std::vector<std::vector<std::string>> indices = bitrev(N);
for (int i = 0; i < N; ++i)
{
// 3rd column contains reversed indices
output[i] = input[std::stoi(indices[i][2])];
}
// Keep left shifting until the bitlength of N is reached
for (int length = 2; length <= N; length <<= 1)
{
double angle = 2 * PI / length;
std::complex<double> wlen(cos(angle), sin(angle));
// Instead of recursion, sum up adjacent entries in output
for (int i = 0; i < N; i += length)
{
std::complex<double> w(1);
// Only need to go through half the length (powers of 2)
for (int j = 0; j < length / 2; ++j)
{
std::complex<double> u = output[i + j];
std::complex<double> v = output[i + j + length / 2] * w;
output[i + j] = u + v;
output[i + j + length / 2] = u - v;
w *= wlen;
}
}
}
// Reverse bits
std::complex<double>* temp = new std::complex<double>[N];
for (int i = 0; i < N; ++i)
{
temp[i] = output[i];
}
for (int i = 1; i < N; ++i)
{
output[i] = temp[N - i];
}
delete [] temp;
}Results
Using the std::chrono header, one can time each of the three implementations and clearly see how much more efficient FFT2 is then the rest:
FUNCTION TIME
====================================
DFT 42 microseconds
FFT 10 microseconds
FFT2 3 microseconds