fix conceptually broken cubic interpolator that failed to handle the i=0 case, and don't rely on inlining for speed there

git-svn-id: svn://localhost/ardour2/branches/3.0@5792 d708f5d6-7413-0410-9779-e7cbd77b26cf

libs/ardour/ardour/interpolation.h

@@ -39,23 +39,11 @@ class Interpolation {
 };
 
 class LinearInterpolation : public Interpolation {
- protected:
-
  public:
      nframes_t interpolate (int channel, nframes_t nframes, Sample* input, Sample* output);
 };
 
 class CubicInterpolation : public Interpolation {
- protected:
-    // 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)
-    {
-        return in + 0.5f * fr * (inp1 - inm1 +
-         fr * (4.0f * inp1 + 2.0f * inm1 - 5.0f * in - inp2 +
-         fr * (3.0f * (in - inp1) - inm1 + inp2)));
-    }
-
  public:
      nframes_t interpolate (int channel, nframes_t nframes, Sample* input, Sample* output);
 };

libs/ardour/interpolation.cc

@@ -61,19 +61,71 @@ CubicInterpolation::interpolate (int channel, nframes_t nframes, Sample *input,
     }
 
     distance = phase[channel];
-    for (nframes_t outsample = 0; outsample < nframes; ++outsample) {
-        i = floor(distance);
-        Sample fractional_phase_part = distance - i;
-        if (fractional_phase_part >= 1.0) {
-            fractional_phase_part -= 1.0;
-            i++;
-        }
 
-        if (input && output) {
+    if (nframes < 3) {
-            // Cubically interpolate into the output buffer
+	    /* no interpolation possible */
-            output[outsample] = cube_interp(fractional_phase_part, input[i-1], input[i], input[i+1], input[i+2]);
+
-        }
+	    for (i = 0; i < nframes; ++i) {
-        distance += _speed + acceleration;
+		    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 (nframes_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;
+	    }
+
+    } else {
+
+	    /* not sure that this is ever utilized - it implies that one of the input/output buffers is missing */
+
+	    for (nframes_t outsample = 0; outsample < nframes; ++outsample) {
+		    distance += _speed + acceleration;
+	    }
     }
 
     i = floor(distance);