}
};
-// 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:
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)
{
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
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)
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;
-}
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;
}
}
/*
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);
+ 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 << "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])); }
- }
- */
+ 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]);
+ }
- // 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;
+ 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);
- /*
- //sine wave
- for (int i = 0; i < NUM_SAMPLES; ++i) {
- input[i] = sin(double(i) * M_2_PI / INTERVAL * 10.0);
- }
+ /* 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);
*/
- 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);
- 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: 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]);
}
-}
-
-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;
- }
- */
+ 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]);
+ }
}
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
ARDOUR::Sample output[NUM_SAMPLES];
ARDOUR::LinearInterpolation linear;
- ARDOUR::SplineInterpolation spline;
- ARDOUR::LibSamplerateInterpolation interpolation;
+ ARDOUR::CubicInterpolation cubic;
public:
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();
};
projects / ardour.git / commitdiff