mmat_16.cc

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

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

// system include files
//
#include <typeinfo>

// isip include files
//
#include "MMatrixMethods.h"
#include "MMatrix.h"
#include <VectorLong.h>

#ifndef ISIP_TEMPLATE_unsigned

// method: eigen
//
// arguments:
//  const MMatrix<TScalar, TIntegral>& this: (input) input matrix
//  MVector<TScalar, TIntegral>& eigvals: (output) eigenvalues
//  MMatrix<TScalar, TIntegral>& eigvects: (output) eigenvectors
//
// return: a boolean value indicating status
//
// this method computes the eigenvalues of a real general matrix.
// we use Hessenberg Reduction + QR to compute eigenvalues and Inverse
// Iteration method to compute eigenvectors.
//
// reference:
//  W. Press, S. Teukolsky, W. Vetterling, B. Flannery,
//  Numerical Recipes in C, Second Edition, Cambridge University Press,
//  New York, New York, USA, pp. 482-495, 1995.
//
// the algorithm used here only supports matrices which have real
// eigenvalues and eigenvectors. if the matrix has complex eigenvalues,
// this method will return the real parts of the eigenvalues and
// eigenvectors.
//
// the eigen vectors are the rows of the output matrix and the eigen values
// are always sorted in decending order.
//
template <class TScalar, class TIntegral>
boolean MMatrixMethods::eigen(const MMatrix<TScalar, TIntegral>& this_a,
			      MVector<TScalar, TIntegral>& eigvals_a,
			      MMatrix<TScalar, TIntegral>& eigvects_a) {

#ifndef ISIP_TEMPLATE_complex

  // check arguments: non-square or singularity
  //
  if (!this_a.isSquare()) {
    this_a.debug(L"this_a");
    return Error::handle(name(), L"eigen()", Error::ARG, __FILE__, __LINE__);
  }

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

  // condition: empty input matrix
  //
  if (nrows <= 0) {
    this_a.debug(L"this_a");    
    return Error::handle(name(), L"eigen()",
			 MMatrix<TScalar, TIntegral>::ERR, 
			 __FILE__, __LINE__, Error::WARNING);
  }

  // initialize the output vector and matrix
  //
  eigvals_a.setLength(nrows);
  eigvects_a.changeType(Integral::FULL);
  eigvects_a.setDimensions(nrows, nrows);

  // if there is only one element in the matrix, the eigenvalue by
  // definition is the element's value.
  //
  if (nrows == 1) {
    eigvals_a(0) = this_a(0, 0);
    eigvects_a.setValue(0, 0, (TIntegral)1);
  }

  // else: do the long calculation
  //
  else {

    // for float and double matrices, use the same current matrix
    //
    if ((typeid(TIntegral) == typeid(float)) ||
	(typeid(TIntegral) == typeid(double)) ||
	(typeid(TIntegral) == typeid(complexfloat))) {

      // declare local variables
      //
      MVector<TScalar, TIntegral> eigimage(nrows);
      MVector<TScalar, TIntegral> tmp_vec(nrows);
      
      // make a copy of current matrix as it is const and we can't convert a
      // const matrix to a balance matrix
      // 
      MMatrix<TScalar, TIntegral> tmp;
      tmp.assign(this_a, Integral::FULL);
    
      // perform three key operations:
      //  - balance the matrix
      //  - use Hessenberg elimination
      //  - use QR to get the eigenvalues
      //
      tmp.eigenBalance();
      tmp.eigenEliminateHessenberg();
      tmp.eigenHessenbergQR(eigvals_a, eigimage);
      eigvals_a.sort(Integral::DESCENDING);

      // compute each eigenvector corresponding to an eigenvalue
      //
      for (long row = 0; row < nrows; row++) {
	tmp.copy(this_a);
	tmp.eigenComputeVector(tmp_vec, (TIntegral)eigvals_a(row));
	eigvects_a.setRow(row, tmp_vec);
      }
    }
    
    // for other matrices, copy the current matrix to a double matrix
    // and call eigen method - this copy is required for numerical
    // precision
    //
    else {

      // declare local variables
      //
      MVector<Double, double> eigvals;
      MMatrix<Double, double> eigvects;
      MMatrix<Double, double> tmp;
      
      // copy the current matrix to a double matrix
      //
      tmp.assign(this_a);
      
      // call eigen method for the matrix
      //
      if (!tmp.eigen(eigvals, eigvects)) {
	return false;
      }
      
      // assign the eigen values and eigen vectors to the output arguments
      //
      eigvals_a.assign(eigvals);
      eigvects_a.assign(eigvects);
    } 
  }
  
  // exit gracefully
  //
  return true;

#else

  // this algorithm is not supported for complex matrix
  //
  return Error::handle(name(), L"eigen: eigenvalues computation is not supported for complex matrix", Error::ARG, __FILE__, __LINE__);
  
#endif
}

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

// these methods are for MMatrix<Float, float> and MMatrix<Double, double>
// only because they use values related to float-point precision. so they
// can be compiled and linked only once
//
#ifdef ISIP_TEMPLATE_fp

// method: eigenBalance
//
// arguments:
//  MMatrix<TScalar, TIntegral>& this: (input/output) input matrix
//
// return: a boolean value indicating status
//
// this method replaces a matrix by a balanced matrix with
// identical eigenvalues. see:
//
//  W. Press, S. Teukolsky, W. Vetterling, B. Flannery,
//  Numerical Recipes in C, Second Edition, Cambridge University Press,
//  New York, New York, USA, pp. 482-484, 1995.
//
// for more details.
//
template<class TScalar, class TIntegral>
boolean MMatrixMethods::eigenBalance(MMatrix<TScalar, TIntegral>& this_a) {

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

  // choose the radix base employed for floating-point arithmetic
  // i.e. 2 for most machines, but 16 for IBM mainframe architectures
  //
  double radix = 2.0;
  double sqrdx = radix * radix;

  // begin the balance operations
  //
  boolean last = false;

  while (!last) {

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

      // initialize accumulators
      //
      double row_norm = 0;
      double col_norm = 0;
      double sum = 0;

      // calculate norms of row and column
      //
      for (long j = 0; j < nrows; j++) {
	if (j != i) {
	  col_norm += Integral::abs((double)this_a(j, i));
	  row_norm += Integral::abs((double)this_a(i, j));
	}
      }

      // when both norms are nonzero
      //
      if ((col_norm != 0) && (row_norm != 0)) {

	// compute a measure of the need to balance the matrix
	//
	sum = col_norm + row_norm;

	// find the integer power of the machine radix that comes
	// closest to balancing the matrix
	//
	double factor = 1;
	double g = row_norm / radix;
	while (col_norm < g) {
	  factor *= radix;
	  col_norm *= sqrdx;
	}
	
	g = row_norm * radix;
	while (col_norm > g) {
	  factor /= radix;
	  col_norm /= sqrdx;
	}

	// determine if a similarity transformation is needed
	//
	if ((MMatrix<TScalar, TIntegral>::THRESH_BALANCE * sum) >
	    (col_norm + row_norm) / factor) {

	  // perform a similarity transformation
	  //
	  last = false;
	  for (long j = 0; j < nrows; j++) {
	    this_a.setValue(i, j, this_a(i, j) / factor);
	  }
	  for (long j = 0; j < nrows; j++) {
	    this_a.setValue(j, i, this_a(j, i) * factor);
	  }
	}
      }
    }
  }

  // exit gracefully
  //
  return true;
}  

// method: eigenEliminateHessenberg
//
// arguments:
//  MMatrix<TScalar, TIntegral>& this: (input) input matrix
//
// return: a boolean value indicating status
//
// this method reduces a matrix to hessenberg's form by the elimination
// method. see:
//
//  W. Press, S. Teukolsky, W. Vetterling, B. Flannery,
//  Numerical Recipes in C, Second Edition, Cambridge University Press,
//  New York, New York, USA, pp. 484-486, 1995.
//
// for more details.
//
template<class TScalar, class TIntegral>
boolean MMatrixMethods::eigenEliminateHessenberg(MMatrix<TScalar,
						 TIntegral>& this_a) {

  // get the dimension of input matrix
  //
  long nrows = this_a.getNumRows();

  // process eliminating operation to the input matrix
  //
  for (long m = 1; m < (nrows - 1); m++) {

    // declare local variables
    //
    double x = 0;
    long i = m;

    // find the pivot
    //
    for (long j = m; j < nrows; j++) {
      if (Integral::abs(this_a(j, m - 1)) > Integral::abs(x)) {
	x = this_a(j, m - 1);
	i = j;
      }
    }

    // if the pivot is found, interchange rows and columns
    //
    if (i != m) {
      for (long j = m - 1; j < nrows; j++) {
	double tmp = this_a(i, j); 
	this_a.setValue(i, j, this_a(m, j)); 
	this_a.setValue(m, j, tmp);
      }
      for (long j = 0; j < nrows; j++) {
	double tmp = this_a(j, i); 
	this_a.setValue(j, i, this_a(j, m)); 
	this_a.setValue(j, m, tmp);
      }
    }

    // carry out the elimination
    //
    if (x != 0) {

      // loop over rows
      //
      for (i = m + 1; i < nrows; i++) {

	double y = this_a(i, m - 1);
	if (y != 0) {
	  y /= x;
	  this_a.setValue(i, m - 1, y);
	  
	  for (long j = m; j < nrows; j++) {
	    this_a.setValue(i, j, this_a(i, j) - y * this_a(m, j));
	  }
	  
	  for (long j = 0; j < nrows; j++) {
	    this_a.setValue(j, m, this_a(j, m) + y * this_a(j, i));
	  }
	}
      }
    }
  }

  // exit gracefully
  //
  return true;
}

// method: eigenHessenbergQR
//
// arguments:
//  MMatrix<TScalar, TIntegral>& this: (input) input matrix
//  MVector<TScalar, TIntegral>& eigval_r: (output) real part of eigenvalues
//  MVector<TScalar, TIntegral>& eigval_i: (output) imag part of eigenvalues
//
// return: a boolean value indicating status
//
// this method finds all eigenvalues of an upper hessenberg matrix. see:
//
//  W. Press, S. Teukolsky, W. Vetterling, B. Flannery,
//  Numerical Recipes in C, Second Edition, Cambridge University Press,
//  New York, New York, USA, pp. 491-493, 1995.
//
// for more details. this is a rather specialized method that has some
// internal constants associated with it. we chose not to make these
// part of the class header file because they are very algorithm specific.
//
// note that this method destroys the input.
//
template<class TScalar, class TIntegral>
boolean
MMatrixMethods::eigenHessenbergQR(MMatrix<TScalar, TIntegral>& this_a,
				  MVector<TScalar, TIntegral>& eigval_r_a,
				  MVector<TScalar, TIntegral>& eigval_i_a) {

  // get the dimension of input matrix
  //
  long nrows = this_a.getNumRows();

  // compute matrix norm for possible use in locating single small
  // subdiagonal element
  //
  double anorm = 0;
  for (long i = 0; i < nrows; i++) {
    for (long j = ((i - 1) < 0 ? 0: i - 1); j < nrows; j++) {
      anorm += Integral::abs(this_a(i, j));
    }
  }
  
  // declare local variables
  //
  double z = 0.0;
  double y = 0.0;
  double x = 0.0;
  double w = 0.0;
  double v = 0.0;
  double u = 0.0;
  double r = 0.0;
  double q = 0.0;
  double p = 0.0; 

  long hqr_max_iterations = 30;

  MVector<Double, double> wr(nrows);
  MVector<Double, double> wi(nrows);
  
  long nn = nrows;
  long l;
  double t = 0.0;
  double s = 0.0;

  // loop over the number of rows
  //
  while (nn >= 1) {

    long its = 0;

    do {

      // look for single small subdiagonal element
      //
      for (l = nn; l >= 2; l--) {
	s = Integral::abs(this_a(l - 2, l - 2)) +
	  Integral::abs(this_a(l - 1, l - 1));
	if (s == 0.0) {
	  s = anorm;
	}
	if ((double)(Integral::abs(this_a(l - 1, l - 2)) + s) == s) {
	  break;
	}
      }

      x = this_a(nn - 1, nn - 1);
      if (l == nn) {

	// one root is found
	//
	wr(nn - 1) = x + t;
	wi(nn - 1) = 0.0;
	nn--;
      }
      else {

	y = this_a(nn - 2, nn - 2);
	w = this_a(nn - 1, nn - 2) * this_a(nn - 2, nn - 1);
	if (l == (nn - 1)) {
	  p = 0.5 * (y - x);
	  q = p * p + w;
	  z = Integral::sqrt((double)Integral::abs(q));
	  x += t;
	  if (q >= 0.0) {
	    if (p > 0) {
	      z = Integral::abs(z);
	    }
	    else {
	      z = -Integral::abs(z);
	    }       
	    z += p;
	    wr(nn - 2) = wr(nn - 1) = x + z;
	    if (z != 0) {
	      wr(nn - 1) = x - w / z;
	    }
	    wi(nn - 2) = wi(nn - 1) = 0.0;
	  }
	  else {
	    wr(nn - 2) = wr(nn - 1) = x + p;
	    wi(nn - 2) = -(wi(nn - 1) = z);
	  }
	  nn -= 2;
	}
	else {
	  if (its >= hqr_max_iterations) {
	    this_a.debug(L"this_a");	    
	    return Error::handle(name(), L"eigenHessenbergQR()",
				 MMatrix<Double, double>::ERR,
	      			 __FILE__, __LINE__, Error::WARNING);
      	  }
	  if (its == 10 || its == 20) {
	    t += x;
	    for (long i = 1; i <= nn; i++) {
	      this_a.setValue(i - 1, i - 1, this_a(i - 1, i - 1) - x);
	    }
	    s = Integral::abs(this_a(nn - 1, nn - 2)) +
	      Integral::abs(this_a(nn - 2, nn - 3));
	    y = x = 0.75 * s;
	    w = -0.4375 * s * s;
	  }
	  its++;

	  long m;
	  for (m = (nn - 2); m >= l; m--) {
	    z = this_a(m - 1, m - 1);
	    r = x - z;
	    s = y - z;
	    p = (r * s - w) / this_a(m, m - 1) + this_a(m - 1, m);
	    q = this_a(m, m) - z - r - s;
	    r = this_a(m + 1, m);
	    s = Integral::abs(p) + Integral::abs(q) + Integral::abs(r);
	    p /= s;
	    q /= s;
	    r /= s;
	    if (m == l) {
	      break;
	    }
	    u = Integral::abs(this_a(m - 1, m - 2)) *
	      (Integral::abs(q) + Integral::abs(r));
	    v = Integral::abs(p) * (Integral::abs(this_a(m - 2, m - 2)) +
				    Integral::abs(z) +
				    Integral::abs(this_a(m, m)));
	    if ((double)(u + v) == v) {
	      break;
	    }
	  }
	  for (long i = m + 2; i <= nn; i++) {
	    this_a.setValue(i - 1, i - 3, 0.0);
	    if  (i != (m + 2)) {
	      this_a.setValue(i - 1, i - 4, 0.0);
	    }
	  }
	  for (long k = m; k <= nn - 1; k++) {
	    if (k != m) {
	      p = this_a(k - 1, k - 2);
	      q = this_a(k, k - 2);
	      r = 0.0;
	      if (k != (nn - 1)) {
		r = this_a(k + 1, k - 2);
	      }
	      x = Integral::abs(p) + Integral::abs(q) + Integral::abs(r);
	      if (x != 0) {
		p /= x;
		q /= x;
		r /= x;
	      }