From bb947b4a53d61bfbd2676ce4d9aee2654ecd2081 Mon Sep 17 00:00:00 2001 From: Robert Date: Fri, 30 Apr 2021 13:25:09 +0200 Subject: [PATCH] added trigonometric approximations --- src/FFT.cpp | 60 ++++++++++++++++++++++++++++++++++++++++++++++------ src/FFT.hpp | 3 ++- src/main.cpp | 6 ++++++ 3 files changed, 62 insertions(+), 7 deletions(-) diff --git a/src/FFT.cpp b/src/FFT.cpp index 6291fb8..f3895ef 100644 --- a/src/FFT.cpp +++ b/src/FFT.cpp @@ -8,11 +8,24 @@ #define POW_OF_TWO(x) (x && !(x & (x - 1))) +constexpr double REC_2_FAC = (double)1.0f / (double)2.0f; +constexpr double REC_3_FAC = (double)1.0f / (double)6.0f; +constexpr double REC_4_FAC = (double)1.0f / (double)24.0f; +constexpr double REC_5_FAC = (double)1.0f / (double)120.0f; +constexpr double REC_6_FAC = (double)1.0f / (double)720.0f; +constexpr double REC_7_FAC = (double)1.0f / (double)5040.0f; +constexpr double REC_8_FAC = (double)1.0f / (double)40320.0f; +constexpr double REC_9_FAC = (double)1.0f / (double)362880.0f; + using namespace std::complex_literals; typedef std::function WindowFunction; +typedef std::function TrigFunction; +typedef std::function(double)> ExpFunction; -static WindowFunction window; +WindowFunction window; +TrigFunction Sin = std::bind((double(*)(double))& std::sin, std::placeholders::_1); +TrigFunction Cos = std::bind((double(*)(double))& std::cos, std::placeholders::_1); inline double WindowRectangle(unsigned int k, unsigned int offset, unsigned int width); inline double WindowVonHann(unsigned int k, unsigned int offset, unsigned int width); @@ -20,6 +33,10 @@ inline double WindowGauss(unsigned int k, unsigned int offset, unsigned int widt inline double WindowTriangle(unsigned int k, unsigned int offset, unsigned int width); inline double WindowBlackman(unsigned int k, unsigned int offset, unsigned int width); +double FastCos(double x); +double FastSin(double x); +std::complex ComplexExp(double x); + std::vector> radix2dit( const std::vector& list, @@ -38,12 +55,12 @@ radix2dit( std::vector> first = radix2dit(list, offset, halfN, s << 1); std::vector> second = radix2dit(list, offset + s, halfN, s << 1); - std::complex coeff = -M_PI * 1.0i / (double)halfN; + double coeff = -M_PI / (double)halfN; for (int k = 0; k < halfN; k++) { std::complex p = first[k]; - std::complex q = std::exp(coeff * (double)k) * second[k]; + std::complex q = ComplexExp(coeff * (double)k) * second[k]; output[k] = p + q; output[halfN + k] = p - q; @@ -76,7 +93,7 @@ FFT(const std::vector::const_iterator& begin, WindowFunction f; - + std::vector> spectrum = radix2dit(signal, 0, N, 1); double freqRes = (double)sampleRate / (double)N; @@ -109,6 +126,12 @@ void SetWindowFunction(WindowFunctions func, unsigned int width) } } +void UseFastFunctions() +{ + Sin = std::bind(FastSin, std::placeholders::_1); + Cos = std::bind(FastCos, std::placeholders::_1); +} + inline double WindowRectangle(unsigned int k, unsigned int offset, unsigned int width) { return ((offset < k) && (k < width)); @@ -116,7 +139,7 @@ inline double WindowRectangle(unsigned int k, unsigned int offset, unsigned int inline double WindowVonHann(unsigned int k, unsigned int offset, unsigned int width) { - return ((offset < k) && (k < width)) ? (0.5f * (1.0f - cos(2.0f * M_PI * k / (width - 1)))) : 0; + return ((offset < k) && (k < width)) ? (0.5f * (1.0f - Cos(2.0f * M_PI * k / (width - 1)))) : 0; } inline double WindowGauss(unsigned int k, unsigned int offset, unsigned int width) @@ -132,5 +155,30 @@ inline double WindowTriangle(unsigned int k, unsigned int offset, unsigned int w inline double WindowBlackman(unsigned int k, unsigned int offset, unsigned int width) { - return 0.5f * (1.0f - 0.16f) - 0.5f * cos(2.0f * M_PI * k / (width - 1)) + 0.5f * 0.16f * cos(4.0f * M_PI * k / (width - 1)); + return (double)0.5f * ((double)1.0f - (double)0.16f) - 0.5f * Cos(2.0f * M_PI * k / (width - 1)) + (double)0.5f * (double)0.16f * Cos(4.0f * M_PI * k / (width - 1)); +} + +double FastCos(double x) +{ + x -= (x > M_PI) * (double)2.0f * M_PI; + x += (x < -M_PI) * (double)2.0f * M_PI; + double xpow2 = x * x; + double xpow4 = xpow2 * x * x; + double xpow6 = xpow4 * x * x; + return (double)1.0f - xpow2 * REC_2_FAC + xpow4 * REC_4_FAC - xpow6 * REC_6_FAC + xpow6 * x * x * REC_8_FAC; +} + +double FastSin(double x) +{ + x -= (x > M_PI) * (double)2.0f * M_PI; + x += (x < -M_PI) * (double)2.0f * M_PI; + double xpow3 = x * x * x; + double xpow5 = xpow3 * x * x; + double xpow7 = xpow5 * x * x; + return (double)x - xpow3 * REC_3_FAC + xpow5 * REC_5_FAC - xpow7 * REC_7_FAC + xpow7 * x * x * REC_9_FAC; +} + +std::complex ComplexExp(double x) +{ + return std::complex(Cos(x), Sin(x)); } diff --git a/src/FFT.hpp b/src/FFT.hpp index 0a04842..6417fd6 100644 --- a/src/FFT.hpp +++ b/src/FFT.hpp @@ -16,4 +16,5 @@ extern std::vector> FFT(const std::vector::con double minFreq, double maxFreq, unsigned int zeropadding); -extern void SetWindowFunction(WindowFunctions func, unsigned int width); \ No newline at end of file +extern void SetWindowFunction(WindowFunctions func, unsigned int width); +extern void UseFastFunctions(); \ No newline at end of file diff --git a/src/main.cpp b/src/main.cpp index 8daff5a..6fcd617 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -26,6 +26,7 @@ struct Settings { double minFreq, maxFreq; unsigned int analyzeChannel; unsigned int zeropadding; + bool approx; WindowFunctions window; }; @@ -36,6 +37,9 @@ int main(int argc, char** argv) Settings setts; setts = Parse(argc, argv); + if (setts.approx) + UseFastFunctions(); + int numFiles = setts.files.size(); for (auto& file : setts.files) { AudioFile audioFile; @@ -141,6 +145,7 @@ Settings Parse(int argc, char** argv) ("p,pad", "Add extra zero-padding. By default, the program will pad the signals with 0s until the number of samples is a power of 2 (this would be equivalent to -p 1). With this option you can tell the program to instead pad until the power of 2 after the next one (-p 2) etc. This increases frequency resolution", cxxopts::value()) ("w,window", "Specify the window function used (rectangle (default), von-hann, gauss, triangle, blackman (3-term))", cxxopts::value()->default_value("rectangle")) ("m,mono", "Analyze only the given channel", cxxopts::value()->default_value("0")) + ("approx", "Use faster, but more inaccurate trigonometric functions instead of the std-functions (EXPERIMENTAL)") ("files", "Files to fourier transform", cxxopts::value>()) ("h,help", "Print usage") ; @@ -176,6 +181,7 @@ Settings Parse(int argc, char** argv) setts.splitInterval = (result.count("interval") ? result["interval"].as() : 0.0f); setts.analyzeChannel = (result.count("mono") ? result["mono"].as() : 0); setts.zeropadding = (result.count("pad") ? result["pad"].as() : 1); + setts.approx = (result.count("approx") ? true : false); if (!result.count("window")) {