mmat_09.cc

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

      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));
      }
    }
  }

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

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

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

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

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

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

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

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

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

template boolean
MMatrixMethods::sub<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::sub<ISIP_TEMPLATE_TARGET, ComplexLong, complexlong>
(MMatrix<ISIP_TEMPLATE_TARGET>&, const MMatrix<ISIP_TEMPLATE_TARGET>&,
 const MMatrix<ComplexLong, complexlong>&);

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

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

// method: mult
//
// arguments:
//  MMatrix<TScalar, TIntegral>& this: (output) class operand
//  const MMatrix<TScalar, TIntegral>& m1: (input) operand 1
//  const MMatrix<TAScalar, TAIntegral>& m2: (input) operand 2
//
// return: a boolean value indicating status
//
// this method multiplies two matrices, then assign result to the
// current matrix
// for example:
//
//       [1 2 3]       [2 0 4]                    [11 12  7]
//  m1 = [4 5 6]  m2 = [0 6 0]  m_out = m1 * m2 = [26 30 22]
//       [7 8 9]       [3 0 1]                    [41 48 37]
//
template<class TScalar, class TIntegral, class TAScalar, class TAIntegral>
boolean MMatrixMethods::mult(MMatrix<TScalar, TIntegral>& this_a, 
			     const MMatrix<TScalar, TIntegral>& m1_a, 
			     const MMatrix<TAScalar, TAIntegral>& m2_a) {

  // if the current matrix is same as one operand, a copy to temporary matrix
  // of that operand is needed
  //
  if (&this_a == &m1_a) {
    MMatrix<TScalar, TIntegral> tmp(this_a);
    return this_a.mult(tmp, m2_a) ;
  }
  else if ((void*)&this_a == (void*)&m2_a) {
    MMatrix<TScalar, TIntegral> tmp(this_a);
    return this_a.mult(m1_a, tmp);
  }

  // check dimensions
  //
  if (m1_a.getNumColumns() != m2_a.getNumRows()) {
    m1_a.debug(L"m1_a");    
    m2_a.debug(L"m2_a");    
    return Error::handle(name(), L"mult", m2_a.ERR_DIM, __FILE__, __LINE__);
  }

  // get the type of two matrices and branch on type
  //
  Integral::MTYPE old_type = this_a.type_d;
  long nrows = m1_a.getNumRows();
  long ncols = m2_a.getNumColumns();
  Integral::MTYPE type1 = m1_a.getType();
  Integral::MTYPE type2 = m2_a.getType();

  // type1: FULL
  //
  if (type1 == Integral::FULL) {

    // type2: DIAGONAL
    //  set the type of output matrix as full, and multiply the elements in
    //  the first matrix with the diagonal elements in the second matrix which
    //  has the same column index
    //
    if (type2 == Integral::DIAGONAL) {
      
      // create space in the output matrix
      //
      this_a.setDimensions(nrows, ncols, false, Integral::FULL);
      
      // loop over elements and multiply
      //
      for (long row = 0; row < nrows; row++) {
        for (long col = 0; col < ncols; col++) {
          this_a.m_d(this_a.index(row, col)) =
	    (TIntegral)m1_a.m_d(m1_a.index(row, col)) *
            (TIntegral)m2_a.m_d(col);
        }
      }
    }

    // type2: FULL, SYMMETRIC, LOWER_TRIAGULAR, UPPER_TRIANGULAR
    //  set the type of output matrix as full, and multiply row of first
    //  matrix by the column of the second matrix - multiply elements in
    //  the given row of the first matrix with the elements in the given
    //  column of the second matrix
    //
    else if ((type2 == Integral::FULL) ||
	     (type2 == Integral::SYMMETRIC) ||
	     (type2 == Integral::LOWER_TRIANGULAR) ||
	     (type2 == Integral::UPPER_TRIANGULAR)) {

      // create space in the output matrix
      //
      this_a.setDimensions(nrows, ncols, false, Integral::FULL);
      
      // loop over all elements
      //
      for (long row = 0; row < nrows; row++) {
        for (long col = 0; col < ncols; col++) {  
          this_a.m_d(this_a.index(row, col)) = 
            this_a.multiplyRowByColumn(m1_a, m2_a, row, col);
        }
      }
    }

    // type2: SPARSE
    //  set the type of output matrix as full. 
    //
    else if (type2 == Integral::SPARSE) {

      // create space in the output matrix
      //
      this_a.setDimensions(nrows, ncols, false, Integral::FULL);
      this_a.clear(Integral::RETAIN);
       
      // declare local variables
      //
      long b_row;
      long b_col;
      TIntegral b_val;
      TIntegral c_val = 0;
      TIntegral sub_val = 0;
      long length = m2_a.m_d.length();

      // loop over the length of the vector of sparse matrix
      //
      for (long b_index = 0; b_index < length; b_index++) {

	// get row index, column index and value from sparse matrix
	//
        b_row = m2_a.row_index_d(b_index);
        b_col = m2_a.col_index_d(b_index);
        b_val = m2_a.m_d(b_index);

	// get the number of rows
	//
        long a_rows = m1_a.getNumRows();

	// loop through the number of rows and multiply
	//
        for (long a_index = 0; a_index < a_rows; a_index++) {
          c_val = this_a.getValue(a_index, b_col);
          sub_val = b_val * m1_a.getValue(a_index, b_row);
          this_a.setValue(a_index, b_col, c_val + sub_val);
        }
      }
    }
  }

  // type1: DIAGONAL
  //
  else if (type1 == Integral::DIAGONAL) {

    // type2: DIAGONAL
    //  if both matrices are diagonal ones, return a diagonal matrix and
    //  assign the products of the diagonal elements of the two input matrices
    //  to the output
    //
    if (type2 == Integral::DIAGONAL) {

      // create space in the output matrix
      //
      this_a.setDimensions(nrows, nrows, false, Integral::DIAGONAL);

      // loop over each diagonal element and multiply
      //
      for (long row = 0; row < nrows; row++) {
        this_a.m_d(row) = (TIntegral)m1_a.m_d(row) * (TIntegral)m2_a.m_d(row);
      }
    }

    // type2: FULL
    //  set the output matrix as full, and multiply the diagonal elements
    //  of the first matrix with the elements which has the same row index.
    //
    else if (type2 == Integral::FULL) {

      // create space in the output matrix
      //
      this_a.setDimensions(nrows, ncols, false, Integral::FULL);

      // loop over the elements and multiply
      //
      for (long row = 0; row < nrows; row++) {
        for (long col = 0; col < ncols; col++) {
          this_a.m_d(this_a.index(row, col)) =
	    (TIntegral)(m1_a.m_d(row)) *
	    (TIntegral)(m2_a.m_d(m2_a.index(row, col)));
        }
      }
    }

    // type2: SYMMETRIC
    //  set output matrix as symmetric, and multiply the diagonal elements
    //  of the first matrix to the symmetric matrix, for elements in the
    //  upper triangle of the symmetric matrix, use its symmetric value.
    //
    else if (type2 == Integral::SYMMETRIC) {

      // create space in the output matrix
      //
      this_a.setDimensions(nrows, ncols, false, Integral::FULL);
      
      // loop over the elements and multiply
      //      
      for (long row = 0; row < nrows; row++) {
	for (long col = 0; col < row; col++) {
          this_a.m_d(this_a.index(row, col)) =
	    (TIntegral)(m1_a.m_d(row)) *
	    (TIntegral)(m2_a.m_d(m2_a.index(row, col)));
        }
        for (long col = row; col < ncols; col++) {
          this_a.m_d(this_a.index(row, col)) =
	    (TIntegral)(m1_a.m_d(row)) *
	    (TIntegral)(m2_a.m_d(m2_a.index(col, row)));
        }
      }
    }

    // type2: LOWER_TRIANGULAR
    //  set the type of output matrix as lower triangular, and then 
    //  multiply the diagonal elements with the corresponding elements in
    //  the lower triangle.
    //
    else if (type2 == Integral::LOWER_TRIANGULAR) {

      // create space in the output matrix
      //
      this_a.setDimensions(nrows, nrows, false, Integral::LOWER_TRIANGULAR);
      
      // loop over the elements and multiply
      //      
      for (long row = 0; row < nrows; row++) {
        for (long col = 0; col <= row; col++) {
          this_a.m_d(this_a.index(row, col)) =
	    (TIntegral)(m1_a.m_d(row)) *
            (TIntegral)(m2_a.m_d(m2_a.index(row, col)));
        }
      }
    }

    // type2: UPPER_TRIANGULAR
    //  set the type of matrix as upper triangular, and multiply the diagonal
    //  elements with the corresponding elements in the upper triangle
    //
    else if (type2 == Integral::UPPER_TRIANGULAR) {

      // create space in the output matrix
      //
      this_a.setDimensions(nrows, nrows, false, Integral::UPPER_TRIANGULAR);
      
      // loop over the elements and multiply
      //      
      for (long row = 0; row < nrows; row++) {
        for (long col = row; col < nrows; col++) {
          this_a.m_d(this_a.index(row, col)) =
	    (TIntegral)(m1_a.m_d(row))*
            (TIntegral)(m2_a.m_d(m2_a.index(row, col)));
        }
      }
    }

    // type2: SPARSE
    //  set the type of output matrix as sparse. the general idea here is
    //  that we change sparse to full matrix, after finishing computation,
    //  we will change it to sparse matrix.
    //
    else if (type2 == Integral::SPARSE) {

      // create space in the output matrix
      //
      this_a.setDimensions(nrows, ncols, false, Integral::FULL);
      this_a.clear(Integral::RETAIN);
      
      // get the number of rows of diagonal matrix and the length of the
      // vector for sparse matrix
      //
      long len1 = m1_a.getNumRows();
      long len2 = m2_a.m_d.length();

      // declare local variables
      //
      long row2 = 0;
      long col2 = 0;
      TAIntegral val2 = 0;
      TIntegral val1 = 0;

      // loop over rows of the diagonal matrix
      //
      for (long row1 = 0; row1 < len1; row1++) {

	// get the element on the diagonal vector
	//
        val1 = m1_a.m_d(row1);

	// loop over the vector of sparse matrix
	//
        for (long index = 0; index < len2; index++) {

	  // get the row index, column index and value
	  //
          row2 = m2_a.row_index_d(index);
          col2 = m2_a.col_index_d(index);
          val2 = m2_a.m_d(index);
          
          // if the row index matches with position of diag vector
          // elements, multiply them
          //
          if (row2 == row1) {
            this_a.setValue(row1, col2, (TIntegral)(val1 * val2));
          }
        }
      }

      // now change type to sparse
      //
      this_a.changeType(Integral::SPARSE);
    }
  }

  // type1: SYMMETRIC
  //
  else if (type1 == Integral::SYMMETRIC) {

    // type2: DIAGONAL
    //  set the type of output matrix as full.
    //
    if (type2 ==  Integral::DIAGONAL) {

      // create space in the output matrix
      //
      this_a.setDimensions(nrows, ncols, false, Integral::FULL);
      
      // loop through the elements and multiply
      //
      for (long row = 0; row < nrows; row++) {
        for (long col = 0; col < row; col++) {
          this_a.m_d(this_a.index(row, col)) =
	    (TIntegral)(m2_a.m_d(col)) *
            (TIntegral)(m1_a.m_d(m1_a.index(row, col)));
        }
        for (long col = row; col < ncols; col++) {