Audacity  2.2.2
Functions
InterpolateAudio.h File Reference
#include "Audacity.h"

Go to the source code of this file.

Functions

void AUDACITY_DLL_API InterpolateAudio (float *buffer, size_t len, size_t firstBad, size_t numBad)
 

Function Documentation

void AUDACITY_DLL_API InterpolateAudio ( float *  buffer,
size_t  len,
size_t  firstBad,
size_t  numBad 
)

Definition at line 82 of file InterpolateAudio.cpp.

References Matrix::Cols(), imax(), imin(), InterpolateAudio(), InvertMatrix(), LinearInterpolateAudio(), MatrixConcatenateCols(), MatrixMultiply(), MatrixSubset(), Matrix::Rows(), TransposeMatrix(), VectorConcatenate(), and VectorSubset().

Referenced by InterpolateAudio(), and EffectRepair::ProcessOne().

84 {
85  const auto N = len;
86 
87  wxASSERT(len > 0 &&
88  firstBad >= 0 &&
89  numBad < len &&
90  firstBad+numBad <= len);
91 
92  if(numBad >= len)
93  return; //should never have been called!
94 
95  if (firstBad == 0) {
96  // The algorithm below has a weird asymmetry in that it
97  // performs poorly when interpolating to the left. If
98  // we're asked to interpolate the left side of a buffer,
99  // we just reverse the problem and try it that way.
100  Floats buffer2{ len };
101  for(size_t i=0; i<len; i++)
102  buffer2[len-1-i] = buffer[i];
103  InterpolateAudio(buffer2.get(), len, len-numBad, numBad);
104  for(size_t i=0; i<len; i++)
105  buffer[len-1-i] = buffer2[i];
106  return;
107  }
108 
109  Vector s(len, buffer);
110 
111  // Choose P, the order of the autoregression equation
112  const int IP =
113  imin(imin(numBad * 3, 50), imax(firstBad - 1, len - (firstBad + numBad) - 1));
114 
115  if (IP < 3 || IP >= (int)N) {
116  LinearInterpolateAudio(buffer, len, firstBad, numBad);
117  return;
118  }
119 
120  size_t P(IP);
121 
122  // Add a tiny amount of random noise to the input signal -
123  // this sounds like a bad idea, but the amount we're adding
124  // is only about 1 bit in 16-bit audio, and it's an extremely
125  // effective way to avoid nearly-singular matrices. If users
126  // run it more than once they get slightly different results;
127  // this is sometimes even advantageous.
128  for(size_t i=0; i<N; i++)
129  s[i] += (rand()-(RAND_MAX/2))/(RAND_MAX*10000.0);
130 
131  // Solve for the best autoregression coefficients
132  // using a least-squares fit to all of the non-bad
133  // data we have in the buffer
134  Matrix X(P, P);
135  Vector b(P);
136 
137  for(size_t i = 0; i + P < len; i++)
138  if (i+P < firstBad || i >= (firstBad + numBad))
139  for(size_t row=0; row<P; row++) {
140  for(size_t col=0; col<P; col++)
141  X[row][col] += (s[i+row] * s[i+col]);
142  b[row] += s[i+P] * s[i+row];
143  }
144 
145  Matrix Xinv(P, P);
146  if (!InvertMatrix(X, Xinv)) {
147  // The matrix is singular! Fall back on linear...
148  // In practice I have never seen this happen if
149  // we add the tiny bit of random noise.
150  LinearInterpolateAudio(buffer, len, firstBad, numBad);
151  return;
152  }
153 
154  // This vector now contains the autoregression coefficients
155  const Vector &a = Xinv * b;
156 
157  // Create a matrix (a "Toeplitz" matrix, as it turns out)
158  // which encodes the autoregressive relationship between
159  // elements of the sequence.
160  Matrix A(N-P, N);
161  for(size_t row=0; row<N-P; row++) {
162  for(size_t col=0; col<P; col++)
163  A[row][row+col] = -a[col];
164  A[row][row+P] = 1;
165  }
166 
167  // Split both the Toeplitz matrix and the signal into
168  // two pieces. Note that this code could be made to
169  // work even in the case where the "bad" samples are
170  // not contiguous, but currently it assumes they are.
171  // "u" is for unknown (bad)
172  // "k" is for known (good)
173  Matrix Au = MatrixSubset(A, 0, N-P, firstBad, numBad);
174  Matrix A_left = MatrixSubset(A, 0, N-P, 0, firstBad);
175  Matrix A_right = MatrixSubset(A, 0, N-P,
176  firstBad+numBad, N-(firstBad+numBad));
177  Matrix Ak = MatrixConcatenateCols(A_left, A_right);
178 
179  const Vector &s_left = VectorSubset(s, 0, firstBad);
180  const Vector &s_right = VectorSubset(s, firstBad+numBad,
181  N-(firstBad+numBad));
182  const Vector &sk = VectorConcatenate(s_left, s_right);
183 
184  // Do some linear algebra to find the best possible
185  // values that fill in the "bad" area
186  Matrix AuT = TransposeMatrix(Au);
187  Matrix X1 = MatrixMultiply(AuT, Au);
188  Matrix X2(X1.Rows(), X1.Cols());
189  if (!InvertMatrix(X1, X2)) {
190  // The matrix is singular! Fall back on linear...
191  LinearInterpolateAudio(buffer, len, firstBad, numBad);
192  return;
193  }
194  Matrix X2b = X2 * -1.0;
195  Matrix X3 = MatrixMultiply(X2b, AuT);
196  Matrix X4 = MatrixMultiply(X3, Ak);
197  // This vector contains our best guess as to the
198  // unknown values
199  const Vector &su = X4 * sk;
200 
201  // Put the results into the return buffer
202  for(size_t i=0; i<numBad; i++)
203  buffer[firstBad+i] = (float)su[i];
204 }
static int imin(int x, int y)
static void LinearInterpolateAudio(float *buffer, int len, int firstBad, int numBad)
Matrix MatrixMultiply(const Matrix &left, const Matrix &right)
Definition: Matrix.cpp:243
Vector VectorSubset(const Vector &other, unsigned start, unsigned len)
Definition: Matrix.cpp:170
Matrix MatrixSubset(const Matrix &input, unsigned startRow, unsigned numRows, unsigned startCol, unsigned numCols)
Definition: Matrix.cpp:256
Matrix MatrixConcatenateCols(const Matrix &left, const Matrix &right)
Definition: Matrix.cpp:267
unsigned Rows() const
Definition: Matrix.h:68
void InterpolateAudio(float *buffer, const size_t len, size_t firstBad, size_t numBad)
static int imax(int x, int y)
Vector VectorConcatenate(const Vector &left, const Vector &right)
Definition: Matrix.cpp:178
Holds a matrix of doubles and supports arithmetic, subsetting, and matrix inversion. Used by InterpolateAudio.
Definition: Matrix.h:57
unsigned Cols() const
Definition: Matrix.h:69
Matrix TransposeMatrix(const Matrix &other)
Definition: Matrix.cpp:280
Holds a matrix of doubles and supports arithmetic operations, including Vector-Matrix operations...
Definition: Matrix.h:33
bool InvertMatrix(const Matrix &input, Matrix &Minv)
Definition: Matrix.cpp:289