mmat_15.cc

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

  long maxmn = (m > n) ? m : n;

  // set the dimensions for the output matrices
  //
  v_a.setDimensions(maxmn, maxmn, false, Integral::FULL);
  v_a.assign(0.0);

  w_a.setDimensions(m, n, false, Integral::FULL);
  w_a.assign(0.0);
  
  // copy the input for manipulation
  //
  u_a.assign(this_a, Integral::FULL);
  u_a.setDimensions(maxmn, maxmn, true);
  
  // Householder reduction to bidiagonal form
  //
  for (i = 0; i < n; i++) {
    l = i + 1;
    rv1(i) = scale * g;
    g = 0.0;
    s = 0.0;
    scale = 0.0;

    if (i < m) {

      for (k = i; k < m; k++) {
	scale = scale + Integral::abs(u_a(k, i));
      }

      if (scale != 0.0) {
	for (k = i; k < m; k++) {
	  u_a.setValue(k, i, u_a(k, i) / scale);
	  s = s + u_a(k, i) * u_a(k, i);
	}  

	f = u_a(i, i);
	g = (f > 0.0) ? -Integral::abs(Integral::sqrt(s)) :
	  Integral::abs(Integral::sqrt(s));
	h = f * g - s;
	u_a.setValue(i, i, f - g);

	for (j = l; j < n; j++) {
	  s = 0.0;
	  for (k = i; k < m; k++) {
	    s = s + u_a(k, i) * u_a(k, j);
	  }
	  f = s / h;

	  for (k = i; k < m; k++) {
	    u_a.setValue(k, j, u_a(k, j) + f * u_a(k, i));
	  }
	}
	for (k = i; k < m; k++) {
	  u_a.setValue(k, i, scale * u_a(k, i));
	}
      }
    }
    
    tmp_w(i) = scale * g;
    g = 0.0;
    s = 0.0;
    scale = 0.0;
    if ((i < m) && (i != n - 1)) {
      for (k = l; k < n; k++) {
	scale = scale + Integral::abs(u_a(i, k));
      }

      if (scale != 0.0) {
	for (k = l; k < n; k++) {
	  u_a.setValue(i, k, u_a(i, k) / scale);
	  s = s + u_a(i, k) * u_a(i, k);
	}

	f = u_a(i, l);
	g = (f > 0.0) ? -Integral::abs(Integral::sqrt(s)) :
	  Integral::abs(Integral::sqrt(s));
	h = f * g - s;
	u_a.setValue(i, l, f - g);

	for (k = l; k < n; k++) {
	  rv1(k) = u_a(i, k) / h;
	}

	for (j = l; j < m; j++) {

	  s = 0.0;

	  for (k = l; k < n; k++) {
	    s = s + u_a(j, k) * u_a(i, k);
	  }

	  for (k = l; k < n; k++) {
	    u_a.setValue(j, k, u_a(j, k) + s * rv1(k));
	  }
	}

	for (k = l; k < n; k++) {
	  u_a.setValue(i, k, scale * u_a(i, k));
	}
      }
    }
    anorm = Integral::max(anorm, (Integral::abs(tmp_w(i)) +
			  Integral::abs(rv1(i))));
  }

  // QR accumulation (compute V)
  //
  for (i = n - 1; i >= 0; i--) {
    if (i < n - 1) {
      if (g != 0.0) {
	for (j = l; j < n; j++) {

	  // double division avoids possible underflow
	  //
	  v_a.setValue(j, i, (u_a(i, j) / u_a(i, l)) / g);
	}

	for (j = l; j < n; j++) {
	  s = 0.0;
	  for (k = l; k < n; k++) {
	    s = s + u_a(i, k) * v_a(k, j);
	  }

	  for (k = l; k < n; k++) {
	    v_a.setValue(k, j, v_a(k, j) + s * v_a(k, i));
	  }
	}
      }
      
      for (j = l; j < n; j++) {
	v_a.setValue(i, j, 0.0);
	v_a.setValue(j, i, 0.0);
      }
    }
    v_a.setValue(i, i, 1.0);
    g = rv1(i);
    l = i;
  }

  // LQ accumulation (compute U)
  //
  for (i = minmn - 1; i >= 0; i--) {
    l = i + 1;
    g = tmp_w(i);

    for (j = l; j < n; j++) {
      u_a.setValue(i, j, 0.0);
    }

    if (g != 0.0) {
      
      for (j = l; j < n; j++) {
	s = 0.0;

	for (k = l; k < m; k++) {
	  s = s + u_a(k, i) * u_a(k, j);
	}

	// double division avoids possible underflow
	//
	f = (s / u_a(i, i)) / g;

	for (k = i; k < m; k++) {
	  u_a.setValue(k, j, u_a(k, j) + f * u_a(k, i));
	}
      }

      for (j = i; j < m; j++) {
	u_a.setValue(j, i, u_a(j, i) / g);
      }
    }
    else {
      for (j = i; j < m; j++) {
	u_a.setValue(j, i, 0.0);
      }
    }
    u_a.setValue(i, i, u_a(i, i) + 1.0);
  }

  // diagonalization of the bidiagonal form (creating the w matrix)
  // and normalizing u and v accordingly
  //
  // loop over singular values
  //
  for (k = n - 1; k >= 0; k--) {

    for (its = 1; its <= 30; its++) {

      boolean flag = true;

      // test for splitting
      //
      for (l = k; l >= 0; l--) {
	nm = l - 1;
	if (((float)anorm + (float)Integral::abs(rv1(l))) == (float)anorm) {
	  flag = false;
	  break;
	}

	// rv1(0) is always zero. thus this section won't be reached so there
	// is no need to worry with nm < 0
	//
	if (((float)anorm + (float)Integral::abs(tmp_w(nm))) == (float)anorm) {
	  break;
	}
      }
      if (flag) {

	// cancellation of rv1(l) if l greater than 1
	//
	c = 0.0;
	s = 1.0;

	for (i = l; i <= k; i++) {
	  f = s * rv1(i);
	  rv1(i) = c * rv1(i);
	  if (((float)anorm + (float)Integral::abs(f)) == (float)anorm) {
	    break;
	  }
	  g = tmp_w(i);

	  // compute (a**2 + b**2) ** (1/2) without destructive underflow
	  // or overflow
	  //
	  absf = Integral::abs(f);
	  absg = Integral::abs(g);
	  if (absf > absg) {
	    h = (absf * Integral::sqrt(1.0 + (absg / absf) * (absg / absf)));
	  }
	  else {
	    if (absg == 0.0) {
	      h = 0.0;
	    }
	    else {
	      h = (absg * Integral::sqrt(1.0 + (absf / absg) *
					 (absf / absg)));
	    }
	  }
	  
	  tmp_w(i) = h;
	  c = g / h;
	  s = -f / h;
	  for (j = 0; j < m; j++) {
	    y = u_a(j, nm);
	    z = u_a(j, i);
	    u_a.setValue(j, nm, y * c + z * s);
	    u_a.setValue(j, i, z * c - y * s);
	  }
	}
      }

      // test for convergence
      //
      z = tmp_w(k);

      if (l == k) {
	    
	// convergence is reached
	//
	if (z < 0.0) {

	  // singular value is made non-negative and we update the sign of the
	  // right singular vectors
	  //
	  tmp_w(k) = -z;
	  for (j = 0; j < n; j++) {
	    v_a.setValue(j, k, -(v_a(j, k)));
	  }
	}
	break;
      }

      // if we've computed the maximum number of iterations and still not
      // converged then give up
      //
      if (its == 30) {
	u_a.debug(L"u_a");	
	this_a.debug(L"this_a");	
	return Error::handle(name(), L"decompositionSVD", u_a.ERR_SINGLR,
			     __FILE__, __LINE__);
      }

      // else we need to iterate again
      //
      // shift from bottom 2 by 2 minor
      x = tmp_w(l);
      nm = k - 1;
      y = tmp_w(nm);
      g = rv1(nm);
      h = rv1(k);
      f = ((y - z) * (y + z) + (g - h) * (g + h)) / (2.0 * h * y);

      // compute (a**2 + b**2) ** (1/2) without destructive underflow
      // or overflow
      //
      absf = Integral::abs(f);
      if (absf > 1.0) {
	g = (absf * Integral::sqrt(1.0 + (1.0 / absf) * (1.0 / absf)));
      }
      else {
	g = Integral::sqrt(1.0 + (absf) * (absf));
      }

      q = (f > 0.0) ? Integral::abs(g) : -Integral::abs(g);
      f = ((x - z) * (x + z) + h * ((y / (f + q)) - h)) / x;

      // next QR transformation
      //
      c = 1.0;
      s = 1.0;

      for (j = l; j <= nm; j++) {
	i = j + 1;
	g = rv1(i);
	y = tmp_w(i);
	h = s * g;
	g = c * g;

	absf = Integral::abs(f);
	absh = Integral::abs(h);
	if (absf > absh) {
	  z = (absf * Integral::sqrt(1.0 + (absh / absf) * (absh / absf)));
	}
	else {
	  if (absh == 0.0) {
	    z = 0.0;
	  }
	  else {
	    z = (absh * Integral::sqrt(1.0 + (absf / absh) *
				       (absf / absh)));
	  }
	}

	rv1(j) = z;
	c = f / z;
	s = h / z;
	f = x * c + g * s;
	g = g * c - x * s;
	h = y * s;
	y = y * c;

	for (jj = 0; jj < n; jj++) {
	  x = v_a(jj, j);
	  z = v_a(jj, i);
	  v_a.setValue(jj, j, x * c + z * s);
	  v_a.setValue(jj, i, z * c - x * s);
	}

	absf = Integral::abs(f);
	absh = Integral::abs(h);
	if (absf > absh) {
	  z = (absf * Integral::sqrt(1.0 + (absh / absf) * (absh / absf)));
	}
	else {
	  if (absh == 0.0) {
	    z = 0.0;
	  }
	  else {
	    z = (absh * Integral::sqrt(1.0 + (absf / absh) *
				       (absf / absh)));
	  }
	}
	tmp_w(j) = z;

	// rotation can be arbitrary if z is zero
	//
	if (z != 0.0) {
	  c = f / z;
	  s = h / z;
	}

	f = c * g + s * y;
	x = c * y - s * g;

	for (jj = 0; jj < m; jj++) {
	  y = u_a(jj, j);
	  z = u_a(jj, i);
	  u_a.setValue(jj, j, y * c + z * s);
	  u_a.setValue(jj, i, z * c - y * s);
	}
      }      

      rv1(l) = 0.0;
      rv1(k) = f;
      tmp_w(k) = x;
    }
  }

  // reorder singular values in descending order. adjust the u and v
  // matrices as well
  //
  // sort the w vector
  //
  VectorLong sort_indices;
  tmp_w.index(sort_indices);

  // reverse the order (the index method sorts in ascending order)
  //
  for (long i = 0; i < n / 2; i++) {
    long tmp = sort_indices(i);
    sort_indices(i) = sort_indices(n - i - 1);
    sort_indices(n - i - 1) = tmp;
  }

  // reorder matrices
  //
  u_a.reorderColumns(sort_indices);
  v_a.reorderColumns(sort_indices);
  tmp_w.reorder(sort_indices);

  // copy data out to the w matrix
  //
  for (i = 0; i < minmn; i++) {
    w_a.setValue(i, i, tmp_w(i));
  }

  // truncate the orthonormal matrices as needed
  //
  v_a.setDimensions(n, n, true);
  u_a.setDimensions(m, m, true);

  // exit gracefully
  //
  return true;
#else
  return Error::handle(name(), L"decompositionSVD", Error::TEMPLATE_TYPE,
		       __FILE__, __LINE__);
#endif
}

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