/*
- Copyright (C) 2012 Paul Davis
+ Copyright (C) 2012 Paul Davis
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
*/
-#include <stdint.h>
+#include <limits>
#include <cstdio>
+#include <stdint.h>
+
#include "ardour/interpolation.h"
#include "ardour/midi_buffer.h"
using namespace ARDOUR;
+using std::cerr;
+using std::endl;
-
-framecnt_t
-LinearInterpolation::interpolate (int channel, framecnt_t nframes, Sample *input, Sample *output)
+CubicInterpolation::CubicInterpolation ()
+ : valid_z_bits (0)
{
- // index in the input buffers
- framecnt_t i = 0;
-
- double acceleration = 0;
-
- if (_speed != _target_speed) {
- acceleration = _target_speed - _speed;
- }
-
- for (framecnt_t outsample = 0; outsample < nframes; ++outsample) {
- double const d = phase[channel] + outsample * (_speed + acceleration);
- i = floor(d);
- Sample fractional_phase_part = d - i;
- if (fractional_phase_part >= 1.0) {
- fractional_phase_part -= 1.0;
- i++;
- }
-
- if (input && output) {
- // Linearly interpolate into the output buffer
- output[outsample] =
- input[i] * (1.0f - fractional_phase_part) +
- input[i+1] * fractional_phase_part;
- }
- }
-
- double const distance = phase[channel] + nframes * (_speed + acceleration);
- i = floor(distance);
- phase[channel] = distance - i;
- return i;
}
-framecnt_t
-CubicInterpolation::interpolate (int channel, framecnt_t nframes, Sample *input, Sample *output)
+samplecnt_t
+CubicInterpolation::interpolate (int channel, samplecnt_t input_samples, Sample *input, samplecnt_t & output_samples, Sample *output)
{
- // index in the input buffers
- framecnt_t i = 0;
-
- double acceleration;
- double distance = 0.0;
+ assert (input_samples > 0);
+ assert (output_samples > 0);
+ assert (input);
+ assert (output);
+ assert (phase.size () > channel);
+
+ _speed = fabs (_speed);
+
+ if (invalid (0)) {
+
+ /* z[0] not set. Two possibilities
+ *
+ * 1) we have just been constructed or ::reset()
+ *
+ * 2) we were only given 1 sample after construction or
+ * ::reset, and stored it in z[1]
+ */
- if (_speed != _target_speed) {
- acceleration = _target_speed - _speed;
- } else {
- acceleration = 0.0;
- }
+ if (invalid (1)) {
+
+ /* first call after construction or after ::reset */
+
+ switch (input_samples) {
+ case 1:
+ /* store one sample for use next time. We don't
+ * have enough points to interpolate or even
+ * compute the first z[0] value, but keep z[1]
+ * around.
+ */
+ z[1] = input[0]; validate (1);
+ output_samples = 0;
+ return 0;
+ case 2:
+ /* store two samples for use next time, and
+ * compute a value for z[0] that will maintain
+ * the slope of the first actual segment. We
+ * still don't have enough samples to interpolate.
+ */
+ z[0] = input[0] - (input[1] - input[0]); validate (0);
+ z[1] = input[0]; validate (1);
+ z[2] = input[1]; validate (2);
+ output_samples = 0;
+ return 0;
+ default:
+ /* We have enough samples to interpolate this time,
+ * but don't have a valid z[0] value because this is the
+ * first call after construction or ::reset.
+ *
+ * First point is based on a requirement to maintain
+ * the slope of the first actual segment
+ */
+ z[0] = input[0] - (input[1] - input[0]); validate (0);
+ break;
+ }
+ } else {
- distance = phase[channel];
+ /* at least one call since construction or
+ * after::reset, since we have z[1] set
+ *
+ * we can now compute z[0] as required
+ */
- if (nframes < 3) {
- /* no interpolation possible */
+ z[0] = z[1] - (input[0] - z[1]); validate (0);
- if (input && output) {
- for (i = 0; i < nframes; ++i) {
- output[i] = input[i];
- }
+ /* we'll check the number of samples we've been given
+ in the next switch() statement below, and either
+ just save some more samples or actual interpolate
+ */
}
- return nframes;
+ assert (is_valid (0));
}
- /* keep this condition out of the inner loop */
-
- if (input && output) {
+ switch (input_samples) {
+ case 1:
+ /* one more sample of input. find the right vX to store
+ it in, and decide if we're ready to interpolate
+ */
+ if (invalid (1)) {
+ z[1] = input[0]; validate (1);
+ /* still not ready to interpolate */
+ output_samples = 0;
+ return 0;
+ } else if (invalid (2)) {
+ /* still not ready to interpolate */
+ z[2] = input[0]; validate (2);
+ output_samples = 0;
+ return 0;
+ } else if (invalid (3)) {
+ z[3] = input[0]; validate (3);
+ /* ready to interpolate */
+ }
+ break;
+ case 2:
+ /* two more samples of input. find the right vX to store
+ them in, and decide if we're ready to interpolate
+ */
+ if (invalid (1)) {
+ z[1] = input[0]; validate (1);
+ z[2] = input[1]; validate (2);
+ /* still not ready to interpolate */
+ output_samples = 0;
+ return 0;
+ } else if (invalid (2)) {
+ z[2] = input[0]; validate (2);
+ z[3] = input[1]; validate (3);
+ /* ready to interpolate */
+ } else if (invalid (3)) {
+ z[3] = input[0]; validate (3);
+ /* ready to interpolate */
+ }
+ break;
+
+ default:
+ /* caller has given us at least enough samples to interpolate a
+ single value.
+ */
+ z[1] = input[0]; validate (1);
+ z[2] = input[1]; validate (2);
+ z[3] = input[2]; validate (3);
+ }
- Sample inm1;
+ /* ready to interpolate using z[0], z[1], z[2] and z[3] */
- if (floor (distance) == 0.0) {
- /* best guess for the fake point we have to add to be able to interpolate at i == 0:
- .... maintain slope of first actual segment ...
- */
- inm1 = input[i] - (input[i+1] - input[i]);
- } else {
- inm1 = input[i-1];
- }
+ assert (is_valid (0));
+ assert (is_valid (1));
+ assert (is_valid (2));
+ assert (is_valid (3));
- for (framecnt_t outsample = 0; outsample < nframes; ++outsample) {
+ /* we can use up to (input_samples - 2) of the input, so compute the
+ * maximum number of output samples that represents.
+ *
+ * Remember that the expected common case here is to be given
+ * input_samples that is substantially larger than output_samples,
+ * thus allowing us to always compute output_samples in one call.
+ */
- float f = floor (distance);
- float fractional_phase_part = distance - f;
+ const samplecnt_t output_from_input = floor ((input_samples - 2) / _speed);
- /* get the index into the input we should start with */
+ /* limit output to either the caller's requested number or the number
+ * determined by the input size.
+ */
- i = lrintf (f);
+ const samplecnt_t limit = std::min (output_samples, output_from_input);
- /* fractional_phase_part only reaches 1.0 thanks to float imprecision. In theory
- it should always be < 1.0. If it ever >= 1.0, then bump the index we use
- and back it off. This is the point where we "skip" an entire sample in the
- input, because the phase part has accumulated so much error that we should
- really be closer to the next sample. or something like that ...
- */
+ samplecnt_t outsample = 0;
+ double distance = phase[channel];
+ samplecnt_t used = floor (distance);
+ samplecnt_t i = 0;
- if (fractional_phase_part >= 1.0) {
- fractional_phase_part -= 1.0;
- ++i;
- }
+ while (outsample < limit) {
- // Cubically interpolate into the output buffer: keep this inlined for speed and rely on compiler
- // optimization to take care of the rest
- // shamelessly ripped from Steve Harris' swh-plugins (ladspa-util.h)
+ i = floor (distance);
- output[outsample] = input[i] + 0.5f * fractional_phase_part * (input[i+1] - inm1 +
- fractional_phase_part * (4.0f * input[i+1] + 2.0f * inm1 - 5.0f * input[i] - input[i+2] +
- fractional_phase_part * (3.0f * (input[i] - input[i+1]) - inm1 + input[i+2])));
+ /* this call may stop the loop from being vectorized */
+ float fractional_phase_part = fmod (distance, 1.0);
- distance += _speed + acceleration;
- inm1 = input[i];
- }
+ /* Cubically interpolate into the output buffer */
+ output[outsample++] = z[1] + 0.5f * fractional_phase_part *
+ (z[2] - z[0] + fractional_phase_part * (4.0f * z[2] + 2.0f * z[0] - 5.0f * z[1] - z[3] +
+ fractional_phase_part * (3.0f * (z[1] - z[2]) - z[0] + z[3])));
- i = floor(distance);
- phase[channel] = distance - floor(distance);
+ distance += _speed;
- } else {
- /* used to calculate play-distance with acceleration (silent roll)
- * (use same algorithm as real playback for identical rounding/floor'ing)
- */
- for (framecnt_t outsample = 0; outsample < nframes; ++outsample) {
- distance += _speed + acceleration;
- }
- i = floor(distance);
+ z[0] = z[1];
+ z[1] = input[i];
+ z[2] = input[i+1];
+ z[3] = input[i+2];
}
- return i;
+ output_samples = outsample;
+ phase[channel] = fmod (distance, 1.0);
+ return i - used;
}
-framecnt_t
-CubicMidiInterpolation::distance (framecnt_t nframes, bool roll)
+void
+CubicInterpolation::reset ()
{
- assert(phase.size() == 1);
-
- framecnt_t i = 0;
-
- double acceleration;
- double distance = 0.0;
-
- if (nframes < 3) {
- return nframes;
- }
-
- if (_speed != _target_speed) {
- acceleration = _target_speed - _speed;
- } else {
- acceleration = 0.0;
- }
-
- distance = phase[0];
-
- for (framecnt_t outsample = 0; outsample < nframes; ++outsample) {
- distance += _speed + acceleration;
- }
-
- if (roll) {
- phase[0] = distance - floor(distance);
- }
-
- i = floor(distance);
+ Interpolation::reset ();
+ valid_z_bits = 0;
+}
- return i;
+samplecnt_t
+CubicInterpolation::distance (samplecnt_t nsamples)
+{
+ assert (phase.size () > 0);
+ return floor (floor (phase[0]) + (_speed * nsamples));
}
projects / ardour.git / blobdiff