libs/ardour/interpolation.cc

   1 /*
   2  * Copyright (C) 2009-2011 David Robillard <d@drobilla.net>
   3  * Copyright (C) 2009-2012 Carl Hetherington <carl@carlh.net>
   4  * Copyright (C) 2009-2017 Paul Davis <paul@linuxaudiosystems.com>
   5  * Copyright (C) 2009 Hans Baier <hansfbaier@googlemail.com>
   6  * Copyright (C) 2013-2017 Robin Gareus <robin@gareus.org>
   7  *
   8  * This program is free software; you can redistribute it and/or modify
   9  * it under the terms of the GNU General Public License as published by
  10  * the Free Software Foundation; either version 2 of the License, or
  11  * (at your option) any later version.
  12  *
  13  * This program is distributed in the hope that it will be useful,
  14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  16  * GNU General Public License for more details.
  17  *
  18  * You should have received a copy of the GNU General Public License along
  19  * with this program; if not, write to the Free Software Foundation, Inc.,
  20  * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
  21  */
  22
  23 #include <limits>
  24 #include <cstdio>
  25
  26 #include <stdint.h>
  27
  28 #include "ardour/interpolation.h"
  29 #include "ardour/midi_buffer.h"
  30
  31 using namespace ARDOUR;
  32 using std::cerr;
  33 using std::endl;
  34
  35 CubicInterpolation::CubicInterpolation ()
  36         : valid_z_bits (0)
  37 {
  38 }
  39
  40 samplecnt_t
  41 CubicInterpolation::interpolate (int channel, samplecnt_t input_samples, Sample *input, samplecnt_t &  output_samples, Sample *output)
  42 {
  43         assert (input_samples > 0);
  44         assert (output_samples > 0);
  45         assert (input);
  46         assert (output);
  47         assert (phase.size () > channel);
  48
  49         _speed = fabs (_speed);
  50
  51         if (invalid (0)) {
  52
  53                 /* z[0] not set. Two possibilities
  54                  *
  55                  * 1) we have just been constructed or ::reset()
  56                  *
  57                  * 2) we were only given 1 sample after construction or
  58                  *    ::reset, and stored it in z[1]
  59                  */
  60
  61                 if (invalid (1)) {
  62
  63                         /* first call after construction or after ::reset */
  64
  65                         switch (input_samples) {
  66                         case 1:
  67                                 /* store one sample for use next time. We don't
  68                                  * have enough points to interpolate or even
  69                                  * compute the first z[0] value, but keep z[1]
  70                                  * around.
  71                                  */
  72                                 z[1] = input[0]; validate (1);
  73                                 output_samples = 0;
  74                                 return 0;
  75                         case 2:
  76                                 /* store two samples for use next time, and
  77                                  * compute a value for z[0] that will maintain
  78                                  * the slope of the first actual segment. We
  79                                  * still don't have enough samples to interpolate.
  80                                  */
  81                                 z[0] = input[0] - (input[1] - input[0]); validate (0);
  82                                 z[1] = input[0]; validate (1);
  83                                 z[2] = input[1]; validate (2);
  84                                 output_samples = 0;
  85                                 return 0;
  86                         default:
  87                                 /* We have enough samples to interpolate this time,
  88                                  * but don't have a valid z[0] value because this is the
  89                                  * first call after construction or ::reset.
  90                                  *
  91                                  * First point is based on a requirement to maintain
  92                                  * the slope of the first actual segment
  93                                  */
  94                                 z[0] = input[0] - (input[1] - input[0]); validate (0);
  95                                 break;
  96                         }
  97                 } else {
  98
  99                         /* at least one call since construction or
 100                          * after::reset, since we have z[1] set
 101                          *
 102                          * we can now compute z[0] as required
 103                          */
 104
 105                         z[0] = z[1] - (input[0] - z[1]); validate (0);
 106
 107                         /* we'll check the number of samples we've been given
 108                            in the next switch() statement below, and either
 109                            just save some more samples or actual interpolate
 110                         */
 111                 }
 112
 113                 assert (is_valid (0));
 114         }
 115
 116         switch (input_samples) {
 117         case 1:
 118                 /* one more sample of input. find the right vX to store
 119                    it in, and decide if we're ready to interpolate
 120                 */
 121                 if (invalid (1)) {
 122                         z[1] = input[0]; validate (1);
 123                         /* still not ready to interpolate */
 124                         output_samples = 0;
 125                         return 0;
 126                 } else if (invalid (2)) {
 127                         /* still not ready to interpolate */
 128                         z[2] = input[0]; validate (2);
 129                         output_samples = 0;
 130                         return 0;
 131                 } else if (invalid (3)) {
 132                         z[3] = input[0]; validate (3);
 133                         /* ready to interpolate */
 134                 }
 135                 break;
 136         case 2:
 137                 /* two more samples of input. find the right vX to store
 138                    them in, and decide if we're ready to interpolate
 139                 */
 140                 if (invalid (1)) {
 141                         z[1] = input[0]; validate (1);
 142                         z[2] = input[1]; validate (2);
 143                         /* still not ready to interpolate */
 144                         output_samples = 0;
 145                         return 0;
 146                 } else if (invalid (2)) {
 147                         z[2] = input[0]; validate (2);
 148                         z[3] = input[1]; validate (3);
 149                         /* ready to interpolate */
 150                 } else if (invalid (3)) {
 151                         z[3] = input[0]; validate (3);
 152                         /* ready to interpolate */
 153                 }
 154                 break;
 155
 156         default:
 157                 /* caller has given us at least enough samples to interpolate a
 158                    single value.
 159                 */
 160                 z[1] = input[0]; validate (1);
 161                 z[2] = input[1]; validate (2);
 162                 z[3] = input[2]; validate (3);
 163         }
 164
 165         /* ready to interpolate using z[0], z[1], z[2] and z[3] */
 166
 167         assert (is_valid (0));
 168         assert (is_valid (1));
 169         assert (is_valid (2));
 170         assert (is_valid (3));
 171
 172         /* we can use up to (input_samples - 2) of the input, so compute the
 173          * maximum number of output samples that represents.
 174          *
 175          * Remember that the expected common case here is to be given
 176          * input_samples that is substantially larger than output_samples,
 177          * thus allowing us to always compute output_samples in one call.
 178          */
 179
 180         const samplecnt_t output_from_input = floor ((input_samples - 2) / _speed);
 181
 182         /* limit output to either the caller's requested number or the number
 183          * determined by the input size.
 184          */
 185
 186         const samplecnt_t limit = std::min (output_samples, output_from_input);
 187
 188         samplecnt_t outsample = 0;
 189         double distance = phase[channel];
 190         samplecnt_t used = floor (distance);
 191         samplecnt_t i = 0;
 192
 193         while (outsample < limit) {
 194
 195                 i = floor (distance);
 196
 197                 /* this call may stop the loop from being vectorized */
 198                 float fractional_phase_part = fmod (distance, 1.0);
 199
 200                 /* Cubically interpolate into the output buffer */
 201                 output[outsample++] = z[1] + 0.5f * fractional_phase_part *
 202                         (z[2] - z[0] + fractional_phase_part * (4.0f * z[2] + 2.0f * z[0] - 5.0f * z[1] - z[3] +
 203                                                               fractional_phase_part * (3.0f * (z[1] - z[2]) - z[0] + z[3])));
 204
 205                 distance += _speed;
 206
 207                 z[0] = z[1];
 208                 z[1] = input[i];
 209                 z[2] = input[i+1];
 210                 z[3] = input[i+2];
 211         }
 212
 213         output_samples = outsample;
 214         phase[channel] = fmod (distance, 1.0);
 215         return i - used;
 216 }
 217
 218 void
 219 CubicInterpolation::reset ()
 220 {
 221         Interpolation::reset ();
 222         valid_z_bits = 0;
 223 }
 224
 225 samplecnt_t
 226 CubicInterpolation::distance (samplecnt_t nsamples)
 227 {
 228         assert (phase.size () > 0);
 229         return floor (floor (phase[0]) + (_speed * nsamples));
 230 }