1
0
mirror of https://github.com/cookiengineer/audacity synced 2025-06-24 16:20:05 +02:00

Remove naked new[] in: linear algebra functions for Repair

This commit is contained in:
Paul Licameli 2017-02-21 19:36:32 -05:00
parent cb05476c45
commit f858d97352
4 changed files with 137 additions and 183 deletions

@ -15,6 +15,7 @@
#include "InterpolateAudio.h"
#include "Matrix.h"
#include "SampleFormat.h"
static inline int imin(int x, int y)
{
@ -78,11 +79,10 @@ static void LinearInterpolateAudio(float *buffer, int len,
// Here's the main interpolate function, using
// Least Squares AutoRegression (LSAR):
void InterpolateAudio(float *buffer, int len,
int firstBad, int numBad)
void InterpolateAudio(float *buffer, const size_t len,
size_t firstBad, size_t numBad)
{
int N = len;
int i, row, col;
const auto N = len;
wxASSERT(len > 0 &&
firstBad >= 0 &&
@ -97,34 +97,35 @@ void InterpolateAudio(float *buffer, int len,
// 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.
float *buffer2 = new float[len];
for(i=0; i<len; i++)
Floats buffer2{ len };
for(size_t i=0; i<len; i++)
buffer2[len-1-i] = buffer[i];
InterpolateAudio(buffer2, len, len-numBad, numBad);
for(i=0; i<len; i++)
InterpolateAudio(buffer2.get(), len, len-numBad, numBad);
for(size_t i=0; i<len; i++)
buffer[len-1-i] = buffer2[i];
delete[] buffer2;
return;
}
Vector s(len, buffer);
// Choose P, the order of the autoregression equation
int P = imin(numBad * 3, 50);
P = imin(P, imax(firstBad - 1, len - (firstBad + numBad) - 1));
const int IP =
imin(imin(numBad * 3, 50), imax(firstBad - 1, len - (firstBad + numBad) - 1));
if (P < 3) {
if (IP < 3 || IP >= 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(i=0; i<N; i++)
for(size_t i=0; i<N; i++)
s[i] += (rand()-(RAND_MAX/2))/(RAND_MAX*10000.0);
// Solve for the best autoregression coefficients
@ -133,10 +134,10 @@ void InterpolateAudio(float *buffer, int len,
Matrix X(P, P);
Vector b(P);
for(i=0; i<len-P; i++)
for(size_t i = 0; i + P < len; i++)
if (i+P < firstBad || i >= (firstBad + numBad))
for(row=0; row<P; row++) {
for(col=0; col<P; col++)
for(size_t row=0; row<P; row++) {
for(size_t col=0; col<P; col++)
X[row][col] += (s[i+row] * s[i+col]);
b[row] += s[i+P] * s[i+row];
}
@ -151,14 +152,14 @@ void InterpolateAudio(float *buffer, int len,
}
// This vector now contains the autoregression coefficients
Vector a = Xinv * b;
const Vector &a = Xinv * b;
// Create a matrix (a "Toeplitz" matrix, as it turns out)
// which encodes the autoregressive relationship between
// elements of the sequence.
Matrix A(N-P, N);
for(row=0; row<N-P; row++) {
for(col=0; col<P; col++)
for(size_t row=0; row<N-P; row++) {
for(size_t col=0; col<P; col++)
A[row][row+col] = -a[col];
A[row][row+P] = 1;
}
@ -175,10 +176,10 @@ void InterpolateAudio(float *buffer, int len,
firstBad+numBad, N-(firstBad+numBad));
Matrix Ak = MatrixConcatenateCols(A_left, A_right);
Vector s_left = VectorSubset(s, 0, firstBad);
Vector s_right = VectorSubset(s, firstBad+numBad,
const Vector &s_left = VectorSubset(s, 0, firstBad);
const Vector &s_right = VectorSubset(s, firstBad+numBad,
N-(firstBad+numBad));
Vector sk = VectorConcatenate(s_left, s_right);
const Vector &sk = VectorConcatenate(s_left, s_right);
// Do some linear algebra to find the best possible
// values that fill in the "bad" area
@ -195,9 +196,9 @@ void InterpolateAudio(float *buffer, int len,
Matrix X4 = MatrixMultiply(X3, Ak);
// This vector contains our best guess as to the
// unknown values
Vector su = X4 * sk;
const Vector &su = X4 * sk;
// Put the results into the return buffer
for(i=0; i<numBad; i++)
for(size_t i=0; i<numBad; i++)
buffer[firstBad+i] = (float)su[i];
}

@ -35,7 +35,7 @@
// side (6x the number of bad samples on either side is great). However,
// it will work with less data, and with the bad samples on one end or
// the other.
void AUDACITY_DLL_API InterpolateAudio(float *buffer, int len,
int firstBad, int numBad);
void AUDACITY_DLL_API InterpolateAudio(float *buffer, size_t len,
size_t firstBad, size_t numBad);
#endif // __AUDACITY_INTERPOLATE_AUDIO__

@ -17,89 +17,74 @@
Vector::Vector()
{
mCopy = false;
mN = 0;
mData = NULL;
}
Vector::Vector(int len, double *data, bool copy)
Vector::Vector(unsigned len, double *data)
: mN{ len }
, mData(len)
{
mN = len;
mCopy = copy;
if (mCopy || !data) {
mCopy = true;
mData = new double[mN];
int i;
for(i=0; i<mN; i++)
if (data)
mData[i] = data[i];
else
mData[i] = 0.0;
}
else {
mCopy = false;
mData = data;
}
if (data)
std::copy(data, data + len, mData.get());
else
std::fill(mData.get(), mData.get() + len, 0.0);
}
Vector::Vector(unsigned len, float *data)
: mN{ len }
, mData{ len }
{
if (data)
std::copy(data, data + len, mData.get());
else
std::fill(mData.get(), mData.get() + len, 0.0);
}
Vector& Vector::operator=(const Vector &other)
{
wxASSERT(Len() == other.Len());
int i;
for(i=0; i<Len(); i++)
mData[i] = other.mData[i];
std::copy(other.mData.get(), other.mData.get() + mN, mData.get());
return *this;
}
Vector::Vector(const Vector &other)
: mN{ other.Len() }
, mData{ mN }
{
CopyFrom(other);
}
void Vector::CopyFrom(const Vector &other)
{
mN = other.Len();
mCopy = true;
mData = new double[mN];
int i;
for(i=0; i<mN; i++)
mData[i] = other.mData[i];
std::copy(other.mData.get(), other.mData.get() + mN, mData.get());
}
Vector::~Vector()
{
if (mCopy)
delete[] mData;
}
Vector::Vector(int len, float *data)
void Vector::Reinit(unsigned len)
{
mCopy = true;
mN = len;
mData = new double[mN];
int i;
for(i=0; i<mN; i++)
mData[i] = (double)data[i];
Vector temp(len);
Swap(temp);
}
void Vector::Swap(Vector &that)
{
std::swap(mN, that.mN);
mData.swap(that.mData);
}
double Vector::Sum() const
{
int i;
double sum = 0.0;
for(i=0; i<Len(); i++)
for(unsigned i = 0; i < Len(); i++)
sum += mData[i];
return sum;
}
Matrix::Matrix(int rows, int cols, double **data)
Matrix::Matrix(unsigned rows, unsigned cols, double **data)
: mRows{ rows }
, mCols{ cols }
, mRowVec{ mRows }
{
mRows = rows;
mCols = cols;
mRowVec = new Vector *[mRows];
int i, j;
for(i=0; i<mRows; i++) {
mRowVec[i] = new Vector(mCols);
for(j=0; j<mCols; j++) {
for(unsigned i = 0; i < mRows; i++) {
mRowVec[i].Reinit( mCols );
for(unsigned j = 0; j < mCols; j++) {
if (data)
(*this)[i][j] = data[i][j];
else
@ -123,44 +108,26 @@ void Matrix::CopyFrom(const Matrix &other)
{
mRows = other.mRows;
mCols = other.mCols;
mRowVec = new Vector *[mRows];
int i;
for(i=0; i<mRows; i++) {
mRowVec[i] = new Vector(mCols);
*mRowVec[i] = *other.mRowVec[i];
mRowVec.reinit(mRows);
for (unsigned i = 0; i < mRows; i++) {
mRowVec[i].Reinit( mCols );
mRowVec[i] = other.mRowVec[i];
}
}
Matrix::~Matrix()
{
int i;
for(i=0; i<mRows; i++)
delete mRowVec[i];
delete[] mRowVec;
}
void Matrix::SwapRows(int i, int j)
void Matrix::SwapRows(unsigned i, unsigned j)
{
Vector *tmp = mRowVec[i];
mRowVec[i] = mRowVec[j];
mRowVec[j] = tmp;
mRowVec[i].Swap(mRowVec[j]);
}
double Matrix::Sum() const
{
int i, j;
double sum = 0.0;
for(i=0; i<Rows(); i++)
for(j=0; j<Cols(); j++)
sum += (*mRowVec[i])[j];
return sum;
}
Matrix IdentityMatrix(int N)
Matrix IdentityMatrix(unsigned N)
{
Matrix M(N, N);
int i;
for(i=0; i<N; i++)
for(unsigned i = 0; i < N; i++)
M[i][i] = 1.0;
return M;
}
@ -169,8 +136,7 @@ Vector operator+(const Vector &left, const Vector &right)
{
wxASSERT(left.Len() == right.Len());
Vector v(left.Len());
int i;
for(i=0; i<left.Len(); i++)
for(unsigned i = 0; i < left.Len(); i++)
v[i] = left[i] + right[i];
return v;
}
@ -179,8 +145,7 @@ Vector operator-(const Vector &left, const Vector &right)
{
wxASSERT(left.Len() == right.Len());
Vector v(left.Len());
int i;
for(i=0; i<left.Len(); i++)
for(unsigned i = 0; i < left.Len(); i++)
v[i] = left[i] - right[i];
return v;
}
@ -189,8 +154,7 @@ Vector operator*(const Vector &left, const Vector &right)
{
wxASSERT(left.Len() == right.Len());
Vector v(left.Len());
int i;
for(i=0; i<left.Len(); i++)
for(unsigned i = 0; i < left.Len(); i++)
v[i] = left[i] * right[i];
return v;
}
@ -198,17 +162,15 @@ Vector operator*(const Vector &left, const Vector &right)
Vector operator*(const Vector &left, double right)
{
Vector v(left.Len());
int i;
for(i=0; i<left.Len(); i++)
for(unsigned i = 0; i < left.Len(); i++)
v[i] = left[i] * right;
return v;
}
Vector VectorSubset(const Vector &other, int start, int len)
Vector VectorSubset(const Vector &other, unsigned start, unsigned len)
{
Vector v(len);
int i;
for(i=0; i<len; i++)
for(unsigned i = 0; i < len; i++)
v[i] = other[start+i];
return v;
}
@ -216,10 +178,9 @@ Vector VectorSubset(const Vector &other, int start, int len)
Vector VectorConcatenate(const Vector& left, const Vector& right)
{
Vector v(left.Len() + right.Len());
int i;
for(i=0; i<left.Len(); i++)
for(unsigned i = 0; i < left.Len(); i++)
v[i] = left[i];
for(i=0; i<right.Len(); i++)
for(unsigned i = 0; i < right.Len(); i++)
v[i + left.Len()] = right[i];
return v;
}
@ -228,10 +189,9 @@ Vector operator*(const Vector &left, const Matrix &right)
{
wxASSERT(left.Len() == right.Rows());
Vector v(right.Cols());
int i, j;
for(i=0; i<right.Cols(); i++) {
for(unsigned i = 0; i < right.Cols(); i++) {
v[i] = 0.0;
for(j=0; j<right.Rows(); j++)
for(unsigned j = 0; j < right.Rows(); j++)
v[i] += left[j] * right[j][i];
}
return v;
@ -241,10 +201,9 @@ Vector operator*(const Matrix &left, const Vector &right)
{
wxASSERT(left.Cols() == right.Len());
Vector v(left.Rows());
int i, j;
for(i=0; i<left.Rows(); i++) {
for(unsigned i = 0; i < left.Rows(); i++) {
v[i] = 0.0;
for(j=0; j<left.Cols(); j++)
for(unsigned j = 0; j < left.Cols(); j++)
v[i] += left[i][j] * right[j];
}
return v;
@ -255,9 +214,8 @@ Matrix operator+(const Matrix &left, const Matrix &right)
wxASSERT(left.Rows() == right.Rows());
wxASSERT(left.Cols() == right.Cols());
Matrix M(left.Rows(), left.Cols());
int i, j;
for(i=0; i<left.Rows(); i++)
for(j=0; j<left.Cols(); j++)
for(unsigned i = 0; i < left.Rows(); i++)
for(unsigned j = 0; j < left.Cols(); j++)
M[i][j] = left[i][j] + right[i][j];
return M;
}
@ -265,9 +223,8 @@ Matrix operator+(const Matrix &left, const Matrix &right)
Matrix operator*(const Matrix &left, const double right)
{
Matrix M(left.Rows(), left.Cols());
int i, j;
for(i=0; i<left.Rows(); i++)
for(j=0; j<left.Cols(); j++)
for(unsigned i = 0; i < left.Rows(); i++)
for(unsigned j = 0; j < left.Cols(); j++)
M[i][j] = left[i][j] * right;
return M;
}
@ -277,9 +234,8 @@ Matrix ScalarMultiply(const Matrix &left, const Matrix &right)
wxASSERT(left.Rows() == right.Rows());
wxASSERT(left.Cols() == right.Cols());
Matrix M(left.Rows(), left.Cols());
int i, j;
for(i=0; i<left.Rows(); i++)
for(j=0; j<left.Cols(); j++)
for(unsigned i = 0; i < left.Rows(); i++)
for(unsigned j = 0; j < left.Cols(); j++)
M[i][j] = left[i][j] * right[i][j];
return M;
}
@ -288,23 +244,22 @@ Matrix MatrixMultiply(const Matrix &left, const Matrix &right)
{
wxASSERT(left.Cols() == right.Rows());
Matrix M(left.Rows(), right.Cols());
int i, j, k;
for(i=0; i<left.Rows(); i++)
for(j=0; j<right.Cols(); j++) {
for(unsigned i = 0; i < left.Rows(); i++)
for(unsigned j = 0; j < right.Cols(); j++) {
M[i][j] = 0.0;
for(k=0; k<left.Cols(); k++)
for(unsigned k = 0; k < left.Cols(); k++)
M[i][j] += left[i][k] * right[k][j];
}
return M;
}
Matrix MatrixSubset(const Matrix &input,
int startRow, int numRows, int startCol, int numCols)
unsigned startRow, unsigned numRows,
unsigned startCol, unsigned numCols)
{
Matrix M(numRows, numCols);
int i, j;
for(i=0; i<numRows; i++)
for(j=0; j<numCols; j++)
for(unsigned i = 0; i < numRows; i++)
for(unsigned j = 0; j < numCols; j++)
M[i][j] = input[startRow+i][startCol+j];
return M;
}
@ -313,11 +268,10 @@ Matrix MatrixConcatenateCols(const Matrix& left, const Matrix& right)
{
wxASSERT(left.Rows() == right.Rows());
Matrix M(left.Rows(), left.Cols() + right.Cols());
int i, j;
for(i=0; i<left.Rows(); i++) {
for(j=0; j<left.Cols(); j++)
for(unsigned i = 0; i < left.Rows(); i++) {
for(unsigned j = 0; j < left.Cols(); j++)
M[i][j] = left[i][j];
for(j=0; j<right.Cols(); j++)
for(unsigned j = 0; j < right.Cols(); j++)
M[i][j+left.Cols()] = right[i][j];
}
return M;
@ -326,9 +280,8 @@ Matrix MatrixConcatenateCols(const Matrix& left, const Matrix& right)
Matrix TransposeMatrix(const Matrix& other)
{
Matrix M(other.Cols(), other.Rows());
int i, j;
for(i=0; i<other.Rows(); i++)
for(j=0; j<other.Cols(); j++)
for(unsigned i = 0; i < other.Rows(); i++)
for(unsigned j = 0; j < other.Cols(); j++)
M[j][i] = other[i][j];
return M;
}
@ -340,20 +293,19 @@ bool InvertMatrix(const Matrix& input, Matrix& Minv)
// Returns true if successful
wxASSERT(input.Rows() == input.Cols());
int N = input.Rows();
int i, j, k;
auto N = input.Rows();
Matrix M = input;
Minv = IdentityMatrix(N);
// Do the elimination one column at a time
for(i=0; i<N; i++) {
for(unsigned i = 0; i < N; i++) {
// Pivot the row with the largest absolute value in
// column i, into row i
double absmax = 0.0;
int argmax=0;
int argmax = 0;
for(j=i; j<N; j++)
for(unsigned j = i; j < N; j++)
if (fabs(M[j][i]) > absmax) {
absmax = fabs(M[j][i]);
argmax = j;
@ -375,13 +327,13 @@ bool InvertMatrix(const Matrix& input, Matrix& Minv)
Minv[i] = Minv[i] * factor;
// Eliminate the rest of the column
for(j=0; j<N; j++) {
if (j==i)
for(unsigned j = 0; j < N; j++) {
if (j == i)
continue;
if (fabs(M[j][i]) > 0) {
// Subtract a multiple of row i from row j
double factor = M[j][i];
for(k=0; k<N; k++) {
for(unsigned k = 0; k < N; k++) {
M[j][k] -= (M[i][k] * factor);
Minv[j][k] -= (Minv[i][k] * factor);
}

@ -26,6 +26,8 @@
#ifndef __AUDACITY_MATRIX__
#define __AUDACITY_MATRIX__
#include "SampleFormat.h"
class Matrix;
class Vector
@ -33,63 +35,61 @@ class Vector
public:
Vector();
Vector(const Vector& copyFrom);
Vector(int len, double *data=NULL, bool copy=true);
Vector(int len, float *data);
Vector(unsigned len, double *data=NULL);
Vector(unsigned len, float *data);
Vector& operator=(const Vector &other);
virtual ~Vector();
~Vector();
inline double& operator[](int i) { return mData[i]; }
inline double operator[](int i) const { return mData[i]; }
inline int Len() const { return mN; }
void Reinit(unsigned len);
void Swap(Vector &that);
inline double& operator[](unsigned i) { return mData[i]; }
inline double operator[](unsigned i) const { return mData[i]; }
inline unsigned Len() const { return mN; }
double Sum() const;
private:
void CopyFrom(const Vector &other);
int mN;
double *mData;
bool mCopy;
unsigned mN{ 0 };
Doubles mData;
};
class Matrix
{
public:
Matrix(const Matrix& copyFrom);
Matrix(int rows, int cols, double **data=NULL);
virtual ~Matrix();
Matrix(unsigned rows, unsigned cols, double **data=NULL);
~Matrix();
Matrix& operator=(const Matrix& other);
inline Vector& operator[](int i) { return *mRowVec[i]; }
inline Vector& operator[](int i) const { return *mRowVec[i]; }
inline int Rows() const { return mRows; }
inline int Cols() const { return mCols; }
inline Vector& operator[](unsigned i) { return mRowVec[i]; }
inline Vector& operator[](unsigned i) const { return mRowVec[i]; }
inline unsigned Rows() const { return mRows; }
inline unsigned Cols() const { return mCols; }
void SwapRows(int i, int j);
double Sum() const;
void SwapRows(unsigned i, unsigned j);
private:
void CopyFrom(const Matrix& other);
int mRows;
int mCols;
Vector **mRowVec;
unsigned mRows;
unsigned mCols;
ArrayOf<Vector> mRowVec;
};
bool InvertMatrix(const Matrix& input, Matrix& Minv);
Matrix TransposeMatrix(const Matrix& M);
Matrix IdentityMatrix(int N);
Matrix IdentityMatrix(unsigned N);
Vector operator+(const Vector &left, const Vector &right);
Vector operator-(const Vector &left, const Vector &right);
Vector operator*(const Vector &left, const Vector &right);
Vector operator*(const Vector &left, double right);
Vector VectorSubset(const Vector &other, int start, int len);
Vector VectorSubset(const Vector &other, unsigned start, unsigned len);
Vector VectorConcatenate(const Vector& left, const Vector& right);
Vector operator*(const Vector &left, const Matrix &right);
@ -103,7 +103,8 @@ Matrix ScalarMultiply(const Matrix &left, const Matrix &right);
Matrix MatrixMultiply(const Matrix &left, const Matrix &right);
Matrix MatrixSubset(const Matrix &M,
int startRow, int numRows, int startCol, int numCols);
unsigned startRow, unsigned numRows,
unsigned startCol, unsigned numCols);
Matrix MatrixConcatenateCols(const Matrix& left, const Matrix& right);