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 #include "ardour/midi_buffer.h"
  25
  26 using namespace ARDOUR;
  27
  28
  29 framecnt_t
  30 LinearInterpolation::interpolate (int channel, framecnt_t nframes, Sample *input, Sample *output)
  31 {
  32         // index in the input buffers
  33         framecnt_t i = 0;
  34
  35         double acceleration = 0;
  36
  37         if (_speed != _target_speed) {
  38                 acceleration = _target_speed - _speed;
  39         }
  40
  41         for (framecnt_t outsample = 0; outsample < nframes; ++outsample) {
  42                 double const d = phase[channel] + outsample * (_speed + acceleration);
  43                 i = floor(d);
  44                 Sample fractional_phase_part = d - i;
  45                 if (fractional_phase_part >= 1.0) {
  46                         fractional_phase_part -= 1.0;
  47                         i++;
  48                 }
  49
  50                 if (input && output) {
  51                         // Linearly interpolate into the output buffer
  52                         output[outsample] =
  53                                 input[i] * (1.0f - fractional_phase_part) +
  54                                 input[i+1] * fractional_phase_part;
  55                 }
  56         }
  57
  58         double const distance = phase[channel] + nframes * (_speed + acceleration);
  59         i = floor(distance);
  60         phase[channel] = distance - i;
  61         return i;
  62 }
  63
  64 framecnt_t
  65 CubicInterpolation::interpolate (int channel, framecnt_t nframes, Sample *input, Sample *output)
  66 {
  67         // index in the input buffers
  68         framecnt_t   i = 0;
  69
  70         double acceleration;
  71         double distance = 0.0;
  72
  73         if (_speed != _target_speed) {
  74                 acceleration = _target_speed - _speed;
  75         } else {
  76                 acceleration = 0.0;
  77         }
  78
  79         distance = phase[channel];
  80
  81         if (nframes < 3) {
  82                 /* no interpolation possible */
  83
  84                 if (input && output) {
  85                         for (i = 0; i < nframes; ++i) {
  86                                 output[i] = input[i];
  87                         }
  88                 }
  89
  90                 return nframes;
  91         }
  92
  93         /* keep this condition out of the inner loop */
  94
  95         if (input && output) {
  96
  97                 Sample inm1;
  98
  99                 if (floor (distance) == 0.0) {
 100                         /* best guess for the fake point we have to add to be able to interpolate at i == 0:
 101                            .... maintain slope of first actual segment ...
 102                            */
 103                         inm1 = input[i] - (input[i+1] - input[i]);
 104                 } else {
 105                         inm1 = input[i-1];
 106                 }
 107
 108                 for (framecnt_t outsample = 0; outsample < nframes; ++outsample) {
 109
 110                         float f = floor (distance);
 111                         float fractional_phase_part = distance - f;
 112
 113                         /* get the index into the input we should start with */
 114
 115                         i = lrintf (f);
 116
 117                         /* fractional_phase_part only reaches 1.0 thanks to float imprecision. In theory
 118                            it should always be < 1.0. If it ever >= 1.0, then bump the index we use
 119                            and back it off. This is the point where we "skip" an entire sample in the
 120                            input, because the phase part has accumulated so much error that we should
 121                            really be closer to the next sample. or something like that ...
 122                            */
 123
 124                         if (fractional_phase_part >= 1.0) {
 125                                 fractional_phase_part -= 1.0;
 126                                 ++i;
 127                         }
 128
 129                         // Cubically interpolate into the output buffer: keep this inlined for speed and rely on compiler
 130                         // optimization to take care of the rest
 131                         // shamelessly ripped from Steve Harris' swh-plugins (ladspa-util.h)
 132
 133                         output[outsample] = input[i] + 0.5f * fractional_phase_part * (input[i+1] - inm1 +
 134                                         fractional_phase_part * (4.0f * input[i+1] + 2.0f * inm1 - 5.0f * input[i] - input[i+2] +
 135                                                 fractional_phase_part * (3.0f * (input[i] - input[i+1]) - inm1 + input[i+2])));
 136
 137                         distance += _speed + acceleration;
 138                         inm1 = input[i];
 139                 }
 140
 141                 i = floor(distance);
 142                 phase[channel] = distance - floor(distance);
 143
 144         } else {
 145                 /* used to calculate play-distance with acceleration (silent roll)
 146                  * (use same algorithm as real playback for identical rounding/floor'ing)
 147                  */
 148                 for (framecnt_t outsample = 0; outsample < nframes; ++outsample) {
 149                         distance += _speed + acceleration;
 150                 }
 151                 i = floor(distance);
 152         }
 153
 154         return i;
 155 }
 156
 157 framecnt_t
 158 CubicMidiInterpolation::distance (framecnt_t nframes, bool roll)
 159 {
 160         assert(phase.size() == 1);
 161
 162         framecnt_t i = 0;
 163
 164         double acceleration;
 165         double distance = 0.0;
 166
 167         if (nframes < 3) {
 168                 return nframes;
 169         }
 170
 171         if (_speed != _target_speed) {
 172                 acceleration = _target_speed - _speed;
 173         } else {
 174                 acceleration = 0.0;
 175         }
 176
 177         distance = phase[0];
 178
 179         for (framecnt_t outsample = 0; outsample < nframes; ++outsample) {
 180                 distance += _speed + acceleration;
 181         }
 182
 183         if (roll) {
 184                 phase[0] = distance - floor(distance);
 185         }
 186
 187         i = floor(distance);
 188
 189         return i;
 190 }