libs/ardour/ardour/interpolation.h

   1 #include <math.h>
   2 #include <samplerate.h>
   3
   4 #include "ardour/types.h"
   5
   6 #ifndef __interpolation_h__
   7 #define __interpolation_h__
   8
   9 namespace ARDOUR {
  10
  11 class Interpolation {
  12  protected:
  13      double   _speed, _target_speed;
  14
  15      // the idea is that when the speed is not 1.0, we have to
  16      // interpolate between samples and then we have to store where we thought we were.
  17      // rather than being at sample N or N+1, we were at N+0.8792922
  18      std::vector<double> phase;
  19
  20
  21  public:
  22      Interpolation ()  { _speed = 1.0; _target_speed = 1.0; }
  23      ~Interpolation () { phase.clear(); }
  24
  25      void set_speed (double new_speed)          { _speed = new_speed; _target_speed = new_speed; }
  26      void set_target_speed (double new_speed)   { _target_speed = new_speed; }
  27
  28      double target_speed()          const { return _target_speed; }
  29      double speed()                 const { return _speed; }
  30
  31      void add_channel_to (int input_buffer_size, int output_buffer_size) { phase.push_back (0.0); }
  32      void remove_channel_from () { phase.pop_back (); }
  33
  34      void reset () {
  35          for (size_t i = 0; i < phase.size(); i++) {
  36               phase[i] = 0.0;
  37           }
  38      }
  39 };
  40
  41 // 40.24 fixpoint math
  42 #define FIXPOINT_ONE 0x1000000
  43
  44 class FixedPointLinearInterpolation : public Interpolation {
  45     protected:
  46     /// speed in fixed point math
  47     uint64_t      phi;
  48
  49     /// target speed in fixed point math
  50     uint64_t      target_phi;
  51
  52     std::vector<uint64_t> last_phase;
  53
  54     // Fixed point is just an integer with an implied scaling factor.
  55     // In 40.24 the scaling factor is 2^24 = 16777216,
  56     // so a value of 10*2^24 (in integer space) is equivalent to 10.0.
  57     //
  58     // The advantage is that addition and modulus [like x = (x + y) % 2^40]
  59     // have no rounding errors and no drift, and just require a single integer add.
  60     // (swh)
  61
  62     static const int64_t fractional_part_mask  = 0xFFFFFF;
  63     static const Sample  binary_scaling_factor = 16777216.0f;
  64
  65     public:
  66
  67         FixedPointLinearInterpolation () : phi (FIXPOINT_ONE), target_phi (FIXPOINT_ONE) {}
  68
  69         void set_speed (double new_speed) {
  70             target_phi = (uint64_t) (FIXPOINT_ONE * fabs(new_speed));
  71             phi = target_phi;
  72         }
  73
  74         uint64_t get_phi() { return phi; }
  75         uint64_t get_target_phi() { return target_phi; }
  76         uint64_t get_last_phase() { assert(last_phase.size()); return last_phase[0]; }
  77         void set_last_phase(uint64_t phase) { assert(last_phase.size()); last_phase[0] = phase; }
  78
  79         void add_channel_to (int input_buffer_size, int output_buffer_size);
  80         void remove_channel_from ();
  81
  82         nframes_t interpolate (int channel, nframes_t nframes, Sample* input, Sample* output);
  83         void reset ();
  84 };
  85
  86 class LinearInterpolation : public Interpolation {
  87  protected:
  88
  89  public:
  90      nframes_t interpolate (int channel, nframes_t nframes, Sample* input, Sample* output);
  91 };
  92
  93 class CubicInterpolation : public Interpolation {
  94  protected:
  95     // shamelessly ripped from Steve Harris' swh-plugins
  96     static inline float cube_interp(const float fr, const float inm1, const float
  97                                     in, const float inp1, const float inp2)
  98     {
  99         return in + 0.5f * fr * (inp1 - inm1 +
 100          fr * (4.0f * inp1 + 2.0f * inm1 - 5.0f * in - inp2 +
 101          fr * (3.0f * (in - inp1) - inm1 + inp2)));
 102     }
 103
 104  public:
 105      nframes_t interpolate (int channel, nframes_t nframes, Sample* input, Sample* output);
 106 };
 107
 108
 109 /**
 110  * @class SplineInterpolation
 111  *
 112  * @brief interpolates using cubic spline interpolation over an input period
 113  *
 114  * Splines are piecewise cubic functions between each samples,
 115  * where the cubic polynomials' values, first and second derivatives are equal
 116  * on each sample point.
 117  *
 118  * Those conditions are equivalent of solving the linear system of equations
 119  * defined by the matrix equation (all indices are zero-based):
 120  *  A * M = d
 121  *
 122  * where A has (n-2) rows and (n-2) columns
 123  *
 124  *  [ 4 1 0 0 ... 0 0 0 0 ]   [ M[1]   ]   [ 6*y[0] - 12*y[1] + 6*y[2] ]
 125  *  [ 1 4 1 0 ... 0 0 0 0 ]   [ M[2]   ]   [ 6*y[1] - 12*y[2] + 6*y[3] ]
 126  *  [ 0 1 4 1 ... 0 0 0 0 ]   [ M[3]   ]   [ 6*y[2] - 12*y[3] + 6*y[4] ]
 127  *  [ 0 0 1 4 ... 0 0 0 0 ]   [ M[4]   ]   [ 6*y[3] - 12*y[4] + 6*y[5] ]
 128  *            ...           *            =            ...
 129  *  [ 0 0 0 0 ... 4 1 0 0 ]   [ M[n-5] ]   [ 6*y[n-6]- 12*y[n-5] + 6*y[n-4] ]
 130  *  [ 0 0 0 0 ... 1 4 1 0 ]   [ M[n-4] ]   [ 6*y[n-5]- 12*y[n-4] + 6*y[n-3] ]
 131  *  [ 0 0 0 0 ... 0 1 4 1 ]   [ M[n-3] ]   [ 6*y[n-4]- 12*y[n-3] + 6*y[n-2] ]
 132  *  [ 0 0 0 0 ... 0 0 1 4 ]   [ M[n-2] ]   [ 6*y[n-3]- 12*y[n-2] + 6*y[n-1] ]
 133  *
 134  *  For our purpose we use natural splines which means the boundary coefficients
 135  *  M[0] = M[n-1] = 0
 136  *
 137  *  The interpolation polynomial in the i-th interval then has the form
 138  *  p_i(x) = a3 (x - i)^3 + a2 (x - i)^2 + a1 (x - i) + a0
 139  *         = ((a3 * (x - i) + a2) * (x - i) + a1) * (x - i) + a0
 140  *     where
 141  *  a3 = (M[i+1] - M[i]) / 6
 142  *  a2 = M[i] / 2
 143  *  a1 = y[i+1] - y[i] - M[i+1]/6 - M[i]/3
 144  *  a0 = y[i]
 145  *
 146  *  We solve the system by LU-factoring the matrix A:
 147  *  A = L * U:
 148  *
 149  *  [ 4 1 0 0 ... 0 0 0 0 ]   [ 1    0    0    0   ... 0      0      0      0 ]   [ m[0] 1    0    0   ... 0      0      0      ]
 150  *  [ 1 4 1 0 ... 0 0 0 0 ]   [ l[0] 1    0    0   ... 0      0      0      0 ]   [ 0    m[1] 1    0   ... 0      0      0      ]
 151  *  [ 0 1 4 1 ... 0 0 0 0 ]   [ 0    l[1] 1    0   ... 0      0      0      0 ]   [ 0    0    m[2] 1   ... 0      0      0      ]
 152  *  [ 0 0 1 4 ... 0 0 0 0 ]   [ 0    0    l[2] 1   ... 0      0      0      0 ]                        ...
 153  *            ...           =                     ...                          *  [ 0    0    0    0   ... 0      0      0      ]
 154  *  [ 0 0 0 0 ... 4 1 0 0 ]   [ 0    0    0    0   ... 1      0      0      0 ]   [ 0    0    0    0   ... 1      0      0      ]
 155  *  [ 0 0 0 0 ... 1 4 1 0 ]   [ 0    0    0    0   ... l[n-6] 1      0      0 ]   [ 0    0    0    0   ... m[n-5] 1      0      ]
 156  *  [ 0 0 0 0 ... 0 1 4 1 ]   [ 0    0    0    0   ... 0      l[n-5] 1      0 ]   [ 0    0    0    0   ... 0      m[n-4] 1      ]
 157  *  [ 0 0 0 0 ... 0 0 1 4 ]   [ 0    0    0    0   ... 0      0      l[n-4] 1 ]   [ 0    0    0    0   ... 0      0      m[n-3] ]
 158  *
 159  *  where the l[i] and m[i] can be precomputed.
 160  *
 161  *  Then we solve the system A * M = L(UM) = d by first solving the system
 162  *    L * t = d
 163  *
 164  *    [ 1    0    0    0   ... 0      0      0      0 ]    [ t[0]   ]    [ 6*y[0] - 12*y[1] + 6*y[2] ]
 165  *    [ l[0] 1    0    0   ... 0      0      0      0 ]    [ t[1]   ]    [ 6*y[1] - 12*y[2] + 6*y[3] ]
 166  *    [ 0    l[1] 1    0   ... 0      0      0      0 ]    [ t[2]   ]    [ 6*y[2] - 12*y[3] + 6*y[4] ]
 167  *    [ 0    0    l[2] 1   ... 0      0      0      0 ]    [ t[3]   ]    [ 6*y[3] - 12*y[4] + 6*y[5] ]
 168  *                        ...                           *             =             ...
 169  *    [ 0    0    0    0   ... 1      0      0      0 ]    [ t[n-6] ]    [ 6*y[n-6]- 12*y[n-5] + 6*y[n-4] ]
 170  *    [ 0    0    0    0   ... l[n-6] 1      0      0 ]    [ t[n-5] ]    [ 6*y[n-5]- 12*y[n-4] + 6*y[n-3] ]
 171  *    [ 0    0    0    0   ... 0      l[n-5] 1      0 ]    [ t[n-4] ]    [ 6*y[n-4]- 12*y[n-3] + 6*y[n-2] ]
 172  *    [ 0    0    0    0   ... 0      0      l[n-4] 1 ]    [ t[n-3] ]    [ 6*y[n-3]- 12*y[n-2] + 6*y[n-1] ]
 173  *
 174  *
 175  *  and then
 176  *    U * M = t
 177  *
 178  *  [ m[0] 1    0    0   ... 0      0      0      ]   [ M[1]   ]    [ t[0]   ]
 179  *  [ 0    m[1] 1    0   ... 0      0      0      ]   [ M[2]   ]    [ t[1]   ]
 180  *  [ 0    0    m[2] 1   ... 0      0      0      ]   [ M[3]   ]    [ t[2]   ]
 181  *                       ...                          [ M[4]   ]    [ t[3]   ]
 182  *  [ 0    0    0    0   ... 0      0      0      ] *            =
 183  *  [ 0    0    0    0   ... 1      0      0      ]   [ M[n-5] ]    [ t[n-6] ]
 184  *  [ 0    0    0    0   ... m[n-5] 1      0      ]   [ M[n-4] ]    [ t[n-5] ]
 185  *  [ 0    0    0    0   ... 0      m[n-4] 1      ]   [ M[n-3] ]    [ t[n-4] ]
 186  *  [ 0    0    0    0   ... 0      0      m[n-3] ]   [ M[n-2] ]    [ t[n-3] ]
 187  *
 188  */
 189 class SplineInterpolation : public Interpolation {
 190  protected:
 191     double _l[19], _m[20];
 192
 193     inline double l(nframes_t i) {  return (i >= 19) ? _l[18] : _l[i]; }
 194     inline double m(nframes_t i) {  return (i >= 20) ? _m[19] : _m[i]; }
 195
 196  public:
 197     SplineInterpolation();
 198     nframes_t interpolate (int channel, nframes_t nframes, Sample* input, Sample* output);
 199 };
 200
 201 class LibSamplerateInterpolation : public Interpolation {
 202  protected:
 203     std::vector<SRC_STATE*>  state;
 204     std::vector<SRC_DATA*>   data;
 205
 206     int        error;
 207
 208     void reset_state ();
 209
 210  public:
 211         LibSamplerateInterpolation ();
 212         ~LibSamplerateInterpolation ();
 213
 214         void   set_speed (double new_speed);
 215         void   set_target_speed (double new_speed)   {}
 216         double speed ()                        const { return _speed;      }
 217
 218         void add_channel_to (int input_buffer_size, int output_buffer_size);
 219         void remove_channel_from ();
 220
 221         nframes_t interpolate (int channel, nframes_t nframes, Sample* input, Sample* output);
 222         void reset() { reset_state (); }
 223 };
 224
 225 } // namespace ARDOUR
 226
 227 #endif