Merge pull request #76 from jamoma/issue/67

Closing Issue/67

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<class T>
@@ -108,6 +108,43 @@ namespace Jamoma {
                 return x1 + delta * (x2-x1);
             }
 		};
+
+        /** 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 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<class T>
+        class Allpass : Base {
+        public:
+            static const int 	delay = 1;
+
+            constexpr T operator()(T x1, T x2, double delta) noexcept {
+                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-mY1);
+                mY1 = out;
+                return out;
+            }
+
+            void reset() {
+                mY1 = T(0.0);
+            }
+
+        private:
+            T mY1 = T(0.0);
+        };
 
 
 		/** Cosine interpolation

test/Interpolator/Interpolator.cpp

@@ -23,6 +23,7 @@ class InterpolatorTest {
         testNone();
         testNearest();
         testLinear();
+        testAllpass();
 		testCosine();
 		testCubic();
         testHermite();
@@ -211,6 +212,121 @@ class InterpolatorTest {
         mTest->TEST_ASSERT("testLinear overloaded operator produced consistent output", badSampleCount == 0);
 	}
 
+
+    void testAllpass() {
+        int		badSampleCount = 0;
+        Jamoma::Interpolator::Allpass<Jamoma::Sample> 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.498039276942416
+        };
+
+        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;
+        my_interp.reset();
+
+        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;

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,22 @@
 	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;
+	% 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);
+	retval = out;
+endfunction
+
 for i = 1:64
 	current_delta = 2.0 + i / 64;
 	output_linear(i) = interp1(x,current_delta);
@@ -102,6 +120,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 +131,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