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 <limits>
  21 #include <cstdio>
  22
  23 #include <stdint.h>
  24
  25 #include "ardour/interpolation.h"
  26 #include "ardour/midi_buffer.h"
  27
  28 using namespace ARDOUR;
  29 using std::cerr;
  30 using std::endl;
  31
  32 CubicInterpolation::CubicInterpolation ()
  33         : valid_z_bits (0)
  34 {
  35 }
  36
  37 samplecnt_t
  38 CubicInterpolation::interpolate (int channel, samplecnt_t input_samples, Sample *input, samplecnt_t &  output_samples, Sample *output)
  39 {
  40         assert (input_samples > 0);
  41         assert (output_samples > 0);
  42         assert (input);
  43         assert (output);
  44         assert (phase.size () > channel);
  45
  46         _speed = fabs (_speed);
  47
  48         if (invalid (0)) {
  49
  50                 /* z[0] not set. Two possibilities
  51                  *
  52                  * 1) we have just been constructed or ::reset()
  53                  *
  54                  * 2) we were only given 1 sample after construction or
  55                  *    ::reset, and stored it in z[1]
  56                  */
  57
  58                 if (invalid (1)) {
  59
  60                         /* first call after construction or after ::reset */
  61
  62                         switch (input_samples) {
  63                         case 1:
  64                                 /* store one sample for use next time. We don't
  65                                  * have enough points to interpolate or even
  66                                  * compute the first z[0] value, but keep z[1]
  67                                  * around.
  68                                  */
  69                                 z[1] = input[0]; validate (1);
  70                                 output_samples = 0;
  71                                 return 0;
  72                         case 2:
  73                                 /* store two samples for use next time, and
  74                                  * compute a value for z[0] that will maintain
  75                                  * the slope of the first actual segment. We
  76                                  * still don't have enough samples to interpolate.
  77                                  */
  78                                 z[0] = input[0] - (input[1] - input[0]); validate (0);
  79                                 z[1] = input[0]; validate (1);
  80                                 z[2] = input[1]; validate (2);
  81                                 output_samples = 0;
  82                                 return 0;
  83                         default:
  84                                 /* We have enough samples to interpolate this time,
  85                                  * but don't have a valid z[0] value because this is the
  86                                  * first call after construction or ::reset.
  87                                  *
  88                                  * First point is based on a requirement to maintain
  89                                  * the slope of the first actual segment
  90                                  */
  91                                 z[0] = input[0] - (input[1] - input[0]); validate (0);
  92                                 break;
  93                         }
  94                 } else {
  95
  96                         /* at least one call since construction or
  97                          * after::reset, since we have z[1] set
  98                          *
  99                          * we can now compute z[0] as required
 100                          */
 101
 102                         z[0] = z[1] - (input[0] - z[1]); validate (0);
 103
 104                         /* we'll check the number of samples we've been given
 105                            in the next switch() statement below, and either
 106                            just save some more samples or actual interpolate
 107                         */
 108                 }
 109
 110                 assert (is_valid (0));
 111         }
 112
 113         switch (input_samples) {
 114         case 1:
 115                 /* one more sample of input. find the right vX to store
 116                    it in, and decide if we're ready to interpolate
 117                 */
 118                 if (invalid (1)) {
 119                         z[1] = input[0]; validate (1);
 120                         /* still not ready to interpolate */
 121                         output_samples = 0;
 122                         return 0;
 123                 } else if (invalid (2)) {
 124                         /* still not ready to interpolate */
 125                         z[2] = input[0]; validate (2);
 126                         output_samples = 0;
 127                         return 0;
 128                 } else if (invalid (3)) {
 129                         z[3] = input[0]; validate (3);
 130                         /* ready to interpolate */
 131                 }
 132                 break;
 133         case 2:
 134                 /* two more samples of input. find the right vX to store
 135                    them in, and decide if we're ready to interpolate
 136                 */
 137                 if (invalid (1)) {
 138                         z[1] = input[0]; validate (1);
 139                         z[2] = input[1]; validate (2);
 140                         /* still not ready to interpolate */
 141                         output_samples = 0;
 142                         return 0;
 143                 } else if (invalid (2)) {
 144                         z[2] = input[0]; validate (2);
 145                         z[3] = input[1]; validate (3);
 146                         /* ready to interpolate */
 147                 } else if (invalid (3)) {
 148                         z[3] = input[0]; validate (3);
 149                         /* ready to interpolate */
 150                 }
 151                 break;
 152
 153         default:
 154                 /* caller has given us at least enough samples to interpolate a
 155                    single value.
 156                 */
 157                 z[1] = input[0]; validate (1);
 158                 z[2] = input[1]; validate (2);
 159                 z[3] = input[2]; validate (3);
 160         }
 161
 162         /* ready to interpolate using z[0], z[1], z[2] and z[3] */
 163
 164         assert (is_valid (0));
 165         assert (is_valid (1));
 166         assert (is_valid (2));
 167         assert (is_valid (3));
 168
 169         /* we can use up to (input_samples - 2) of the input, so compute the
 170          * maximum number of output samples that represents.
 171          *
 172          * Remember that the expected common case here is to be given
 173          * input_samples that is substantially larger than output_samples,
 174          * thus allowing us to always compute output_samples in one call.
 175          */
 176
 177         const samplecnt_t output_from_input = floor ((input_samples - 2) / _speed);
 178
 179         /* limit output to either the caller's requested number or the number
 180          * determined by the input size.
 181          */
 182
 183         const samplecnt_t limit = std::min (output_samples, output_from_input);
 184
 185         samplecnt_t outsample = 0;
 186         double distance = phase[channel];
 187         samplecnt_t used = floor (distance);
 188         samplecnt_t i = 0;
 189
 190         while (outsample < limit) {
 191
 192                 i = floor (distance);
 193
 194                 /* this call may stop the loop from being vectorized */
 195                 float fractional_phase_part = fmod (distance, 1.0);
 196
 197                 /* Cubically interpolate into the output buffer */
 198                 output[outsample++] = z[1] + 0.5f * fractional_phase_part *
 199                         (z[2] - z[0] + fractional_phase_part * (4.0f * z[2] + 2.0f * z[0] - 5.0f * z[1] - z[3] +
 200                                                               fractional_phase_part * (3.0f * (z[1] - z[2]) - z[0] + z[3])));
 201
 202                 distance += _speed;
 203
 204                 z[0] = z[1];
 205                 z[1] = input[i];
 206                 z[2] = input[i+1];
 207                 z[3] = input[i+2];
 208         }
 209
 210         output_samples = outsample;
 211         phase[channel] = fmod (distance, 1.0);
 212         return i - used;
 213 }
 214
 215 void
 216 CubicInterpolation::reset ()
 217 {
 218         Interpolation::reset ();
 219         valid_z_bits = 0;
 220 }
 221
 222 samplecnt_t
 223 CubicInterpolation::distance (samplecnt_t nsamples)
 224 {
 225         assert (phase.size () > 0);
 226         return floor (floor (phase[0]) + (_speed * nsamples));
 227 }