Audacity 3.2.0
InterpolateAudio.cpp
Go to the documentation of this file.
1/**********************************************************************
2
3 Audacity: A Digital Audio Editor
4
5 InterpolateAudio.cpp
6
7 Dominic Mazzoni
8
9**********************************************************************/
10
11#include "InterpolateAudio.h"
12
13#include <math.h>
14#include <stdlib.h>
15
16#include <wx/defs.h>
17
18#include "Matrix.h"
19
20static inline int imin(int x, int y)
21{
22 return x<y? x: y;
23}
24
25static inline int imax(int x, int y)
26{
27 return x>y? x: y;
28}
29
30// This function is a really dumb, simple way to interpolate audio,
31// if the more general InterpolateAudio function below doesn't have
32// enough data to work with. If the bad samples are in the middle,
33// it's literally linear. If it's on either edge, we add some decay
34// back to zero.
35static void LinearInterpolateAudio(float *buffer, int len,
37{
38 int i;
39
40 float decay = 0.9f;
41
42 if (firstBad==0) {
43 float delta = buffer[numBad] - buffer[numBad+1];
44 float value = buffer[numBad];
45 i = numBad - 1;
46 while (i >= 0) {
47 value += delta;
48 buffer[i] = value;
49 value *= decay;
50 delta *= decay;
51 i--;
52 }
53 }
54 else if (firstBad + numBad == len) {
55 float delta = buffer[firstBad-1] - buffer[firstBad-2];
56 float value = buffer[firstBad-1];
57 i = firstBad;
58 while (i < firstBad + numBad) {
59 value += delta;
60 buffer[i] = value;
61 value *= decay;
62 delta *= decay;
63 i++;
64 }
65 }
66 else {
67 float v1 = buffer[firstBad-1];
69 float value = v1;
70 float delta = (v2 - v1) / (numBad+1);
71 i = firstBad;
72 while (i < firstBad + numBad) {
73 value += delta;
74 buffer[i] = value;
75 i++;
76 }
77 }
78}
79
80// Here's the main interpolate function, using
81// Least Squares AutoRegression (LSAR):
82void InterpolateAudio(float *buffer, const size_t len,
84{
85 const auto N = len;
86
87 wxASSERT(len > 0 &&
88 firstBad >= 0 &&
89 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];
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) {
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.
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,
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,
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...
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)
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
Matrix MatrixSubset(const Matrix &input, unsigned startRow, unsigned numRows, unsigned startCol, unsigned numCols)
Definition: Matrix.cpp:256
Matrix TransposeMatrix(const Matrix &other)
Definition: Matrix.cpp:280
bool InvertMatrix(const Matrix &input, Matrix &Minv)
Definition: Matrix.cpp:289
Matrix MatrixConcatenateCols(const Matrix &left, const Matrix &right)
Definition: Matrix.cpp:267
Vector VectorSubset(const Vector &other, unsigned start, unsigned len)
Definition: Matrix.cpp:170
Matrix MatrixMultiply(const Matrix &left, const Matrix &right)
Definition: Matrix.cpp:243
General routine to interpolate (or even extrapolate small amounts) audio when a few of the samples ar...
#define P(T)
Definition: ToChars.cpp:56
#define A(N)
Definition: ToChars.cpp:62
Holds a matrix of doubles and supports arithmetic, subsetting, and matrix inversion....
Definition: Matrix.h:58
unsigned Cols() const
Definition: Matrix.h:69
unsigned Rows() const
Definition: Matrix.h:68
Holds a matrix of doubles and supports arithmetic operations, including Vector-Matrix operations....
Definition: Matrix.h:34