From ba9179c7991eeaf1b8d59a0db3975f049b2735b7 Mon Sep 17 00:00:00 2001 From: Roland Stigge Date: Sun, 27 Jan 2019 13:32:40 +0100 Subject: Added half FFT and magnitudes --- main.cpp | 246 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 246 insertions(+) create mode 100644 main.cpp (limited to 'main.cpp') diff --git a/main.cpp b/main.cpp new file mode 100644 index 0000000..42cd466 --- /dev/null +++ b/main.cpp @@ -0,0 +1,246 @@ +// FFT Test + +#include "fft.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +const std::complex 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> dft(const std::vector> &v) { + std::vector> 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 freqs{ 2, 5, 11, 17, 29 }; // known freqs for testing + +void generate(std::vector>& 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>::iterator a, int n) { + std::vector> b(n/2); + for(int i=0; i>::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 e = X[k ]; // even + std::complex o = X[k+N/2]; // odd + // w is the "twiddle-factor" + std::complex w = exp( std::complex(0,-2.*M_PI*k/N) ); + X[k ] = e + w * o; + X[k+N/2] = e - w * o; + } + } +} + +// Cooley-Tukey +std::vector> fft(const std::vector> &v) { + std::vector> result(v); + + fft_recursive(std::begin(result), result.size()); + + return result; +} +}; // namespace CooleyTukey + +class Measure { + Timer mTimer; + +public: + using Data = std::vector>; + +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& x, const std::complex&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> 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; +} + + -- cgit v1.2.3