/* -*- MaTX -*-
 *
 *  NAME
 *      detrend_col() - Remove a linear trend for each column
 *      detrend_row() - Remove a linear trend for each row
 *      detrend()     - Same as detrend_row()
 *
 *  SYNOPSIS
 *      y = detrend_col(x)
 *         Array y;
 *         Array x;
 *
 *      y = detrend_col(x,tr)
 *         Array y;
 *         Array x;
 *         Integer tr;
 *
 *      y = detrend_row(x)
 *         Array y;
 *         Array x;
 *
 *      y = detrend_row(x,tr)
 *         Array y;
 *         Array x;
 *         Integer tr;
 *
 *      y = detrend(x)
 *         Array y;
 *         Array x;
 *
 *      y = detrend(x,tr)
 *         Array y;
 *         Array x;
 *         Integer o;
 *
 *  DESCRIPTION
 *      detrend() returns an array obtained by removing the linear trend 
 *      from the data in X.
 *
 *      detrend(X,0) removes just the mean value from X.
 *
 *      If X is a matrix, 
 *
 *      detrend_col() removes the trend from each column of the matrix.
 *      detrend_row() removes the trend from each row of the matrix.
 *      detrend() works in the same way as detrend_row().
 */

Func Array detrend(X, tr, ...)
	Array X;
	Integer tr;
{
	error(nargchk(1, 2, nargs, "detrend"));

	if (nargs == 1) {
		detrend_row(X);
	} else {
		detrend_row(X, tr);
	}
}

Func Array detrend_row(X, tr, ...)
	Array X;
	Integer tr;
{
	error(nargchk(1, 2, nargs, "detrend_row"));
	
	if (nargs == 1) {
		trans(detrend_col(trans(X)));
	} else {
		trans(detrend_col(trans(X), tr));
	}
}

Func Array detrend_col(X_, tr, ...)
	Array X_;
	Integer tr;
{
	Matrix A,X,Y;
	Integer m;

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

	if (nargs == 1) {
		tr = 1;
	}
	
	X = makecolv(X_);
	m = Rows(X);

	if (tr == 0) { // Remove the mean
		Y = X - ONE(m,1)*mean_row(trans(X));
	} else {      // Remove line fit
		A = [trans([1:m])/m, ONE(m,1)];
		Y = X - A*(A \ X);
	}

	if (Rows(X_) == 1) {
		Y = trans(Y);
	}
	
	return Array(Y);
}