Audacity  2.2.2
Matrix.cpp
Go to the documentation of this file.
1 /**********************************************************************
2 
3  Audacity: A Digital Audio Editor
4 
5  Matrix.cpp
6 
7  Dominic Mazzoni
8 
9 **********************************************************************/
10 
11 #include <stdlib.h>
12 #include <math.h>
13 
14 #include <wx/defs.h>
15 
16 #include "Matrix.h"
17 
19 {
20 }
21 
22 Vector::Vector(unsigned len, double *data)
23  : mN{ len }
24  , mData(len)
25 {
26  if (data)
27  std::copy(data, data + len, mData.get());
28  else
29  std::fill(mData.get(), mData.get() + len, 0.0);
30 }
31 
32 Vector::Vector(unsigned len, float *data)
33  : mN{ len }
34  , mData{ len }
35 {
36  if (data)
37  std::copy(data, data + len, mData.get());
38  else
39  std::fill(mData.get(), mData.get() + len, 0.0);
40 }
41 
43 {
44  wxASSERT(Len() == other.Len());
45  std::copy(other.mData.get(), other.mData.get() + mN, mData.get());
46  return *this;
47 }
48 
49 Vector::Vector(const Vector &other)
50  : mN{ other.Len() }
51  , mData{ mN }
52 {
53  std::copy(other.mData.get(), other.mData.get() + mN, mData.get());
54 }
55 
57 {
58 }
59 
60 void Vector::Reinit(unsigned len)
61 {
62  Vector temp(len);
63  Swap(temp);
64 }
65 
66 void Vector::Swap(Vector &that)
67 {
68  std::swap(mN, that.mN);
69  mData.swap(that.mData);
70 }
71 
72 double Vector::Sum() const
73 {
74  double sum = 0.0;
75  for(unsigned i = 0; i < Len(); i++)
76  sum += mData[i];
77  return sum;
78 }
79 
80 Matrix::Matrix(unsigned rows, unsigned cols, double **data)
81  : mRows{ rows }
82  , mCols{ cols }
83  , mRowVec{ mRows }
84 {
85  for(unsigned i = 0; i < mRows; i++) {
86  mRowVec[i].Reinit( mCols );
87  for(unsigned j = 0; j < mCols; j++) {
88  if (data)
89  (*this)[i][j] = data[i][j];
90  else
91  (*this)[i][j] = 0.0;
92  }
93  }
94 }
95 
97 {
98  CopyFrom(other);
99  return *this;
100 }
101 
102 Matrix::Matrix(const Matrix &other)
103 {
104  CopyFrom(other);
105 }
106 
107 void Matrix::CopyFrom(const Matrix &other)
108 {
109  mRows = other.mRows;
110  mCols = other.mCols;
111  mRowVec.reinit(mRows);
112  for (unsigned i = 0; i < mRows; i++) {
113  mRowVec[i].Reinit( mCols );
114  mRowVec[i] = other.mRowVec[i];
115  }
116 }
117 
119 {
120 }
121 
122 void Matrix::SwapRows(unsigned i, unsigned j)
123 {
124  mRowVec[i].Swap(mRowVec[j]);
125 }
126 
128 {
129  Matrix M(N, N);
130  for(unsigned i = 0; i < N; i++)
131  M[i][i] = 1.0;
132  return M;
133 }
134 
135 Vector operator+(const Vector &left, const Vector &right)
136 {
137  wxASSERT(left.Len() == right.Len());
138  Vector v(left.Len());
139  for(unsigned i = 0; i < left.Len(); i++)
140  v[i] = left[i] + right[i];
141  return v;
142 }
143 
144 Vector operator-(const Vector &left, const Vector &right)
145 {
146  wxASSERT(left.Len() == right.Len());
147  Vector v(left.Len());
148  for(unsigned i = 0; i < left.Len(); i++)
149  v[i] = left[i] - right[i];
150  return v;
151 }
152 
153 Vector operator*(const Vector &left, const Vector &right)
154 {
155  wxASSERT(left.Len() == right.Len());
156  Vector v(left.Len());
157  for(unsigned i = 0; i < left.Len(); i++)
158  v[i] = left[i] * right[i];
159  return v;
160 }
161 
162 Vector operator*(const Vector &left, double right)
163 {
164  Vector v(left.Len());
165  for(unsigned i = 0; i < left.Len(); i++)
166  v[i] = left[i] * right;
167  return v;
168 }
169 
170 Vector VectorSubset(const Vector &other, unsigned start, unsigned len)
171 {
172  Vector v(len);
173  for(unsigned i = 0; i < len; i++)
174  v[i] = other[start+i];
175  return v;
176 }
177 
178 Vector VectorConcatenate(const Vector& left, const Vector& right)
179 {
180  Vector v(left.Len() + right.Len());
181  for(unsigned i = 0; i < left.Len(); i++)
182  v[i] = left[i];
183  for(unsigned i = 0; i < right.Len(); i++)
184  v[i + left.Len()] = right[i];
185  return v;
186 }
187 
188 Vector operator*(const Vector &left, const Matrix &right)
189 {
190  wxASSERT(left.Len() == right.Rows());
191  Vector v(right.Cols());
192  for(unsigned i = 0; i < right.Cols(); i++) {
193  v[i] = 0.0;
194  for(unsigned j = 0; j < right.Rows(); j++)
195  v[i] += left[j] * right[j][i];
196  }
197  return v;
198 }
199 
200 Vector operator*(const Matrix &left, const Vector &right)
201 {
202  wxASSERT(left.Cols() == right.Len());
203  Vector v(left.Rows());
204  for(unsigned i = 0; i < left.Rows(); i++) {
205  v[i] = 0.0;
206  for(unsigned j = 0; j < left.Cols(); j++)
207  v[i] += left[i][j] * right[j];
208  }
209  return v;
210 }
211 
212 Matrix operator+(const Matrix &left, const Matrix &right)
213 {
214  wxASSERT(left.Rows() == right.Rows());
215  wxASSERT(left.Cols() == right.Cols());
216  Matrix M(left.Rows(), left.Cols());
217  for(unsigned i = 0; i < left.Rows(); i++)
218  for(unsigned j = 0; j < left.Cols(); j++)
219  M[i][j] = left[i][j] + right[i][j];
220  return M;
221 }
222 
223 Matrix operator*(const Matrix &left, const double right)
224 {
225  Matrix M(left.Rows(), left.Cols());
226  for(unsigned i = 0; i < left.Rows(); i++)
227  for(unsigned j = 0; j < left.Cols(); j++)
228  M[i][j] = left[i][j] * right;
229  return M;
230 }
231 
232 Matrix ScalarMultiply(const Matrix &left, const Matrix &right)
233 {
234  wxASSERT(left.Rows() == right.Rows());
235  wxASSERT(left.Cols() == right.Cols());
236  Matrix M(left.Rows(), left.Cols());
237  for(unsigned i = 0; i < left.Rows(); i++)
238  for(unsigned j = 0; j < left.Cols(); j++)
239  M[i][j] = left[i][j] * right[i][j];
240  return M;
241 }
242 
243 Matrix MatrixMultiply(const Matrix &left, const Matrix &right)
244 {
245  wxASSERT(left.Cols() == right.Rows());
246  Matrix M(left.Rows(), right.Cols());
247  for(unsigned i = 0; i < left.Rows(); i++)
248  for(unsigned j = 0; j < right.Cols(); j++) {
249  M[i][j] = 0.0;
250  for(unsigned k = 0; k < left.Cols(); k++)
251  M[i][j] += left[i][k] * right[k][j];
252  }
253  return M;
254 }
255 
257  unsigned startRow, unsigned numRows,
258  unsigned startCol, unsigned numCols)
259 {
260  Matrix M(numRows, numCols);
261  for(unsigned i = 0; i < numRows; i++)
262  for(unsigned j = 0; j < numCols; j++)
263  M[i][j] = input[startRow+i][startCol+j];
264  return M;
265 }
266 
267 Matrix MatrixConcatenateCols(const Matrix& left, const Matrix& right)
268 {
269  wxASSERT(left.Rows() == right.Rows());
270  Matrix M(left.Rows(), left.Cols() + right.Cols());
271  for(unsigned i = 0; i < left.Rows(); i++) {
272  for(unsigned j = 0; j < left.Cols(); j++)
273  M[i][j] = left[i][j];
274  for(unsigned j = 0; j < right.Cols(); j++)
275  M[i][j+left.Cols()] = right[i][j];
276  }
277  return M;
278 }
279 
281 {
282  Matrix M(other.Cols(), other.Rows());
283  for(unsigned i = 0; i < other.Rows(); i++)
284  for(unsigned j = 0; j < other.Cols(); j++)
285  M[j][i] = other[i][j];
286  return M;
287 }
288 
289 bool InvertMatrix(const Matrix& input, Matrix& Minv)
290 {
291  // Very straightforward implementation of
292  // Gauss-Jordan elimination to invert a matrix.
293  // Returns true if successful
294 
295  wxASSERT(input.Rows() == input.Cols());
296  auto N = input.Rows();
297 
298  Matrix M = input;
299  Minv = IdentityMatrix(N);
300 
301  // Do the elimination one column at a time
302  for(unsigned i = 0; i < N; i++) {
303  // Pivot the row with the largest absolute value in
304  // column i, into row i
305  double absmax = 0.0;
306  unsigned int argmax = 0;
307 
308  for(unsigned j = i; j < N; j++)
309  if (fabs(M[j][i]) > absmax) {
310  absmax = fabs(M[j][i]);
311  argmax = j;
312  }
313 
314  // If no row has a nonzero value in that column,
315  // the matrix is singular and we have to give up.
316  if (absmax == 0)
317  return false;
318 
319  if (i != argmax) {
320  M.SwapRows(i, argmax);
321  Minv.SwapRows(i, argmax);
322  }
323 
324  // Divide this row by the value of M[i][i]
325  double factor = 1.0 / M[i][i];
326  M[i] = M[i] * factor;
327  Minv[i] = Minv[i] * factor;
328 
329  // Eliminate the rest of the column
330  for(unsigned j = 0; j < N; j++) {
331  if (j == i)
332  continue;
333  if (fabs(M[j][i]) > 0) {
334  // Subtract a multiple of row i from row j
335  double factor = M[j][i];
336  for(unsigned k = 0; k < N; k++) {
337  M[j][k] -= (M[i][k] * factor);
338  Minv[j][k] -= (Minv[i][k] * factor);
339  }
340  }
341  }
342  }
343 
344  return true;
345 }
unsigned mN
Definition: Matrix.h:53
Matrix & operator=(const Matrix &other)
Definition: Matrix.cpp:96
void SwapRows(unsigned i, unsigned j)
Definition: Matrix.cpp:122
unsigned Len() const
Definition: Matrix.h:48
void reinit(Integral count, bool initialize=false)
Definition: MemoryX.h:117
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
Vector operator*(const Vector &left, const Vector &right)
Definition: Matrix.cpp:153
unsigned mRows
Definition: Matrix.h:76
Vector operator-(const Vector &left, const Vector &right)
Definition: Matrix.cpp:144
Matrix ScalarMultiply(const Matrix &left, const Matrix &right)
Definition: Matrix.cpp:232
General routine to interpolate (or even extrapolate small amounts) audio when a few of the samples ar...
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
Vector operator+(const Vector &left, const Vector &right)
Definition: Matrix.cpp:135
Vector()
Definition: Matrix.cpp:18
Matrix IdentityMatrix(unsigned N)
Definition: Matrix.cpp:127
double Sum() const
Definition: Matrix.cpp:72
Vector VectorConcatenate(const Vector &left, const Vector &right)
Definition: Matrix.cpp:178
~Matrix()
Definition: Matrix.cpp:118
void Swap(Vector &that)
Definition: Matrix.cpp:66
Holds a matrix of doubles and supports arithmetic, subsetting, and matrix inversion. Used by InterpolateAudio.
Definition: Matrix.h:57
void Reinit(unsigned len)
Definition: Matrix.cpp:60
Matrix(const Matrix &copyFrom)
Definition: Matrix.cpp:102
Doubles mData
Definition: Matrix.h:54
~Vector()
Definition: Matrix.cpp:56
ArrayOf< Vector > mRowVec
Definition: Matrix.h:78
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
Vector & operator=(const Vector &other)
Definition: Matrix.cpp:42
unsigned mCols
Definition: Matrix.h:77
void CopyFrom(const Matrix &other)
Definition: Matrix.cpp:107