From dba25a967e97d9726fa8972f7c007f3723456f93 Mon Sep 17 00:00:00 2001 From: nwolek Date: Wed, 30 Dec 2015 12:36:51 -0500 Subject: [PATCH 1/6] Interpolator: first draft of Allpass option. needs testing. issue #67 --- include/core/JamomaInterpolator.h | 28 ++++++++++++++++++++++++++++ 1 file changed, 28 insertions(+) diff --git a/include/core/JamomaInterpolator.h b/include/core/JamomaInterpolator.h index 2bbc2c5..3751736 100644 --- a/include/core/JamomaInterpolator.h +++ b/include/core/JamomaInterpolator.h @@ -108,6 +108,34 @@ namespace Jamoma { return x1 + delta * (x2-x1); } }; + + /** Allpass interpolation + @param x0 Unused sample value + @param x1 Sample value at prior integer index + @param x2 Sample value at next integer index + @param x3 Unused sample value + @param delta Fractional location between x1 (delta=0) and x1 (delta=1) + @return The interpolated value + */ + template + class Allpass : Base { + public: + static const int delay = 1; + T last_out = T(0.0); + + constexpr T operator()(T x1, T x2, double delta) noexcept { + T out = x1 + delta * (x2-last_out); + last_out = out; + return out; + } + + constexpr T operator()(T x0, T x1, T x2, T x3, double delta) noexcept { + // NW: ideally we would call the operator above to remain DRY, but I could not get syntax right + T out = x1 + delta * (x2-last_out); + last_out = out; + return out; + } + }; /** Cosine interpolation From 1da0b247d6b58c11ced61549720e4dab39095982 Mon Sep 17 00:00:00 2001 From: nwolek Date: Wed, 30 Dec 2015 22:16:48 -0500 Subject: [PATCH 2/6] Interpolator: setting up test output in Octave for Allpass, not happy with results. issue #67 --- test/Interpolator/InterpolatorTargetOutput.m | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/test/Interpolator/InterpolatorTargetOutput.m b/test/Interpolator/InterpolatorTargetOutput.m index 688c4ed..03df712 100644 --- a/test/Interpolator/InterpolatorTargetOutput.m +++ b/test/Interpolator/InterpolatorTargetOutput.m @@ -24,6 +24,8 @@ output_splinegen = double (1 : 64); output_cubicgen = double (1 : 64); output_cosinegen = double (1 : 64); +output_allpassgen = double (1 : 64); +lastout_allpass = 0.0; % the following function is adapted from gen~.interpolation example from Max 7.1 function retval = interp_hermitegen(v,delta) @@ -93,6 +95,17 @@ retval = x + a2*(y-x); endfunction +% reference: https://ccrma.stanford.edu/~jos/pasp/First_Order_Allpass_Interpolation.html +function retval = interp_allpassgen(v,delta,history) + retval = 0.0; + delta_int = fix(delta); + a = delta - delta_int; + x1 = v(delta_int); + x2 = v(delta_int+1); + out = x1 + a*(x2-history); + retval = out; +endfunction + for i = 1:64 current_delta = 2.0 + i / 64; output_linear(i) = interp1(x,current_delta); @@ -102,6 +115,8 @@ output_splinegen(i) = interp_splinegen(x,current_delta); output_cubicgen(i) = interp_cubicgen(x,current_delta); output_cosinegen(i) = interp_cosinegen(x,current_delta); + output_allpassgen(i) = interp_allpassgen(x,current_delta,lastout_allpass); + lastout_allpass = output_allpassgen(i); endfor save expectedOutput.mat output_linear @@ -111,3 +126,4 @@ save -append expectedOutput.mat output_splinegen save -append expectedOutput.mat output_cubicgen save -append expectedOutput.mat output_cosinegen +save -append expectedOutput.mat output_allpassgen From 9d4e13233fb7c2a49925174184094d83942dcd35 Mon Sep 17 00:00:00 2001 From: nwolek Date: Thu, 31 Dec 2015 16:04:29 -0500 Subject: [PATCH 3/6] Interpolator: adding tests of the Allpass option, fails comparison w Octave output. issue #67 --- test/Interpolator/Interpolator.cpp | 115 +++++++++++++++++++++++++++++ 1 file changed, 115 insertions(+) diff --git a/test/Interpolator/Interpolator.cpp b/test/Interpolator/Interpolator.cpp index fec42b3..ba9fe31 100644 --- a/test/Interpolator/Interpolator.cpp +++ b/test/Interpolator/Interpolator.cpp @@ -23,6 +23,7 @@ class InterpolatorTest { testNone(); testNearest(); testLinear(); + testAllpass(); testCosine(); testCubic(); testHermite(); @@ -211,6 +212,120 @@ class InterpolatorTest { mTest->TEST_ASSERT("testLinear overloaded operator produced consistent output", badSampleCount == 0); } + + void testAllpass() { + int badSampleCount = 0; + Jamoma::Interpolator::Allpass my_interp; + + auto x0 = -1.0; + auto x1 = 2.0; + auto x2 = 1.0; + auto x3 = 4.0; + + // The following output was generated using the Octave code + // in InterpolatorTargetOutput.m by NW + Jamoma::SampleVector expectedOutputAllpass = { + 2.015625, + 1.96826171875, + 1.954612731933594, + 1.94033670425415, + 1.926536194980145, + 1.913137231720611, + 1.900125615280558, + 1.88748429808993, + 1.875197520581104, + 1.863250387409203, + 1.851628839664043, + 1.840319592562992, + 1.829310082760642, + 1.818588419396109, + 1.808143339204037, + 1.797964165198991, + 1.788040768619018, + 1.778363533825901, + 1.768923325895436, + 1.759711460657676, + 1.7507196769717, + 1.741940111040978, + 1.733365272594648, + 1.724988022777007, + 1.716801553602732, + 1.70879936884889, + 1.700975266266874, + 1.693323321008243, + 1.68583787016814, + 1.678513498358684, + 1.671345024232512, + 1.664327487883744, + 1.657456139059945, + 1.650726426124404, + 1.644133985713216, + 1.637674633036316, + 1.63134435277588, + 1.625139290539321, + 1.619055744827601, + 1.613090159482749, + 1.607239116581364, + 1.60149932974348, + 1.595867637828599, + 1.590340998992838, + 1.584916485083161, + 1.579591276346478, + 1.574362656433055, + 1.569228007675209, + 1.564184806623668, + 1.559230619825259, + 1.554363099826747, + 1.549579981390768, + 1.54487907791077, + 1.540258278012788, + 1.535715542332761, + 1.531248900458834, + 1.526856448028851, + 1.522536343973854, + 1.518286807899103, + 1.51410611759459, + 1.509992606667656, + 1.505944662290708, + 1.501960723057584, + 1 + }; + + Jamoma::Sample temp = 0.0; + Jamoma::Sample tempExpected = 0.0; + double delta = 0.0; + + for (int i = 0; i < expectedOutputAllpass.size(); i++) { + delta = (i + 1.0) / 64.0; + temp = my_interp(x1,x2,delta); + tempExpected = expectedOutputAllpass[i]; + if (! mTest->compare(temp, tempExpected, true, 8) ) { + badSampleCount++; + std::cout << "sample " << i << " had a difference of " << std::fabs(temp - tempExpected) << std::endl; + } + } + + mTest->TEST_ASSERT("testAllpass produced correct interpolation output", badSampleCount == 0); + + // reset varaiables + badSampleCount = 0; + temp = 0.0; + tempExpected = 0.0; + delta = 0.0; + + for (int i = 0; i < expectedOutputAllpass.size(); i++) { + delta = (i + 1.0) / 64.0; + temp = my_interp(x0,x1,x2,x3,delta); + tempExpected = expectedOutputAllpass[i]; + if (! mTest->compare(temp, tempExpected, true, 8) ) { + badSampleCount++; + std::cout << "sample " << i << " had a difference of " << std::fabs(temp - tempExpected) << std::endl; + } + } + + mTest->TEST_ASSERT("testAllpass overloaded operator produced consistent output", badSampleCount == 0); + } + void testCosine() { int badSampleCount = 0; From 417a367f43c659005bfcf221ee0259acf15503d6 Mon Sep 17 00:00:00 2001 From: nwolek Date: Thu, 31 Dec 2015 16:08:20 -0500 Subject: [PATCH 4/6] Interpolator: adding reset() function to Allpass, fixes inconsistency in overload operator. issue #67 --- include/core/JamomaInterpolator.h | 4 ++++ test/Interpolator/Interpolator.cpp | 1 + 2 files changed, 5 insertions(+) diff --git a/include/core/JamomaInterpolator.h b/include/core/JamomaInterpolator.h index 3751736..e911232 100644 --- a/include/core/JamomaInterpolator.h +++ b/include/core/JamomaInterpolator.h @@ -135,6 +135,10 @@ namespace Jamoma { last_out = out; return out; } + + void reset() { + last_out = T(0.0); + } }; diff --git a/test/Interpolator/Interpolator.cpp b/test/Interpolator/Interpolator.cpp index ba9fe31..2efb49a 100644 --- a/test/Interpolator/Interpolator.cpp +++ b/test/Interpolator/Interpolator.cpp @@ -312,6 +312,7 @@ class InterpolatorTest { temp = 0.0; tempExpected = 0.0; delta = 0.0; + my_interp.reset(); for (int i = 0; i < expectedOutputAllpass.size(); i++) { delta = (i + 1.0) / 64.0; From aeb80f04e0f0f42956e4200c1ebfed161ea18a2e Mon Sep 17 00:00:00 2001 From: nwolek Date: Thu, 31 Dec 2015 17:15:05 -0500 Subject: [PATCH 5/6] Interpolator: correcting for difference in the final test sample, now passes test. leaving multiple comments about the limitations of this algorithm. see issue #67 --- include/core/JamomaInterpolator.h | 7 +++++-- test/Interpolator/Interpolator.cpp | 2 +- test/Interpolator/InterpolatorTargetOutput.m | 5 +++++ 3 files changed, 11 insertions(+), 3 deletions(-) diff --git a/include/core/JamomaInterpolator.h b/include/core/JamomaInterpolator.h index e911232..837925c 100644 --- a/include/core/JamomaInterpolator.h +++ b/include/core/JamomaInterpolator.h @@ -91,7 +91,7 @@ namespace Jamoma { @param x1 Sample value at prior integer index @param x2 Sample value at next integer index @param x3 Unused sample value - @param delta Fractional location between x1 (delta=0) and x1 (delta=1) + @param delta Fractional location between x1 (delta=0) and x2 (delta=1) @return The interpolated value */ template @@ -110,11 +110,14 @@ namespace Jamoma { }; /** Allpass interpolation + Testing shows this algorithm will become less accurate the more points it computes between two known samples. + Also, because it uses an internal history, the reset() function should be used when switching between non-continuous segments of sampled audio data. @param x0 Unused sample value @param x1 Sample value at prior integer index @param x2 Sample value at next integer index @param x3 Unused sample value - @param delta Fractional location between x1 (delta=0) and x1 (delta=1) + @param delta Fractional location between x1 (delta=0) and x2 (delta=1) @n + Be aware that delta=1.0 may not return the exact value at x2 given the nature of this algorithm. @return The interpolated value */ template diff --git a/test/Interpolator/Interpolator.cpp b/test/Interpolator/Interpolator.cpp index 2efb49a..80b31ac 100644 --- a/test/Interpolator/Interpolator.cpp +++ b/test/Interpolator/Interpolator.cpp @@ -288,7 +288,7 @@ class InterpolatorTest { 1.509992606667656, 1.505944662290708, 1.501960723057584, - 1 + 1.498039276942416 }; Jamoma::Sample temp = 0.0; diff --git a/test/Interpolator/InterpolatorTargetOutput.m b/test/Interpolator/InterpolatorTargetOutput.m index 03df712..33f0b65 100644 --- a/test/Interpolator/InterpolatorTargetOutput.m +++ b/test/Interpolator/InterpolatorTargetOutput.m @@ -100,6 +100,11 @@ retval = 0.0; delta_int = fix(delta); a = delta - delta_int; + % the following if statement corrects for a difference between deltas of 0.0 and 1.0 in this algorithm + if (a == 0.0) + delta_int = delta_int - 1; + a = 1.0; + endif x1 = v(delta_int); x2 = v(delta_int+1); out = x1 + a*(x2-history); From b9d2d87732fccc4ba54f0154e71a01507f6b6cdb Mon Sep 17 00:00:00 2001 From: nwolek Date: Fri, 1 Jan 2016 13:36:58 -0500 Subject: [PATCH 6/6] Interpolator:Allpass changing last_out variable to private mY1. issue #67 --- include/core/JamomaInterpolator.h | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/include/core/JamomaInterpolator.h b/include/core/JamomaInterpolator.h index 837925c..d564cee 100644 --- a/include/core/JamomaInterpolator.h +++ b/include/core/JamomaInterpolator.h @@ -124,24 +124,26 @@ namespace Jamoma { class Allpass : Base { public: static const int delay = 1; - T last_out = T(0.0); constexpr T operator()(T x1, T x2, double delta) noexcept { - T out = x1 + delta * (x2-last_out); - last_out = out; + T out = x1 + delta * (x2-mY1); + mY1 = out; return out; } constexpr T operator()(T x0, T x1, T x2, T x3, double delta) noexcept { // NW: ideally we would call the operator above to remain DRY, but I could not get syntax right - T out = x1 + delta * (x2-last_out); - last_out = out; + T out = x1 + delta * (x2-mY1); + mY1 = out; return out; } void reset() { - last_out = T(0.0); + mY1 = T(0.0); } + + private: + T mY1 = T(0.0); };