mmat_07.cc

这是一个从音频信号里提取特征参量的程序

// file: $isip/class/math/matrix/MMatrix/mmat_07.cc
// version: $Id: mmat_07.cc,v 1.42 2002/03/12 22:24:30 zheng Exp $
//

// isip include files
//
#include "MMatrixMethods.h"
#include "MMatrix.h"
#include <typeinfo>

// method: inverse
//
// arguments:
//  MMatrix<TScalar, TIntegral>& this: (output) class operand
//  const MMatrix<TScalar, TIntegral>& m: (input) matrix operand
//
// return: a boolean value indicating status
//
// this method gets the inverse of the input matrix and assign it to
// the current matrix
//
template<class TScalar, class TIntegral>
boolean MMatrixMethods::inverse(MMatrix<TScalar, TIntegral>& this_a, 
				const MMatrix<TScalar, TIntegral>& m_a)  {

  // this variant uses double internally, it is defined only for float
  // and double matrices
  //
#ifdef ISIP_TEMPLATE_fp
  
  // declare local variables
  //
  long nrows = m_a.nrows_d;
  long ncols = m_a.ncols_d;
  
  // condition: the two input matrices are same
  //  if the current matrix is same as input matrix, assign the current
  //  matrix to a temporary matrix and call inverse
  //
  if (&this_a == &m_a) {
    MMatrix<TScalar, TIntegral> tmp(m_a);
    return this_a.inverse(tmp);
  }
  
  // condition: non-square matrix
  //
  if (!m_a.isSquare()) {
    m_a.debug(L"m_a");
    this_a.debug(L"this_a");    
    return Error::handle(name(), L"inverse", Error::ARG, __FILE__, __LINE__);
  }
  
  // get the type of the input matrix
  //
  long type = m_a.getType();

  // get the type of the output matrix
  //
  Integral::MTYPE out_type = this_a.getType();

  // create space in the current matrix
  //
  this_a.setDimensions(m_a, false, m_a.type_d);
    
  // if the determinant is zero then error as the matrix is a
  // singular matrix
  //  
  if (m_a.isSingular()) {
    m_a.debug(L"m_a");    
    this_a.debug(L"this_a");    
    return Error::handle(name(), L"inverse", 
			 MMatrix<TScalar, TIntegral>::ERR_SINGLR, 
			 __FILE__, __LINE__, Error::WARNING);
  }
  
  // type: DIAGONAL
  //
  if (type == Integral::DIAGONAL) {
    
    // compute the number of rows of the input matrix
    //
    long nrows = m_a.getNumRows();

    // create space in the current matrix
    //
    this_a.setDimensions(m_a, false, m_a.type_d);
    
    // compute inverse - this can be computed simply by taking the
    // inverse of the elements of the input matrix. these elements are
    // not zero, otherwise the matrix would be singular.
    //
    for (long row = 0; row < nrows; row++) {
      this_a.m_d(row) = 1 / (m_a.m_d(row));
    }
  }
  
  // type: FULL
  //
  else if (type == Integral::FULL) {
    
    // create a temporary double matrix to hold the input matrix
    //
    MMatrix<Double, double> a;
    a.assign(m_a);
    
    // initially store the identity matrix in current matrix. after the last
    // pass this_a will be replaced by the inverse of a.
    //
    this_a.makeIdentity(nrows, Integral::FULL);
    
    // initialize the value of the determinant
    //
    double det = 1;
    
    // start computing inverse matrix - the current pivot row is
    // ipass, for each pass, first find the maximum element in the
    // pivot column
    //
    for (long ipass = 0; ipass < nrows; ipass++) {
      long imx = ipass;
      for (long irow = ipass; irow < nrows; irow++) {
	if (fabs(a.getValue(irow, ipass)) > fabs(a.getValue(imx, ipass))) {
	  imx = irow;
	}
      }
      
      // interchange the elements of row ipass and row imx in both a and this_a
      //
      if (imx != ipass) {
	for (long icol = 0; icol < nrows; icol++) {
	  
	  double temp = this_a.getValue(ipass, icol);
	  this_a.setValue(ipass, icol, this_a.getValue(imx, icol));
	  this_a.setValue(imx, icol, temp);
	  if (icol >= ipass) {
	    temp = a.getValue(ipass, icol);
	    a.setValue(ipass, icol, a.getValue(imx, icol));
	    a.setValue(imx, icol, temp);
	  }
	}
      }
      
      // the current pivot is now a.getValue(ipass, ipass) - the determinant
      // is the product of the pivot elements
      //
      double pivot = a.getValue(ipass, ipass);
      
      det = det * pivot;
      if (det == 0) {
	this_a.debug(L"this_a");	
	return Error::handle(name(), L"inverse", 
			     MMatrix<TScalar, TIntegral>::ERR_SINGLR, 
			     __FILE__, __LINE__, Error::WARNING);
      }

      // normalize the pivot row by dividing across by the pivot element
      //
      for (long icol = 0; icol < nrows; icol++) { 
	this_a.setValue(ipass, icol, this_a.getValue(ipass, icol) / pivot);
	if (icol >= ipass) {
	  a.setValue(ipass, icol, a.getValue(ipass, icol) / pivot);
	}
      }
      
      // now replace each row by the row plus a multiple of the
      // pivot row with a factor chosen so that the element of a on
      // the pivot column is 0
      //
      double factor = 0;
      for (long irow = 0; irow < nrows; irow++) {
	
	if (irow != ipass) {
	  factor = a.getValue(irow, ipass);
	}
	
	for (long icol = 0; icol < nrows; icol++) {
	  if (irow != ipass) {
	    this_a.setValue(irow, icol, this_a.getValue(irow, icol) -
			  factor * this_a.getValue(ipass, icol));
	    a.setValue(irow, icol, a.getValue(irow, icol) -
		       factor * a.getValue(ipass, icol));
	  }
	}
      }
    }
  }
  
  // type: SYMMETRIC, LOWER_TRIANGULAR, UPPER_TRIANGULAR
  //  for symmetric and triangular matrix, change them into a full matrix
  //  and get the inverse
  //
  else if ((type == Integral::SYMMETRIC) ||
	   (type == Integral::LOWER_TRIANGULAR) ||
	   (type == Integral::UPPER_TRIANGULAR)) {
    this_a.assign(m_a, Integral::FULL);
    this_a.inverse();
  }
  
  // type: SPARSE
  //
  else {
    
    // declare temporary sparse matrix to hold the minor, and use this_a
    // to perform cofactor matrix
    //
    MMatrix<TScalar, TIntegral> mat_minor(nrows - 1, ncols - 1,
					  Integral::SPARSE);
    this_a.clear(Integral::RETAIN);
    
    // declare local variable
    //
    long length = 0;
    
    // compute determinant of input matrix
    //
    double det_in = m_a.determinant();
      
    // compute matrix of co-factor for input matrix
    //
    // loop over all the elements of matrix
    //
    for (long i = 0; i < nrows; i++) {
      for (long j = 0; j < ncols; j++) {
	
	// call minor of the element at this (i, j) position
	//
	m_a.getMinor(mat_minor, i, j);
	
	// calculate the determinant of minor matrix
	//
	double det = mat_minor.determinant();

	// check if the matrix is non-singular, if it is then compute inverse
	//
	if (det != 0) {
	  
	  // increasing the number of elements
	  //
	  length++;
	  
	  // checking sign of the elements and assigning it to
	  // their determinant to get the cofactor matrix
	  //
	  if ((i + j) % 2 == 0) {
	    
	    // assign positive sign if sum of row index and column
	    // index is even
	    // 
	    this_a.setValue(i, j, det);
	  }
	  
	  else {
	    
	    // assign negative sign if sum of row index and column
	    // index is odd
	    // 
	    this_a.setValue(i, j, -det);
	  }
	}
      }
    }
    
    // gives adjoint matrix, which is, transpose of co-factor matrix
    //
    this_a.transpose();
    this_a.div(det_in);
  }

  // solve the round-off problem during inverse computation
  //  1. for symmtric, copy the lower triangular part to upper triangular
  //     part, making it completely symmetric
  //  2. for lower triangular, set zeros in complete upper part, in case
  //     some non-zero but very close to zero elements exist;
  //  3. for upper triangular, set zeros in complete lower part, in case
  //     some non-zero but very close to zero elements exist;
  //
  if (out_type == Integral::SYMMETRIC) {
    for (long row = 0; row < nrows; row++) {
      for (long col = row + 1; col < ncols; col++) {
	this_a.setValue(row, col, this_a(col, row));
      }
    }
  }
  else if (out_type == Integral::LOWER_TRIANGULAR) {
    for (long row = 0; row < nrows; row++) {
      for (long col = row + 1; col < ncols; col++) {
	this_a.setValue(row, col, 0);
      }
    }
  }
  else if (out_type == Integral::UPPER_TRIANGULAR) {
    for (long row = 0; row < nrows; row++) {
      for (long col = 0; col < row; col++) {
	this_a.setValue(row, col, 0);
      }
    }
  }
  
  // change back to the oringinal type
  //
  if (!this_a.changeType(out_type)) {
    this_a.debug(L"this_a");    
    return Error::handle(name(), L"inverse",
			 MMatrix<TScalar, TIntegral>::ERR_OPTYPE,
			 __FILE__, __LINE__);
  }
  
  // exit gracefully
  //
  return true;

#endif

  // this variant uses complexdouble internally, it is defined only
  // for complexfloat and complexdouble matrices
  //
#ifdef ISIP_TEMPLATE_complex
#ifndef ISIP_TEMPLATE_ComplexLong_complexlong

  // declare local variables
  //
  long nrows = m_a.nrows_d;
  long ncols = m_a.ncols_d;
  
  // condition: the two input matrices are same
  //  if the current matrix is same as input matrix, assign the current
  //  matrix to a temporary matrix and call inverse
  //
  if (&this_a == &m_a) {
    MMatrix<TScalar, TIntegral> tmp(m_a);
    return this_a.inverse(tmp);
  }
  
  // condition: non-square matrix
  //
  if (!m_a.isSquare()) {
    m_a.debug(L"m_a");    
    this_a.debug(L"this_a");    
    return Error::handle(name(), L"inverse", Error::ARG, __FILE__, __LINE__);
  }
  
  // get the type of the input matrix
  //
  long type = m_a.getType();

  // get the type of the output matrix
  //
  Integral::MTYPE out_type = this_a.getType();

  // create space in the current matrix
  //
  this_a.setDimensions(m_a, false, m_a.type_d);
    
  // if the determinant is zero then error as the matrix is a
  // singular matrix
  //  
  if (m_a.isSingular()) {
    m_a.debug(L"m_a");    
    this_a.debug(L"this_a");    
    return Error::handle(name(), L"inverse", 
			 MMatrix<TScalar, TIntegral>::ERR_SINGLR, 
			 __FILE__, __LINE__, Error::WARNING);
  }
  
  // type: DIAGONAL
  //
  if (type == Integral::DIAGONAL) {
    
    // compute the number of rows of the input matrix
    //
    long nrows = m_a.getNumRows();

    // create space in the current matrix
    //
    this_a.setDimensions(m_a, false, m_a.type_d);
    
    // compute inverse - this can be computed simply by taking the
    // inverse of the elements of the input matrix. these elements are
    // not zero, otherwise the matrix would be singular.
    //
    for (long row = 0; row < nrows; row++) {
      this_a.m_d(row) = (TIntegral)1 / (m_a.m_d(row));
    }
  }
  
  // type: FULL
  //
  else if (type == Integral::FULL) {
    
    // create a temporary double matrix to hold the input matrix
    //
    MMatrix<ComplexDouble, complexdouble> a;
    a.assign(m_a);
    
    // initially store the identity matrix in current matrix. after the last
    // pass this_a will be replaced by the inverse of a.
    //
    this_a.makeIdentity(nrows, Integral::FULL);
    
    // initialize the value of the determinant
    //
    complexdouble det = 1;
    
    // start computing inverse matrix - the current pivot row is
    // ipass, for each pass, first find the maximum element in the
    // pivot column
    //
    for (long ipass = 0; ipass < nrows; ipass++) {
      long imx = ipass;
      for (long irow = ipass; irow < nrows; irow++) {
	
	ComplexDouble abs_irow, abs_imx;
	ComplexDouble temp_irow, temp_imx;	

	temp_irow = a.getValue(irow, ipass);
	temp_imx = a.getValue(imx, ipass);
	abs_irow = temp_irow.abs();
	abs_imx = temp_imx.abs();
	
	if (abs_irow > abs_imx) {
	  imx = irow;
	}
      }
      
      // interchange the elements of row ipass and row imx in both a and this_a
      //
      if (imx != ipass) {
	for (long icol = 0; icol < nrows; icol++) {
	  
	  complexdouble temp = this_a.getValue(ipass, icol);
	  this_a.setValue(ipass, icol, this_a.getValue(imx, icol));
	  this_a.setValue(imx, icol, (TIntegral)temp);
	  if (icol >= ipass) {
	    temp = a.getValue(ipass, icol);
	    a.setValue(ipass, icol, a.getValue(imx, icol));
	    a.setValue(imx, icol, temp);
	  }
	}
      }
      
      // the current pivot is now a.getValue(ipass, ipass) - the determinant
      // is the product of the pivot elements
      //
      complexdouble pivot = a.getValue(ipass, ipass);
      
      det = det * pivot;
      if (det == 0) {
	this_a.debug(L"this_a");	
	return Error::handle(name(), L"inverse", 
			     MMatrix<TScalar, TIntegral>::ERR_SINGLR, 
			     __FILE__, __LINE__, Error::WARNING);
      }

      // normalize the pivot row by dividing across by the pivot element
      //
      for (long icol = 0; icol < nrows; icol++) { 
	this_a.setValue(ipass, icol, this_a.getValue(ipass, icol) / pivot);
	if (icol >= ipass) {