Add fft.hpp, clean up math and dsp headers

include/common.hpp

@@ -15,6 +15,9 @@
 #include <string>
 
 
+namespace rack {
+
+
 /** Deprecation notice for functions
 E.g.
 	DEPRECATED void foo();
@@ -130,4 +133,7 @@ DeferWrapper<F> deferWrapper(F f) {
 	return DeferWrapper<F>(f);
 }
 
-#define DEFER(code) auto CONCAT(_defer_, __COUNTER__) = deferWrapper([&]() code)
+#define DEFER(code) auto CONCAT(_defer_, __COUNTER__) = rack::deferWrapper([&]() code)
+
+
+} // namespace rack

include/dsp/common.hpp

@@ -6,11 +6,15 @@ namespace rack {
 namespace dsp {
 
 
+// Useful constants
+
 static const float FREQ_C4 = 261.6256f;
 static const float FREQ_A4 = 440.0000f;
 
 
-/** The normalized sinc function. https://en.wikipedia.org/wiki/Sinc_function */
+/** The normalized sinc function
+See https://en.wikipedia.org/wiki/Sinc_function
+*/
 inline float sinc(float x) {
 	if (x == 0.f)
 		return 1.f;
@@ -18,6 +22,8 @@ inline float sinc(float x) {
 	return std::sin(x) / x;
 }
 
+// Functions for parameter scaling
+
 inline float quadraticBipolar(float x) {
 	float x2 = x*x;
 	return (x >= 0.f) ? x2 : -x2;
@@ -33,7 +39,7 @@ inline float quarticBipolar(float x) {
 }
 
 inline float quintic(float x) {
-	// optimal with --fast-math
+	// optimal with -fassociative-math
 	return x*x*x*x*x;
 }
 
@@ -47,6 +53,8 @@ inline float exponentialBipolar(float b, float x) {
 	return (std::pow(b, x) - std::pow(b, -x)) / a;
 }
 
+// Useful conversion functions
+
 inline float amplitudeToDb(float amp) {
 	return std::log10(amp) * 20.f;
 }

include/dsp/fft.hpp

@@ -0,0 +1,134 @@
+#pragma once
+#include "dsp/common.hpp"
+#include <pffft.h>
+
+
+namespace rack {
+namespace dsp {
+
+
+template<typename T>
+T *alignedNew(size_t len) {
+	return pffft_aligned_malloc(len * sizeof(T));
+}
+
+template<typename T>
+void alignedDelete(T *p) {
+	pffft_aligned_free(p);
+}
+
+
+/** Real-valued FFT context
+Wraps PFFFT (https://bitbucket.org/jpommier/pffft/)
+*/
+struct RealFFT {
+	PFFFT_Setup *setup;
+	int length;
+
+	RealFFT(size_t length) {
+		this->length = length;
+		setup = pffft_new_setup(length, PFFFT_REAL);
+	}
+
+	~RealFFT() {
+		pffft_destroy_setup(setup);
+	}
+
+	/** Performs the real FFT.
+	Input and output must be aligned using the above align*() functions.
+	Input is `length` elements. Output is `2*length` elements.
+	Output is arbitrarily ordered for performance reasons.
+	However, this ordering is consistent, so element-wise multiplication with line up with other results, and the inverse FFT will return a correctly ordered result.
+	*/
+	void rfftUnordered(const float *input, float *output) {
+		pffft_transform(setup, input, output, NULL, PFFFT_FORWARD);
+	}
+
+	/** Performs the inverse real FFT.
+	Input is `2*length` elements. Output is `length` elements.
+	Scaling is such that IRFFT(RFFT(x)) = N*x.
+	*/
+	void irfftUnordered(const float *input, float *output) {
+		pffft_transform(setup, input, output, NULL, PFFFT_BACKWARD);
+	}
+
+	/** Slower than the above methods, but returns results in the "canonical" FFT order as follows.
+		output[0] = F(0)
+		output[1] = F(n/2)
+		output[2] = real(F(1))
+		output[3] = imag(F(1))
+		output[4] = real(F(2))
+		output[5] = imag(F(2))
+		...
+		output[length - 2] = real(F(n/2 - 1))
+		output[length - 1] = imag(F(n/2 - 1))
+	*/
+	void rfft(const float *input, float *output) {
+		pffft_transform_ordered(setup, input, output, NULL, PFFFT_FORWARD);
+	}
+
+	void irfft(const float *input, float *output) {
+		pffft_transform_ordered(setup, input, output, NULL, PFFFT_BACKWARD);
+	}
+
+	/** Scales the RFFT so that
+	scale(IFFT(FFT(x))) = x
+	*/
+	void scale(float *x) {
+		float a = 1.f / length;
+		for (int i = 0; i < length; i++) {
+			x[i] *= a;
+		}
+	}
+};
+
+
+struct ComplexFFT {
+	PFFFT_Setup *setup;
+	int length;
+
+	ComplexFFT(size_t length) {
+		this->length = length;
+		setup = pffft_new_setup(length, PFFFT_COMPLEX);
+	}
+
+	~ComplexFFT() {
+		pffft_destroy_setup(setup);
+	}
+
+	/** Performs the complex FFT.
+	Input and output must be aligned using the above align*() functions.
+	Input is `2*length` elements. Output is `2*length` elements.
+	*/
+	void fftUnordered(const float *input, float *output) {
+		pffft_transform(setup, input, output, NULL, PFFFT_FORWARD);
+	}
+
+	/** Performs the inverse complex FFT.
+	Input is `2*length` elements. Output is `2*length` elements.
+	Scaling is such that FFT(IFFT(x)) = N*x.
+	*/
+	void ifftUnordered(const float *input, float *output) {
+		pffft_transform(setup, input, output, NULL, PFFFT_BACKWARD);
+	}
+
+	void fft(const float *input, float *output) {
+		pffft_transform_ordered(setup, input, output, NULL, PFFFT_FORWARD);
+	}
+
+	void ifft(const float *input, float *output) {
+		pffft_transform_ordered(setup, input, output, NULL, PFFFT_BACKWARD);
+	}
+
+	void scale(float *x) {
+		float a = 1.f / length;
+		for (int i = 0; i < length; i++) {
+			x[2*i+0] *= a;
+			x[2*i+1] *= a;
+		}
+	}
+};
+
+
+} // namespace dsp
+} // namespace rack

include/dsp/fir.hpp

@@ -24,7 +24,7 @@ inline void boxcarLowpassIR(float *out, int len, float cutoff = 0.5f) {
 	}
 }
 
-inline void blackmanHarrisWindow(float *x, int len) {
+inline void blackmanHarris(float *x, int len) {
 	// Constants from https://en.wikipedia.org/wiki/Window_function#Blackman%E2%80%93Harris_window
 	const float a0 = 0.35875f;
 	const float a1 = 0.48829f;
@@ -34,9 +34,9 @@ inline void blackmanHarrisWindow(float *x, int len) {
 	for (int i = 0; i < len; i++) {
 		x[i] *=
 			+ a0
-			- a1 * cosf(1*factor * i)
-			+ a2 * cosf(2*factor * i)
-			- a3 * cosf(3*factor * i);
+			- a1 * std::cos(1*factor * i)
+			+ a2 * std::cos(2*factor * i)
+			- a3 * std::cos(3*factor * i);
 	}
 }
 

include/dsp/minblep.hpp

@@ -6,31 +6,28 @@ namespace rack {
 namespace dsp {
 
 
-// Pre-made minBLEP samples in minBLEP.cpp
-extern const float minblep_16_32[];
-
-
 template<int ZERO_CROSSINGS>
 struct MinBLEP {
-	float buf[2*ZERO_CROSSINGS] = {};
+	float buf[2 * ZERO_CROSSINGS] = {};
 	int pos = 0;
 	const float *minblep;
 	int oversample;
 
-	/** Places a discontinuity with magnitude dx at -1 < p <= 0 relative to the current frame */
-	void jump(float p, float dx) {
-		if (p <= -1 || 0 < p)
+	/** Places a discontinuity with magnitude `x` at -1 < p <= 0 relative to the current frame */
+	void insertDiscontinuity(float p, float x) {
+		if (!(-1 < p && p <= 0))
 			return;
-		for (int j = 0; j < 2*ZERO_CROSSINGS; j++) {
+		for (int j = 0; j < 2 * ZERO_CROSSINGS; j++) {
 			float minblepIndex = ((float)j - p) * oversample;
-			int index = (pos + j) % (2*ZERO_CROSSINGS);
-			buf[index] += dx * (-1.0 + math::interpolateLinear(minblep, minblepIndex));
+			int index = (pos + j) % (2 * ZERO_CROSSINGS);
+			buf[index] += x * (-1.f + math::interpolateLinear(minblep, minblepIndex));
 		}
 	}
-	float shift() {
+
+	float process() {
 		float v = buf[pos];
-		buf[pos] = 0.0;
-		pos = (pos + 1) % (2*ZERO_CROSSINGS);
+		buf[pos] = 0.f;
+		pos = (pos + 1) % (2 * ZERO_CROSSINGS);
 		return v;
 	}
 };

include/dsp/resampler.hpp

@@ -109,7 +109,7 @@ struct Decimator {
 
 	Decimator(float cutoff = 0.9f) {
 		boxcarLowpassIR(kernel, OVERSAMPLE*QUALITY, cutoff * 0.5f / OVERSAMPLE);
-		blackmanHarrisWindow(kernel, OVERSAMPLE*QUALITY);
+		blackmanHarris(kernel, OVERSAMPLE*QUALITY);
 		reset();
 	}
 	void reset() {

include/logger.hpp

@@ -1,17 +1,6 @@
 #pragma once
 
 
-/** Example usage:
-	DEBUG("error: %d", errno);
-will print something like
-	[0.123 debug myfile.cpp:45] error: 67
-*/
-#define DEBUG(format, ...) rack::logger::log(rack::logger::DEBUG_LEVEL, __FILE__, __LINE__, format, ##__VA_ARGS__)
-#define INFO(format, ...) rack::logger::log(rack::logger::INFO_LEVEL, __FILE__, __LINE__, format, ##__VA_ARGS__)
-#define WARN(format, ...) rack::logger::log(rack::logger::WARN_LEVEL, __FILE__, __LINE__, format, ##__VA_ARGS__)
-#define FATAL(format, ...) rack::logger::log(rack::logger::FATAL_LEVEL, __FILE__, __LINE__, format, ##__VA_ARGS__)
-
-
 namespace rack {
 namespace logger {
 
@@ -29,5 +18,16 @@ void destroy();
 void log(Level level, const char *filename, int line, const char *format, ...);
 
 
+/** Example usage:
+	DEBUG("error: %d", errno);
+will print something like
+	[0.123 debug myfile.cpp:45] error: 67
+*/
+#define DEBUG(format, ...) rack::logger::log(rack::logger::DEBUG_LEVEL, __FILE__, __LINE__, format, ##__VA_ARGS__)
+#define INFO(format, ...) rack::logger::log(rack::logger::INFO_LEVEL, __FILE__, __LINE__, format, ##__VA_ARGS__)
+#define WARN(format, ...) rack::logger::log(rack::logger::WARN_LEVEL, __FILE__, __LINE__, format, ##__VA_ARGS__)
+#define FATAL(format, ...) rack::logger::log(rack::logger::FATAL_LEVEL, __FILE__, __LINE__, format, ##__VA_ARGS__)
+
+
 } // namespace logger
 } // namespace rack

include/math.hpp

@@ -119,15 +119,15 @@ See https://en.wikipedia.org/wiki/Euclidean_division
 */
 inline float eucMod(float a, float base) {
 	float mod = std::fmod(a, base);
-	return (mod >= 0.0f) ? mod : mod + base;
+	return (mod >= 0.f) ? mod : mod + base;
 }
 
-inline bool isNear(float a, float b, float epsilon = 1.0e-6f) {
+inline bool isNear(float a, float b, float epsilon = 1e-6f) {
 	return std::abs(a - b) <= epsilon;
 }
 
 /** If the magnitude of x if less than epsilon, return 0 */
-inline float chop(float x, float epsilon = 1.0e-6f) {
+inline float chop(float x, float epsilon = 1e-6f) {
 	return isNear(x, 0.f, epsilon) ? 0.f : x;
 }
 
@@ -157,7 +157,6 @@ inline void complexMult(float *cr, float *ci, float ar, float ai, float br, floa
 	*ci = ar * bi + ai * br;
 }
 
-
 ////////////////////
 // 2D vector and rectangle
 ////////////////////
@@ -232,7 +231,7 @@ struct Vec {
 		return x == b.x && y == b.y;
 	}
 	bool isZero() const {
-		return x == 0.0f && y == 0.0f;
+		return x == 0.f && y == 0.f;
 	}
 	bool isFinite() const {
 		return std::isfinite(x) && std::isfinite(y);

include/rack.hpp

@@ -79,6 +79,7 @@
 
 #include "dsp/common.hpp"
 #include "dsp/digital.hpp"
+#include "dsp/fft.hpp"
 #include "dsp/filter.hpp"
 #include "dsp/fir.hpp"
 #include "dsp/frame.hpp"