mmat_15.cc

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

// file: $isip/class/math/matrix/MMatrix/mmat_15.cc
// version: $Id: mmat_15.cc,v 1.32 2002/02/27 20:54:30 alphonso Exp $
//

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

// method: decompositionLU
//
// arguments:
//  const MMatrix<TScalar, TIntegral>& this: (input) input matrix
//  MMatrix<TScalar, TIntegral>& l: (output) lower triangle matrix
//  MMatrix<TScalar, TIntegral>& u: (output) upper triangle matrix
//  MVector<Long, long>& index: (output) index of pivoting
//  long& sign: (output) indicate the number of row interchanges even or odd
//  double stabilize: (input) a stabilization factor
//
// return: a boolean value indicating status
//
// this method constructs the LU decomposition of the input matrix. the
// algorithm used is the right-looking algorithm:
//
//  W. Press, S. Teukolsky, W. Vetterling, B. Flannery,
//  Numerical Recipes in C, Second Edition, Cambridge University Press,
//  New York, New York, USA, pp. 46-47, 1995.
//
// the algorithm used in the above reference uses "implicit pivoting".
// we slightly modified the algorithm to use "partial pivoting". in
// implicit pivoting, elements of the matrix are scaled by the largest
// element of that row before comparisons are made for pivoting, which is
// computationally less efficient. in partial pivoting, comparisons
// are made on elements without any scaling, which is simpler and
// computationally efficient.
//
template <class TScalar, class TIntegral>
boolean
MMatrixMethods::decompositionLU(const MMatrix<TScalar, TIntegral>& this_a,
				MMatrix<TScalar, TIntegral>& l_a,
				MMatrix<TScalar, TIntegral>& u_a,
				MVector<Long, long>& index_a,
				long& sign_a,
				double stabilize_a) {

  // this method is defined only for float and double matrices only
  //
#ifdef ISIP_TEMPLATE_fp
    
  // check arguments: matrices must be square
  //
  if (!this_a.isSquare()) {
    this_a.debug(L"this_a");    
    return Error::handle(name(), L"decompositionLU()", Error::ARG,
			 __FILE__, __LINE__);
  }

  // initialize sign
  //
  sign_a = 1;
  
  // initialize pivoting index
  //
  index_a.setLength(this_a.getNumRows());
  index_a.ramp(0, 1);
  
  // type: DIAGONAL
  //
  if (this_a.type_d == Integral::DIAGONAL) {

    // the lower triangular matrix is an identity matrix
    //
    l_a.makeIdentity(this_a.getNumRows(), Integral::LOWER_TRIANGULAR);

    // the upper triangular matrix is the same as the current matrix
    //
    u_a.assign(this_a, Integral::UPPER_TRIANGULAR);
  }

  // type: UPPER_TRIANGULAR
  //
  else if (this_a.type_d == Integral::UPPER_TRIANGULAR) {

    // the lower triangular matrix is identity matrix
    //    
    l_a.makeIdentity(this_a.getNumRows(), Integral::LOWER_TRIANGULAR);

    // upper triangular matrix is the same as the current matrix
    //    
    u_a.assign(this_a);
  }

  // type: FULL, SYMMETRIC, LOWER_TRIANGULAR, SPARSE
  //
  else {

    // get the number of rows of the matrix
    //
    long nrows = this_a.getNumRows();

    // make a copy of the current matrix
    //
    MMatrix<TScalar, TIntegral> a;
    a.assign(this_a, Integral::FULL);

    // process each row in the decomposition
    //
    for (long j = 0; j < nrows; j++) {

      // declare local variables
      //
      double max = 0;
      long imax = j;

      // calculate the elements in the upper triangular matrix
      // iterating over columns. these elements either have the same row
      // indices but smaller column indices, or the same column indices but
      // smaller row indices:
      //
      //    u(i, j) = a(i, j) - sum[l(i, k) * u(k, j)],
      //   				     for (k=1,...,k-1) when i <= j
      //
      for (long i = 0; i < j; i++) {
	double sum = a(i, j);
	for (long k = 0; k < i; k++) {
	  sum -= a(i, k) * a(k, j);
	}
	a.setValue(i, j, sum);
      }

      // calculate the elements in the lower triangular matrix
      // iterating over columns. these elements either have same row indices
      // but smaller column indices, or the same column indices but smaller
      // row indices:
      //
      //    l(i, j) = a(i, j) - sum[l(i, k) * u(k, j)] / u(j, j),
      //					for (k=1,...,k-1) when i > j
      //
      for (long i = j; i < nrows; i++) {

	// accumulate the sum
	//
	double sum = a(i, j);
	for (long k = 0; k < j; k++) {
	  sum -= a(i, k) * a(k, j);
	}
	a.setValue(i, j, sum);

	// find the biggest number in column j and record the row index -
	// this is required for pivoting.
	//
	double value = Integral::abs(sum);
	if (value > max) {
	  max = value;
	  imax = i;
	}
      }

      // pivoting: check whether we need to interchange rows. we swap rows
      // if the diagonal element for a particular column is not the maximum
      // magnitude element. record the index of the pivot in the output array.
      //
      if (j != imax) {
	a.swapRows(imax, j);
	long k = index_a(j);
	index_a(j) = index_a(imax);
	index_a(imax) = k;
	sign_a = -sign_a;
      }

      // divide by the pivot element: the pivot element will not be zero
      // because the matrix is not singular. or it should be substituted
      // with stabilization factor
      //
      max = a(j, j);
      if ((max == 0) && (stabilize_a != 0)) {
	max = stabilize_a;
	a.setValue(j, j, (TIntegral)stabilize_a);
      }

      if (max != 0) {
	if (j != (nrows - 1)) {
	  TIntegral dum = 1.0 / a(j, j);
	  for (long i = j + 1; i < nrows; i++) {
	    a.setValue(i, j, a(i, j) * dum);
	  }
	}
      }
    }

    // copy the lower triangular part of the result matrix into l_a and
    // set the diagonal elements of l_a matrix to (TIntegral)1
    //
    l_a.setDimensions(nrows, nrows, false, Integral::LOWER_TRIANGULAR);
    l_a.setLower(a);
    l_a.setDiagonal((TIntegral)1);

    // get the upper triangle part of result matrix and assign to u_a
    //
    u_a.setDimensions(nrows, nrows, false, Integral::UPPER_TRIANGULAR);
    u_a.setUpper(a);
  }

  // exit gracefully
  //
  return true;

  // for matrices other than float and double, return error
  //
#else
  return Error::handle(name(), L"decompositionLU", Error::TEMPLATE_TYPE,
		       __FILE__, __LINE__);
#endif
}

// explicit instantiations
//
template boolean
MMatrixMethods::decompositionLU(const MMatrix<ISIP_TEMPLATE_TARGET>&,
				MMatrix<ISIP_TEMPLATE_TARGET>&,
			        MMatrix<ISIP_TEMPLATE_TARGET>&,
				MVector<Long, long>&, long&,
				double);

// method: decompositionCholesky
//
// arguments:
//  const MMatrix<TScalar, TIntegral>& this: (input) input matrix
//  MMatrix<TScalar, TIntegral>& l: (output) lower triangular matrix
//
// return: a boolean value indicating status
//
// this method constructs the Cholesky decomposition of an input matrix:
//
//  W. Press, S. Teukolsky, W. Vetterling, B. Flannery,
//  Numerical Recipes in C, Second Edition, Cambridge University Press,
//  New York, New York, USA, pp. 96-97, 1995.
//
// upon output, l' * l = this
// note that this method only operates on symmetric matrices.
//
template <class TScalar, class TIntegral>
boolean
MMatrixMethods::decompositionCholesky(const
				      MMatrix<TScalar, TIntegral>& this_a,
				      MMatrix<TScalar, TIntegral>& l_a) {

  // this method is defined only for float and double matrices only
  //
#ifdef ISIP_TEMPLATE_fp
      
  // condition: non-symmetric
  //
  if (!this_a.isSymmetric()) {
    this_a.debug(L"this_a");    
    return Error::handle(name(), L"decompositionCholesky()", Error::ARG,
			 __FILE__, __LINE__);
  }

  // type: DIAGONAL
  //
  if (this_a.type_d == Integral::DIAGONAL) {

    // the diagonal elements should be positive or zero
    //
    if ((double)this_a.m_d.min() < 0) {
      this_a.debug(L"this_a");      
      return Error::handle(name(), L"decompositionCholesky()", Error::ARG,
			   __FILE__, __LINE__, Error::WARNING);
    }

    // assign the input matrix to lower triangular matrix
    //
    l_a.assign(this_a);
    l_a.m_d.sqrt();
    l_a.changeType(Integral::LOWER_TRIANGULAR);
  }

  // type: FULL, SYMMETRIC, LOWER_TRIANGULAR, UPPER_TRIANGULAR, SPARSE
  //
  else {

    // get number of rows of current matrix
    //
    long nrows = this_a.getNumRows();

    // make a copy of the current matrix
    //
    MMatrix<TScalar, TIntegral> a;
    a.assign(this_a, Integral::FULL);

    // loop over all rows
    //
    for (long i = 0; i < nrows; i++) {

      // constructs a lower triangular matrix directly iterating over columns.
      // remember an upper triangular matrix is the transpose matrix of a
      // lower triangular matrix in our storage format.
      //
      for (long j = i; j < nrows; j++) {

	// accumulate the sum
	//
	double sum = a(i, j);
	for (long k = i - 1; k >= 0; k--) {
	  sum -= a(i, k) * a(j, k);
	}

	// handle diagonal elements separately
	//      
	if (i == j) {

	  // if the sum is not greater than zero, the matrix is
	  // not positive definite
	  //
	  if (sum <= 0.0) {
	    this_a.debug(L"this_a");	    
	    return Error::handle(name(), L"decompositionCholesky()", 
				 MMatrix<TScalar, TIntegral>::ERR_POSDEF,
				 __FILE__, __LINE__, Error::WARNING);
	  }

	  // assign the value
	  //
	  a.setValue(i, i, Integral::sqrt(sum));
	}

	// other elements are handled with the normal calculation
	//      
	else {
	  a.setValue(j, i, sum / a(i, i));
	}
      }
    }
    
    // copy the lower triangular part of the result matrix into l_a and
    // set the diagonal elements of l_a matrix to (TIntegral)1
    //
    l_a.setDimensions(nrows, nrows, false, Integral::LOWER_TRIANGULAR);
    l_a.setLower(a);
  }

  return true;

#else
  return Error::handle(name(), L"decompositionCholesky",
		       Error::TEMPLATE_TYPE,
		       __FILE__, __LINE__);
#endif
}

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

// method: decompositionSVD
//
// arguments:
//  const MMatrix<TScalar, TIntegral>& this: (input) input matrix
//  MMatrix<TScalar, TIntegral>& u: (output) row orthonormal
//  MMatrix<TScalar, TIntegral>& w: (output) singular value matrix
//  MMatrix<TScalar, TIntegral>& v: (output) column orthonormal
//
// return: a boolean value indicating status
//
// this method constructs the singular value decomposition of the input
// matrix such that u * w * transpose(v) = this
//
// this method relies heavily on concepts from the EISPACK math libraries
//   Burton S. Garbow, Mathematics and Computer Science Div, Argonne National
//   Laboratory, August 1983.
//
// Major modifications have been made to incorporate this routine into the
// ISIP IFC libraries
// 
// Householder bidiagonalization and a variant of the qr algorithm are used.
//
// An error is reported if 30 iterations do not produce convergence in the
// singular value decomposition. The number 30 is arbitrary but seems to be
// common to all implementations in math libraries such as LaPack, etc.
//
template <class TScalar, class TIntegral>
boolean
MMatrixMethods::decompositionSVD(const MMatrix<TScalar, TIntegral>& this_a,
				 MMatrix<TScalar, TIntegral>& u_a,
				 MMatrix<TScalar, TIntegral>& w_a,
				 MMatrix<TScalar, TIntegral>& v_a) {

  // this method is defined only for float and double matrices
  //
#ifdef ISIP_TEMPLATE_fp

  // check for matrix shape: SVD is defined for matrices where the number
  // of rows is greater than or equal to the number of columns. In the
  // converse case, we transpose A and solve.
  //
  long m = this_a.getNumRows();
  long n = this_a.getNumColumns();
  if (m < n) {
    MMatrix<TScalar, TIntegral> tmp(this_a);
    tmp.transpose();
    tmp.decompositionSVD(v_a, w_a, u_a);
    u_a.transpose();
    v_a.transpose();
    w_a.transpose();
  }

  // declare local variables
  //
  MVector<TScalar, TIntegral> rv1(n);
  rv1.assign(0.0);
  MVector<TScalar, TIntegral> tmp_w(n);
  tmp_w.assign(0.0);
  
  TIntegral anorm = 0.0;
  TIntegral c = 0.0;
  TIntegral f = 0.0;
  TIntegral g = 0.0;
  TIntegral h = 0.0;
  TIntegral s = 0.0;
  TIntegral scale = 0.0;
  TIntegral x = 0.0;
  TIntegral y = 0.0;
  TIntegral z = 0.0;
  TIntegral q = 0.0;
  
  TIntegral absf = 0.0;
  TIntegral absg = 0.0;
  TIntegral absh = 0.0;

  long i = 0;
  long its = 0;
  long j = 0;
  long jj = 0;
  long k = 0;
  long l = 0;
  long nm = 0;
  long minmn = (m < n) ? m : n;