Simple Discrete Fourier Transform

simple-dft.cpp

// Based on: http://nayuki.eigenstate.org/page/how-to-implement-the-discrete-fourier-transform
#include <iostream>
#include <vector>
#include <complex>
#include <cmath>
#include <cstdlib>
using namespace std;

vector<complex<double> > compute_dft(vector<complex<double> > in) {
	int n = in.size();
	vector<complex<double> > out(n);
	// k-th harmonic
	for (int k = 0; k < n; k++) {
		double sumreal = 0;
		double sumimag = 0;
		// t-th sample
		for (int t = 0; t < n; t++) {
			// From Eulera's formula e^ix = cos x + i sin x
			double arg = 2 * M_PI * t * k / n;
			double cosinus = cos(arg);
			double sinus = sin(arg);
			sumreal += in[t].real() * cosinus + in[t].imag() * sinus;
			sumimag += -in[t].real() * sinus + in[t].imag() * cosinus;
		}
		out[k] = complex<double>(sumreal, sumimag);
	}
	return out;
}

template<class T>
void getInput(T &input) {
	bool inputOK;
	cin.exceptions(ios_base::failbit); // Throw exceptions on fail
	do {
		try {
			inputOK = true;
			cin >> input;
		} catch (const ios_base::failure& e) {
			cout << "Bad input. Please retry." << endl;
			inputOK = false;
			cin.clear(); // Clear failbit
			cin.ignore(); // Clear input
		}
	} while (!inputOK);
}

// You can check correctness in octave:
// $ octave
// > x = [ 1 + 3i, 2 + 2i, 3 + 1i ]
// > fft(x)
int main() {
	int n; // Input samples count
	vector<complex<double> > in, out;

	cout << "How many samples?" << endl;
	getInput(n);

	in.resize(n);
	for (int i = 0; i < n; i++) {
		cout << "Reading sample " << i + 1 << "." << endl;
		// Input in format: 1, (1,0), (1,3), etc.
		getInput(in[i]);
	}

  out = compute_dft(in);
	cout << scientific;
	cout << "n = " << n << endl;
	for (int i = 0; i < n; i++) {
		cout << i + 1 << ". z = " << out[i] << ", |z| = " << abs(out[i])
				<< endl; // Prints harmonic and module
	}

	return EXIT_SUCCESS;
}