path: root/main.cpp

// FFT Test

#include "fft.h"

#include <complex>
#include <chrono>
#include <cmath>
#include <ctime>
#include <exception>
#include <functional>
#include <iostream>
#include <memory>
#include <numeric>
#include <string>
#include <vector>

const std::complex<double> i(0, 1);

const double tolerance = 0.000001;

using namespace std::complex_literals;

/// CPU time clock
class Timer {
public:
	Timer(): mStart(std::clock()){}
	void start(){mStart = std::clock();}
	double elapsed() { return (double(std::clock()) - mStart) / CLOCKS_PER_SEC; }
	static double precision() {
		static double result{0};
		
		if (result != 0)
			return result;
		
		std::clock_t start = std::clock();
		while (start == std::clock()) {}
		start = std::clock();

		std::clock_t end = start;
		while (end == std::clock()) {}
		end = std::clock();

		result = (double(end) - start) / CLOCKS_PER_SEC;
		return result;
	}
private:
	std::clock_t mStart;
};

// (Slow) DFT
std::vector<std::complex<double>> dft(const std::vector<std::complex<double>> &v) {
	std::vector<std::complex<double>> result(v.size(), 0);

	const double N = v.size();

	for (int k = 0; k < v.size(); k++) {
		for (int j = 0; j < result.size(); j++) {
			result[k] += v[j] * std::exp(-i * 2. * M_PI * double(k * j) / N);
		}
	}

	return result;
}

std::vector<double> freqs{ 2, 5, 11, 17, 29 }; // known freqs for testing

void generate(std::vector<std::complex<double>>& v) {
    for (int i = 0; i < v.size(); i++) {
        v[i] = 0.;
        // sum several known sinusoids into v[]
        for(int j = 0; j < freqs.size(); j++)
            v[i] += sin(2 * M_PI * freqs[j] * i / v.size() );
    }
}

namespace CooleyTukey {

// separate even/odd elements to lower/upper halves of array respectively.
// Due to Butterfly combinations, this turns out to be the simplest way 
// to get the job done without clobbering the wrong elements.
void separate(std::vector<std::complex<double>>::iterator a, int n) {
	std::vector<std::complex<double>> b(n/2);
	for(int i=0; i<n/2; i++)    // copy all odd elements to heap storage
		b[i] = a[i*2+1];
	for(int i=0; i<n/2; i++)    // copy all even elements to lower-half of a[]
		a[i] = a[i*2];
	for(int i=0; i<n/2; i++)    // copy all odd (from heap) to upper-half of a[]
		a[i+n/2] = b[i];
}

// N must be a power-of-2, or bad things will happen.
// Currently no check for this condition.
//
// N input samples in X[] are FFT'd and results left in X[].
// Because of Nyquist theorem, N samples means 
// only first N/2 FFT results in X[] are the answer.
// (upper half of X[] is a reflection with no new information).
void fft_recursive(std::vector<std::complex<double>>::iterator X, int N) {
    if(N < 2) {
        // bottom of recursion.
        // Do nothing here, because already X[0] = x[0]
    } else {
        separate(X,N);      // all evens to lower half, all odds to upper half
        fft_recursive(X,     N/2);   // recurse even items
        fft_recursive(X+N/2, N/2);   // recurse odd  items
        // combine results of two half recursions
        for(int k=0; k<N/2; k++) {
            std::complex<double> e = X[k    ];   // even
            std::complex<double> o = X[k+N/2];   // odd
                         // w is the "twiddle-factor"
            std::complex<double> w = exp( std::complex<double>(0,-2.*M_PI*k/N) );
            X[k    ] = e + w * o;
            X[k+N/2] = e - w * o;
        }
    }
}

// Cooley-Tukey
std::vector<std::complex<double>> fft(const std::vector<std::complex<double>> &v) {
	std::vector<std::complex<double>> result(v);

	fft_recursive(std::begin(result), result.size());
	
	return result;
}
}; // namespace CooleyTukey

class Measure {
	Timer mTimer;

public:
	using Data = std::vector<std::complex<double>>;
	
protected:
	Measure* mReference; // Reference measurement: we should calculate similar to this one and be faster than this one
	
	const Data& mIn;
	Data mResult;

	std::string mName;
	double mElapsed;

public:
	Measure(const Data& in): mReference(nullptr), mIn(in), mResult(in.size()) {}

	virtual void run_impl() = 0; // implemented by subclasses

	void run() {
		mTimer.start();
		run_impl();
		mElapsed = mTimer.elapsed();
		std::cout << mName << ": " << mElapsed << "s\n";

		if (mReference) {
			if (mResult != mReference->mResult) {
				double diff = std::transform_reduce(std::begin(mResult), std::end(mResult), std::begin(mReference->mResult), double(0.0),
								    [](const double& x, const double& y) -> double
								    { return x + y;},
								    [](const std::complex<double>& x, const std::complex<double>&y) -> double
								    { return abs(y - x);}
								    );
				if (diff > tolerance)
					std::cout << "Error: Results diff: " << diff << "\n";
			}

			if (mElapsed > mReference->mElapsed) {
				std::cout << "Error: " << mName << " too slow!\n";
			}
		}
	}

	double elapsed() { return mElapsed; }

	Data& result() { return mResult; }
};

class MeasureDFT: public Measure {
public:
	MeasureDFT(const Data& in): Measure(in){ mName = "DFT";}
	void run_impl() override {
		mResult = dft(mIn);
	}
};

class MeasureFFT: public Measure {
public:
	MeasureFFT(const Data& in, Measure& reference): Measure(in) { mName = "FFT Cooley-Tukey"; mReference = &reference;}
	void run_impl() override {
		mResult = CooleyTukey::fft(mIn);
	}
};

class MeasureFFT_RR: public Measure {
	RIT::FFT mRR;
public:
	MeasureFFT_RR(const Data& in, Measure& reference): Measure(in), mRR(in.size()){ mName = "FFT RR"; mReference = &reference;}
	void run_impl() override {
		mResult = mRR(mIn);
	}
};

class MeasureFFT_RR_half: public Measure {
	RIT::FFT mRR;
public:
	MeasureFFT_RR_half(const Data& in, Measure& reference): Measure(in), mRR(in.size(), true){ mName = "FFT RR half"; mReference = &reference; }
	void run_impl() override {
		mResult = mRR(mIn);
	}
};

class MeasureFFT_RR_half_magnitudes: public Measure {
	RIT::FFT mRR;
public:
	MeasureFFT_RR_half_magnitudes(const Data& in, Measure& reference): Measure(in), mRR(in.size(), true){ mName = "FFT RR half magnitudes"; mReference = &reference; }
	void run_impl() override {
		mResult = mRR(mIn);
		RIT::magnitudes(mResult);
	}
};

int main(int argc, char* argv[]) {
	std::vector<std::complex<double>> v(4096, 0);

	generate(v);

	std::cout << "Timer precision: " << Timer::precision() << "s\n";

	MeasureDFT measureDFT(v);
	measureDFT.run();

	MeasureFFT measureFFT(v, measureDFT);
	measureFFT.run();

	MeasureFFT_RR measureFFT_RR(v, measureFFT);
	measureFFT_RR.run();

	MeasureFFT_RR_half measureFFT_RR_half(v, measureFFT_RR);
	measureFFT_RR_half.run();
	
	MeasureFFT_RR_half_magnitudes measureFFT_RR_half_magnitudes(v, measureDFT);
	measureFFT_RR_half_magnitudes.run();
	
	return 0;
}