mmat_07.cc

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

  //
  if (this_a.type_d == Integral::FULL) {
    for (long row = 0, index = 0; row < nrows; row++) {
      sum += this_a.m_d(index);
      index += ncols + 1;
    }
  }

  // type: DIAGONAL
  //  do this efficiently using a vector function
  //
  else if (this_a.type_d == Integral::DIAGONAL) {
    sum = this_a.sum();
  }

  // type: SYMMETRIC, LOWER_TRIANGULAR, UPPER_TRIANGULAR
  //  do this efficiently using direct indexing
  //
  else if ((this_a.type_d == Integral::SYMMETRIC) ||
	   (this_a.type_d == Integral::LOWER_TRIANGULAR) ||
	   (this_a.type_d == Integral::UPPER_TRIANGULAR)) {
    for (long row = 0, index = 0; row < nrows; row++) {
      sum += this_a.m_d(index);
      index += 2 + row;
    }
  }

  // type: SPARSE
  //  use indirect addressing through the index function
  //  
  else {
    
    // loop over rows and add the diagonal elements
    //
    for (long row = 0; row < this_a.nrows_d; row++) {
      sum += this_a.getValue(row, row);
    }
  }

  // exit gracefully
  //
  return sum;
}

// explicit instantiations
//
template ISIP_TEMPLATE_T1   
MMatrixMethods::trace<ISIP_TEMPLATE_TARGET>
(const MMatrix<ISIP_TEMPLATE_TARGET>&);

// method: rank
//
// arguments:
//  const MMatrix<TScalar, TIntegral>& this: (input) class operand
//
// return: a value containing the rank of the matrix
//
// this method calculates the rank of a matrix (someday soon...)
//
template<class TScalar, class TIntegral>
long MMatrixMethods::rank(const MMatrix<TScalar, TIntegral>& this_a) {

  // this method is defined only for float and double matrices 
  //
#ifdef ISIP_TEMPLATE_fp
  
  MMatrix<TScalar, TIntegral> u;  // row orthonormal
  MMatrix<TScalar, TIntegral> w;  // singular value matrix
  MMatrix<TScalar, TIntegral> v;  // column orthonormal

  this_a.decompositionSVD(u, w, v);
  
  long m = w.getNumRows(); // w has same row and column 
  long i;

  // set up near-singular values.
  //
  TIntegral wmin = 0.0;
  TIntegral wmax = this_a.max();

  if (typeid(TIntegral) == typeid(double)) {
    wmin = wmax * m * w.THRESH_SINGULAR_DOUBLE;
  }
  else if (typeid(TIntegral) == typeid(float)) {
    wmin = wmax * m * w.THRESH_SINGULAR_FLOAT;
  }

  // check the data in w matrix
  //
  for (i = 0; i < m; i++) {
    if (w.getValue(i, i) < wmin) break;
  }

  return i;
  
#else
  return Error::handle(name(), L"rank", Error::TEMPLATE_TYPE,
		       __FILE__, __LINE__);
#endif

  this_a.debug(L"this_a");

  return Error::handle(name(), L"rank", Error::NOT_IMPLEM, __FILE__, __LINE__);

}

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

// method: swapRows
//
// arguments:
//  MMatrix<TScalar, TIntegral>& this: (output) class operand
//  const MMatrix<TScalar, TIntegral>& matrix: (input) source matrix
//  long row1: (input) first row for swap
//  long row2: (input) second row for swap
//
// return: a boolean value indicating status
//
// this method swaps two rows of the input matrix. the input matrix can only
// be full or sparse, since this operation will change the type of
// a diagonal, symmetric, or triangular matrix.
//
template<class TScalar, class TIntegral>
boolean MMatrixMethods::swapRows(MMatrix<TScalar, TIntegral>& this_a, 
				 const MMatrix<TScalar, TIntegral>& matrix_a, 
				 long row1_a, long row2_a) {

  // check the arguments
  //
  if (row1_a == row2_a) {
    return true;
  }
  else if ((row1_a < 0) || (row2_a < 0) ||
	   (row1_a >= matrix_a.nrows_d) || (row2_a >= matrix_a.nrows_d)) {
    this_a.debug(L"this_a");
    matrix_a.debug(L"matrix_a");    
    return Error::handle(name(), L"swapRows", Error::ARG, __FILE__, __LINE__);
  }

  // copy input matrix to output matrix
  //
  Integral::MTYPE old_type = this_a.type_d;
  if (&this_a != &matrix_a) {
    this_a.assign(matrix_a);
  }
  
  // implement this using getRow:
  //  create copies of the rows
  //
  MVector<TScalar, TIntegral> row1;
  MVector<TScalar, TIntegral> row2;
  this_a.getRow(row1, row1_a);
  this_a.getRow(row2, row2_a);

  // swap them
  //
  this_a.setRow(row1_a, row2);
  this_a.setRow(row2_a, row1);

  // restore the type
  //
  return this_a.changeType(old_type);
}

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

// method: swapColumns
//
// arguments:
//  MMatrix<TScalar, TIntegral>& this: (output) class operand
//  const MMatrix<TScalar, TIntegral>& matrix: (input) class operand
//  long col1: (input) first column for swap
//  long col2: (input) second column for swap
//
// return: a boolean value indicating status
//
// this method swaps two columns of the input matrix and stores it
// into the output matrix. the input matrix can only be full or sparse,
// since this operation will change the type of a diagonal, symmetric,
// or triangular matrix.
//
template<class TScalar, class TIntegral>
boolean
MMatrixMethods::swapColumns(MMatrix<TScalar, TIntegral>& this_a, 
			    const MMatrix<TScalar, TIntegral>& matrix_a, 
			    long col1_a, long col2_a) {

  // check the arguments
  //
  if (col1_a == col2_a) {
    return true;
  }
  else if ((col1_a < 0) || (col2_a < 0) ||
	   (col1_a >= matrix_a.ncols_d) || (col2_a >= matrix_a.ncols_d)) {
    this_a.debug(L"this_a");
    matrix_a.debug(L"matrix_a");  
    return Error::handle(name(), L"swapColumns",
			 Error::ARG, __FILE__, __LINE__);
  }
  
  // copy input matrix to output matrix, save the type
  //
  Integral::MTYPE old_type = this_a.type_d;
  this_a.assign(matrix_a);
  
  // implement this using getColumn: create copies of the columns
  //
  MVector<TScalar, TIntegral> col1;
  MVector<TScalar, TIntegral> col2;
  this_a.getColumn(col1, col1_a);
  this_a.getColumn(col2, col2_a);

  // swap them
  //
  this_a.setColumn(col1_a, col2);
  this_a.setColumn(col2_a, col1);

  // exit gracefully
  //
  return this_a.changeType(old_type);
}

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

// method: determinantLU
//
// arguments:
//  const MMatrix<TScalar, TIntegral>& this: (output) class operand
//
// return: a TIntegral value containing the determinant
//
// this method calculates the determinant of the matrix, attempting to
// do this by copying the matrix into a MMatrixDouble matrix and use
// LU decomposition as the fallback algorithm
//
// through the combination of division and subtraction in the LU
// decomposition method, we may potentially encounter precision problems
// in the LU calculations. elements on the diagonal that should be
// zero may be returned by the LU computation as very small (on the
// order of 1e-8 for floats and 1e-16 for doubles) non-zero numbers
// thus, we will floor these numbers.
//
// unfortunately, this means we will have to loop over the
// upper-triangular matrix's diagonal elements twice - once to find
// the dynamic range and once to zero out the elements which fall
// outside of an epsilon of the dynamic range. finding a way to move
// this flooring directly to the LU decomposition would be a better
// way to handle this.
//
template<class TScalar, class TIntegral>
TIntegral MMatrixMethods::determinantLU(const MMatrix<TScalar, TIntegral>& this_a) {

  // declare slack variables to control round-off error. The machine epsilon
  // for 32-bit floats is typically 3e-8. We back off of that by an order
  // of magnitude to give some breathing room.
  //
  static double det_slack_var = 0.0;

  if (typeid(TIntegral) == typeid(float)) {
    det_slack_var = 3e-7;
  }
  else if (typeid(TIntegral) == typeid(double)) {
    det_slack_var = 3e-15;
  }
  else {
    det_slack_var = 3e-7;
  }
  
  // declare local variables
  //
  MMatrix<TScalar, TIntegral> l;
  MMatrix<TScalar, TIntegral> u;
  MVector<Long, long> index;
  long sign;

  this_a.decompositionLU(l, u, index, sign, 0);

  // get the mean of the elements on the diagonal of the upper-triangular
  //
  double mean = 0;
  for (long i = 0; i < this_a.nrows_d; i++) {
    mean += Integral::abs(u(i, i));
  }
  mean /= (double)this_a.nrows_d;

  // get the minimum magnitude that we will accept - make sure that the
  // floor is reasonably close to zero
  //
  double floor = Integral::abs(mean) * det_slack_var;
  if (floor > det_slack_var) {
    floor = det_slack_var;
  }

  // loop over the diagonal elements again and zero out those that
  // are below the floor so long as the floor is reasonably close to zero
  //
  for (long i = 0; i < this_a.nrows_d; i++) {
    if (Integral::abs(u(i, i)) < floor) {
      return 0;
    }
  }
  
  // since the lower triangular component has 1's along the diagonal,
  // we can just calculate the determinant of the upper triangular
  // matrix. We have to make a special case for complex numbers so the
  // compiler does not complain
  //
#ifdef ISIP_TEMPLATE_complex
  return ((TIntegral)sign * u.determinant());
#else
  return (sign * u.determinant());
#endif
}

// explicit instantiations
//
template ISIP_TEMPLATE_T1
MMatrixMethods::determinantLU<ISIP_TEMPLATE_TARGET>
(const MMatrix<ISIP_TEMPLATE_TARGET>&);

// method: determinantMinor
//
// arguments:
//  const MMatrix<TScalar, TIntegral>& this: (output) class operand
//
// return: a TIntegral value containing the determinant
//
// this method calculates the determinant of a sparse matrix, attempting to
// do this using recursive getMinor method
//
template<class TScalar, class TIntegral>
TIntegral MMatrixMethods::determinantMinor(const MMatrix<TScalar, TIntegral>& this_a) {

  // declare local variables
  //
  TIntegral det = 0;

  // type: FULL
  //
  if (this_a.type_d == Integral::FULL) {

    // define full matrix to hold the minor matrix
    //
    MMatrix<TScalar, TIntegral> mat_minor(this_a.nrows_d, this_a.ncols_d, 
					  Integral::FULL);

    // get the minor matrix and use recursive method to compute determinant
    //    
    long sign = 1;
    for (int i = 0; i < this_a.ncols_d; i++) {
      
      // get the minor matrix
      //
      this_a.getMinor(mat_minor, 0, i);
      if (sign == 1) {
	det += this_a.m_d(i) * mat_minor.determinant();
      }
      else {
	det -= this_a.m_d(i) * mat_minor.determinant();
      }
      sign = -sign;
    }
  }

  // type: SYMMETRIC
  //  for symmetric matrix, change it into a full matrix before calculating
  //  the determinant
  //
  else if (this_a.type_d == Integral::SYMMETRIC) {
    MMatrix<TScalar, TIntegral> temp;
    temp.assign(this_a, Integral::FULL);
    return (TIntegral)temp.determinant(); 
  }

  // type: SPARSE
  //
  else {
    
    // define row and column vectors
    //
    MVector<Long, long> rows(this_a.nrows_d);
    MVector<Long, long> cols(this_a.ncols_d);
    
    // define sparse matrix to hold the minor matrix
    //
    MMatrix<TScalar, TIntegral> mat_minor(this_a.nrows_d, this_a.ncols_d, 
					  Integral::SPARSE);
    
    // assigning number to every row and column of the real matrix
    //
    for (long i = 0; i < this_a.m_d.length(); i++) {
      rows(this_a.row_index_d(i)) += 1;
      cols(this_a.col_index_d(i)) += 1;
    }
    
    // initializing row_to_use (the row having lowest number of
    // non-zero elements) to zero
    //
    long row_to_use = 0;
    
    // if any row has all zeros then determinant is zero
    //
    if (rows.min(row_to_use) == 0) {
      return 0;
    }
    
    // if all rows have at least one non-zero element
    //
    else {
      for (long i = 0; i < this_a.ncols_d; i++) {
	
	// calculate only if the value is non-zero
	//
	TIntegral tmp_val = this_a.getValue(row_to_use, i);
	if (tmp_val != 0) {
	  
	  // compute the minor of the matrix
	  //
	  this_a.getMinor(mat_minor, row_to_use, i);
	  
	  // get the sign for a particular value and multiply value by
	  // the determinant of the minor to get the determinant
	  //
	  if (((row_to_use + i) % 2) == 0) {
	    det += mat_minor.determinant() * tmp_val;
	  }
	  else {
	    det -= mat_minor.determinant() * tmp_val;
	  }
	}
      }
    }
  }

  // exit gracefully: return the value
  //
  return (TIntegral)det;
}

// explicit instantiations
//
template ISIP_TEMPLATE_T1
MMatrixMethods::determinantMinor<ISIP_TEMPLATE_TARGET>
(const MMatrix<ISIP_TEMPLATE_TARGET>&);