+/*
+ 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
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program; if not, write to the Free Software
+ Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+
+*/
+
+#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;
-nframes_t
-LinearInterpolation::interpolate (nframes_t nframes, Sample *input, Sample *output)
+CubicInterpolation::CubicInterpolation ()
+ : valid_z_bits (0)
{
- // the idea 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
-
- // index in the input buffers
- nframes_t i = 0;
-
- double acceleration;
- double distance = 0.0;
-
- if (_speed != _target_speed) {
- acceleration = _target_speed - _speed;
- } else {
- acceleration = 0.0;
+}
+
+samplecnt_t
+CubicInterpolation::interpolate (int channel, samplecnt_t input_samples, Sample *input, samplecnt_t & output_samples, Sample *output)
+{
+ 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 (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 {
+
+ /* at least one call since construction or
+ * after::reset, since we have z[1] set
+ *
+ * we can now compute z[0] as required
+ */
+
+ z[0] = z[1] - (input[0] - z[1]); validate (0);
+
+ /* 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
+ */
+ }
+
+ assert (is_valid (0));
}
- for (nframes_t outsample = 0; outsample < nframes; ++outsample) {
- i = distance;
- Sample fractional_phase_part = distance - 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;
+ 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 */
}
- distance += _speed + acceleration;
+ 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);
+ }
+
+ /* ready to interpolate using z[0], z[1], z[2] and z[3] */
+
+ assert (is_valid (0));
+ assert (is_valid (1));
+ assert (is_valid (2));
+ assert (is_valid (3));
+
+ /* 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.
+ */
+
+ const samplecnt_t output_from_input = floor ((input_samples - 2) / _speed);
+
+ /* limit output to either the caller's requested number or the number
+ * determined by the input size.
+ */
+
+ const samplecnt_t limit = std::min (output_samples, output_from_input);
+
+ samplecnt_t outsample = 0;
+ double distance = phase[channel];
+ samplecnt_t used = floor (distance);
+ samplecnt_t i = 0;
+
+ while (outsample < limit) {
+
+ i = floor (distance);
+
+ /* this call may stop the loop from being vectorized */
+ float fractional_phase_part = fmod (distance, 1.0);
+
+ /* 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])));
+
+ distance += _speed;
+
+ z[0] = z[1];
+ z[1] = input[i];
+ z[2] = input[i+1];
+ z[3] = input[i+2];
}
-
- i = (distance + 0.5L);
- // playback distance
- return i;
+
+ output_samples = outsample;
+ phase[channel] = fmod (distance, 1.0);
+ return i - used;
+}
+
+void
+CubicInterpolation::reset ()
+{
+ Interpolation::reset ();
+ valid_z_bits = 0;
+}
+
+samplecnt_t
+CubicInterpolation::distance (samplecnt_t nsamples)
+{
+ assert (phase.size () > 0);
+ return floor (floor (phase[0]) + (_speed * nsamples));
}