added trigonometric approximations

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<double(unsigned int)> WindowFunction;
+typedef std::function<double(double)> TrigFunction;
+typedef std::function<std::complex<double>(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<double> ComplexExp(double x);
+
 std::vector<std::complex<double>> 
 radix2dit(
 	const std::vector<double>& list,
@@ -38,12 +55,12 @@ radix2dit(
 		std::vector<std::complex<double>> first = radix2dit(list, offset, halfN, s << 1);
 		std::vector<std::complex<double>> second = radix2dit(list, offset + s, halfN, s << 1);
 
-		std::complex<double> coeff = -M_PI * 1.0i / (double)halfN;
+		double coeff = -M_PI / (double)halfN;
 
 		for (int k = 0; k < halfN; k++)
 		{
 			std::complex<double> p = first[k];
-			std::complex<double> q = std::exp(coeff * (double)k) * second[k];
+			std::complex<double> q = ComplexExp(coeff * (double)k) * second[k];
 
 			output[k] = p + q;
 			output[halfN + k] = p - q;
@@ -76,7 +93,7 @@ FFT(const std::vector<double>::const_iterator& begin,
 
 	WindowFunction f;
 	
-
+	
 
 	std::vector<std::complex<double>> 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<double> ComplexExp(double x)
+{
+	return std::complex<double>(Cos(x), Sin(x));
 }

src/FFT.hpp

@@ -16,4 +16,5 @@ extern std::vector<std::pair<double, double>> FFT(const std::vector<double>::con
 	double minFreq, double maxFreq,
 	unsigned int zeropadding);
 
-extern void SetWindowFunction(WindowFunctions func, unsigned int width);
+extern void SetWindowFunction(WindowFunctions func, unsigned int width);
+extern void UseFastFunctions();

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<double> 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<unsigned int>())
 			("w,window", "Specify the window function used (rectangle (default), von-hann, gauss, triangle, blackman (3-term))", cxxopts::value<std::string>()->default_value("rectangle"))
 			("m,mono", "Analyze only the given channel", cxxopts::value<unsigned int>()->default_value("0"))
+			("approx", "Use faster, but more inaccurate trigonometric functions instead of the std-functions (EXPERIMENTAL)")
 			("files", "Files to fourier transform", cxxopts::value<std::vector<std::filesystem::path>>())
 			("h,help", "Print usage")
 			;
@@ -176,6 +181,7 @@ Settings Parse(int argc, char** argv)
 		setts.splitInterval = (result.count("interval") ? result["interval"].as<float>() : 0.0f);
 		setts.analyzeChannel = (result.count("mono") ? result["mono"].as<unsigned int>() : 0);
 		setts.zeropadding = (result.count("pad") ? result["pad"].as<unsigned int>() : 1);
+		setts.approx = (result.count("approx") ? true : false);
 
 		if (!result.count("window"))
 		{