mmat_08.cc

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

	   (this_a.type_d == Integral::SPARSE)) {

    // declare temporary vector
    //
    MVector<TScalar, TIntegral> temp;

    // do vector matrix multiplication and store the result into a
    // temporary vector
    //
    this_a.vmult(temp, v_in_a);

    // compute the dot product of the two vectors
    // note that dotProduct for complex numbers conjugates elements of
    // the second vector
    //
    output_a = v_in_a.dotProduct(temp);
  }

  // exit gracefully
  //  
  return true;
}

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


// method: quadratic
//
// arguments:
//  MMatrix<TScalar, TIntegral>& output: (output) matrix output
//  const MMatrix<TScalar, TIntegral>& mat1: (input) input matrix
//  const MMatrix<TScalar, TIntegral>& mat2: (input) input matrix
//
// return: a boolean value indicating status
//
// this method computes mat2 * mat1 * transpose(mat2)
// the result is stored in the output matrix. Note that mat1
// must be square.
//
//         [2 0 1]        [0 1 2] 
//  mat2 = [4 5 7] mat1 = [3 4 5] 
//                        [6 7 8]
//
//                         [0 1 2]   [2 4]             
//      output = [2 0 1] * [3 4 5] * [0 5] = [ 24  153]
//               [4 5 7]   [6 7 8]   [1 7]   [203 1216]
//
template<class TScalar, class TIntegral>
boolean
MMatrixMethods::quadratic(MMatrix<TScalar, TIntegral>& output_a,
			  const MMatrix<TScalar, TIntegral>& mat1_a,
			  const MMatrix<TScalar, TIntegral>& mat2_a) {

  // sanity checks
  //
  if (!mat1_a.isSquare() || (mat2_a.getNumColumns() != mat1_a.getNumRows())) {
    mat1_a.debug(L"mat1_a");
    mat2_a.debug(L"mat2_a");    
    return Error::handle(name(), L"quadratic", output_a.ERR_DIM,
			 __FILE__, __LINE__);
  }

  // make sure that none of the inputs is equal to the output
  //
  if (&mat1_a == &output_a) {
    MMatrix<TScalar, TIntegral> tmp(mat1_a);
    return quadratic(output_a, tmp, mat2_a) ;
  }
  else if (&output_a == &mat2_a) {
    MMatrix<TScalar, TIntegral> tmp(mat2_a);
    return quadratic(output_a, mat1_a, tmp);
  }
  
  // for the moment, unless the center matrix is diagonal, we do a simple
  // (but resource-intensive for larger matrices) temporary matrix computation
  //
  if (mat1_a.getType() != Integral::DIAGONAL) {

    MMatrix<TScalar, TIntegral> tmp_matrix;

    // compute the first multiplication: mat1 * mat2
    // and then the second: (mat1 * mat2) * transpose(mat1)
    //
    return (tmp_matrix.mult(mat2_a, mat1_a) &&
	    output_a.outerProduct(tmp_matrix, mat2_a));
  }

  // if the center matrix is diagonal then the output will be
  // symmetric and we need only calculate half of the output.
  // Further, there is a simple loop structure to the calculation
  //
  else {

    Integral::MTYPE old_type  = output_a.getType();
    long nrows = mat2_a.getNumRows();
    long ncols = mat2_a.getNumColumns();

    // if the output matrix is symmetric then we need not make any
    // change to the type. otherwise we will set it to be full temporarily
    // and determine later if the type can be changed back. Only
    // a symmetric, full, sparse, or (perhaps) diagonal are appropriate
    //
    if (output_a.getType() != Integral::SYMMETRIC) {
      output_a.setDimensions(nrows, nrows, false, Integral::FULL);
    }
    else {
      output_a.setDimensions(nrows, nrows, false);
    }
    
    // now loop and compute the output
    //
    for (long row = 0; row < nrows; row++) {
      for (long col = 0; col <= row; col++) {
	TIntegral sum = (TIntegral)0;
	for (long row2 = 0; row2 < ncols; row2++) {
	  sum += mat2_a(row,row2) * mat2_a(col,row2) * mat1_a.m_d(row2);
	}
	output_a.setValue(row,col,sum);
	output_a.setValue(col,row,sum);
      }
    }

    // try to change the output type to the desired type
    // restore the type if possible
    //
    if (output_a.isTypePossible(old_type)) {
      return output_a.changeType(old_type);
    }
    else {
      return output_a.changeType(Integral::FULL);
    }
  }
}

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


// method: outerProduct
//
// 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 computes the outer product of two matrices. The
// column dimension of both matrices must agree.
// for example:
//
//       [1 4]       [2 0]                               [2 24  4 42]
//  m1 = [2 5]  m2 = [0 6]  m_out = m1 * transpose(m2) = [4 30  8 57]
//       [3 6]       [4 0]                               [6 36 12 72]
//                   [6 9]
//
template<class TScalar, class TIntegral, class TAScalar, class TAIntegral>
boolean MMatrixMethods::outerProduct(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.outerProduct(tmp, m2_a) ;
  }
  else if ((void*)&this_a == (void*)&m2_a) {
    MMatrix<TScalar, TIntegral> tmp(this_a);
    return this_a.outerProduct(m1_a, tmp);
  }

  // check dimensions
  //
  if (m1_a.getNumColumns() != m2_a.getNumColumns()) {
    m1_a.debug(L"m1_a");
    m2_a.debug(L"m2_a");    
    return Error::handle(name(), L"outerProduct", 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.getNumRows();
  Integral::MTYPE type1 = m1_a.getType();
  Integral::MTYPE type2 = m2_a.getType();

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

    // create space in the output matrix
    //
    this_a.setDimensions(nrows, ncols, false, 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) {
      
      // 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)((TAIntegral)m1_a.m_d(m1_a.index(row, col)) *
	       (TAIntegral)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)) {

      // 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)) = 
            (TIntegral)this_a.multiplyRowByRow(m1_a, m2_a, row, col);
        }
      }
    }

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

      // clear the output matrix
      //
      this_a.clear(Integral::RETAIN);
       
      // declare local variables
      //
      long b_row;
      long b_col;
      TAIntegral 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_row);
          sub_val = (TIntegral)((TAIntegral)b_val *
				m1_a.getValue(a_index, b_col));
          this_a.setValue(a_index, b_row, 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)((TAIntegral)m1_a.m_d(row) * (TAIntegral)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.
    //  since a diagonal matrix is equal to its transpose, we can do normal
    //  matrix multiplication here
    //
    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)((TAIntegral)(m1_a.m_d(row)) *
	    (TAIntegral)(m2_a.m_d(m2_a.index(col, row))));
        }
      }
    }

    // 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)((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)((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::UPPER_TRIANGULAR);
      
      // loop over the elements and multiply
      //      
      for (long row1 = 0; row1 < nrows; row1++) {
        for (long row2 = row1; row2 < nrows; row2++) {
          this_a.m_d(this_a.index(row1, row2)) =
	    (TIntegral)((TAIntegral)(m1_a.m_d(row1)) *
            (TAIntegral)(m2_a.m_d(m2_a.index(row2, row1))));
        }
      }
    }

    // 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::LOWER_TRIANGULAR);
      
      // loop over the elements and multiply
      //
      for (long row1 = 0; row1 < nrows; row1++) {
        for (long row2 = 0; row2 <= row1; row2++) {
          this_a.m_d(this_a.index(row1, row2)) =
	    (TIntegral)((TAIntegral)(m1_a.m_d(row1))*
            (TAIntegral)(m2_a.m_d(m2_a.index(row2, row1))));
        }
      }
    }

    // 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.getNumColumns();
      long len2 = m2_a.m_d.length();

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

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

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

	// 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 (col2 == col1) {
            this_a.setValue(col1, row2, (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)) =
	    (TScalar)((TIntegral)(m2_a.m_d(col)) *
		      (TIntegral)(m1_a.m_d(m1_a.index(row, col))));