/*  -*- MaTX -*-
 *
 *  NAME
 *      xcorr() - Cross-correlation function
 *
 *  SYNOPSIS
 *      c = xcorr(x)
 *         Array c;
 *         Array x;
 *
 *      c = xcorr(x,y)
 *         Array c;
 *         Array x,y;
 *
 *  DESCRIPTION
 *      xcorr(x,y) returns the cross-correlation of x and y in a 
 *      (2*m-1)-point vector, where x and y are m-point vectors.
 *
 *      xcorr(x), when x is a vector, is the auto-correlation.
 *
 *      xcorr(x), when x is an m-by-n matrix, is a matrix with
 *      2*m-1 rows whose n^2 columns contain the cross-correlation
 *      sequences for all combinations of the columns of x.
 *
 *	SEE ALSO
 *      xcov and conv
 */

Func Array xcorr(x_, y_, ...)
	Array x_, y_;
{
	Integer i,j;
	Integer nx,ny,xm,xn,m,k1,k2;
	Array X,x,y,c,cc,tmp;

	error(nargchk(1, 2, nargs, "xcorr"));

	x = x_;

	if (nargs == 1) {
		if (min(Cols(x), Rows(x)) == 1) {
			X = [makecolv(x) makecolv(x)];
		} else {
			X = x;
		}
	} else {
		y = y_;

		if (min(Cols(x), Rows(x)) != 1 && min(Cols(y), Rows(y)) != 1) {
			error("xcorr(): Arguments must be vectors.\n");
		}

		nx = length(x);
		ny = length(y);
		if (nx > ny) {
			y(nx) = 0;
		} else if (nx < ny) {
			x(ny) = 0;
		}
		X = [makecolv(x) makecolv(y)];
	}

	xm = Rows(X);
	xn = Cols(X);
	m = 2*xm - 1;
	c = Z(m,xn^2);

	for (i = 1; i <= xn; i++) {
		tmp = conj(X(xm:-1:1,i));
		tmp(m) = 0;
		for (j = 1; j <= i; j++) {
			{cc} = filter(X(:,j),[1],tmp);
			k1 = (i-1)*xn + j;
			k2 = (j-1)*xn + i;
			c(:,k2) = trans(cc);
			if (i != j) {
				c(:,k1) = conj(c(m:-1:1,k2));
			}
		}
	}
	
	c = c(:,nargs);

	if (Rows(x_) == 1) {
		c = trans(c);
	}

	return c;
}