From 86a4d3cdd784061ad8211bd039e44c35ac42207d Mon Sep 17 00:00:00 2001 From: Steven Atkinson Date: Thu, 4 Apr 2024 01:22:34 -0700 Subject: [PATCH 1/2] Integer math in resampler --- .../Dependencies/LanczosResampler.h | 66 +++++++++++++++---- 1 file changed, 52 insertions(+), 14 deletions(-) diff --git a/dsp/ResamplingContainer/Dependencies/LanczosResampler.h b/dsp/ResamplingContainer/Dependencies/LanczosResampler.h index 9849506..0d8af78 100644 --- a/dsp/ResamplingContainer/Dependencies/LanczosResampler.h +++ b/dsp/ResamplingContainer/Dependencies/LanczosResampler.h @@ -149,11 +149,12 @@ class LanczosResampler */ LanczosResampler(float inputRate, float outputRate) : mInputSampleRate(inputRate) - , mOutputSamplerate(outputRate) - , mPhaseOutIncr(mInputSampleRate / mOutputSamplerate) + , mOutputSampleRate(outputRate) { + SetPhases(); ClearBuffer(); + auto kernel = [](double x) { if (std::fabs(x) < 1e-7) return T(1.0); @@ -200,7 +201,7 @@ class LanczosResampler * Use the fact that mPhaseInIncr = mInputSampleRate and find * res > (A+1) - (mPhaseIn - mPhaseOut + mPhaseOutIncr * desiredOutputs) * sri */ - auto res = A + 1.0 - (mPhaseIn - mPhaseOut - mPhaseOutIncr * nOutputSamples); + double res = A + 1.0 - ((double) (mPhaseInNumerator - mPhaseOutNumerator - mPhaseOutIncrNumerator * (long) nOutputSamples)) / ((double) mPhaseDenominator); return static_cast(std::max(res + 1.0, 0.0)); } @@ -216,17 +217,17 @@ class LanczosResampler } mWritePos = (mWritePos + 1) & (kBufferSize - 1); - mPhaseIn += mPhaseInIncr; + mPhaseInNumerator += mPhaseInIncrNumerator; } } size_t PopBlock(T** outputs, size_t max) { int populated = 0; - while (populated < max && (mPhaseIn - mPhaseOut) > A + 1) + while (populated < max && (mPhaseInNumerator - mPhaseOutNumerator) > mPhaseDenominator * (A + 1)) { - ReadSamples((mPhaseIn - mPhaseOut), outputs, populated); - mPhaseOut += mPhaseOutIncr; + ReadSamples(((double)(mPhaseInNumerator - mPhaseOutNumerator)) / ((double) mPhaseDenominator), outputs, populated); + mPhaseOutNumerator += mPhaseOutIncrNumerator; populated++; } return populated; @@ -234,8 +235,8 @@ class LanczosResampler inline void RenormalizePhases() { - mPhaseIn -= mPhaseOut; - mPhaseOut = 0; + mPhaseInNumerator -= mPhaseOutNumerator; + mPhaseOutNumerator = 0; } void Reset() { ClearBuffer(); } @@ -341,6 +342,40 @@ class LanczosResampler } } #endif + void SetPhases() { + // This is going to assume I can treat the sample rates as longs... + // But if they're not, then things will sound just a little wrong and honestly I'm fine with that. + // It's your fault for not using something normal like 44.1k, 48k, or their multiples. + // (Looking at you, VST3PluginTestHost!) + auto AssertLongLikeSampleRate = [](double x) { + if ((double)((long)x) != x) + { + std::cerr << "Expected long-like sample rate; got " << x << " instead! Truncating..." << std::endl; + } + return (long)x; + }; + + // Greatest common denominator + auto gcd = [](long a, long b) -> long { + while (b != 0) + { + long temp = b; + b = a % b; + a = temp; + } + return a; + }; + + const long inputSampleRate = AssertLongLikeSampleRate(mInputSampleRate); + const long outputSampleRate = AssertLongLikeSampleRate(mOutputSampleRate); + const long g = gcd(inputSampleRate, outputSampleRate); + + // mPhaseInIncr = 1.0 + // mPhaseOutIncr = mInputSampleRate / mOutputSampleRate + mPhaseInIncrNumerator = outputSampleRate / g; + mPhaseDenominator = mPhaseInIncrNumerator; // val / val = 1 + mPhaseOutIncrNumerator = inputSampleRate / g; + }; static T sTable alignas(16)[kTablePoints + 1][kFilterWidth]; static T sDeltaTable alignas(16)[kTablePoints + 1][kFilterWidth]; @@ -349,11 +384,14 @@ class LanczosResampler T mInputBuffer[NCHANS][kBufferSize * 2]; int mWritePos = 0; const float mInputSampleRate; - const float mOutputSamplerate; - double mPhaseIn = 0.0; - double mPhaseOut = 0.0; - double mPhaseInIncr = 1.0; - double mPhaseOutIncr = 0.0; + const float mOutputSampleRate; + // Phase is treated as rational numbers to ensure floating point errors don't accumulate and we stay exactly on. + // (Issue 15) + long mPhaseInNumerator = 0; + long mPhaseOutNumerator = 0; + long mPhaseInIncrNumerator = 1; + long mPhaseOutIncrNumerator = 1; + long mPhaseDenominator = 1; }; template From cd40877793dc696e24c6d6f61fcc62589707e5a1 Mon Sep 17 00:00:00 2001 From: Steven Atkinson Date: Thu, 4 Apr 2024 08:36:05 -0700 Subject: [PATCH 2/2] Formatting --- .../Dependencies/LanczosResampler.h | 73 ++++++++++--------- dsp/ResamplingContainer/ResamplingContainer.h | 8 +- 2 files changed, 43 insertions(+), 38 deletions(-) diff --git a/dsp/ResamplingContainer/Dependencies/LanczosResampler.h b/dsp/ResamplingContainer/Dependencies/LanczosResampler.h index 0d8af78..eaab23c 100644 --- a/dsp/ResamplingContainer/Dependencies/LanczosResampler.h +++ b/dsp/ResamplingContainer/Dependencies/LanczosResampler.h @@ -135,7 +135,7 @@ class LanczosResampler // WARNING: hard-coded to accommodate 8192 samples, from 44.1 to 192k! // If someone is past this, then maybe they know what they're doing ;) // (size_t)(8192 * 3 * 192000 / 44100)+1 = 106998 - static constexpr size_t kBufferSize = 131072; // Round up to pow2. I don't know why, but it doesn't work otherwise. + static constexpr size_t kBufferSize = 131072; // Round up to pow2. I don't know why, but it doesn't work otherwise. // The filter width. 2x because the filter goes from -A to A static constexpr size_t kFilterWidth = A * 2; // The discretization resolution for the filter table. @@ -201,7 +201,9 @@ class LanczosResampler * Use the fact that mPhaseInIncr = mInputSampleRate and find * res > (A+1) - (mPhaseIn - mPhaseOut + mPhaseOutIncr * desiredOutputs) * sri */ - double res = A + 1.0 - ((double) (mPhaseInNumerator - mPhaseOutNumerator - mPhaseOutIncrNumerator * (long) nOutputSamples)) / ((double) mPhaseDenominator); + double res = A + 1.0 + - ((double)(mPhaseInNumerator - mPhaseOutNumerator - mPhaseOutIncrNumerator * (long)nOutputSamples)) + / ((double)mPhaseDenominator); return static_cast(std::max(res + 1.0, 0.0)); } @@ -226,7 +228,7 @@ class LanczosResampler int populated = 0; while (populated < max && (mPhaseInNumerator - mPhaseOutNumerator) > mPhaseDenominator * (A + 1)) { - ReadSamples(((double)(mPhaseInNumerator - mPhaseOutNumerator)) / ((double) mPhaseDenominator), outputs, populated); + ReadSamples(((double)(mPhaseInNumerator - mPhaseOutNumerator)) / ((double)mPhaseDenominator), outputs, populated); mPhaseOutNumerator += mPhaseOutIncrNumerator; populated++; } @@ -342,41 +344,42 @@ class LanczosResampler } } #endif - void SetPhases() { - // This is going to assume I can treat the sample rates as longs... - // But if they're not, then things will sound just a little wrong and honestly I'm fine with that. - // It's your fault for not using something normal like 44.1k, 48k, or their multiples. - // (Looking at you, VST3PluginTestHost!) - auto AssertLongLikeSampleRate = [](double x) { - if ((double)((long)x) != x) - { - std::cerr << "Expected long-like sample rate; got " << x << " instead! Truncating..." << std::endl; - } - return (long)x; - }; + void SetPhases() + { + // This is going to assume I can treat the sample rates as longs... + // But if they're not, then things will sound just a little wrong and honestly I'm fine with that. + // It's your fault for not using something normal like 44.1k, 48k, or their multiples. + // (Looking at you, VST3PluginTestHost!) + auto AssertLongLikeSampleRate = [](double x) { + if ((double)((long)x) != x) + { + std::cerr << "Expected long-like sample rate; got " << x << " instead! Truncating..." << std::endl; + } + return (long)x; + }; - // Greatest common denominator - auto gcd = [](long a, long b) -> long { - while (b != 0) - { - long temp = b; - b = a % b; - a = temp; - } - return a; - }; - - const long inputSampleRate = AssertLongLikeSampleRate(mInputSampleRate); - const long outputSampleRate = AssertLongLikeSampleRate(mOutputSampleRate); - const long g = gcd(inputSampleRate, outputSampleRate); - - // mPhaseInIncr = 1.0 - // mPhaseOutIncr = mInputSampleRate / mOutputSampleRate - mPhaseInIncrNumerator = outputSampleRate / g; - mPhaseDenominator = mPhaseInIncrNumerator; // val / val = 1 - mPhaseOutIncrNumerator = inputSampleRate / g; + // Greatest common denominator + auto gcd = [](long a, long b) -> long { + while (b != 0) + { + long temp = b; + b = a % b; + a = temp; + } + return a; }; + const long inputSampleRate = AssertLongLikeSampleRate(mInputSampleRate); + const long outputSampleRate = AssertLongLikeSampleRate(mOutputSampleRate); + const long g = gcd(inputSampleRate, outputSampleRate); + + // mPhaseInIncr = 1.0 + // mPhaseOutIncr = mInputSampleRate / mOutputSampleRate + mPhaseInIncrNumerator = outputSampleRate / g; + mPhaseDenominator = mPhaseInIncrNumerator; // val / val = 1 + mPhaseOutIncrNumerator = inputSampleRate / g; + }; + static T sTable alignas(16)[kTablePoints + 1][kFilterWidth]; static T sDeltaTable alignas(16)[kTablePoints + 1][kFilterWidth]; static bool sTablesInitialized; diff --git a/dsp/ResamplingContainer/ResamplingContainer.h b/dsp/ResamplingContainer/ResamplingContainer.h index 4c4c359..3099e42 100644 --- a/dsp/ResamplingContainer/ResamplingContainer.h +++ b/dsp/ResamplingContainer/ResamplingContainer.h @@ -193,10 +193,12 @@ class ResamplingContainer { std::cerr << "Did not yield enough samples (" << populated2 << ") to provide the required output buffer (expected" << nFrames << ")! Filling with last sample..." << std::endl; - for (int c = 0; c < NCHANS; c++) { + for (int c = 0; c < NCHANS; c++) + { const T lastSample = populated2 > 0 ? outputs[c][populated2 - 1] : 0.0; - for (int i = populated2; i < nFrames; i++) { - outputs[c][i] = lastSample; + for (int i = populated2; i < nFrames; i++) + { + outputs[c][i] = lastSample; } } }