From db0d20a6c322211593fecf209898f2435468e813 Mon Sep 17 00:00:00 2001 From: Roland Stigge Date: Sat, 26 Jan 2019 23:19:38 +0100 Subject: Improved Cooley-Tukey --- fft.cpp | 247 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++------- 1 file changed, 221 insertions(+), 26 deletions(-) (limited to 'fft.cpp') diff --git a/fft.cpp b/fft.cpp index 8613dd9..e7d5de6 100644 --- a/fft.cpp +++ b/fft.cpp @@ -1,15 +1,51 @@ +// FFT + +#include +#include +#include +#include +#include +#include #include +#include +#include #include #include -#include -#include -#include const double PI = std::acos(-1); 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); @@ -36,19 +72,19 @@ void generate(std::vector>& v) { } } +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::complex* a, int n) { - std::complex* b = new std::complex[n/2]; // get temp heap storage - for(int i=0; i>::iterator a, int n) { + std::vector> b(n/2); + for(int i=0; i* a, int n) { // 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 fft2 (std::complex* X, int N) { +void fft_recursive(std::vector>::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 - fft2(X, N/2); // recurse even items - fft2(X+N/2, N/2); // recurse odd items + 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 @@ -79,27 +115,186 @@ void fft2 (std::complex* X, int N) { } // Cooley-Tukey -// TODO: Use fft2 via iterators! -std::vector> fft1(const std::vector> &v) { - std::vector> result(v.size(), 0); +std::vector> fft(const std::vector> &v) { + std::vector> result(v); + + fft_recursive(std::begin(result), result.size()); return result; } +}; // namespace CooleyTukey -int main(int argc, char* argv[]) { - std::vector> v(1024, 0); +// FFT Roland Reichwein +class RR { +private: + int mSize; + std::vector order; + std::vector> expLUT; - generate(v); +public: + RR(int size): mSize(size), order(size), expLUT(size/2) { + // reorder LUT + for (int i = 0; i < size; ++i) { + order[i] = bitreverse(i); + } + + // exp LUT + for (int i = 0; i < size / 2; ++i) { + expLUT[i] = exp(std::complex(0,-2.*M_PI*i/size)); + } + } + + std::vector> fft(const std::vector> &v) { + if (v.size() != mSize) + throw std::length_error("Bad input size"); + + std::vector> result; + + reorder(v, result); + fft_recursive(std::begin(result), mSize); + + return result; + } - std::vector> v_dft = dft(v); +// TODO: option for only half FFT due to real redundancy - std::vector> v_fft = fft1(v); +private: + int bitreverse(int i) { + int size{mSize}; + int result{0}; + + while (size > 1) { + result <<= 1; + result |= i & 1; + i >>= 1; + size >>= 1; + } + + return result; + } + + void reorder(const std::vector>& src, std::vector>& dst) { + dst.resize(src.size()); + for (int i = 0; i < mSize; ++i) { + dst[order[i]] = src[i]; + } + } - if (v_dft != v_fft) { - std::cout << "Error: Results differ!\n"; - return 1; + // 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>::iterator X, int N) { + if(N > 2) { + 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 = expLUT[k * mSize / N]; + X[k ] = e + w * o; + X[k+N/2] = e - w * o; + } } + +}; // class RR + +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(mReference->mResult), std::end(mReference->mResult), std::begin(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 { + RR 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.fft(mIn); + } +}; + +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(); + return 0; } -- cgit v1.2.3