407 lines
8.2 KiB
C++
407 lines
8.2 KiB
C++
/**********************************************************************
|
|
|
|
Audacity: A Digital Audio Editor
|
|
|
|
Matrix.cpp
|
|
|
|
Dominic Mazzoni
|
|
|
|
**********************************************************************/
|
|
|
|
#include <stdlib.h>
|
|
#include <math.h>
|
|
|
|
#include <wx/defs.h>
|
|
|
|
#include "Matrix.h"
|
|
|
|
Vector::Vector()
|
|
{
|
|
mCopy = false;
|
|
mN = 0;
|
|
mData = NULL;
|
|
}
|
|
|
|
Vector::Vector(int len, double *data, bool copy)
|
|
{
|
|
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;
|
|
}
|
|
}
|
|
|
|
Vector& Vector::operator=(const Vector &other)
|
|
{
|
|
wxASSERT(Len() == other.Len());
|
|
int i;
|
|
for(i=0; i<Len(); i++)
|
|
mData[i] = other.mData[i];
|
|
return *this;
|
|
}
|
|
|
|
Vector::Vector(const Vector &other)
|
|
{
|
|
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];
|
|
}
|
|
|
|
Vector::~Vector()
|
|
{
|
|
if (mCopy)
|
|
delete[] mData;
|
|
}
|
|
|
|
Vector::Vector(int len, float *data)
|
|
{
|
|
mCopy = true;
|
|
mN = len;
|
|
mData = new double[mN];
|
|
int i;
|
|
for(i=0; i<mN; i++)
|
|
mData[i] = (double)data[i];
|
|
}
|
|
|
|
double Vector::Sum() const
|
|
{
|
|
int i;
|
|
double sum = 0.0;
|
|
for(i=0; i<Len(); i++)
|
|
sum += mData[i];
|
|
return sum;
|
|
}
|
|
|
|
Matrix::Matrix(int rows, int cols, double **data)
|
|
{
|
|
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++) {
|
|
if (data)
|
|
(*this)[i][j] = data[i][j];
|
|
else
|
|
(*this)[i][j] = 0.0;
|
|
}
|
|
}
|
|
}
|
|
|
|
Matrix& Matrix::operator=(const Matrix &other)
|
|
{
|
|
CopyFrom(other);
|
|
return *this;
|
|
}
|
|
|
|
Matrix::Matrix(const Matrix &other)
|
|
{
|
|
CopyFrom(other);
|
|
}
|
|
|
|
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];
|
|
}
|
|
}
|
|
|
|
Matrix::~Matrix()
|
|
{
|
|
int i;
|
|
for(i=0; i<mRows; i++)
|
|
delete mRowVec[i];
|
|
delete[] mRowVec;
|
|
}
|
|
|
|
void Matrix::SwapRows(int i, int j)
|
|
{
|
|
Vector *tmp = mRowVec[i];
|
|
mRowVec[i] = mRowVec[j];
|
|
mRowVec[j] = tmp;
|
|
}
|
|
|
|
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 M(N, N);
|
|
int i;
|
|
for(i=0; i<N; i++)
|
|
M[i][i] = 1.0;
|
|
return M;
|
|
}
|
|
|
|
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++)
|
|
v[i] = left[i] + right[i];
|
|
return v;
|
|
}
|
|
|
|
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++)
|
|
v[i] = left[i] - right[i];
|
|
return v;
|
|
}
|
|
|
|
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++)
|
|
v[i] = left[i] * right[i];
|
|
return v;
|
|
}
|
|
|
|
Vector operator*(const Vector &left, double right)
|
|
{
|
|
Vector v(left.Len());
|
|
int i;
|
|
for(i=0; i<left.Len(); i++)
|
|
v[i] = left[i] * right;
|
|
return v;
|
|
}
|
|
|
|
Vector VectorSubset(const Vector &other, int start, int len)
|
|
{
|
|
Vector v(len);
|
|
int i;
|
|
for(i=0; i<len; i++)
|
|
v[i] = other[start+i];
|
|
return v;
|
|
}
|
|
|
|
Vector VectorConcatenate(const Vector& left, const Vector& right)
|
|
{
|
|
Vector v(left.Len() + right.Len());
|
|
int i;
|
|
for(i=0; i<left.Len(); i++)
|
|
v[i] = left[i];
|
|
for(i=0; i<right.Len(); i++)
|
|
v[i + left.Len()] = right[i];
|
|
return v;
|
|
}
|
|
|
|
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++) {
|
|
v[i] = 0.0;
|
|
for(j=0; j<right.Rows(); j++)
|
|
v[i] += left[j] * right[j][i];
|
|
}
|
|
return v;
|
|
}
|
|
|
|
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++) {
|
|
v[i] = 0.0;
|
|
for(j=0; j<left.Cols(); j++)
|
|
v[i] += left[i][j] * right[j];
|
|
}
|
|
return v;
|
|
}
|
|
|
|
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++)
|
|
M[i][j] = left[i][j] + right[i][j];
|
|
return M;
|
|
}
|
|
|
|
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++)
|
|
M[i][j] = left[i][j] * right;
|
|
return M;
|
|
}
|
|
|
|
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++)
|
|
M[i][j] = left[i][j] * right[i][j];
|
|
return M;
|
|
}
|
|
|
|
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++) {
|
|
M[i][j] = 0.0;
|
|
for(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)
|
|
{
|
|
Matrix M(numRows, numCols);
|
|
int i, j;
|
|
for(i=0; i<numRows; i++)
|
|
for(j=0; j<numCols; j++)
|
|
M[i][j] = input[startRow+i][startCol+j];
|
|
return M;
|
|
}
|
|
|
|
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++)
|
|
M[i][j] = left[i][j];
|
|
for(j=0; j<right.Cols(); j++)
|
|
M[i][j+left.Cols()] = right[i][j];
|
|
}
|
|
return M;
|
|
}
|
|
|
|
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++)
|
|
M[j][i] = other[i][j];
|
|
return M;
|
|
}
|
|
|
|
bool InvertMatrix(const Matrix& input, Matrix& Minv)
|
|
{
|
|
// Very straightforward implementation of
|
|
// Gauss-Jordan elimination to invert a matrix.
|
|
// Returns true if successful
|
|
|
|
wxASSERT(input.Rows() == input.Cols());
|
|
int N = input.Rows();
|
|
int i, j, k;
|
|
|
|
Matrix M = input;
|
|
Minv = IdentityMatrix(N);
|
|
|
|
// Do the elimination one column at a time
|
|
for(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;
|
|
|
|
for(j=i; j<N; j++)
|
|
if (fabs(M[j][i]) > absmax) {
|
|
absmax = fabs(M[j][i]);
|
|
argmax = j;
|
|
}
|
|
|
|
// If no row has a nonzero value in that column,
|
|
// the matrix is singular and we have to give up.
|
|
if (absmax == 0)
|
|
return false;
|
|
|
|
if (i != argmax) {
|
|
M.SwapRows(i, argmax);
|
|
Minv.SwapRows(i, argmax);
|
|
}
|
|
|
|
// Divide this row by the value of M[i][i]
|
|
double factor = 1.0 / M[i][i];
|
|
M[i] = M[i] * factor;
|
|
Minv[i] = Minv[i] * factor;
|
|
|
|
// Eliminate the rest of the column
|
|
for(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++) {
|
|
M[j][k] -= (M[i][k] * factor);
|
|
Minv[j][k] -= (Minv[i][k] * factor);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
return true;
|
|
}
|
|
|
|
// Indentation settings for Vim and Emacs.
|
|
// Please do not modify past this point.
|
|
//
|
|
// Local Variables:
|
|
// c-basic-offset: 3
|
|
// indent-tabs-mode: nil
|
|
// End:
|
|
//
|
|
// vim: et sts=3 sw=3
|
|
//
|
|
|
|
|