cancoh = function(xmat, i1, i2) {
#
# Canonical coherency analysis using Moore-Penrose square root
#
# xmat	cross spectral density matrix
# i1	indices of first block of series
# i2	indices of second block of series
#
# returns an extended svd object:
# $d	canonical coherencies (absolute, not squared)
# $u.orig	coefficients of canonical components of first block
# $v.orig	coefficients of canonical components of second block
#
	e1 = eigen(xmat[i1, i1]);
	e2 = eigen(xmat[i2, i2]);
	a = xmat[i1, i2];
	n = dimnames(a);
	H = function(x) t(Conj(x));
	a = H(e1$vectors) %*% a %*% e2$vectors;
	a = a / sqrt(e1$values);
	a = t(t(a) / sqrt(e2$values));
	dimnames(a) = n;
	b = svd(a);
	b$u.orig = b$u / sqrt(e1$values);
	b$u.orig = e1$vectors %*% b$u.orig;
	dimnames(b$u.orig) = list(n[[1]], character(0));
	b$v.orig = b$v / sqrt(e2$values);
	b$v.orig = e2$vectors %*% b$v.orig;
	dimnames(b$v.orig) = list(n[[2]], character(0));
	for (j in 1:ncol(b$u)) {
		# make largest element real and positive
                i = which.max(Mod(b$u.orig[ , j]));
		z = b$u.orig[i, j];
		z = Conj(z / Mod(z));
		b$u.orig[ , j] = b$u.orig[ , j] * z;
		b$v.orig[ , j] = b$v.orig[ , j] * z;
	}
	b;
}