Add window functions to dsp. Implement minBlepImpulse().

include/dsp/common.hpp

@@ -6,11 +6,12 @@ namespace rack {
 namespace dsp {
 
 
-// Useful constants
+// Constants
 
 static const float FREQ_C4 = 261.6256f;
 static const float FREQ_A4 = 440.0000f;
 
+// Mathematical functions
 
 /** The normalized sinc function
 See https://en.wikipedia.org/wiki/Sinc_function
@@ -22,6 +23,84 @@ inline float sinc(float x) {
 	return std::sin(x) / x;
 }
 
+// Window functions
+
+/** Hann window function
+p: proportion from [0, 1], usually `i / (len - 1)`
+https://en.wikipedia.org/wiki/Window_function#Hann_and_Hamming_windows
+*/
+inline float hann(float p) {
+	return 0.5f * (1.f - std::cos(2*M_PI * p));
+}
+
+inline void hannWindow(float *x, int len) {
+	for (int i = 0; i < len; i++) {
+		x[i] *= hann((float) i / (len - 1));
+	}
+}
+
+/** Blackman window function
+https://en.wikipedia.org/wiki/Window_function#Blackman_window
+A typical alpha value is 0.16.
+*/
+inline float blackman(float alpha, float p) {
+	return
+		+ (1 - alpha) / 2.f
+		- 1 / 2.f * std::cos(2*M_PI * p)
+		+ alpha / 2.f * std::cos(4*M_PI * p);
+}
+
+inline void blackmanWindow(float alpha, float *x, int len) {
+	for (int i = 0; i < len; i++) {
+		x[i] *= blackman(alpha, (float) i / (len - 1));
+	}
+}
+
+
+/** Blackman-Nuttall window function
+https://en.wikipedia.org/wiki/Window_function#Blackman%E2%80%93Nuttall_window
+*/
+inline float blackmanNuttall(float p) {
+	return
+		+ 0.3635819f
+		- 0.4891775f * std::cos(2*M_PI * p)
+		+ 0.1365995f * std::cos(4*M_PI * p)
+		- 0.0106411f * std::cos(6*M_PI * p);
+}
+
+inline void blackmanNuttallWindow(float *x, int len) {
+	for (int i = 0; i < len; i++) {
+		x[i] *= blackmanNuttall((float) i / (len - 1));
+	}
+}
+
+/** Blackman-Harris window function
+https://en.wikipedia.org/wiki/Window_function#Blackman%E2%80%93Harris_window
+*/
+inline float blackmanHarris(float p) {
+	return
+		+ 0.35875f
+		- 0.48829f * std::cos(2*M_PI * p)
+		+ 0.14128f * std::cos(4*M_PI * p)
+		- 0.01168f * std::cos(6*M_PI * p);
+}
+
+inline void blackmanHarrisWindow(float *x, int len) {
+	for (int i = 0; i < len; i++) {
+		x[i] *= blackmanHarris((float) i / (len - 1));
+	}
+}
+
+// Conversion functions
+
+inline float amplitudeToDb(float amp) {
+	return std::log10(amp) * 20.f;
+}
+
+inline float dbToAmplitude(float db) {
+	return std::pow(10.f, db / 20.f);
+}
+
 // Functions for parameter scaling
 
 inline float quadraticBipolar(float x) {
@@ -53,16 +132,6 @@ 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;
-}
-
-inline float dbToAmplitude(float db) {
-	return std::pow(10.f, db / 20.f);
-}
-
 
 } // namespace dsp
 } // namespace rack

include/dsp/fft.hpp

@@ -9,7 +9,7 @@ namespace dsp {
 
 template<typename T>
 T *alignedNew(size_t len) {
-	return pffft_aligned_malloc(len * sizeof(T));
+	return (T*) pffft_aligned_malloc(len * sizeof(T));
 }
 
 template<typename T>

include/dsp/fir.hpp

@@ -24,22 +24,6 @@ inline void boxcarLowpassIR(float *out, int len, float cutoff = 0.5f) {
 	}
 }
 
-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;
-	const float a2 = 0.14128f;
-	const float a3 = 0.01168f;
-	float factor = 2*M_PI / (len - 1);
-	for (int i = 0; i < len; i++) {
-		x[i] *=
-			+ a0
-			- a1 * std::cos(1*factor * i)
-			+ a2 * std::cos(2*factor * i)
-			- a3 * std::cos(3*factor * i);
-	}
-}
-
 
 struct RealTimeConvolver {
 	// `kernelBlocks` number of contiguous FFT blocks of size `blockSize`

include/dsp/minblep.hpp

@@ -6,28 +6,41 @@ namespace rack {
 namespace dsp {
 
 
-template<int ZERO_CROSSINGS>
-struct MinBLEP {
-	float buf[2 * ZERO_CROSSINGS] = {};
+/** Computes the minimum-phase bandlimited step (MinBLEP)
+z: number of zero-crossings
+o: oversample factor
+output: must be length `2 * z * o`.
+https://www.cs.cmu.edu/~eli/papers/icmc01-hardsync.pdf
+*/
+void minBlepImpulse(int z, int o, float *output);
+
+
+template<int Z, int O>
+struct MinBlepGenerator {
+	float buf[2 * Z] = {};
 	int pos = 0;
-	const float *minblep;
-	int oversample;
+	float impulse[2 * Z * O + 1];
+
+	MinBlepGenerator() {
+		minBlepImpulse(Z, O, impulse);
+		impulse[2 * Z * O] = 1.f;
+	}
 
 	/** 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++) {
-			float minblepIndex = ((float)j - p) * oversample;
-			int index = (pos + j) % (2 * ZERO_CROSSINGS);
-			buf[index] += x * (-1.f + math::interpolateLinear(minblep, minblepIndex));
+		for (int j = 0; j < 2 * Z; j++) {
+			float minBlepIndex = ((float)j - p) * O;
+			int index = (pos + j) % (2 * Z);
+			buf[index] += x * (-1.f + math::interpolateLinear(impulse, minBlepIndex));
 		}
 	}
 
 	float process() {
 		float v = buf[pos];
 		buf[pos] = 0.f;
-		pos = (pos + 1) % (2 * ZERO_CROSSINGS);
+		pos = (pos + 1) % (2 * Z);
 		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);
-		blackmanHarris(kernel, OVERSAMPLE*QUALITY);
+		blackmanHarrisWindow(kernel, OVERSAMPLE*QUALITY);
 		reset();
 	}
 	void reset() {