mirror of
https://github.com/Ardour/ardour.git
synced 2026-01-06 13:45:43 +01:00
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
This commit is contained in:
Hans Baier
parent
272cad6241
commit
b0bb5342dd
4 changed files with 119 additions and 28 deletions
|
|
@ -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();
|
||||
|
|
|
|||
|
|
@ -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;
|
||||
}
|
||||
|
|
|
|||
|
|
@ -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 ()
|
||||
{
|
||||
|
|
|
|||
|
|
@ -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();
|
||||
|
||||
};
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue