mmat_09.cc

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

// file: $isip/class/math/matrix/MMatrix/mmat_09.cc
// version: $Id: mmat_09.cc,v 1.41 2002/03/15 23:29:58 zheng Exp $
//

// isip include files
//
#include "MMatrixMethods.h"
#include "MMatrix.h"

// method: add
//
// 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::add(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 add the value to the storage vector
  //
  if (!this_a.m_d.add(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::add<ISIP_TEMPLATE_TARGET>
(MMatrix<ISIP_TEMPLATE_TARGET>&, ISIP_TEMPLATE_T1);

// method: add
//
// arguments:
//  MMatrix<TScalar, TIntegral>& this: (output) class operand
//  const MMatrix<TAScalar, TAIntegral>& m2: (input) operand 1
//
// return: a boolean value indicating status
//
// this method adds a matrix to "this". we use the following flowchart:
//
// (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
//
template<class TScalar, class TIntegral, class TAScalar, class TAIntegral>
boolean MMatrixMethods::add(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"add", 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.add(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), then
  // add 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
  // add 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) {
      this_a.changeType(Integral::SPARSE);
    }
    else {
      this_a.changeType(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)));
    }

    // restore the type
    //
    return this_a.changeType(type1);
  }
    
  // if either type is FULL, we assign the output type to FULL,
  // then add elements of the second
  //
  else if ((type1 == Integral::FULL) || (type2 == Integral::FULL)) {
    this_a.changeType(Integral::FULL);
  }

  //---------------------------------------------------------------------------
  // 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.
  // 
  else if (type1 == Integral::DIAGONAL) {
    this_a.changeType(type2);
  }
    
  // type1: SYMMETRIC, *_TRIANGULAR
  //  the remaining types are DIAGONAL, SYMMETRIC, and *_TRIANGULAR. 
  //  the output types vary with the input.
  // 
  else {
    
    // set up the output matrix: a diagonal matrix preserves the type
    // of the input; all other types force a full matrix.
    //
    if (type2 != Integral::DIAGONAL) {
      this_a.changeType(Integral::FULL);
    }
  }

  // now, we can add elements of the second matrix to "this"
  //
  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.m_d(this_a.index(row, col)) += 
	(TIntegral)m2_a.m_d(m2_a.index(row, col));
    }
  }

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

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

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

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

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

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

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

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

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

template boolean
MMatrixMethods::add<ISIP_TEMPLATE_TARGET, Double, double>
(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<ComplexLong, complexlong>&);

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

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

// method: add
//
// 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
//
// add the two argument matrices and store result in the current matrix
// for example:
//
//       [1 2 3]       [2 0 4]                    [3  2  7]
//  m1 = [4 5 6]  m2 = [0 6 0]  m_out = m1 + m2 = [4  11 6]
//       [7 8 9]       [8 0 1]                    [15 8 10]
//
// this method has evolved somewhat as our ability to efficiently
// copy and assign matrices has improved. this version implements
// some complicated logic to do things as efficiently as possible.
// 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 c = 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
//
// addition has an advantage over subtraction in that we can interchange
// the matrices such that the least number of operations are performed.
//
template<class TScalar, class TIntegral, class TAScalar, class TAIntegral>
boolean MMatrixMethods::add(MMatrix<TScalar, TIntegral>& this_a, 
			    const MMatrix<TScalar, TIntegral>& m1_a, 
			    const MMatrix<TAScalar, TAIntegral>& m2_a) {

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

  // declare local variables
  //
  long nrows = m1_a.nrows_d;
  long ncols = m2_a.ncols_d;
  Integral::MTYPE old_type = this_a.type_d;
  Integral::MTYPE type1 = m1_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::DIAGONAL) ||
	(type1 == Integral::SYMMETRIC) ||
	(type1 == Integral::LOWER_TRIANGULAR) ||
	(type1 == Integral::UPPER_TRIANGULAR)) {
      this_a.setDimensions(nrows, ncols, false, m1_a.type_d);
      this_a.m_d.add(m1_a.m_d, m2_a.m_d);
    }

    // type: SPARSE
    //
    else {

      // set the output matrix to SPARSE
      //
      this_a.assign(m1_a, Integral::SPARSE);

      // 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)(m1_a(row, col) + m2_a.m_d(index)));
      }
    }
    
    // restore the type and return
    //
    return this_a.changeType(old_type);
  }

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

    // create a FULL matrix for all types except DIAGONAL, and zero it
    //
    if (type2 == Integral::DIAGONAL) {
      this_a.setDimensions(nrows, ncols, false, Integral::SPARSE);