mmat_17.cc

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

// file: $isip/class/math/matrix/MMatrix/mmat_17.cc
// version: $Id: mmat_17.cc,v 1.17 2002/02/27 20:54:30 alphonso Exp $
//

// isip include files
//
#include "MMatrixMethods.h"
#include "MMatrix.h"

// method: getDiagonal
//
// arguments:
//  const MMatrix<TScalar, TIntegral>& this: (input) class operand
//  MVector<TScalar, TIntegral>& vector: (output) diagonal values
//
// return: a boolean value indicating status
//
template<class TScalar, class TIntegral>
boolean MMatrixMethods::getDiagonal(const MMatrix<TScalar, TIntegral>& this_a,
				    MVector<TScalar, TIntegral>& vector_a) {

  // the matrix must be square
  //
  if (!this_a.isSquare()) {
    this_a.debug(L"this_a");
    return Error::handle(name(), L"getDiagonal", Error::ARG,
			 __FILE__, __LINE__);
  }
    
  // type: DIAGONAL
  //  simply copy the output matrix
  //
  if (this_a.type_d == Integral::DIAGONAL) {
    return vector_a.assign(this_a.m_d);
  }

  // type: non-DIAGONAL
  //
  else {

    // create output space
    //
    vector_a.setLength(this_a.nrows_d, false);

    // loop over all elements
    //
    for (long row = 0; row < this_a.nrows_d; row++) {
      vector_a(row) = this_a.getValue(row, row);
    }
  }
  
  // exit gracefully
  //
  return true;
}

// explicit instantiations
//
template boolean
MMatrixMethods::getDiagonal<ISIP_TEMPLATE_TARGET>
(const MMatrix<ISIP_TEMPLATE_TARGET>&, MVector<ISIP_TEMPLATE_TARGET>&);

// method: getLower
//
// arguments:
//  const MMatrix<TScalar, TIntegral>& this: (input) class operand
//  MMatrix<TScalar, TIntegral>& arg: (output) lower triangular values
//
// return: a boolean value indicating status
//
// note that this method clears the output matrix
//
template<class TScalar, class TIntegral>
boolean MMatrixMethods::getLower(const MMatrix<TScalar, TIntegral>& this_a,
				 MMatrix<TScalar, TIntegral>& arg_a) {

  // check for compatibility
  //
  if (((arg_a.type_d == Integral::DIAGONAL) && (!this_a.isDiagonal())) ||
      ((arg_a.type_d == Integral::UPPER_TRIANGULAR) &&
       (!this_a.isDiagonal())) ||
      ((arg_a.type_d == Integral::SYMMETRIC) && (!this_a.isDiagonal()))) {
    arg_a.debug(L"arg_a");    
    this_a.debug(L"this_a");    
    return Error::handle(name(), L"getLower", Error::ARG, __FILE__, __LINE__);
  }

  // check for both types being LOWER_TRIANGULAR (in which case this
  // method shouldn't even be called)
  //
  else if ((this_a.type_d == Integral::LOWER_TRIANGULAR) &&
	   (arg_a.type_d == Integral::LOWER_TRIANGULAR)) {
    arg_a.nrows_d = this_a.nrows_d;
    arg_a.ncols_d = this_a.ncols_d;
    return arg_a.m_d.assign(this_a.m_d);
  }

  // all other cases: loop over all elements and manually copy data
  //
  arg_a.setDimensions(this_a.nrows_d, this_a.ncols_d, false);
  arg_a.clear(Integral::RETAIN);

  for (long row = 0; row < this_a.nrows_d; row++) {
    for (long col = 0; col <= row; col++) {
      arg_a.setValue(row, col, this_a.getValue(row, col));
    }
  }
  
  // exit gracefully
  //
  return true;
}

// explicit instantiations
//
template boolean
MMatrixMethods::getLower<ISIP_TEMPLATE_TARGET>
(const MMatrix<ISIP_TEMPLATE_TARGET>&, MMatrix<ISIP_TEMPLATE_TARGET>&);

// method: getUpper
//
// arguments:
//  const MMatrix<TScalar, TIntegral>& this: (input) class operand
//  MMatrix<TScalar, TIntegral>& arg: (output) upper triangular values
//
// return: a boolean value indicating status
//
// note that this method clears the output matrix
//
template<class TScalar, class TIntegral>
boolean MMatrixMethods::getUpper(const MMatrix<TScalar, TIntegral>& this_a,
				 MMatrix<TScalar, TIntegral>& arg_a) {

  // check for compatibility:
  //  there are some constraints on when this method can be invoked.
  //  the destination matrix must be square,
  //  and not be the wrong type (UPPER_TRIANGULAR).
  //
  if (((arg_a.type_d == Integral::DIAGONAL) && (!this_a.isDiagonal())) ||
      ((arg_a.type_d == Integral::LOWER_TRIANGULAR) &&
       (!this_a.isDiagonal())) ||
      ((arg_a.type_d == Integral::SYMMETRIC) && (!this_a.isDiagonal()))) {
    arg_a.debug(L"arg_a");    
    this_a.debug(L"this_a");    
    return Error::handle(name(), L"getUpper", Error::ARG, __FILE__, __LINE__);
  }

  // check for both types being UPPER_TRIANGULAR (in which case this
  // method shouldn't even be called)
  //
  else if ((this_a.type_d == Integral::UPPER_TRIANGULAR) &&
	   (arg_a.type_d == Integral::UPPER_TRIANGULAR)) {
    arg_a.nrows_d = this_a.nrows_d;
    arg_a.ncols_d = this_a.ncols_d;
    return arg_a.m_d.assign(this_a.m_d);
  }

  // all other cases: loop over all elements and manually copy data
  //
  arg_a.setDimensions(this_a.nrows_d, this_a.ncols_d, false);
  arg_a.clear(Integral::RETAIN);

  for (long row = 0; row < this_a.nrows_d; row++) {
    for (long col = row; col < this_a.ncols_d; col++) {
      arg_a.setValue(row, col, this_a.getValue(row, col));
    }
  }
  
  // exit gracefully
  //
  return true;
}

// explicit instantiations
//
template boolean
MMatrixMethods::getUpper<ISIP_TEMPLATE_TARGET>
(const MMatrix<ISIP_TEMPLATE_TARGET>&, MMatrix<ISIP_TEMPLATE_TARGET>&);

// method: getMinor
//
// arguments:
//  const MMatrix<TScalar, TIntegral>& this: (output) class operand
//  MMatrix<TScalar, TIntegral>& arg: (output) contains result
//  long row_index: (input) row for which minor is to be calculated
//  long col_index: (input) column for which minor is to be calculated
//
// return: a boolean value indicating status
//
// this method computes the minor matrix for a particular element of
// the input matrix. the output will be stored in the current matrix
//
template<class TScalar, class TIntegral>
boolean MMatrixMethods::getMinor(const MMatrix<TScalar, TIntegral>& this_a, 
				 MMatrix<TScalar, TIntegral>& arg_a, 
				 long row_index_a, long col_index_a) {

  // calculate the dimensions for the output
  //
  long nrows_new = (long)(this_a.nrows_d) - 1;
  long ncols_new = (long)(this_a.ncols_d) - 1;

  // check the arguments and dimensions
  //
  if ((nrows_new <= 0) || (ncols_new <= 0) ||
      (row_index_a < 0) || (col_index_a < 0) ||
      (row_index_a >= this_a.nrows_d) || (col_index_a >= this_a.ncols_d)) {
    this_a.debug(L"this_a");    
    return Error::handle(name(), L"getMinor", Error::ARG,
			 __FILE__, __LINE__);
  }
  else if ((nrows_new == 0) && (ncols_new == 0)) {
    return arg_a.clear(Integral::RESET);
  }

  // create the space in the output matrix
  //
  arg_a.setDimensions(nrows_new, ncols_new, false, Integral::UNCHANGED);

  // type: non-SPARSE input; non-SPARSE output
  //
  if ((this_a.type_d != Integral::SPARSE) &&
      (arg_a.type_d != Integral::SPARSE)) {
  
    // loop over the rows of the current matrix and copy elements:
    //  be sure to skip the excluded row
    //
    for (long row = 0, orow = 0; row < this_a.nrows_d; row++) {
      
      // exclude the specified row
      //
      if (row != row_index_a) {

	// loop over columns in this row
	//
	for (long col = 0, ocol = 0; col < this_a.ncols_d; col++) {
	  
	  // exclude the specified column
	  //
	  if (col != col_index_a) {
	    arg_a.setValue(orow, ocol++, this_a.getValue(row, col));
	  }
	}

	// increment the output row counter
	//
	orow++;
      }
    }
  }  
  
  // type: non-SPARSE input; SPARSE output
  //  
  else if ((this_a.type_d != Integral::SPARSE) &&
	   (arg_a.type_d == Integral::SPARSE)) {
  
    // compute the number of non-zero elements and make (extra) space
    //
    long num_elements = this_a.numNotEqual(0);
    arg_a.m_d.setCapacity(num_elements);
    arg_a.row_index_d.setCapacity(num_elements);
    arg_a.col_index_d.setCapacity(num_elements);

    // loop over the rows of the current matrix and copy elements:
    //  be sure to skip the excluded row
    //
    for (long row = 0, index = 0; row < this_a.nrows_d; row++) {
      
      // exclude the specified row
      //
      if (row != row_index_a) {

	// loop over columns in this row
	//
	for (long col = 0; col < this_a.ncols_d; col++) {
	  
	  // exclude the specified column
	  //
	  if (col != col_index_a) {

	    // get the current value
	    //
	    TIntegral tmp = this_a.getValue(row, col);

	    // add this value if it is not zero
	    //
	    long new_row = row;
	    long new_col = col;
	    
	    if (tmp != 0) {

	      // compute the output row and column index: we must subtract one
	      // if this element occurs after the row and column to be deleted
	      //
	      if (new_row > row_index_a) {
		new_row--;
	      }
	      if (new_col > col_index_a) {
		new_col--;
	      }

	      // create space
	      //
	      arg_a.m_d.setLength(index + 1);
	      arg_a.row_index_d.setLength(index + 1);
	      arg_a.col_index_d.setLength(index + 1);
	    
	      // add this to the matrix
	      //
	      arg_a.m_d(index) = tmp;
	      arg_a.row_index_d(index) = new_row;
	      arg_a.col_index_d(index) = new_col;

	      // increment the counter
	      //
	      index++;
	    }
	  }
	}
      }
    }
  }

  // type: SPARSE input; SPARSE output
  //  
  else {
  
    // compute the number of non-zero elements and make space
    //
    long num_elements = this_a.m_d.length();
    arg_a.m_d.setCapacity(num_elements);
    arg_a.row_index_d.setCapacity(num_elements);
    arg_a.col_index_d.setCapacity(num_elements);

    // loop over the rows of the current matrix and copy elements:
    //  be sure to skip the excluded row
    //
    for (long index = 0, oindex = 0; index < num_elements; index++) {
      
      // exclude the specified row
      //
      long new_row = this_a.row_index_d(index);
      long new_col = this_a.col_index_d(index);

      if ((new_row != row_index_a) &&
	  (new_col != col_index_a)) {

	// get the current value
	//
	TIntegral tmp = this_a.m_d(index);

	// add this value if it is not zero
	//
	if (tmp != 0) {

	  // compute the output row and column index: we must subtract one
	  // if this element occurs after the row and column to be deleted
	  //
	  if (new_row > row_index_a) {
	    new_row--;
	  }
	  if (new_col > col_index_a) {
	    new_col--;
	  }

	  // create space
	  //
	  arg_a.m_d.setLength(oindex + 1);
	  arg_a.row_index_d.setLength(oindex + 1);
	  arg_a.col_index_d.setLength(oindex + 1);
	    
	  // add this to the matrix
	  //
	  arg_a.m_d(oindex) = tmp;
	  arg_a.row_index_d(oindex) = new_row;
	  arg_a.col_index_d(oindex) = new_col;

	  // increment the counter
	  //
	  oindex++;
	}
      }
    }
  }

  // exit gracefully
  //
  return true;
}

// explicit instantiations
//
template boolean
MMatrixMethods::getMinor<ISIP_TEMPLATE_TARGET>
(const MMatrix<ISIP_TEMPLATE_TARGET>&, MMatrix<ISIP_TEMPLATE_TARGET>&, long, long);

// method: setDiagonal
//
// arguments:
//  MMatrix<TScalar, TIntegral>& this: (output) class operand
//  const MVector<TScalar, TIntegral>& arg: (input) new values
//
// return: a boolean value indicating status
//
// this method sets the diagonal elements of the current matrix to be
// the same as the elements of the input vector
//
template<class TScalar, class TIntegral>
boolean MMatrixMethods::setDiagonal(MMatrix<TScalar, TIntegral>& this_a,
				    const MVector<TScalar, TIntegral>& arg_a) {

  // declare local variables
  //
  long nrows = arg_a.length();

  // check if the current matrix is square and the dimension of the
  // vector and matrix match
  //
  if ((!this_a.isSquare()) || (this_a.nrows_d != nrows)) {
    this_a.debug(L"this_a");    
    return Error::handle(name(), L"setDiagonal", Error::ARG, 
			 __FILE__, __LINE__);
  }

  // type: non-SPARSE
  //
  if ((this_a.type_d == Integral::FULL) ||
      (this_a.type_d == Integral::DIAGONAL) ||
      (this_a.type_d == Integral::SYMMETRIC) ||
      (this_a.type_d == Integral::LOWER_TRIANGULAR) ||
      (this_a.type_d == Integral::UPPER_TRIANGULAR)) {

    // only copy the diagonal elements from this matrix
    //
    for (long row = 0; row < nrows; row++) {
      this_a.setValue(row, row, arg_a(row));
    }
  }

  // type: SPARSE
  //
  else {

    // the only really clean way to do this is to manually copy all
    // the non-zero elements. for a truly sparse matrix, this won't be
    // prohibitively expensive.
    //
    for (long row = 0; row < nrows; row++) {

      // since looking up this value is expensive, we want to save it
      //
      TIntegral tmp = arg_a(row);
      
      // add the value if it is non-zero
      //
      if (tmp != 0) {
	this_a.setValue(row, row, tmp);
      }
    }
  }

  // exit gracefully
  //
  return true;
}

// explicit instantiations
//
template boolean
MMatrixMethods::setDiagonal<ISIP_TEMPLATE_TARGET>
(MMatrix<ISIP_TEMPLATE_TARGET>&, const MVector<ISIP_TEMPLATE_TARGET>&);

// method: setLower
//
// arguments:
//  MMatrix<TScalar, TIntegral>& this: (output) class operand
//  const MMatrix<TScalar, TIntegral>& arg: (input) new values
//
// return: a boolean value indicating status
//
// this method sets the lower triangular elements of the matrix to be
// the same as the lower triangular elements of the source matrix
//
template<class TScalar, class TIntegral>
boolean MMatrixMethods::setLower(MMatrix<TScalar, TIntegral>& this_a,
				 const MMatrix<TScalar, TIntegral>& arg_a) {