13
0

interpolation.cc/h: Remove all failed and obsolete attempts, leave linear and cubic

git-svn-id: svn://localhost/ardour2/branches/3.0@5424 d708f5d6-7413-0410-9779-e7cbd77b26cf
This commit is contained in:
Hans Baier 2009-07-24 05:27:49 +00:00
parent 3e88c8aa25
commit 16b964020f
4 changed files with 78 additions and 486 deletions

View File

@ -38,51 +38,6 @@ class Interpolation {
}
};
// 40.24 fixpoint math
#define FIXPOINT_ONE 0x1000000
class FixedPointLinearInterpolation : public Interpolation {
protected:
/// speed in fixed point math
uint64_t phi;
/// target speed in fixed point math
uint64_t target_phi;
std::vector<uint64_t> last_phase;
// Fixed point is just an integer with an implied scaling factor.
// In 40.24 the scaling factor is 2^24 = 16777216,
// so a value of 10*2^24 (in integer space) is equivalent to 10.0.
//
// The advantage is that addition and modulus [like x = (x + y) % 2^40]
// have no rounding errors and no drift, and just require a single integer add.
// (swh)
static const int64_t fractional_part_mask = 0xFFFFFF;
static const Sample binary_scaling_factor = 16777216.0f;
public:
FixedPointLinearInterpolation () : phi (FIXPOINT_ONE), target_phi (FIXPOINT_ONE) {}
void set_speed (double new_speed) {
target_phi = (uint64_t) (FIXPOINT_ONE * fabs(new_speed));
phi = target_phi;
}
uint64_t get_phi() { return phi; }
uint64_t get_target_phi() { return target_phi; }
uint64_t get_last_phase() { assert(last_phase.size()); return last_phase[0]; }
void set_last_phase(uint64_t phase) { assert(last_phase.size()); last_phase[0] = phase; }
void add_channel_to (int input_buffer_size, int output_buffer_size);
void remove_channel_from ();
nframes_t interpolate (int channel, nframes_t nframes, Sample* input, Sample* output);
void reset ();
};
class LinearInterpolation : public Interpolation {
protected:
@ -92,7 +47,7 @@ class LinearInterpolation : public Interpolation {
class CubicInterpolation : public Interpolation {
protected:
// shamelessly ripped from Steve Harris' swh-plugins
// shamelessly ripped from Steve Harris' swh-plugins (ladspa-util.h)
static inline float cube_interp(const float fr, const float inm1, const float
in, const float inp1, const float inp2)
{
@ -105,63 +60,6 @@ class CubicInterpolation : public Interpolation {
nframes_t interpolate (int channel, nframes_t nframes, Sample* input, Sample* output);
};
/**
* @class SplineInterpolation
*
* @brief interpolates using cubic spline interpolation over an input period
*
* Splines are piecewise cubic functions between each samples,
* where the cubic polynomials' values, first and second derivatives are equal
* on each sample point.
*
* The interpolation polynomial in the i-th interval then has the form
* p_i(x) = a3 (x - i)^3 + a2 (x - i)^2 + a1 (x - i) + a0
* = ((a3 * (x - i) + a2) * (x - i) + a1) * (x - i) + a0
* where
* a3 = (M[i+1] - M[i]) / 6
* a2 = M[i] / 2
* a1 = y[i+1] - y[i] - M[i+1]/6 - M[i]/3
* a0 = y[i]
*
* The M's are calculated recursively:
* M[i+2] = 6.0 * (y[i] - 2y[i+1] + y[i+2]) - 4M[i+1] - M[i]
*
*/
class SplineInterpolation : public Interpolation {
protected:
double M[2];
public:
void reset ();
SplineInterpolation();
nframes_t interpolate (int channel, nframes_t nframes, Sample* input, Sample* output);
};
class LibSamplerateInterpolation : public Interpolation {
protected:
std::vector<SRC_STATE*> state;
std::vector<SRC_DATA*> data;
int error;
void reset_state ();
public:
LibSamplerateInterpolation ();
~LibSamplerateInterpolation ();
void set_speed (double new_speed);
void set_target_speed (double new_speed) {}
double speed () const { return _speed; }
void add_channel_to (int input_buffer_size, int output_buffer_size);
void remove_channel_from ();
nframes_t interpolate (int channel, nframes_t nframes, Sample* input, Sample* output);
void reset() { reset_state (); }
};
} // namespace ARDOUR
#endif

View File

@ -5,69 +5,6 @@
using namespace ARDOUR;
nframes_t
FixedPointLinearInterpolation::interpolate (int channel, nframes_t nframes, Sample *input, Sample *output)
{
// the idea behind phase is that when the speed is not 1.0, we have to
// interpolate between samples and then we have to store where we thought we were.
// rather than being at sample N or N+1, we were at N+0.8792922
// so the "phase" element, if you want to think about this way,
// varies from 0 to 1, representing the "offset" between samples
uint64_t the_phase = last_phase[channel];
// acceleration
int64_t phi_delta;
// phi = fixed point speed
if (phi != target_phi) {
phi_delta = ((int64_t)(target_phi - phi)) / nframes;
} else {
phi_delta = 0;
}
// index in the input buffers
nframes_t i = 0;
for (nframes_t outsample = 0; outsample < nframes; ++outsample) {
i = the_phase >> 24;
Sample fractional_phase_part = (the_phase & fractional_part_mask) / binary_scaling_factor;
if (input && output) {
// Linearly interpolate into the output buffer
output[outsample] =
input[i] * (1.0f - fractional_phase_part) +
input[i+1] * fractional_phase_part;
}
the_phase += phi + phi_delta;
}
last_phase[channel] = (the_phase & fractional_part_mask);
// playback distance
return i;
}
void
FixedPointLinearInterpolation::add_channel_to (int /*input_buffer_size*/, int /*output_buffer_size*/)
{
last_phase.push_back (0);
}
void
FixedPointLinearInterpolation::remove_channel_from ()
{
last_phase.pop_back ();
}
void
FixedPointLinearInterpolation::reset()
{
for (size_t i = 0; i <= last_phase.size(); i++) {
last_phase[i] = 0;
}
}
nframes_t
LinearInterpolation::interpolate (int channel, nframes_t nframes, Sample *input, Sample *output)
@ -144,179 +81,3 @@ CubicInterpolation::interpolate (int channel, nframes_t nframes, Sample *input,
return i;
}
SplineInterpolation::SplineInterpolation()
{
reset ();
}
void SplineInterpolation::reset()
{
Interpolation::reset();
M[0] = 0.0;
M[1] = 0.0;
M[2] = 0.0;
}
nframes_t
SplineInterpolation::interpolate (int channel, nframes_t nframes, Sample *input, Sample *output)
{
// now interpolate
// index in the input buffers
nframes_t i = 0, delta_i = 0;
double acceleration;
double distance = 0.0;
if (_speed != _target_speed) {
acceleration = _target_speed - _speed;
} else {
acceleration = 0.0;
}
distance = phase[channel];
assert(distance >= 0.0 && distance < 1.0);
for (nframes_t outsample = 0; outsample < nframes; outsample++) {
i = floor(distance);
double x = double(distance) - double(i);
// if distance is something like 0.999999999999
// it will get rounded to 1 in the conversion to float above
while (x >= 1.0) {
x -= 1.0;
i++;
}
assert(x >= 0.0 && x < 1.0);
if (input && output) {
// if i changed, recalculate coefficients
if (delta_i == 1) {
// if i changed, rotate the M's
M[0] = M[1];
M[1] = M[2];
M[2] = 6.0 * (input[i] - 2.0*input[i+1] + input[i+2]) - 4.0*M[1] - M[0];
printf("\ny[%d] = %lf\n", i, input[i]);
printf("y[%d] = %lf\n", i+1, input[i+1]);
printf("y[%d] = %lf\n\n", i+2, input[i+2]);
printf("M[2] = %lf M[1] = %lf M[0] = %lf y-term: %lf M-term: %lf\n",
M[2], M[1], M[0], 6.0 * (input[i] - 2.0*input[i+1] + input[i+2]),
- 4.0*M[1] - M[0]);
}
double a3 = (M[1] - M[0]) / 6.0;
double a2 = M[0] / 2.0;
double a1 = input[i+1] - input[i] - (M[1] + 2.0*M[0]) / 6.0;
double a0 = input[i];
// interpolate into the output buffer
output[outsample] = ((a3*x + a2)*x + a1)*x + a0;
//printf( "input[%d/%d] = %lf/%lf distance: %lf output[%d] = %lf\n", i, i+1, input[i], input[i+1], distance, outsample, output[outsample]);
}
distance += _speed + acceleration;
delta_i = floor(distance) - i;
}
i = floor(distance);
phase[channel] = distance - floor(distance);
assert (phase[channel] >= 0.0 && phase[channel] < 1.0);
return i;
}
LibSamplerateInterpolation::LibSamplerateInterpolation() : state (0)
{
_speed = 1.0;
}
LibSamplerateInterpolation::~LibSamplerateInterpolation()
{
for (size_t i = 0; i < state.size(); i++) {
state[i] = src_delete (state[i]);
}
}
void
LibSamplerateInterpolation::set_speed (double new_speed)
{
_speed = new_speed;
for (size_t i = 0; i < state.size(); i++) {
src_set_ratio (state[i], 1.0/_speed);
}
}
void
LibSamplerateInterpolation::reset_state ()
{
printf("INTERPOLATION: reset_state()\n");
for (size_t i = 0; i < state.size(); i++) {
if (state[i]) {
src_reset (state[i]);
} else {
state[i] = src_new (SRC_SINC_FASTEST, 1, &error);
}
}
}
void
LibSamplerateInterpolation::add_channel_to (int input_buffer_size, int output_buffer_size)
{
SRC_DATA* newdata = new SRC_DATA;
/* Set up sample rate converter info. */
newdata->end_of_input = 0 ;
newdata->input_frames = input_buffer_size;
newdata->output_frames = output_buffer_size;
newdata->input_frames_used = 0 ;
newdata->output_frames_gen = 0 ;
newdata->src_ratio = 1.0/_speed;
data.push_back (newdata);
state.push_back (0);
reset_state ();
}
void
LibSamplerateInterpolation::remove_channel_from ()
{
SRC_DATA* d = data.back ();
delete d;
data.pop_back ();
if (state.back ()) {
src_delete (state.back ());
}
state.pop_back ();
reset_state ();
}
nframes_t
LibSamplerateInterpolation::interpolate (int channel, nframes_t nframes, Sample *input, Sample *output)
{
if (!data.size ()) {
printf ("ERROR: trying to interpolate with no channels\n");
return 0;
}
data[channel]->data_in = input;
data[channel]->data_out = output;
data[channel]->input_frames = nframes * _speed;
data[channel]->output_frames = nframes;
data[channel]->src_ratio = 1.0/_speed;
if ((error = src_process (state[channel], data[channel]))) {
printf ("\nError : %s\n\n", src_strerror (error));
exit (1);
}
//printf("INTERPOLATION: channel %d input_frames_used: %d\n", channel, data[channel]->input_frames_used);
return data[channel]->input_frames_used;
}

View File

@ -16,10 +16,7 @@ InterpolationTest::linearInterpolationTest ()
for (int i = 0; 3*i < NUM_SAMPLES - 1024;) {
linear.set_speed (double(1.0)/double(3.0));
linear.set_target_speed (double(1.0)/double(3.0));
//printf ("Interpolate: input: %d, output: %d, i: %d\n", input + i, output + i, i);
result = linear.interpolate (0, 1024, input + i, output + i*3);
//printf ("Result: %d\n", result);
//CPPUNIT_ASSERT_EQUAL ((uint32_t)((NUM_SAMPLES - 100) * interpolation.speed()), result);
i += result;
}
@ -87,146 +84,87 @@ InterpolationTest::linearInterpolationTest ()
}
/*
for (int i=0; i < NUM_SAMPLES; ++i) {
cout << "input[" << i << "] = " << input[i] << " output[" << i << "] = " << output[i] << endl;
}
*/
cout << i << " " << output[i] << endl;
}
*/
}
void
InterpolationTest::splineInterpolationTest ()
InterpolationTest::cubicInterpolationTest ()
{
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 + 2*i);
i += result;
}
for (int i=0; i < NUM_SAMPLES - one_period; ++i) {
//cout << "input[" << i << "] = " << input[i] << " 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 wave
for (int i = 0; i < NUM_SAMPLES; ++i) {
if (i % (INTERVAL/2) < INTERVAL/4 ) {
input[i] = 1.0f;
} else {
input[i] = 0.0f;
}
output[i] = 0.0f;
}
/*
//sine wave
for (int i = 0; i < NUM_SAMPLES; ++i) {
input[i] = sin(double(i) * M_2_PI / INTERVAL * 10.0);
}
*/
one_period = 512;
cout << "\nSpeed: 1/60" << endl;
spline.reset();
spline.set_speed (1.0/90.0);
for (int i = 0, o = 0; 90 * i < NUM_SAMPLES - one_period; o++) {
result = spline.interpolate (0, one_period, input + i, output + o * one_period);
//printf ("Result: %d\n", result);
cout << "\nCubic Interpolation Test\n";
cout << "\nSpeed: 1/3";
for (int i = 0; 3*i < NUM_SAMPLES - 1024;) {
cubic.set_speed (double(1.0)/double(3.0));
cubic.set_target_speed (double(1.0)/double(3.0));
result = cubic.interpolate (0, 1024, input + i, output + i*3);
i += result;
}
for (int i=0; i < NUM_SAMPLES - one_period; ++i) {
cout << 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])); }
cout << "\nSpeed: 1.0";
cubic.reset();
cubic.set_speed (1.0);
cubic.set_target_speed (cubic.speed());
result = cubic.interpolate (0, NUM_SAMPLES, input, output);
CPPUNIT_ASSERT_EQUAL ((uint32_t)(NUM_SAMPLES * cubic.speed()), result);
for (int i = 0; i < NUM_SAMPLES; i += INTERVAL) {
CPPUNIT_ASSERT_EQUAL (1.0f, output[i]);
}
cout << "\nSpeed: 0.5";
cubic.reset();
cubic.set_speed (0.5);
cubic.set_target_speed (cubic.speed());
result = cubic.interpolate (0, NUM_SAMPLES, input, output);
CPPUNIT_ASSERT_EQUAL ((uint32_t)(NUM_SAMPLES * cubic.speed()), result);
for (int i = 0; i < NUM_SAMPLES; i += (INTERVAL / cubic.speed() +0.5)) {
CPPUNIT_ASSERT_EQUAL (1.0f, output[i]);
}
cout << "\nSpeed: 0.2";
cubic.reset();
cubic.set_speed (0.2);
cubic.set_target_speed (cubic.speed());
result = cubic.interpolate (0, NUM_SAMPLES, input, output);
CPPUNIT_ASSERT_EQUAL ((uint32_t)(NUM_SAMPLES * cubic.speed()), result);
cout << "\nSpeed: 0.02";
cubic.reset();
cubic.set_speed (0.02);
cubic.set_target_speed (cubic.speed());
result = cubic.interpolate (0, NUM_SAMPLES, input, output);
CPPUNIT_ASSERT_EQUAL ((uint32_t)(NUM_SAMPLES * cubic.speed()), result);
/* This one fails due too error accumulation
cout << "\nSpeed: 0.002";
cubic.reset();
cubic.set_speed (0.002);
cubic.set_target_speed (cubic.speed());
result = cubic.interpolate (0, NUM_SAMPLES, input, output);
cubic.speed();
CPPUNIT_ASSERT_EQUAL ((nframes_t)(NUM_SAMPLES * cubic.speed()), result);
*/
cout << "\nSpeed: 2.0";
cubic.reset();
cubic.set_speed (2.0);
cubic.set_target_speed (cubic.speed());
result = cubic.interpolate (0, NUM_SAMPLES / 2, input, output);
CPPUNIT_ASSERT_EQUAL ((uint32_t)(NUM_SAMPLES / 2 * cubic.speed()), result);
for (int i = 0; i < NUM_SAMPLES / 2; i += (INTERVAL / cubic.speed() +0.5)) {
CPPUNIT_ASSERT_EQUAL (1.0f, output[i]);
}
cout << "\nSpeed: 10.0";
cubic.set_speed (10.0);
cubic.set_target_speed (cubic.speed());
result = cubic.interpolate (0, NUM_SAMPLES / 10, input, output);
CPPUNIT_ASSERT_EQUAL ((uint32_t)(NUM_SAMPLES / 10 * cubic.speed()), result);
for (int i = 0; i < NUM_SAMPLES / 10; i += (INTERVAL / cubic.speed() +0.5)) {
CPPUNIT_ASSERT_EQUAL (1.0f, output[i]);
}
}
void
InterpolationTest::libSamplerateInterpolationTest ()
{
nframes_t result;
cout << "\nLibSamplerate Interpolation Test\n";
/*
cout << "\nSpeed: 1.0";
interpolation.set_speed (1.0);
for (int i = 0; i < NUM_SAMPLES;) {
interpolation.set_speed (1.0);
result = interpolation.interpolate (0, INTERVAL/10, input + i, output + i);
CPPUNIT_ASSERT_EQUAL ((uint32_t)(INTERVAL/10 * interpolation.speed()), result);
i += result;
}
for (int i = 0; i < NUM_SAMPLES; i += INTERVAL) {
CPPUNIT_ASSERT_EQUAL (1.0f, output[i+1]);
}
*/
cout << "\nSpeed: 0.5";
for (int i = 0; i < NUM_SAMPLES;) {
interpolation.set_speed (0.5);
//printf ("Interpolate: input: %d, output: %d, i: %d\n", input + i, output + i, i);
result = interpolation.interpolate (0, NUM_SAMPLES - 100, input + i, output + i);
printf ("Result: %d\n", result);
//CPPUNIT_ASSERT_EQUAL ((uint32_t)((NUM_SAMPLES - 100) * interpolation.speed()), result);
//i += result;
break;
}
for (int i=0; i < NUM_SAMPLES; ++i) {
cout << "input[" << i << "] = " << input[i] << " output[" << i << "] = " << output[i] << endl;
}
cout << "\nSpeed: 0.2";
interpolation.set_speed (0.2);
result = interpolation.interpolate (0, NUM_SAMPLES, input, output);
CPPUNIT_ASSERT_EQUAL ((uint32_t)(NUM_SAMPLES * interpolation.speed()), result);
cout << "\nSpeed: 0.02";
interpolation.set_speed (0.02);
result = interpolation.interpolate (0, NUM_SAMPLES, input, output);
CPPUNIT_ASSERT_EQUAL ((uint32_t)(NUM_SAMPLES * interpolation.speed()), result);
cout << "\nSpeed: 0.002";
interpolation.set_speed (0.002);
result = interpolation.interpolate (0, NUM_SAMPLES, input, output);
CPPUNIT_ASSERT_EQUAL ((uint32_t)(NUM_SAMPLES * interpolation.speed()), result);
cout << "\nSpeed: 2.0";
interpolation.set_speed (2.0);
result = interpolation.interpolate (0, NUM_SAMPLES / 2, input, output);
CPPUNIT_ASSERT_EQUAL ((uint32_t)(NUM_SAMPLES / 2 * interpolation.speed()), result);
for (int i = 0; i < NUM_SAMPLES / 2; i += (INTERVAL / interpolation.speed() +0.5)) {
CPPUNIT_ASSERT_EQUAL (1.0f, output[i]);
}
cout << "\nSpeed: 10.0";
interpolation.set_speed (10.0);
result = interpolation.interpolate (0, NUM_SAMPLES / 10, input, output);
CPPUNIT_ASSERT_EQUAL ((uint32_t)(NUM_SAMPLES / 10 * interpolation.speed()), result);
for (int i = 0; i < NUM_SAMPLES / 10; i += (INTERVAL / interpolation.speed() +0.5)) {
CPPUNIT_ASSERT_EQUAL (1.0f, output[i]);
}
/*
for (int i=0; i < NUM_SAMPLES; ++i) {
cout << "input[" << i << "] = " << input[i] << " output[" << i << "] = " << output[i] << endl;
}
*/
}

View File

@ -26,9 +26,8 @@
class InterpolationTest : public CppUnit::TestFixture
{
CPPUNIT_TEST_SUITE(InterpolationTest);
CPPUNIT_TEST(splineInterpolationTest);
//CPPUNIT_TEST(linearInterpolationTest);
//CPPUNIT_TEST(libSamplerateInterpolationTest);
CPPUNIT_TEST(cubicInterpolationTest);
CPPUNIT_TEST(linearInterpolationTest);
CPPUNIT_TEST_SUITE_END();
#define NUM_SAMPLES 1000000
@ -38,8 +37,7 @@ class InterpolationTest : public CppUnit::TestFixture
ARDOUR::Sample output[NUM_SAMPLES];
ARDOUR::LinearInterpolation linear;
ARDOUR::SplineInterpolation spline;
ARDOUR::LibSamplerateInterpolation interpolation;
ARDOUR::CubicInterpolation cubic;
public:
@ -53,15 +51,12 @@ 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);
cubic.add_channel_to (NUM_SAMPLES, NUM_SAMPLES);
}
void tearDown() {
}
void linearInterpolationTest();
void splineInterpolationTest();
void libSamplerateInterpolationTest();
void cubicInterpolationTest();
};