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