One fix.
[ardour.git] / libs / ardour / interpolation.cc
index 7aece6453cf2623261d70b9f945e7bf99cd64190..286030c26db487cb8d825d0bd712119d3be9116a 100644 (file)
+/*
+    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 <stdint.h>
+#include <cstdio>
+
 #include "ardour/interpolation.h"
+#include "ardour/midi_buffer.h"
+
+using namespace ARDOUR;
 
-nframes_t
-LinearInterpolation::interpolate (nframes_t nframes, Sample *input, Sample *output)
+
+framecnt_t
+LinearInterpolation::interpolate (int channel, framecnt_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    phase = last_phase;
-       
-       // 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
+       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)
+{
        // index in the input buffers
-       nframes_t   i = 0;
-
-       for (nframes_t outsample = 0; outsample < nframes; ++outsample) {
-               i = phase >> 24;
-               Sample fractional_phase_part = (phase & fractional_part_mask) / binary_scaling_factor;
-               
-        // Linearly interpolate into the output buffer
-        // using fixed point math
-               output[outsample] = 
-                       input[i] * (1.0f - fractional_phase_part) +
-                       input[i+1] * fractional_phase_part;
-               phase += phi + phi_delta;
-       }
-
-       last_phase = (phase & fractional_part_mask);
-       
-       // playback distance
+       framecnt_t   i = 0;
+
+       double acceleration;
+       double distance = 0.0;
+
+       if (_speed != _target_speed) {
+               acceleration = _target_speed - _speed;
+       } else {
+               acceleration = 0.0;
+       }
+
+       distance = phase[channel];
+
+       if (nframes < 3) {
+               /* no interpolation possible */
+
+               if (input && output) {
+                       for (i = 0; i < nframes; ++i) {
+                               output[i] = input[i];
+                       }
+               }
+
+               return nframes;
+       }
+
+       /* keep this condition out of the inner loop */
+
+       if (input && output) {
+
+               Sample inm1;
+
+               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];
+               }
+
+               for (framecnt_t outsample = 0; outsample < nframes; ++outsample) {
+
+                       float f = floor (distance);
+                       float fractional_phase_part = distance - f;
+
+                       /* get the index into the input we should start with */
+
+                       i = lrintf (f);
+
+                       /* 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 ...
+                          */
+
+                       if (fractional_phase_part >= 1.0) {
+                               fractional_phase_part -= 1.0;
+                               ++i;
+                       }
+
+                       // 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)
+
+                       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])));
+
+                       distance += _speed + acceleration;
+                       inm1 = input[i];
+               }
+
+               i = floor(distance);
+               phase[channel] = distance - floor(distance);
+
+       } 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);
+       }
+
+       return i;
+}
+
+framecnt_t
+CubicMidiInterpolation::distance (framecnt_t nframes, bool roll)
+{
+       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);
+
        return i;
 }