Function InterpolateAudio

libraries/lib-math/InterpolateAudio.cpp:82–204

Here's the main interpolate function, using Least Squares AutoRegression (LSAR):

// Here's the main interpolate function, using
// Least Squares AutoRegression (LSAR):
void InterpolateAudio(float *buffer, const size_t len,
                      size_t firstBad, size_t numBad)
{
   const auto N = len;

   wxASSERT(len > 0 &&
            firstBad >= 0 &&
            numBad < len &&
            firstBad+numBad <= len);

   if(numBad >= len)
      return;  //should never have been called!

   if (firstBad == 0) {
      // The algorithm below has a weird asymmetry in that it
      // performs poorly when interpolating to the left.  If
      // we're asked to interpolate the left side of a buffer,
      // we just reverse the problem and try it that way.
      Floats buffer2{ len };
      for(size_t i=0; i<len; i++)
         buffer2[len-1-i] = buffer[i];
      InterpolateAudio(buffer2.get(), len, len-numBad, numBad);
      for(size_t i=0; i<len; i++)
         buffer[len-1-i] = buffer2[i];
      return;
   }

   Vector s(len, buffer);

   // Choose P, the order of the autoregression equation
   const int IP =
      imin(imin(numBad * 3, 50), imax(firstBad - 1, len - (firstBad + numBad) - 1));

   if (IP < 3 || IP >= (int)N) {
      LinearInterpolateAudio(buffer, len, firstBad, numBad);
      return;
   }

   size_t P(IP);

   // Add a tiny amount of random noise to the input signal -
   // this sounds like a bad idea, but the amount we're adding
   // is only about 1 bit in 16-bit audio, and it's an extremely
   // effective way to avoid nearly-singular matrices.  If users
   // run it more than once they get slightly different results;
   // this is sometimes even advantageous.
   for(size_t i=0; i<N; i++)
      s[i] += (rand()-(RAND_MAX/2))/(RAND_MAX*10000.0);

   // Solve for the best autoregression coefficients
   // using a least-squares fit to all of the non-bad
   // data we have in the buffer
   Matrix X(P, P);
   Vector b(P);

   for(size_t i = 0; i + P < len; i++)
      if (i+P < firstBad || i >= (firstBad + numBad))
         for(size_t row=0; row<P; row++) {