mmat_08.cc

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

//                [7]
//
// if length(m1) = r and length(m2) = c, then m_out will be a r by c matrix.
//
template<class TScalar, class TIntegral, class TAScalar, class TAIntegral>
boolean
MMatrixMethods::outerProduct(MMatrix<TScalar, TIntegral>& this_a, 
			     const MVector<TScalar, TIntegral>& v1_a, 
			     const MVector<TAScalar, TAIntegral>& v2_a) {
  
  // determine the new dimensions
  //
  long nrows = v1_a.length();
  long ncols = v2_a.length();

  // get the current type of the output matrix. if it is SPARSE then leave it
  // as SPARSE, otherwise change it to FULL
  //
  Integral::MTYPE output_type = this_a.type_d;

  if (output_type != Integral::SPARSE) {
    output_type = Integral::FULL;
  }

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

  // loop over all elements of both vectors and create the outer
  // product matrix
  //
  for (long i = 0; i < nrows; i++) {
    for (long j = 0; j < ncols; j++) {
      this_a.setValue(i, j, (TIntegral)(v1_a(i) * (TIntegral)v2_a(j)));
    }
  }

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

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

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

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

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

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

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

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

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

template boolean
MMatrixMethods::outerProduct<ISIP_TEMPLATE_TARGET, Double, double>
(MMatrix<ISIP_TEMPLATE_TARGET>&, const MVector<ISIP_TEMPLATE_TARGET>&,
 const MVector<Double, double>&);
#endif
#ifdef ISIP_TEMPLATE_complex
template boolean
MMatrixMethods::outerProduct<ISIP_TEMPLATE_TARGET, ComplexLong, complexlong>
(MMatrix<ISIP_TEMPLATE_TARGET>&, const MVector<ISIP_TEMPLATE_TARGET>&,
 const MVector<ComplexLong, complexlong>&);

template boolean
MMatrixMethods::outerProduct<ISIP_TEMPLATE_TARGET, ComplexFloat, complexfloat>
(MMatrix<ISIP_TEMPLATE_TARGET>&, const MVector<ISIP_TEMPLATE_TARGET>&,
 const MVector<ComplexFloat, complexfloat>&);
 
template boolean
MMatrixMethods::outerProduct<ISIP_TEMPLATE_TARGET, ComplexDouble, complexdouble>
(MMatrix<ISIP_TEMPLATE_TARGET>&, const MVector<ISIP_TEMPLATE_TARGET>&,
 const MVector<ComplexDouble, complexdouble>&);
#endif

// method: multiplyRowByRow
//
// arguments:
//  const MMatrix<TScalar, TIntegral>& this: (output) class operand
//  const MMatrix<TScalar, TIntegral>& m1: (input) operand 1
//  const MMatrix<TAScalar, TAIntegral>& m2: (input) operand 2
//  long row_m1: (input) row index in m1
//  long row_m2: (input) row index in m2
//
// return: a TIntegral value that results from multiplying the specified row
// in m1 by the specified row in m2
//
// this method is designed for the outer-product multiplication of full,
// symmetric, lower and upper triangular matrices
//
template<class TScalar, class TIntegral, class TAScalar, class TAIntegral>
TIntegral
MMatrixMethods::multiplyRowByRow(const MMatrix<TScalar, TIntegral>& this_a, 
				 const MMatrix<TScalar, TIntegral>& m1_a, 
				 const MMatrix<TAScalar, TAIntegral>& m2_a,
				 long row_m1_a, long row_m2_a) {

  // declare local variables
  //
  TIntegral sum = 0;
  
  // get the number of columns of the first matrix
  // 
  long last_index = m1_a.getNumRows();
  long ncols = m1_a.getNumColumns();

  // get the type of the two matrices and branch on type
  //
  Integral::MTYPE type1 = m1_a.getType();
  Integral::MTYPE type2 = m2_a.getType();

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

    // type2: FULL
    //
    if (type2 == Integral::FULL) {
      
      // compute the dot product
      //  
      for (long index = 0; index < ncols; index++) {
        sum += (TIntegral)((TIntegral)m1_a.m_d(m1_a.index(row_m1_a, index)) *
			  (TAIntegral)m2_a.m_d(m2_a.index(row_m2_a, index)));
      }      
    }

    // type2: SYMMETRIC
    //  multiply each element in the given column of the first matrix
    //  with the elements in the corresponding column of the second matrix
    //
    else if (type2 == Integral::SYMMETRIC) {
      
      // loop through the upper triangular elements of the second matrix
      //
      for (long index = 0; index < row_m2_a; index++) {
        sum += (TIntegral)((TIntegral)m1_a.m_d(m1_a.index(row_m1_a, index)) *
			   (TAIntegral)m2_a.m_d(m2_a.index(index, row_m2_a)));
      }

      // loop through the lower triangular elements of the second matrix
      //
      for (long index = row_m2_a; index < ncols; index++) {
        sum += (TIntegral)((TIntegral)m1_a.m_d(m1_a.index(row_m1_a, index)) *
			   (TAIntegral)m2_a.m_d(m2_a.index(row_m2_a, index)));
      }
    }

    // type2: LOWER_TRIANGULAR
    //  multiply the lower triangular elements in the given column of
    //  the second matrix with the corresponding elements in the
    //  given column of the first matrix
    //
    else if (type2 == Integral::LOWER_TRIANGULAR) {
      for (long index = 0; index <= row_m2_a; index++) {
        sum += (TIntegral)((TIntegral)m1_a.m_d(m1_a.index(row_m1_a, index)) *
			   (TAIntegral)m2_a.m_d(m2_a.index(row_m2_a, index)));
      }    
    }

    // type2: UPPER_TRIANGULAR
    //  multiply the upper triangular elements in the given column of
    //  the second matrix with the corresponding elements in the
    //  given column of the first matrix
    //      
    else if (type2 == Integral::UPPER_TRIANGULAR) {
      for (long index = row_m2_a; index < ncols; index++) {
        sum += (TIntegral)((TIntegral)m1_a.m_d(m1_a.index(row_m1_a, index)) *
			   (TAIntegral)m2_a.m_d(m2_a.index(row_m2_a, index)));
      }    
    }
  }

  // type1: SYMMETRIC
  //  a symmetric matrix is equal to its transpose so we can follow normal
  //  multiplication rules here
  //
  else if (type1 == Integral::SYMMETRIC) {

    // type2: FULL
    //  multiply each elements in the given row of the first matrix with
    //  the elements in the corresponding column of the second matrix
    // 
    if (type2 == Integral::FULL) {  
      for (long index = 0; index < row_m1_a; index++) {
        sum += (TIntegral)((TIntegral)m1_a.m_d(m1_a.index(index, row_m1_a)) *
			   (TAIntegral)m2_a.m_d(m2_a.index(row_m2_a, index)));
      }      
      for (long index = row_m1_a; index < last_index; index++) {
        sum += (TIntegral)((TIntegral)m1_a.m_d(m1_a.index(row_m1_a, index)) *
			   (TAIntegral)m2_a.m_d(m2_a.index(row_m2_a, index)));
      }    
    }

    // type2: SYMMETRIC
    //  multiply each elements in the given row of the first matrix with
    //  the elements in the corresponding column of the second matrix.
    //  the code needs to be separated into several loops to make sure when
    //  a upper triangle element is needed, its symmetric value will be used
    // 
    else if (type2 == Integral::SYMMETRIC) {

      if (row_m1_a > row_m2_a) {
	for (long index = 0; index < row_m2_a; index++) {
          sum += (TIntegral)((TIntegral)m1_a.m_d(m1_a.index(index, row_m1_a)) *
			     (TAIntegral)m2_a.m_d(m2_a.index(index, row_m2_a)));
        }
        for (long index = row_m2_a; index < row_m1_a; index++) {
          sum += (TIntegral)((TIntegral)m1_a.m_d(m1_a.index(index, row_m1_a)) *
			     (TAIntegral)m2_a.m_d(m2_a.index(row_m2_a, index)));
        }       
        for (long index = row_m1_a; index < last_index; index++) {
          sum += (TIntegral)((TIntegral)m1_a.m_d(m1_a.index(row_m1_a, index)) *
			     (TAIntegral)m2_a.m_d(m2_a.index(row_m2_a, index)));
        }
      }

      else {
        for (long index = 0; index < row_m1_a; index++) {
          sum += (TIntegral)((TIntegral)m1_a.m_d(m1_a.index(index, row_m1_a)) *
			     (TAIntegral)m2_a.m_d(m2_a.index(index, row_m2_a)));
        }
        for (long index = row_m1_a; index < row_m2_a; index++) {
          sum += (TIntegral)((TIntegral)m1_a.m_d(m1_a.index(row_m1_a, index)) *
			     (TAIntegral)m2_a.m_d(m2_a.index(index, row_m2_a)));
        }               
        for (long index = row_m2_a; index < last_index; index++) {
          sum += (TIntegral)((TIntegral)m1_a.m_d(m1_a.index(row_m1_a, index)) *
			     (TAIntegral)m2_a.m_d(m2_a.index(row_m2_a, index)));
        }
      }
    }

    // type2: LOWER_TRIANGULAR
    //  multiply the elements in the lower triangle of the second matrix with
    //  its corresponding elements in the first matrix
    //
    else if (type2 == Integral::LOWER_TRIANGULAR) {

      for (long index = 0; index <= row_m2_a; index++) {
	sum += (TIntegral)((TIntegral)m1_a.m_d(m1_a.index(row_m1_a, index)) *
			   (TAIntegral)m2_a.m_d(m2_a.index(row_m2_a, index)));
      }
    }

    // type2: UPPER_TRIANGULAR
    //  multiply the elements in the upper triangle of the second matrix with
    //  its corresponding elements in the first matrix
    //
    else if (type2 == Integral::UPPER_TRIANGULAR) {

      for (long index = row_m2_a; index < ncols; index++) {
	sum += (TIntegral)((TIntegral)m1_a.m_d(m1_a.index(row_m1_a, index)) *
			   (TAIntegral)m2_a.m_d(m2_a.index(row_m2_a, index)));
      }
    }
  }

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

    // type2: FULL
    //  multiply each element in the given column (and in the lower
    //  triangle) of the first matrix with its corresponding one in
    //  the given column of the second matrix, return the summation
    //  of the products
    //        
    if (type2 == Integral::FULL) {
      for (long index = 0; index <= row_m1_a; index++) {
        sum += (TIntegral)((TIntegral)m1_a.m_d(m1_a.index(row_m1_a, index)) *
			   (TAIntegral)m2_a.m_d(m2_a.index(row_m2_a, index)));
      }      
    }

    // type2: SYMMETRIC
    //  multiply each elements in the given column (and in the lower
    //  triangle) of the first matrix with its corresponding one (if
    //  it's in the upper triangle, use its symmetric value) in the
    //  given column of the second matrix, return the summation of
    //  the products
    //
    else if (type2 == Integral::SYMMETRIC) {

      for (long index = 0; index <= row_m1_a; index++) {
	sum += (TIntegral)((TIntegral)m1_a.m_d(m1_a.index(row_m1_a, index)) *
			   (TAIntegral)m2_a.m_d(m2_a.index(row_m2_a, index)));
      }
    }

    // type2: LOWER_TRIANGULAR
    //  multiply the elements in lower triangle
    //
    else if (type2 == Integral::LOWER_TRIANGULAR) {

      // find the start index
      //
      long last_index = (row_m1_a > row_m2_a) ? row_m2_a : row_m1_a;
      for (long index = 0; index <= last_index; index++) {
	sum += (TIntegral)((TIntegral)m1_a.m_d(m1_a.index(row_m1_a, index)) *
			   (TAIntegral)m2_a.m_d(m2_a.index(row_m2_a, index)));
      }
    }

    // type: UPPER_TRIANGULAR
    //
    else if (type2 == Integral::UPPER_TRIANGULAR) {

      long start_index = row_m2_a;
      last_index = row_m1_a;
      for (long index = start_index; index <= last_index; index++) {
	sum += (TIntegral)((TIntegral)m1_a.m_d(m1_a.index(row_m1_a, index)) *
			   (TAIntegral)m2_a.m_d(m2_a.index(row_m2_a, index)));
      }
    }
  }

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

    // type2: FULL
    //  multiply each element in the given column (and in the upper triangle)
    //  of the first matrix with its corresponding one in the given column
    //  of the second matrix, return the summation of the products
    //
    if (type2 == Integral::FULL) {
      for (long index = row_m1_a; index < ncols; index++) {
        sum += (TIntegral)((TIntegral)m1_a.m_d(m1_a.index(row_m1_a, index)) *
			   (TAIntegral)m2_a.m_d(m2_a.index(row_m2_a, index)));
      }      
    }

    // type2: SYMMETRIC
    //  multiply each element in the given column (and in the lower triangle)
    //  of the first matrix with its corresponding one (if it's in the upper
    //  triangle, use its symmetric value) in the given column
    //  of the second matrix, return the summation of the products
    // 
    else if (type2 == Integral::SYMMETRIC) {

      for (long index = row_m1_a; index < ncols; index++) {
	sum += (TIntegral)((TIntegral)m1_a.m_d(m1_a.index(row_m1_a, index)) *
			   (TAIntegral)m2_a.m_d(m2_a.index(row_m2_a, index)));
      }
    }

    // type2: LOWER_TRIANGULAR
    //
    else if (type2 == Integral::LOWER_TRIANGULAR) {

      // find the start index
      //
      long start_index = row_m1_a;
      last_index = row_m2_a;
      for (long index = start_index; index <= last_index; index++) {
	sum += (TIntegral)((TIntegral)m1_a.m_d(m1_a.index(row_m1_a, index)) *
			   (TAIntegral)m2_a.m_d(m2_a.index(row_m2_a, index)));
      }
    }

    // type2: UPPER_TRIANGULAR
    //
    else if (type2 == Integral::UPPER_TRIANGULAR) {

      // multiply upper triangular elements
      //
      long start_index = (row_m1_a < row_m2_a) ? row_m2_a : row_m1_a;
      for (long index = start_index; index < ncols; index++) {
	sum += (TIntegral)((TIntegral)m1_a.m_d(m1_a.index(row_m1_a, index)) *
			   (TAIntegral)m2_a.m_d(m2_a.index(row_m2_a, index)));
      }
    }
  }

  // exit gracefully
  //
  return sum;
}

// explicit instantiations
//
// we only allow complex matrices to operate with other complex matrices
//
#ifndef ISIP_TEMPLATE_complex
template ISIP_TEMPLATE_T1
MMatrixMethods::multiplyRowByRow<ISIP_TEMPLATE_TARGET, Byte, byte>
(const MMatrix<ISIP_TEMPLATE_TARGET>&, const MMatrix<ISIP_TEMPLATE_TARGET>&,
 const MMatrix<Byte, byte>&, long, long);

template ISIP_TEMPLATE_T1
MMatrixMethods::multiplyRowByRow<ISIP_TEMPLATE_TARGET, Ushort, ushort>
(const MMatrix<ISIP_TEMPLATE_TARGET>&, const MMatrix<ISIP_TEMPLATE_TARGET>&,
 const MMatrix<Ushort, ushort>&, long, long);

template ISIP_TEMPLATE_T1
MMatrixMethods::multiplyRowByRow<ISIP_TEMPLATE_TARGET, Ulong, ulong>
(const MMatrix<ISIP_TEMPLATE_TARGET>&, const MMatrix<ISIP_TEMPLATE_TARGET>&,
 const MMatrix<Ulong, ulong>&, long, long);

template ISIP_TEMPLATE_T1
MMatrixMethods::multiplyRowByRow<ISIP_TEMPLATE_TARGET, Ullong, ullong>
(const MMatrix<ISIP_TEMPLATE_TARGET>&, const MMatrix<ISIP_TEMPLATE_TARGET>&,
 const MMatrix<Ullong, ullong>&, long, long);

template ISIP_TEMPLATE_T1
MMatrixMethods::multiplyRowByRow<ISIP_TEMPLATE_TARGET, Short, short>
(const MMatrix<ISIP_TEMPLATE_TARGET>&, const MMatrix<ISIP_TEMPLATE_TARGET>&,
 const MMatrix<Short, short>&, long, long);

template ISIP_TEMPLATE_T1
MMatrixMethods::multiplyRowByRow<ISIP_TEMPLATE_TARGET, Long, long>
(const MMatrix<ISIP_TEMPLATE_TARGET>&, const MMatrix<ISIP_TEMPLATE_TARGET>&,
 const MMatrix<Long, long>&, long, long);

template ISIP_TEMPLATE_T1
MMatrixMethods::multiplyRowByRow<ISIP_TEMPLATE_TARGET, Llong, llong>
(const MMatrix<ISIP_TEMPLATE_TARGET>&, const MMatrix<ISIP_TEMPLATE_TARGET>&,
 const MMatrix<Llong, llong>&, long, long);

template ISIP_TEMPLATE_T1
MMatrixMethods::multiplyRowByRow<ISIP_TEMPLATE_TARGET, Float, float>
(const MMatrix<ISIP_TEMPLATE_TARGET>&, const MMatrix<ISIP_TEMPLATE_TARGET>&,
 const MMatrix<Float, float>&, long, long);

template ISIP_TEMPLATE_T1
MMatrixMethods::multiplyRowByRow<ISIP_TEMPLATE_TARGET, Double, double>
(const MMatrix<ISIP_TEMPLATE_TARGET>&, const MMatrix<ISIP_TEMPLATE_TARGET>&,
 const MMatrix<Double, double>&, long, long);
#endif

#ifdef ISIP_TEMPLATE_complex
template ISIP_TEMPLATE_T1
MMatrixMethods::multiplyRowByRow<ISIP_TEMPLATE_TARGET, ComplexLong, complexlong>
(const MMatrix<ISIP_TEMPLATE_TARGET>&, const MMatrix<ISIP_TEMPLATE_TARGET>&,
 const MMatrix<ComplexLong, complexlong>&, long, long);

template ISIP_TEMPLATE_T1
MMatrixMethods::multiplyRowByRow<ISIP_TEMPLATE_TARGET, ComplexFloat, complexfloat>
(const MMatrix<ISIP_TEMPLATE_TARGET>&, const MMatrix<ISIP_TEMPLATE_TARGET>&,
 const MMatrix<ComplexFloat, complexfloat>&, long, long);

template ISIP_TEMPLATE_T1
MMatrixMethods::multiplyRowByRow<ISIP_TEMPLATE_TARGET, ComplexDouble, complexdouble>
(const MMatrix<ISIP_TEMPLATE_TARGET>&, const MMatrix<ISIP_TEMPLATE_TARGET>&,
 const MMatrix<ComplexDouble, complexdouble>&, long, long);
#endif