Audacity 3.2.0
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 "Matrix.h"
12
13#include <stdlib.h>
14#include <math.h>
15
16#include <wx/defs.h>
17
19{
20}
21
22Vector::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
32Vector::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
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
60void Vector::Reinit(unsigned len)
61{
62 Vector temp(len);
63 Swap(temp);
64}
65
67{
68 std::swap(mN, that.mN);
69 mData.swap(that.mData);
70}
71
72double 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
80Matrix::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
103{
104 CopyFrom(other);
105}
106
107void Matrix::CopyFrom(const Matrix &other)
108{
109 mRows = other.mRows;
110 mCols = other.mCols;
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
122void 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
135Vector 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
144Vector 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
153Vector 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
162Vector 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
170Vector 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
178Vector 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
188Vector 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
200Vector 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
212Matrix 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
223Matrix 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
232Matrix 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
243Matrix 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
267Matrix 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
289bool 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 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}
Matrix IdentityMatrix(unsigned N)
Definition: Matrix.cpp:127
Vector VectorConcatenate(const Vector &left, const Vector &right)
Definition: Matrix.cpp:178
Vector operator-(const Vector &left, const Vector &right)
Definition: Matrix.cpp:144
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
Vector operator+(const Vector &left, const Vector &right)
Definition: Matrix.cpp:135
bool InvertMatrix(const Matrix &input, Matrix &Minv)
Definition: Matrix.cpp:289
Matrix MatrixConcatenateCols(const Matrix &left, const Matrix &right)
Definition: Matrix.cpp:267
Vector operator*(const Vector &left, const Vector &right)
Definition: Matrix.cpp:153
Vector VectorSubset(const Vector &other, unsigned start, unsigned len)
Definition: Matrix.cpp:170
Matrix ScalarMultiply(const Matrix &left, const Matrix &right)
Definition: Matrix.cpp:232
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...
void reinit(Integral count, bool initialize=false)
Definition: MemoryX.h:57
Holds a matrix of doubles and supports arithmetic, subsetting, and matrix inversion....
Definition: Matrix.h:58
unsigned Cols() const
Definition: Matrix.h:69
Matrix(const Matrix &copyFrom)
Definition: Matrix.cpp:102
Matrix & operator=(const Matrix &other)
Definition: Matrix.cpp:96
unsigned Rows() const
Definition: Matrix.h:68
ArrayOf< Vector > mRowVec
Definition: Matrix.h:78
void SwapRows(unsigned i, unsigned j)
Definition: Matrix.cpp:122
~Matrix()
Definition: Matrix.cpp:118
unsigned mRows
Definition: Matrix.h:76
void CopyFrom(const Matrix &other)
Definition: Matrix.cpp:107
unsigned mCols
Definition: Matrix.h:77
Holds a matrix of doubles and supports arithmetic operations, including Vector-Matrix operations....
Definition: Matrix.h:34
void Reinit(unsigned len)
Definition: Matrix.cpp:60
~Vector()
Definition: Matrix.cpp:56
Vector & operator=(const Vector &other)
Definition: Matrix.cpp:42
Doubles mData
Definition: Matrix.h:54
unsigned Len() const
Definition: Matrix.h:48
void Swap(Vector &that)
Definition: Matrix.cpp:66
double Sum() const
Definition: Matrix.cpp:72
Vector()
Definition: Matrix.cpp:18
unsigned mN
Definition: Matrix.h:53
void swap(std::unique_ptr< Alg_seq > &a, std::unique_ptr< Alg_seq > &b)
Definition: NoteTrack.cpp:753