interpolation.cc/.h: Spline-Bugfixes: Crash bug at tempos close to 0, wrong calculation of M, unbounded precalculated L/U Matrices

git-svn-id: svn://localhost/ardour2/branches/3.0@5410 d708f5d6-7413-0410-9779-e7cbd77b26cf

4 files changed, 119 insertions, 28 deletions

diff --git a/libs/ardour/ardour/interpolation.h b/libs/ardour/ardour/interpolation.h
index 6ceb63e527..f417c8c46e 100644
--- a/libs/ardour/ardour/interpolation.h
+++ b/libs/ardour/ardour/interpolation.h
@@ -19,7 +19,8 @@ class Interpolation {
 
              
  public:
-     Interpolation () { _speed = 1.0; _target_speed = 1.0; }
+     Interpolation ()  { _speed = 1.0; _target_speed = 1.0; }
+     ~Interpolation () { phase.clear(); }
  
      void set_speed (double new_speed)          { _speed = new_speed; _target_speed = new_speed; }
      void set_target_speed (double new_speed)   { _target_speed = new_speed; }
@@ -90,7 +91,6 @@ class LinearInterpolation : public Interpolation {
 };
  
 
-#define MAX_PERIOD_SIZE 4096
 /**
  * @class SplineInterpolation
  * 
@@ -143,14 +143,40 @@ class LinearInterpolation : public Interpolation {
  *
  *  where the l[i] and m[i] can be precomputed.
  * 
- *  Then we solve the system A * M = d by first solving the system
+ *  Then we solve the system A * M = L(UM) = d by first solving the system
  *    L * t = d 
+ *    
+ *    [ 1    0    0    0   ... 0      0      0      0 ]    [ t[0]   ]    [ 6*y[0] - 12*y[1] + 6*y[2] ]
+ *    [ l[0] 1    0    0   ... 0      0      0      0 ]    [ t[1]   ]    [ 6*y[1] - 12*y[2] + 6*y[3] ]
+ *    [ 0    l[1] 1    0   ... 0      0      0      0 ]    [ t[2]   ]    [ 6*y[2] - 12*y[3] + 6*y[4] ]
+ *    [ 0    0    l[2] 1   ... 0      0      0      0 ]    [ t[3]   ]    [ 6*y[3] - 12*y[4] + 6*y[5] ]
+ *                        ...                           *             =             ...            
+ *    [ 0    0    0    0   ... 1      0      0      0 ]    [ t[n-6] ]    [ 6*y[n-6]- 12*y[n-5] + 6*y[n-4] ]
+ *    [ 0    0    0    0   ... l[n-6] 1      0      0 ]    [ t[n-5] ]    [ 6*y[n-5]- 12*y[n-4] + 6*y[n-3] ]
+ *    [ 0    0    0    0   ... 0      l[n-5] 1      0 ]    [ t[n-4] ]    [ 6*y[n-4]- 12*y[n-3] + 6*y[n-2] ]
+ *    [ 0    0    0    0   ... 0      0      l[n-4] 1 ]    [ t[n-3] ]    [ 6*y[n-3]- 12*y[n-2] + 6*y[n-1] ]
+ *    
+ *    
  *  and then
- *    R * M = t
+ *    U * M = t
+ *  
+ *  [ m[0] 1    0    0   ... 0      0      0      ]   [ M[1]   ]    [ t[0]   ]
+ *  [ 0    m[1] 1    0   ... 0      0      0      ]   [ M[2]   ]    [ t[1]   ]
+ *  [ 0    0    m[2] 1   ... 0      0      0      ]   [ M[3]   ]    [ t[2]   ]
+ *                       ...                          [ M[4]   ]    [ t[3]   ]
+ *  [ 0    0    0    0   ... 0      0      0      ] *            =            
+ *  [ 0    0    0    0   ... 1      0      0      ]   [ M[n-5] ]    [ t[n-6] ]
+ *  [ 0    0    0    0   ... m[n-5] 1      0      ]   [ M[n-4] ]    [ t[n-5] ]
+ *  [ 0    0    0    0   ... 0      m[n-4] 1      ]   [ M[n-3] ]    [ t[n-4] ]
+ *  [ 0    0    0    0   ... 0      0      m[n-3] ]   [ M[n-2] ]    [ t[n-3] ]
+ *  
  */
 class SplineInterpolation : public Interpolation {
  protected:
-    double l[MAX_PERIOD_SIZE], m[MAX_PERIOD_SIZE];
+    double _l[19], _m[20];
+    
+    inline double l(nframes_t i) {  return (i >= 19) ? _l[18] : _l[i]; }
+    inline double m(nframes_t i) {  return (i >= 20) ? _m[19] : _m[i]; }
     
  public:
     SplineInterpolation();
diff --git a/libs/ardour/interpolation.cc b/libs/ardour/interpolation.cc
index 8b4bb862ed..716ebd4139 100644
--- a/libs/ardour/interpolation.cc
+++ b/libs/ardour/interpolation.cc
@@ -120,10 +120,12 @@ SplineInterpolation::SplineInterpolation()
 {
     // precompute LU-factorization of matrix A
     // see "Teubner Taschenbuch der Mathematik", p. 1105
-    m[0] = 4.0;
-    for (int i = 0; i <= MAX_PERIOD_SIZE - 2; i++) {
-        l[i] = 1.0 / m[i];
-        m[i+1] = 4.0 - l[i];
+    // We only need to calculate up to 20, because they
+    // won't change any more above that
+    _m[0] = 4.0;
+    for (int i = 0; i <= 20 - 2; i++) {
+        _l[i] = 1.0 / _m[i];
+        _m[i+1] = 4.0 - _l[i];
     }
 }
 
@@ -131,9 +133,12 @@ nframes_t
 SplineInterpolation::interpolate (int channel, nframes_t nframes, Sample *input, Sample *output)
 {
     // How many input samples we need
-    nframes_t n = ceil (double(nframes) * _speed) + 2;
-    //   |------------------------------------------^
-    // this won't be here in the debugged version.
+    nframes_t n = ceil (double(nframes) * _speed + phase[channel]) + 1;
+    //printf("n = %d\n", n);
+
+    if (n <= 3) {
+        return 0;
+    }
     
     double M[n], t[n-2];
     
@@ -142,20 +147,19 @@ SplineInterpolation::interpolate (int channel, nframes_t nframes, Sample *input,
     M[n - 1] = 0.0;
     
     // solve L * t = d
-    // see "Teubner Taschenbuch der Mathematik", p. 1105
     t[0] = 6.0 * (input[0] - 2*input[1] + input[2]); 
     for (nframes_t i = 1; i <= n - 3; i++) {
         t[i] = 6.0 * (input[i] - 2*input[i+1] + input[i+2])
-               - l[i-1] * t[i-1];
+               - l(i-1) * t[i-1];
     }
     
-    // solve R * M = t
-    // see "Teubner Taschenbuch der Mathematik", p. 1105
-    M[n-2] = -t[n-3] / m[n-3];
+    // solve U * M = t
+    M[n-2] = t[n-3] / m(n-3);
     for (nframes_t i = n-4;; i--) {
-        M[i+1] = -(t[i] + M[i+2]) / m[i];
+        M[i+1] = (t[i]-M[i+2])/m(i);
         if ( i == 0 ) break;
     }
+    assert (M[0] == 0.0 && M[n-1] == 0.0);
     
     // now interpolate
     // index in the input buffers
@@ -174,29 +178,32 @@ SplineInterpolation::interpolate (int channel, nframes_t nframes, Sample *input,
     for (nframes_t outsample = 0; outsample < nframes; outsample++) {
         i = floor(distance);
         
-        Sample x = distance - i;
+        Sample x = double(distance) - double(i);
         
-        /* this would break the assertion below
+        // if distance is something like 0.999999999999
+        // it will get rounded to 1 in the conversion to float above
         if (x >= 1.0) {
-            x -= 1.0;
+            x = 0.0;
             i++;
         }
-        */
+        
+        assert(x >= 0.0 && x < 1.0);
         
         if (input && output) {
             assert (i <= n-1);
-            double a0 = input[i];
-            double a1 = input[i+1] - input[i] - M[i+1]/6.0 - M[i]/3.0;
-            double a2 = M[i] / 2.0;
             double a3 = (M[i+1] - M[i]) / 6.0;
+            double a2 = M[i] / 2.0;
+            double a1 = input[i+1] - input[i] - (M[i+1] + 2.0*M[i])/6.0;
+            double a0 = input[i];
             // interpolate into the output buffer
-            output[outsample] = ((a3*x +a2)*x +a1)*x + a0;
+            output[outsample] = ((a3*x + a2)*x + a1)*x + a0;
         }
         distance += _speed + acceleration;
     }
     
     i = floor(distance);
     phase[channel] = distance - floor(distance);
+    assert (phase[channel] >= 0.0 && phase[channel] < 1.0);
     
     return i;
 }
diff --git a/libs/ardour/tests/interpolation-test.cc b/libs/ardour/tests/interpolation-test.cc
index ec14abd8d3..6df46ad194 100644
--- a/libs/ardour/tests/interpolation-test.cc
+++ b/libs/ardour/tests/interpolation-test.cc
@@ -92,6 +92,60 @@ InterpolationTest::linearInterpolationTest ()
 }
 
 void
+InterpolationTest::splineInterpolationTest ()
+{
+        nframes_t result = 0;
+         cout << "\nspline Interpolation Test\n";
+         
+         cout << "\nSpeed: 1/2" << endl;
+         spline.reset();
+         spline.set_speed (0.5);
+         int one_period = 1024;
+         /*
+         
+         for (int i = 0; 2 * i < NUM_SAMPLES - one_period;) {
+             result = spline.interpolate (0, one_period, input + i, output + int(2*i));
+             i += result;
+         }
+         for (int i=0; i < NUM_SAMPLES - one_period; ++i) {
+             //cout << "output[" << i << "] = " << output[i] << endl;    
+             if (i % 200 == 0) { CPPUNIT_ASSERT_EQUAL (double(1.0), double(output[i])); }
+             else if (i % 2 == 0) { CPPUNIT_ASSERT_EQUAL (double(0.0), double(output[i])); }
+         }
+         */
+         
+         /*
+         // square function
+         
+         for (int i = 0; i < NUM_SAMPLES; ++i) {
+             if (i % INTERVAL/8 < INTERVAL/16 ) {
+                 input[i] = 1.0f;
+             } else {
+                 input[i] = 0.0f;
+             }
+             output[i] = 0.0f;
+         }
+         */
+         
+         cout << "\nSpeed: 1/60" << endl;
+         spline.reset();
+         spline.set_speed (1.0/60.0);
+         
+         one_period = 8192;
+         
+         for (int i = 0; 60 * i < NUM_SAMPLES - one_period;) {
+             result = spline.interpolate (0, one_period, input + i, output + int(60*i));
+             printf ("Result: %d\n", result);
+             i += result;
+         }
+         for (int i=0; i < NUM_SAMPLES - one_period; ++i) {
+             cout << "input[" << i << "] = " << input[i] << "  output[" << i << "] = " << output[i] << endl; 
+             //if (i % 333 == 0) { CPPUNIT_ASSERT_EQUAL (double(1.0), double(output[i])); }
+             //else if (i % 2 == 0) { CPPUNIT_ASSERT_EQUAL (double(0.0), double(output[i])); }
+         }
+}
+
+void
 InterpolationTest::libSamplerateInterpolationTest ()
 {
     nframes_t result;
diff --git a/libs/ardour/tests/interpolation-test.h b/libs/ardour/tests/interpolation-test.h
index 34ec0daaae..c5cc3b67b1 100644
--- a/libs/ardour/tests/interpolation-test.h
+++ b/libs/ardour/tests/interpolation-test.h
@@ -26,7 +26,8 @@
 class InterpolationTest : public CppUnit::TestFixture
 {
     CPPUNIT_TEST_SUITE(InterpolationTest);
-    CPPUNIT_TEST(linearInterpolationTest);
+    //CPPUNIT_TEST(linearInterpolationTest);
+    CPPUNIT_TEST(splineInterpolationTest);
     //CPPUNIT_TEST(libSamplerateInterpolationTest);
     CPPUNIT_TEST_SUITE_END();
     
@@ -37,13 +38,14 @@ class InterpolationTest : public CppUnit::TestFixture
     ARDOUR::Sample output[NUM_SAMPLES];
     
     ARDOUR::LinearInterpolation linear;
+    ARDOUR::SplineInterpolation spline;
     ARDOUR::LibSamplerateInterpolation interpolation;
 
     public:
        	
         void setUp() {
             for (int i = 0; i < NUM_SAMPLES; ++i) {
-                if (i % INTERVAL == 0) {
+                if (i % INTERVAL == 50) {
                     input[i] = 1.0f;
                 } else {
                     input[i] = 0.0f;
@@ -51,6 +53,7 @@ class InterpolationTest : public CppUnit::TestFixture
                 output[i] = 0.0f;
             }
             linear.add_channel_to (NUM_SAMPLES, NUM_SAMPLES);
+            spline.add_channel_to (NUM_SAMPLES, NUM_SAMPLES);
             interpolation.add_channel_to (NUM_SAMPLES, NUM_SAMPLES);
         }
         
@@ -58,6 +61,7 @@ class InterpolationTest : public CppUnit::TestFixture
         }
 
         void linearInterpolationTest();
+        void splineInterpolationTest();
         void libSamplerateInterpolationTest();
 
 };