mmat_09.cc

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

    }
    else {
      this_a.setDimensions(nrows, ncols, false, Integral::FULL);
    }
    this_a.copy(m2_a);

    // add elements of the first matrix
    //
    for (long index = 0; index < m1_a.m_d.length(); index++) {
      long row = m1_a.row_index_d(index);
      long col = m1_a.col_index_d(index);
      this_a.setValue(row, col, this_a(row, col) + m1_a(row, col));
    }
  }

  // if the second type is sparse, and the first type isn't sparse,
  // we assign the output type to FULL (except for DIAGONAL), copy the
  // elements of the first, and then add elements of the second
  //
  else if ((type1 != Integral::SPARSE) && (type2 == Integral::SPARSE)) {

    // create a FULL matrix for all types except DIAGONAL, and zero it
    //
    if (m1_a.type_d == Integral::DIAGONAL) {
      this_a.assign(m1_a, Integral::SPARSE);
    }
    else {
      this_a.assign(m1_a, Integral::FULL);
    }
    
    // add elements of the second sparse matrix
    //
    for (long index = 0; index < m2_a.m_d.length(); index++) {
      long row = m2_a.row_index_d(index);
      long col = m2_a.col_index_d(index);
      this_a.setValue(row, col, (TIntegral)(this_a(row, col) + m2_a(row, col)));
    }
  }

  // if either type is FULL, we assign the output type to FULL,
  // copy the elements from the first, then add elements of the second
  //
  else if (type1 == Integral::FULL) {

    // create a FULL matrix through assignment
    //
    this_a.assign(m1_a, Integral::FULL);
    
    // add elements of the second matrix
    //
    for (long row = 0; row < nrows; row++) {

      // add the values
      //
      long start_col = m2_a.startColumn(row, this_a.type_d, type2);
      long stop_col = m2_a.stopColumn(row, this_a.type_d, type2);

      for (long col = start_col; col <= stop_col; col++) {
	this_a.setValue(row, col, (TIntegral)(this_a(row, col) + m2_a(row, col)));
      }
    }
  }

  // if either type is FULL, we assign the output type to FULL,
  // copy the elements from the first, then add elements of the second
  //
  else if (type2 == Integral::FULL) {

    // create a FULL matrix through assignment
    //
    this_a.assign(m2_a, Integral::FULL);
    
    // add elements of the second matrix
    //
    for (long row = 0; row < nrows; row++) {

      // add the values
      //
      long start_col = m1_a.startColumn(row, type1, this_a.type_d);
      long stop_col = m1_a.stopColumn(row, type1, this_a.type_d);

      for (long col = start_col; col <= stop_col; col++) {
	this_a.setValue(row, col, this_a(row, col) + m1_a(row, col));
      }
    }
  }
  
  //---------------------------------------------------------------------------
  // in this section of the code, we deal with the remaining cases where
  // the matrices are of different types, and some intelligent conversion
  // must be done.
  //---------------------------------------------------------------------------
  
  // type1: DIAGONAL
  //  the remaining types are SYMMETRIC and *_TRIANGULAR. 
  //  the output types follow the type of the second operand.
  //  the approach is to create a matrix, copy the second operand,
  //  and add the first operand.
  // 
  else if (type1 == Integral::DIAGONAL) {
    
    // set up the output matrix and assign values
    //
    this_a.assign(m2_a, type2);

    // add the diagonal
    //
    for (long row = 0; row < nrows; row++) {  
      this_a.m_d(this_a.index(row, row)) += (TIntegral)m1_a.m_d(row);
    }
  }

  else if (type2 == Integral::SYMMETRIC) {

    // set up the output matrix and assign values
    //
    this_a.assign(m2_a, Integral::FULL);

    // add elements of the first matrix
    //
    for (long row = 0; row < nrows; row++) {

      // add the values
      //
      long start_col = m1_a.startColumn(row, type1, this_a.type_d);
      long stop_col = m1_a.stopColumn(row, type1, this_a.type_d);

      for (long col = start_col; col <= stop_col; col++) {
	this_a.m_d(this_a.index(row, col)) += 
	  (TIntegral)m1_a.m_d(m1_a.index(row, col));
      }
    }
  }
    
  // type1: SYMMETRIC, *_TRIANGULAR
  //  the remaining types are DIAGONAL, SYMMETRIC, and *_TRIANGULAR. 
  //  the output types vary with the input.
  //  the approach is to copy the values of the first operand,
  //  and add the values of the second operand.
  // 
  else {
    
    // set up the output matrix: a diagonal matrix preserves the type
    // of the input; all other types force a full matrix.
    //
    if ((type1 == Integral::SYMMETRIC) &&
	(type2 == Integral::DIAGONAL)) {
      this_a.assign(m1_a, Integral::SYMMETRIC);
    }
    if ((type1 == Integral::LOWER_TRIANGULAR) &&
	(type2 == Integral::DIAGONAL)) {
      this_a.assign(m1_a, Integral::LOWER_TRIANGULAR);
    }
    else if ((type1 == Integral::UPPER_TRIANGULAR) &&
	     (type2 == Integral::DIAGONAL)) {
      this_a.assign(m1_a, Integral::UPPER_TRIANGULAR);
    }
    else {
      this_a.assign(m1_a, Integral::FULL);
    }

    // add elements of the second matrix
    //
    for (long row = 0; row < nrows; row++) {

      // add the values
      //
      long start_col = m2_a.startColumn(row, this_a.type_d, type2);
      long stop_col = m2_a.stopColumn(row, this_a.type_d, type2);

      for (long col = start_col; col <= stop_col; col++) {
	this_a.m_d(this_a.index(row, col)) += 
	  (TIntegral)m2_a.m_d(m2_a.index(row, col));
      }
    }
  }
  
  // exit gracefully
  //
  return this_a.changeType(old_type);
}

// explicit instantiations
//
#ifndef ISIP_TEMPLATE_complex
template boolean
MMatrixMethods::add<ISIP_TEMPLATE_TARGET, Byte, byte>
(MMatrix<ISIP_TEMPLATE_TARGET>&, const MMatrix<ISIP_TEMPLATE_TARGET>&,
 const MMatrix<Byte, byte>&);

template boolean
MMatrixMethods::add<ISIP_TEMPLATE_TARGET, Ushort, ushort>
(MMatrix<ISIP_TEMPLATE_TARGET>&, const MMatrix<ISIP_TEMPLATE_TARGET>&,
 const MMatrix<Ushort, ushort>&);

template boolean
MMatrixMethods::add<ISIP_TEMPLATE_TARGET, Ulong, ulong>
(MMatrix<ISIP_TEMPLATE_TARGET>&, const MMatrix<ISIP_TEMPLATE_TARGET>&,
 const MMatrix<Ulong, ulong>&);

template boolean
MMatrixMethods::add<ISIP_TEMPLATE_TARGET, Ullong, ullong>
(MMatrix<ISIP_TEMPLATE_TARGET>&, const MMatrix<ISIP_TEMPLATE_TARGET>&,
 const MMatrix<Ullong, ullong>&);

template boolean
MMatrixMethods::add<ISIP_TEMPLATE_TARGET, Short, short>
(MMatrix<ISIP_TEMPLATE_TARGET>&, const MMatrix<ISIP_TEMPLATE_TARGET>&,
 const MMatrix<Short, short>&);

template boolean
MMatrixMethods::add<ISIP_TEMPLATE_TARGET, Long, long>
(MMatrix<ISIP_TEMPLATE_TARGET>&, const MMatrix<ISIP_TEMPLATE_TARGET>&,
 const MMatrix<Long, long>&);

template boolean
MMatrixMethods::add<ISIP_TEMPLATE_TARGET, Llong, llong>
(MMatrix<ISIP_TEMPLATE_TARGET>&, const MMatrix<ISIP_TEMPLATE_TARGET>&,
 const MMatrix<Llong, llong>&);

template boolean
MMatrixMethods::add<ISIP_TEMPLATE_TARGET, Float, float>
(MMatrix<ISIP_TEMPLATE_TARGET>&, const MMatrix<ISIP_TEMPLATE_TARGET>&,
 const MMatrix<Float, float>&);

template boolean
MMatrixMethods::add<ISIP_TEMPLATE_TARGET, Double, double>
(MMatrix<ISIP_TEMPLATE_TARGET>&, const MMatrix<ISIP_TEMPLATE_TARGET>&,
 const MMatrix<Double, double>&);
#endif

#ifdef ISIP_TEMPLATE_complex
template boolean
MMatrixMethods::add<ISIP_TEMPLATE_TARGET, ComplexLong, complexlong>
(MMatrix<ISIP_TEMPLATE_TARGET>&, const MMatrix<ISIP_TEMPLATE_TARGET>&,
 const MMatrix<ComplexLong, complexlong>&);

template boolean
MMatrixMethods::add<ISIP_TEMPLATE_TARGET, ComplexFloat, complexfloat>
(MMatrix<ISIP_TEMPLATE_TARGET>&, const MMatrix<ISIP_TEMPLATE_TARGET>&,
 const MMatrix<ComplexFloat, complexfloat>&);

template boolean
MMatrixMethods::add<ISIP_TEMPLATE_TARGET, ComplexDouble, complexdouble>
(MMatrix<ISIP_TEMPLATE_TARGET>&, const MMatrix<ISIP_TEMPLATE_TARGET>&,
 const MMatrix<ComplexDouble, complexdouble>&);
#endif

// method: sub
//
// arguments:
//  MMatrix<TScalar, TIntegral>& this: (output) class operand
//  TIntegral value: (input) scalar value
//
// return: a boolean value indicating status
//
template<class TScalar, class TIntegral>
boolean MMatrixMethods::sub(MMatrix<TScalar, TIntegral>& this_a, 
			    TIntegral value_a) {

  // get the type of the matrix and branch on type
  //
  Integral::MTYPE type = this_a.getType();
  
  // if the current matrix is a sparse, change it to full for the
  // computation. this is admittedly inefficient, but there is really
  // no way for this computation to be efficient anyway.
  //
  if (type == Integral::SPARSE) {
    this_a.changeType(Integral::FULL);
  }

  //  now we simply need to subtract the value to the storage vector
  //
  if (!this_a.m_d.sub(value_a)) {
    return false;
  }

  // for sparse we need to change the type back
  //
  if (type == Integral::SPARSE) {
    this_a.changeType(Integral::SPARSE);
  }

  // exit gracefully
  //
  return true;
}

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

// method: sub
//
// arguments:
//  MMatrix<TScalar, TIntegral>& this: (output) class operand
//  const MMatrix<TAScalar, TAIntegral>& m2: (input) operand 2
//
// return: a boolean value indicating status
//
// the approach can best be summarized by the following flowchart and matrix:
//
// (1) special cases:
//      types are equal: process internal data directly for all types
//      one of the types is sparse: process internal data directly
//      one of the types is full: the result is always a full matrix
//      
// (2) other cases:
//      this matrix describes the type of the output matrix as a function
//      of the input types. For the operation a -= b,
//      type1 is the type of a; type 2 is the type of b.
//
//    TYPE2 --------> FULL  DIAG  SYMM  LTRI  UTRI  SPRS
//    TYPE1 ->  FULL: full  full  full  full  full  full
//              DIAG: full  diag  symm  ltri  utri  sprs
//              SYMM: full  symm  symm  full  full  full
//              LTRI: full  ltri  full  ltri  full  full
//              UTRI: full  utri  full  full  utri  full
//              SPRS: full  sprs  full  full  full  sprs
//
// the code that implements each combination attempts to exploit the
// properties of the matrices as much as possible.
//
template<class TScalar, class TIntegral, class TAScalar, class TAIntegral>
boolean MMatrixMethods::sub(MMatrix<TScalar, TIntegral>& this_a, 
			    const MMatrix<TAScalar, TAIntegral>& m2_a) {

  // check if the two input matrices have the dimensions
  //
  if (!this_a.checkDimensions(m2_a))  {
    m2_a.debug(L"m2_a");    
    this_a.debug(L"this_a");    
    return Error::handle(name(), L"sub", Error::ARG, __FILE__, __LINE__);
  }

  // declare local variables
  //
  long nrows = this_a.nrows_d;
  Integral::MTYPE type1 = this_a.type_d;
  Integral::MTYPE type2 = m2_a.type_d;
  
  //---------------------------------------------------------------------------
  // special cases:
  //  in this section, we deal with special cases that share the same
  //  code structure.
  //---------------------------------------------------------------------------

  // if the two types are the same, we can simply operate on the storage
  // vectors directly
  //
  if (type1 == type2) {

    // type: FULL, DIAGONAL, SYMMETRIC, LOWER_TRIANGULAR, UPPER_TRIANGULAR
    //
    if ((type1 == Integral::FULL) ||
	(type1 == Integral::SYMMETRIC) ||
	(type1 == Integral::DIAGONAL) ||
	(type1 == Integral::LOWER_TRIANGULAR) ||
	(type1 == Integral::UPPER_TRIANGULAR)) {
      return this_a.m_d.sub(m2_a.m_d);
    }

    // type: SPARSE
    //
    else {

      // loop through each element of the second argument
      //
      for (long index = 0; index < m2_a.m_d.length(); index++) {
	long row = m2_a.row_index_d(index);
	long col = m2_a.col_index_d(index);
        this_a.setValue(row, col, (TIntegral)(this_a(row, col) - m2_a.m_d(index)));
      }
      return true;
    }
  }

  // if the first type is sparse, and the second type isn't sparse,
  // we assign the output type to FULL (except for DIAG),
  // and then subtract elements of the second.
  //
  else if ((type1 == Integral::SPARSE) && (type2 != Integral::SPARSE)) {

    // create a FULL matrix for all types except DIAGONAL
    //
    if (type2 != Integral::DIAGONAL) {
      this_a.changeType(Integral::FULL);
    }

    // add elements of the second matrix
    //
    for (long row = 0; row < nrows; row++) {
      long start_col = m2_a.startColumn(row, this_a.type_d, type2);
      long stop_col = m2_a.stopColumn(row, this_a.type_d, type2);
      for (long col = start_col; col <= stop_col; col++) {
	this_a.setValue(row, col, (TIntegral)(this_a(row, col) - m2_a(row, col)));
      }
    }

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

  // if the second type is sparse, and the first type isn't sparse,
  // we assign the output type to FULL (except for DIAGONAL), and
  // subtract the elements of the second.
  //
  else if ((type1 != Integral::SPARSE) && (type2 == Integral::SPARSE)) {

    // create a FULL matrix for all types except DIAGONAL
    //
    if (this_a.type_d == Integral::DIAGONAL) {