mmat_08.cc

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

// file: $isip/class/math/matrix/MMatrix/mmat_08.cc
// version: $Id: mmat_08.cc,v 1.31 2002/07/29 14:53:18 hamaker Exp $
//

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

// method: multv
//
// arguments:
//  const MMatrix<TScalar, TIntegral>& this: (output) class operand
//  MVector<TScalar, TIntegral>& v_out: (output) output vector
//  const MVector<TScalar, TIntegral>& v_in: (input) input vector
//
// return: a boolean value indicating status
//
// this method multiplies the current matrix with the input vector and
// assign result to the output vector, for example:
//
//           [0 1 2]         [2]          [0 1 2]   [2]   [2 ]
//  matrix = [3 4 5]  v_in = [0]  v_out = [3 4 5] * [0] = [11]
//           [6 7 8]	     [1]   	  [6 7 8]   [1]   [20]
//
// for this multiplication, the vector should be a column vector and
// the dimensions should agree
//
template<class TScalar, class TIntegral>
boolean MMatrixMethods::multv(const MMatrix<TScalar, TIntegral>& this_a, 
			      MVector<TScalar, TIntegral>& v_out_a, 
			      const MVector<TScalar, TIntegral>& v_in_a) {

  // if the input and output are the same, then a temporary copy is needed
  //
  if (&v_in_a == &v_out_a) {
    MVector<TScalar, TIntegral> tmp(v_in_a);
    return multv(this_a, v_out_a, tmp);
  }
  
  // check the arguments:
  //  the columns of the matrix should be equal to the length of input vector
  //
  if (this_a.getNumColumns() != v_in_a.length()) {
    v_in_a.debug(L"v_in_a");    
    this_a.debug(L"this_a");    
    return Error::handle(name(), L"multv", Error::ARG, __FILE__, __LINE__);
  }

  // get the number of rows and columns of the current matrix
  //
  long nrows = (long)this_a.nrows_d;
  long ncols = (long)this_a.ncols_d;

  // create space in the output vector
  // 
  if (!v_out_a.setLength(nrows, false)) {
    v_out_a.debug(L"v_out_a");
    this_a.debug(L"this_a");    
    return Error::handle(name(), L"multv", Error::MEM, __FILE__, __LINE__);
  }

  // type: FULL
  //  the general idea is that we loop through each element of matrix and get
  //  the dot product
  //
  if (this_a.type_d == Integral::FULL) {

    // loop over all elements
    //
    for (long row = 0, index = 0; row < nrows; row++) {

      // initialize the accumulator for the dot product
      //
      TIntegral sum = 0;

      // compute the dot product
      //
      for (long col = 0; col < ncols; col++) {
	sum += this_a.m_d(index++) * v_in_a(col);
      }
      
      // save the result
      //
      v_out_a(row) = sum;
    }
  }
  
  // type: DIAGONAL
  //  the general idea here is that we just multiply corresponding elements of
  //  diagonal and input vector
  //
  else if (this_a.type_d == Integral::DIAGONAL) {

    // only multiply the diagonal values to the corresponding elements in the
    // input vector
    //
    for (long row = 0; row < nrows; row++) {
      v_out_a(row) = this_a.m_d(row) * v_in_a(row);
    }
  }

  // type: SYMMETRIC
  //  the general idea is that we loop through each row of the matrix.
  //  for elements in the upper triangular part, we use the values from
  //  the lower triangular part.
  //
  else if (this_a.type_d == Integral::SYMMETRIC) {

    // loop over rows of the matrix
    //
    for (long row = 0; row < nrows; row++) {

      // initialize the sum of dot products for current row
      //      
      TIntegral sum = 0;

      // loop over the lower triangular elements of the matrix
      //      
      for (long col = 0; col < row; col++) {
        sum += this_a.m_d(this_a.index(row, col)) * v_in_a(col);
      }

      // loop over the upper triangular elements of the matrix
      //            
      for (long col = row; col < nrows; col++) {
        sum += this_a.m_d(this_a.index(col, row)) * v_in_a(col);
      }

      // save the result
      //
      v_out_a(row) = sum;
    }
  }
  
  // type: LOWER_TRIANGULAR
  //  the general idea is that we only loop through the elements
  //  in the lower triangular part of the matrix
  //
  else if (this_a.type_d == Integral::LOWER_TRIANGULAR) {

    // loop over rows of the matrix
    //
    for (long row = 0; row < nrows; row++) {

      // initialize the sum of dot products for current row
      //      
      TIntegral sum = 0;

      // loop over the lower triangular elements of the matrix
      //      
      for (long col = 0; col <= row; col++) {
        sum += this_a.m_d(this_a.index(row, col)) * v_in_a(col);
      }

      // save the result
      //
      v_out_a(row) = sum;
    }
  }
  
  // type: UPPER_TRIANGULAR
  //
  else if (this_a.type_d == Integral::UPPER_TRIANGULAR) {

    // loop over rows of the matrix
    //
    for (long row = 0; row < nrows; row++) {

      // initialize the sum of dot products for current row
      //      
      TIntegral sum = 0;

      // loop over the lower triangular elements of the matrix
      //      
      for (long col = row; col < ncols; col++) {
        sum += this_a.m_d(this_a.index(row, col)) * v_in_a(col);
      }

      // save the result
      //
      v_out_a(row) = sum;
    }
  }
  
  // type: SPARSE
  //  only use the non-zero elements
  //
  else if (this_a.type_d == Integral::SPARSE) {

    // clear the output vector since we will use it as an accumulator
    //
    v_out_a.clear(Integral::RETAIN);
    
    // loop through non-zero elements of matrix, multiply and store in
    // output vector
    //
    long last_index = this_a.m_d.length();

    for (long index = 0; index < last_index; index++) {
      long row = this_a.row_index_d(index);
      long col = this_a.col_index_d(index);
      v_out_a(row) += this_a.m_d(index) * v_in_a(col);
    }
  }

  // exit gracefully
  //  
  return true;
}

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

// method: vmult
//
// arguments:
//  const MMatrix<TScalar, TIntegral>& this: (output) class operand
//  MVector<TScalar, TIntegral>& v_out: (output) output vector
//  const MVector<TScalar, TIntegral>& v_in: (input) input vector
//
// return: a boolean value indicating status
//
// this method multiply the input vector with the current matrix and assign
// result to the output vector, for example:
//
//                           [0 1 2]                    [0 1 2]   
//  v_in = [2 0 1]  matrix = [3 4 5]  v_out = [2 0 1] * [3 4 5] = [6 9 12]
//                           [6 7 8]	        	[6 7 8]   
//
// for this multiplication, the vector should be a row vector and
// dimensions should agree
//
template<class TScalar, class TIntegral>
boolean MMatrixMethods::vmult(const MMatrix<TScalar, TIntegral>& this_a, 
			      MVector<TScalar, TIntegral>& v_out_a, 
			      const MVector<TScalar, TIntegral>& v_in_a) {

  // if the input and output are the same, then a temporary copy is needed
  //
  if (&v_in_a == &v_out_a) {
    MVector<TScalar, TIntegral> tmp(v_in_a);
    return multv(this_a, v_out_a, tmp);
  }
  
  // get the number of rows and columns of the matrix
  //
  long nrows = this_a.getNumRows();
  long ncols = this_a.getNumColumns();
  
  // check the dimensions
  //  the row of matrix should be equal to the length of input vector
  //
  if (nrows != v_in_a.length()) {
    v_in_a.debug(L"v_in_a");    
    this_a.debug(L"this_a");    
    return Error::handle(name(), L"vmult()", Error::ARG, __FILE__, __LINE__);
  }

  // create space in the output vector
  //
  if (!v_out_a.setLength(ncols, false)) {
    v_out_a.debug(L"v_out_a");    
    this_a.debug(L"this_a");    
    return Error::handle(name(), L"vmult", Error::MEM, __FILE__, __LINE__);
  }

  // type: FULL
  //  compute the dot product directly
  //
  if (this_a.type_d == Integral::FULL) {

    // loop over the number of columns
    //
    for (long col = 0; col < ncols; col++) {

      // initialize the sum of dot products for current column
      //            
      TIntegral sum = 0;

      // calculate the dot product of the row
      //
      for (long row = 0; row < nrows; row++) {
	sum += this_a.m_d(this_a.index(row, col)) * v_in_a(row);
      }
      
      // save the result
      //
      v_out_a(col) = sum;
    }
  }

  // type: DIAGONAL
  //  multiply corresponding elements of diagonal and input vector
  //
  else if (this_a.type_d == Integral::DIAGONAL) {
    for (long row = 0; row < nrows; row++) {
      v_out_a(row) = v_in_a(row) * this_a.m_d(row);
    }
  }

  // type: SYMMETRIC
  //  loop through each column of the matrix and compute the dot product
  //  for elements in the upper triangle, we use the values from
  //  lower triangle 
  //
  else if (this_a.type_d == Integral::SYMMETRIC) {

    // loop over the rows of the matrix
    //          
    for (long col = 0; col < ncols; col++) {

      // loop over the upper triangular elements of the matrix
      //            
      TIntegral sum = 0;

      for (long row = 0; row < col; row++) {
        sum += v_in_a(row) * this_a.m_d(this_a.index(col, row));
      }

      // loop over the lower triangular elements of the matrix
      //
      for (long row = col; row < nrows; row++) {
        sum  += v_in_a(row) * this_a.m_d(this_a.index(row, col));
      }

      // output the results
      //
      v_out_a(col) = sum;
    }
  }

  // type: LOWER_TRIANGULAR
  //  same as SYMMETRIC, but only the lower triangular values are used
  //
  else if (this_a.type_d == Integral::LOWER_TRIANGULAR) {

    // loop over the rows of the matrix
    //          
    for (long col = 0; col < ncols; col++) {

      // loop over the upper triangular elements of the matrix
      //            
      TIntegral sum = 0;

      for (long row = col; row < nrows; row++) {
        sum  += v_in_a(row) * this_a.m_d(this_a.index(row, col));
      }

      // output the results
      //
      v_out_a(col) = sum;
    }
  }

  // type: UPPER_TRIANGULAR
  //  same as SYMMETRIC, but only the upper triangle values are used
  //
  else if (this_a.type_d == Integral::UPPER_TRIANGULAR) {

    // loop over the rows of the matrix
    //          
    for (long col = 0; col < ncols; col++) {

      // loop over the upper triangular elements of the matrix
      //            
      TIntegral sum = 0;

      for (long row = 0; row <= col; row++) {
        sum += v_in_a(row) * this_a.m_d(this_a.index(row, col));
      }

      // output the results
      //
      v_out_a(col) = sum;
    }
  }

  // type: SPARSE
  //  the general idea is that we loop through the non-zero elements only
  //  and get the dot product
  //
  else if (this_a.type_d == Integral::SPARSE) {

    // clear the contents of output vector
    //
    v_out_a.clear(Integral::RETAIN);
    
    // loop through non-zero elements of matrix, multiply and store in
    // output vector
    //
    long last_index = this_a.m_d.length();

    for (long index = 0; index < last_index; index++) {
      long row = this_a.row_index_d(index);
      long col = this_a.col_index_d(index);
      v_out_a(col) += this_a.m_d(index) * v_in_a(row);
    }
  }

  // exit gracefully
  //  
  return true;
}

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

// method: quadratic
//
// arguments:
//  TIntegral& output: (output) scalar return value
//  const MMatrix<TScalar, TIntegral>& this: (input) class operand
//  const MVector<TScalar, TIntegral>& v_in: (input) input vector
//
// return: a boolean value indicating status
//   v_in * this_a * transpose(v_in)
//
// this method computes v_in * this_a * transpose(v_in) and returns a
// scalar value, for example:
//
//                           [0 1 2]                     [0 1 2]   [2]
//  v_in = [2 0 1]  matrix = [3 4 5]  output = [2 0 1] * [3 4 5] * [0] = 24
//                           [6 7 8]	          	 [6 7 8]   [1]
//
template<class TScalar, class TIntegral>
boolean
MMatrixMethods::quadratic(TIntegral& output_a,
			  const MMatrix<TScalar, TIntegral>& this_a, 
			  const MVector<TScalar, TIntegral>& v_in_a) {

  // initialize the output
  //
  output_a = 0;
  
  // get the number of rows and columns of the matrix
  //
  long nrows = this_a.getNumRows();
  long ncols = this_a.getNumColumns();

  // check dimensions:
  //  (1) the number of rows and number of columns of the matrix can't be
  //      less than 1
  //  (2) the rows and columns of matrix must be equal to the length of
  //      input vector
  //
  if ((nrows < 1) || (ncols < 1) ||
      (nrows != v_in_a.length()) || (ncols != v_in_a.length())) {
    v_in_a.debug(L"v_in_a");
    this_a.debug(L"this_a");    
    return Error::handle(name(), L"quadratic", this_a.ERR_DIM,
			 __FILE__, __LINE__);
  }
  
  // type: DIAGONAL
  //  the general idea is that we only loop through the elements of diagonal
  //
  if (this_a.type_d == Integral::DIAGONAL) {

    // loop over diagonal elements of the matrix and multiply
    //
    for (long i = 0; i < nrows; i++) {
      output_a += (v_in_a(i) * v_in_a(i)) * this_a.m_d(i);
    }
  }

  // type: NON-DIAGONAL
  //  the general idea here is that we are going to build a vector to do
  //  vector * matrix first, and then do dot product of vector
  //
  else if ((this_a.type_d == Integral::FULL) ||
	   (this_a.type_d == Integral::SYMMETRIC) ||
	   (this_a.type_d == Integral::LOWER_TRIANGULAR) ||
	   (this_a.type_d == Integral::UPPER_TRIANGULAR) ||