libs/ardour/interpolation.cc

   1 /*
   2     Copyright (C) 2012 Paul Davis
   3
   4     This program is free software; you can redistribute it and/or modify
   5     it under the terms of the GNU General Public License as published by
   6     the Free Software Foundation; either version 2 of the License, or
   7     (at your option) any later version.
   8
   9     This program is distributed in the hope that it will be useful,
  10     but WITHOUT ANY WARRANTY; without even the implied warranty of
  11     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  12     GNU General Public License for more details.
  13
  14     You should have received a copy of the GNU General Public License
  15     along with this program; if not, write to the Free Software
  16     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  17
  18 */
  19
  20 #include <stdint.h>
  21 #include <cstdio>
  22
  23 #include "ardour/interpolation.h"
  24
  25 using namespace ARDOUR;
  26
  27
  28 framecnt_t
  29 LinearInterpolation::interpolate (int channel, framecnt_t nframes, Sample *input, Sample *output)
  30 {
  31         // index in the input buffers
  32         framecnt_t i = 0;
  33
  34         double acceleration = 0;
  35
  36         if (_speed != _target_speed) {
  37                 acceleration = _target_speed - _speed;
  38         }
  39
  40         for (framecnt_t outsample = 0; outsample < nframes; ++outsample) {
  41                 double const d = phase[channel] + outsample * (_speed + acceleration);
  42                 i = floor(d);
  43                 Sample fractional_phase_part = d - i;
  44                 if (fractional_phase_part >= 1.0) {
  45                         fractional_phase_part -= 1.0;
  46                         i++;
  47                 }
  48
  49                 if (input && output) {
  50                 // Linearly interpolate into the output buffer
  51                         output[outsample] =
  52                                 input[i] * (1.0f - fractional_phase_part) +
  53                                 input[i+1] * fractional_phase_part;
  54                 }
  55         }
  56
  57         double const distance = phase[channel] + nframes * (_speed + acceleration);
  58         i = floor(distance);
  59         phase[channel] = distance - i;
  60         return i;
  61 }
  62
  63 framecnt_t
  64 CubicInterpolation::interpolate (int channel, framecnt_t nframes, Sample *input, Sample *output)
  65 {
  66     // index in the input buffers
  67     framecnt_t   i = 0;
  68
  69     double acceleration;
  70     double distance = 0.0;
  71
  72     if (_speed != _target_speed) {
  73         acceleration = _target_speed - _speed;
  74     } else {
  75             acceleration = 0.0;
  76     }
  77
  78     distance = phase[channel];
  79
  80     if (nframes < 3) {
  81             /* no interpolation possible */
  82
  83             for (i = 0; i < nframes; ++i) {
  84                     output[i] = input[i];
  85             }
  86
  87             return nframes;
  88     }
  89
  90     /* keep this condition out of the inner loop */
  91
  92     if (input && output) {
  93
  94             Sample inm1;
  95
  96             if (floor (distance) == 0.0) {
  97                     /* best guess for the fake point we have to add to be able to interpolate at i == 0:
  98                        .... maintain slope of first actual segment ...
  99                     */
 100                     inm1 = input[i] - (input[i+1] - input[i]);
 101             } else {
 102                     inm1 = input[i-1];
 103             }
 104
 105             for (framecnt_t outsample = 0; outsample < nframes; ++outsample) {
 106
 107                     float f = floor (distance);
 108                     float fractional_phase_part = distance - f;
 109
 110                     /* get the index into the input we should start with */
 111
 112                     i = lrintf (f);
 113
 114                     /* fractional_phase_part only reaches 1.0 thanks to float imprecision. In theory
 115                        it should always be < 1.0. If it ever >= 1.0, then bump the index we use
 116                        and back it off. This is the point where we "skip" an entire sample in the
 117                        input, because the phase part has accumulated so much error that we should
 118                        really be closer to the next sample. or something like that ...
 119                     */
 120
 121                     if (fractional_phase_part >= 1.0) {
 122                             fractional_phase_part -= 1.0;
 123                             ++i;
 124                     }
 125
 126                     // Cubically interpolate into the output buffer: keep this inlined for speed and rely on compiler
 127                     // optimization to take care of the rest
 128                     // shamelessly ripped from Steve Harris' swh-plugins (ladspa-util.h)
 129
 130                     output[outsample] = input[i] + 0.5f * fractional_phase_part * (input[i+1] - inm1 +
 131                                                           fractional_phase_part * (4.0f * input[i+1] + 2.0f * inm1 - 5.0f * input[i] - input[i+2] +
 132                                                                 fractional_phase_part * (3.0f * (input[i] - input[i+1]) - inm1 + input[i+2])));
 133
 134                     distance += _speed + acceleration;
 135                     inm1 = input[i];
 136             }
 137
 138             i = floor(distance);
 139             phase[channel] = distance - floor(distance);
 140
 141     } else {
 142             /* used to calculate play-distance with acceleration (silent roll)
 143              * (use same algorithm as real playback for identical rounding/floor'ing)
 144              */
 145             for (framecnt_t outsample = 0; outsample < nframes; ++outsample) {
 146                     distance += _speed + acceleration;
 147             }
 148             i = floor(distance);
 149     }
 150
 151     return i;
 152 }