mmat_00.cc

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

  //
  // note that since the matrix is sparse, and values are not preserved,
  // we simply clear the SPARSE vectors.
  //
  else {

    // free space
    //
    m_d.clear(Integral::RELEASE);

    // allocate space
    //
    m_d.setLength(0);
    row_index_d.setLength(0);
    col_index_d.setLength(0);
  }

  // set the dimensions
  //
  nrows_d = nrows_a;
  ncols_d = ncols_a;

  // install the new type
  //
  type_d = type_a;

  // exit gracefully
  //
  return true;
}

// method: vecResetCapacity
//
// arguments:
//  long nrows: (input) number of rows
//  long ncols: (input) number of columns
//  Integral::MTYPE type: (input) desired type
//
// return: a boolean value containing the status
//  
// this method is a very specific method used to facilitate setCapacity.
// it deletes resizes old data, creating new capacity according to the type.
// old data is, by definition, destroyed.
//
template<class TScalar, class TIntegral>
boolean
MMatrix<TScalar, TIntegral>::vecResetCapacity(long nrows_a, long ncols_a,
					      Integral::MTYPE type_a) {

  // note that note that we only need to check for
  // a crossconnect between SPARSE and all other types, because
  // SPARSE uses a unique storage format. setLength will take care
  // of the rest.
  //
  // condition: new matrix is not SPARSE
  // action: delete space used by the SPARSE representation
  //
  if (type_a != Integral::SPARSE) {
    
    // free space
    //
    row_index_d.clear(Integral::RELEASE);
    col_index_d.clear(Integral::RELEASE);
    
    // allocate space
    //
    long new_size = vecLength(nrows_a, ncols_a, type_a);
    if (m_d.length() > new_size) {
      m_d.setLength(new_size, false);
    }
    m_d.setCapacity(new_size, false);
  }
  
  // condition: new matrix is SPARSE
  // action: delete space used by the non-SPARSE representation
  //
  // note that since the matrix is sparse, and values are not preserved,
  // we simply clear the SPARSE vectors.
  //
  else {
    
    // free space
    //
    m_d.clear(Integral::RELEASE);
    
    // allocate space
    //
    long num_elements = nrows_a * ncols_a;

    if (m_d.length() > num_elements) {
      m_d.setLength(num_elements, false);
    }
    m_d.setCapacity(num_elements, false);
    if (row_index_d.length() > num_elements) {
      row_index_d.setLength(num_elements, false);
    }
    row_index_d.setCapacity(num_elements, false);
    if (col_index_d.length() > num_elements) {
      col_index_d.setLength(num_elements, false);
    }
    col_index_d.setCapacity(num_elements, false);
  }
  
  // install the new type: but we still have to set the dimensions
  //
  type_d = type_a;

  // exit gracefully
  //
  return true;
}

// method: vecResizeCapacity
//
// arguments:
//  long nrows: (input) number of rows
//  long ncols: (input) number of columns
//
// return: a boolean value containing the status
//  
// this method is a very specific method used to facilitate setCapacity.
// it is only called when the type of the matrix need not be changed.
// it invokes setCapacity for the internal data.
// old data is, by definition, retained.
//
template<class TScalar, class TIntegral>
boolean MMatrix<TScalar, TIntegral>::
vecResizeCapacity(long nrows_a, long ncols_a) {

  // note that note that we only need to check for
  // a crossconnect between SPARSE and all other types, because
  // SPARSE uses a unique storage format.
  //
  // condition: new matrix is not SPARSE
  //
  if (type_d != Integral::SPARSE) {
    long new_size = vecLength(nrows_a, ncols_a, type_d);
    if (m_d.length() > new_size) {
      m_d.setLength(new_size, true);
    }
    m_d.setCapacity(new_size, true);
  }

  // condition: new matrix is SPARSE
  //
  else {

    // allocate space
    //
    long num_elements = nrows_a * ncols_a;
    if (m_d.length() > num_elements) {
      m_d.setLength(num_elements, true);
    }
    m_d.setCapacity(num_elements, true);
    if (row_index_d.length() > num_elements) {
      row_index_d.setLength(num_elements, true);
    }
    row_index_d.setCapacity(num_elements, true);
    if (col_index_d.length() > num_elements) {
      col_index_d.setLength(num_elements, true);
    }
    col_index_d.setCapacity(num_elements, true);
  }

  // exit gracefully
  //
  return true;
}

// method: vecDimensionLength
//
// arguments:
//  long nrows: (input) number of rows
//  long ncols: (input) number of columns
//  boolean preserve_values: (input) flag to save existing values
//
// return: a boolean value containing the status
//  
// this method is a very specific method used to facilitate setDimensions.
// it is only called when the type of the matrix need not be changed.
// it invokes a change of length for the internal data.
// old data is, by definition, retained.
//
template<class TScalar, class TIntegral>
boolean MMatrix<TScalar, TIntegral>::
vecDimensionLength(long nrows_a, long ncols_a,boolean preserve_values_a) {

  // condition: SPARSE
  //  this is really easy since sparse only saves non-zero values element
  //  by element. there is one catch: we might need to delete elements if
  //  the new dimensions are smaller.
  //
  if (type_d == Integral::SPARSE) {

    // check for a request to set the lengths to zero: it is best to do
    // this efficiently
    //
    if ((nrows_a == 0) && (ncols_a == 0)) {
      row_index_d.clear(Integral::RESET);
      col_index_d.clear(Integral::RESET);
      m_d.clear(Integral::RESET);
    }

    // else: delete element by element
    //
    else {

      // loop over all elements
      //
      long num_elements = row_index_d.length();
      long i = 0;

      while (i < num_elements) {
	if ((row_index_d(i) >= nrows_a) || (col_index_d(i) >= ncols_a)) {
	  row_index_d.deleteRange(i, 1);
	  col_index_d.deleteRange(i, 1);
	  m_d.deleteRange(i, 1);
	  num_elements--;
	}
	else {
	  i++;
	}
      }
    }

    // set the dimensions
    //
    nrows_d = nrows_a;
    ncols_d = ncols_a;
    return true;
  }

  // condition: preserve_values = false
  //  this is also easy. we simply need to set the vector length.
  //  the vector length method is smart enough to figure out what it
  //  has to do with the capacity.
  //
  else if (!preserve_values_a) {

    // set the dimensions
    //
    nrows_d = nrows_a;
    ncols_d = ncols_a;

    // set the length
    //
    return m_d.setLength(vecLength(nrows_a, ncols_a, type_d), false);
  }

  // condition: preserve_values = true
  //  at this point in the code, we can only have a non-sparse matrix.
  //  that makes things a little easier. but we still have to do some
  //  type-specific grungy stuff.
  //
  else {

    // copy the data
    //
    TVector m_old(m_d);

    // increase/decrease the length of the data vector. be sure to
    // clear it, since the contents are unpredictable.
    //
    m_d.setLength(vecLength(nrows_a, ncols_a, type_d), false);    
    m_d.clear(Integral::RETAIN);

    // copy only the intersection of the old and new matrices
    //
    long nrows_new = (long)Integral::min(nrows_a, nrows_d);
    long ncols_new = (long)Integral::min(ncols_a, ncols_d);

    // type: FULL
    //  loop over all elements (not the most efficient way, but simple)
    //
    if (type_d == Integral::FULL) {
      for (long row = 0; row < nrows_new; row++) {
	for (long col = 0; col < ncols_new; col++) {
	  m_d(index(row, col, nrows_a, ncols_a)) =
	    m_old(index(row, col, nrows_d, ncols_d));
	}
      }
    }

    // type: DIAGONAL
    //  copy the vector directly
    //
    else if (type_d == Integral::DIAGONAL) {
      m_d.move(m_old, nrows_new, 0, 0);
    }

    // type: SYMMETRIC or LOWER_TRIANGULAR
    //  loop over all elements
    //
    else if ((type_d == Integral::SYMMETRIC) ||
	     (type_d == Integral::LOWER_TRIANGULAR)) {
      for (long row = 0; row < nrows_new; row++) {
	for (long col = 0; col <= row; col++) {
	  m_d(index(row, col, nrows_a, ncols_a)) =
	    m_old(index(row, col, nrows_d, ncols_d));
	}
      }
    }

    // type: UPPER_TRIANGULAR
    //  loop over all elements
    //
    else if ((type_d == Integral::SYMMETRIC) ||
	     (type_d == Integral::UPPER_TRIANGULAR)) {
      for (long row = 0; row < nrows_new; row++) {
	for (long col = row; col < ncols_new; col++) {
	  m_d(index(row, col, nrows_a, ncols_a)) =
	    m_old(index(row, col, nrows_d, ncols_d));
	}
      }
    }

    // update the dimensions
    //
    nrows_d = nrows_a;
    ncols_d = ncols_a;

    // exit gracefully
    //
    return true;
  }
}

// method: sortSparse
//
// arguments: none
//
// return: a boolean value containing the status
//  
// this method is a very specific method used to facilitate sparse matrix
// manipulations. it takes an existing sparse matrix and sorts the row
// and column indices to make sure they are in the proper order.
//
// this code uses a radix sort. the first pass sorts the columns, the second
// sorts the rows using a stable approach. the algorithm is described in:
//
//  Thomas H. Cormen, Charles E. Leiserson, Ronald L. Rivest,
//  Introduction to Algorithms, The MIT Press, Cambridge, Massachusetts, USA,
//  McGraw-Hill Book Company, New York, USA, pp. 175-179, 1997.
//
template<class TScalar, class TIntegral>
boolean MMatrix<TScalar, TIntegral>::sortSparse() {

  // check the type
  //
  if (type_d != Integral::SPARSE) {
    this->debug(L"this");    
    Error::handle(name(), L"sortSparse", Error::ARG, __FILE__, __LINE__);
  }
  
  // declare a buffer: a counting sort needs scratch space for counts
  //
  MVector<Long, long> counts((long)Integral::max(nrows_d, ncols_d));

  // this algorithm does not sort in place, so we need temporary
  // buffers for the three vectors that comprise a sparse matrix.
  //
  long nelem = m_d.length();
  TVector tmp_m(nelem);
  MVector<Long, long> tmp_ri(nelem);
  MVector<Long, long> tmp_ci(nelem);

  // first sort the columns with a counting sort:
  //   A: col_index_d
  //   B: tmp_ci (note that the other two vectors are reordered simultaneously)
  //   C: counts  
  // after this loop counts(i) will contain the number of elements equal to i
  //
  for (long i = 0; i < nelem; i++) {
    counts(col_index_d(i)) += 1;
  }
  
  // after this loop counts(i) will contain the number of elements
  // less than or equal to index
  //
  for (long i = 1; i < ncols_d; i++) {
    counts(i) += counts(i - 1);
  }

  // copy the data vectors into the locations specified by counts()
  //
  for (long j = nelem - 1; j >= 0; j--) {
    long k = (long)counts(col_index_d(j)) - 1;
    tmp_ri(k) = row_index_d(j);
    tmp_ci(k) = col_index_d(j);
    tmp_m(k) = m_d(j);
    counts(col_index_d(j)) -= 1;
  }

  // we now have the three temporary vectors sorted by columns, so
  // perform the second pass of counting sort to sort by rows.
  //   A: tmp_ri
  //   B: row_index_d
  //   C: counts
  //
  counts.clear(Integral::RETAIN);
  
  // after this loop, counts(i) will contain the number of elements equal to i
  //
  for (long i = 0; i < nelem; i++) {
    counts(tmp_ri(i)) += 1;
  }

  // after this loop, counts(i) will contain the number of elements
  // less than or equal to i
  //
  for (long i = 1; i < nrows_d; i++) {
    counts(i) += counts(i - 1);
  }

  // now copy the data vectors into the locations specified by counts()
  //
  for (long j = nelem - 1; j >= 0; j--) {
    long k = (long)counts(tmp_ri(j)) - 1;
    row_index_d(k) = tmp_ri(j);
    col_index_d(k) = tmp_ci(j);
    m_d(k) = tmp_m(j);
    counts(tmp_ri(j)) -= 1;
  }

  // exit gracefully
  //
  return true;
}

// method: isSingular
//
// arguments:
//  const MMatrix& matrix: (input) the matrix to compute singularity on
//  double thresh: (input) if the determinant is less than this value then the
//                         matrix is considered singular
//
// return: a boolean value containing the status
//
// this method checks if the matrix is a singular matrix or not. this
// method doesn't work for unsigned matrices as the determinant for
// unsigned matrices is not defined
//
template<class TScalar, class TIntegral>
boolean MMatrix<TScalar, TIntegral>::isSingular(const MMatrix& matrix_a, double thresh_a) const {

#ifdef ISIP_TEMPLATE_unsigned
  this->debug(L"this");  
  return Error::handle(name(), L"isSingular", Error::TEMPLATE_TYPE, __FILE__, __LINE__);
#else

  // for diagonal or triangular matrices, the matrix is non-singular
  // so long as there are no zero-entries on the diagonal and we are not
  // thresholding
  //
  if ((thresh_a == 0) && ((matrix_a.type_d == Integral::DIAGONAL) ||
			  (matrix_a.type_d == Integral::LOWER_TRIANGULAR) ||
			  (matrix_a.type_d == Integral::LOWER_TRIANGULAR))) { 
    long nrows = matrix_a.nrows_d;
    for (long row = 0; row < nrows; row++) {
      if (matrix_a.m_d(matrix_a.index(row, row)) == (TIntegral)0.0) {
	return true;
      }
    }
    return false;
  }
  else {
    return (((thresh_a == 0) && (matrix_a.determinant() == 0)) ||
	    ((thresh_a != 0) && (Integral::abs(matrix_a.determinant()) < thresh_a)));
  }
#endif
}

// explicit instantiations
//
template class MMatrix<ISIP_TEMPLATE_TARGET>;