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 --- fft.cpp | 336 ++++++++++++++++------------------------------------------------ 1 file changed, 80 insertions(+), 256 deletions(-) (limited to 'fft.cpp') diff --git a/fft.cpp b/fft.cpp index e7d5de6..29fd405 100644 --- a/fft.cpp +++ b/fft.cpp @@ -1,6 +1,8 @@ // FFT -#include +#include "fft.h" + +#include #include #include #include @@ -10,291 +12,113 @@ #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; +namespace { // Helper functions + bool is_power_of_two(unsigned int n) { + return n != 0 && (n & (n - 1)) == 0; } -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. * PI * double(k * j) / N); - } - } +} +std::vector RIT::magnitudes(std::vector>& v) { + std::vector result(v.size()); + std::transform(std::begin(v), std::end(v), std::begin(result), std::abs); return result; } -std::vector freqs{ 2, 5, 11, 17, 29 }; // known freqs for testing +RIT::FFT::FFT(int size, bool halfOnly): mSize(size), order(size), expLUT(size/2), mFlagHalfOnly(halfOnly) { -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() ); - } -} + if (!is_power_of_two(size)) + throw std::invalid_argument("Size must be a power of two"); -namespace CooleyTukey { + // reorder LUT + for (int i = 0; i < size; ++i) { + order[i] = bitreverse(i); + } -// 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(0,-2.*M_PI*i/size)); + } } -// 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) { - // 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); +std::vector> RIT::FFT::operator()(const std::vector> &v) { + if (v.size() != mSize) + throw std::length_error("Bad input size"); + + std::vector> result; - fft_recursive(std::begin(result), result.size()); + reorder(v, result); + if (mFlagHalfOnly) { + fft_half(std::begin(result), mSize); + result.resize(mSize/2); + } else + fft_recursive(std::begin(result), mSize); return result; } -}; // namespace CooleyTukey - -// FFT Roland Reichwein -class RR { -private: - int mSize; - std::vector order; - std::vector> expLUT; - -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"); +RIT::FFT& RIT::FFT::SetHalfOnly(bool enable) { + mFlagHalfOnly = enable; + return *this; +} - std::vector> result; +int RIT::FFT::bitreverse(int i) { + int size{mSize}; + int result{0}; - reorder(v, result); - fft_recursive(std::begin(result), mSize); - - return result; + while (size > 1) { + result <<= 1; + result |= i & 1; + i >>= 1; + size >>= 1; } -// TODO: option for only half FFT due to real redundancy - -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; - } + 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]; - } - } +void RIT::FFT::reorder(const std::vector>& src, std::vector>& dst) { + int size = src.size(); - // 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; - } + dst.resize(size); + for (int i = 0; i < size; ++i) { + dst[order[i]] = src[i]; } +} -}; // 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"; - } - } +// 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 RIT::FFT::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 } - - 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); + // 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 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); +// Same as fft_recursive, but results in only half the result due to symmetry in real input +void RIT::FFT::fft_half(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 } -}; - -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); + // combine results of two half recursions + for(int k=0; k> 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