diff --git a/include/dsp/common.hpp b/include/dsp/common.hpp index fb96b68c..c54adb4f 100644 --- a/include/dsp/common.hpp +++ b/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 diff --git a/include/dsp/fft.hpp b/include/dsp/fft.hpp index 513b2d39..89de5aa0 100644 --- a/include/dsp/fft.hpp +++ b/include/dsp/fft.hpp @@ -9,7 +9,7 @@ namespace dsp { template T *alignedNew(size_t len) { - return pffft_aligned_malloc(len * sizeof(T)); + return (T*) pffft_aligned_malloc(len * sizeof(T)); } template diff --git a/include/dsp/fir.hpp b/include/dsp/fir.hpp index 6e3a5e92..77c7e883 100644 --- a/include/dsp/fir.hpp +++ b/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` diff --git a/include/dsp/minblep.hpp b/include/dsp/minblep.hpp index 884a37c1..ea4a126d 100644 --- a/include/dsp/minblep.hpp +++ b/include/dsp/minblep.hpp @@ -6,28 +6,41 @@ namespace rack { namespace dsp { -template -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 +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; } }; diff --git a/include/dsp/resampler.hpp b/include/dsp/resampler.hpp index bc3b8bf0..03c27798 100644 --- a/include/dsp/resampler.hpp +++ b/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() {