From b0bb5342dd0fddf867208ea1e4c4a87a32a7b214 Mon Sep 17 00:00:00 2001 From: Hans Baier Date: Wed, 22 Jul 2009 08:42:33 +0000 Subject: 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 --- libs/ardour/ardour/interpolation.h | 36 +++++++++++++++++++--- libs/ardour/interpolation.cc | 49 +++++++++++++++++------------- libs/ardour/tests/interpolation-test.cc | 54 +++++++++++++++++++++++++++++++++ libs/ardour/tests/interpolation-test.h | 8 +++-- 4 files changed, 119 insertions(+), 28 deletions(-) (limited to 'libs') 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 @@ -91,6 +91,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 () { 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(); }; -- cgit v1.2.3