/* -*- MaTX -*-
 *
 *  NAME
 *      filter() - Digital filter
 *
 *  SYNOPSIS
 *      {y,z} = filter(b,a,x)
 *         Array y,z;
 *         Matrix b,a;
 *         Array x;
 *
 *      {y,z} = filter(b,a,x,zi)
 *         Array y,z;
 *         Matrix b,a;
 *         Array x;
 *         Matrix zi;
 *
 *  DESCRIPTION
 *      filter(b,a,x) filters the data in x with the filter given by
 *
 *         y(n) = b(1)*x(n) + b(2)*x(n-1) + ... + b(nb+1)*x(n-nb)
 *                          - a(2)*y(n-1) - ... - a(na+1)*y(n-na)
 *
 *      filter(...,zi) gives access to initial and final conditions.
 *
 *  SEE ALSO
 *
 */

Func List filter(b_, a_, x_, zi_, ...)
	Matrix b_;
	Matrix a_;
	Array x_;
	Matrix zi_;
{
	Integer i,na,nb,nx;
	Matrix a,b,a2t,b2t,x_i,y_i,zi,z,zab;
	Array x,y;

	error(nargchk(3, 4, nargs, "filter"));

	b = makerowv(b_);
	a = makerowv(a_);
	x = makerowv(x_);

	na = length(a) - 1;
	nb = length(b) - 1;
	nx = length(x);

	if (nargs == 3) {
		zi = Z(max(na,nb),1);
	} else {
		zi = makecolv(zi_);
	}

	if (na > 0) {
		if (abs(a(1)) <= EPS) {
			error("filter(): 1st coefficient of denominator must be nonzero.\n");
		}
		a = a/a(1);
		b = b/a(1);
	}
	
	if (na > 0) {
		a2t = trans(a(2:));
	}
	if (nb > 0) {
		b2t = trans(b(2:));
	}
	
	if (iscomplex(a) || iscomplex(b) || iscomplex(x) || iscomplex(zi)) {
		y = (Z(nx,1),:);
	} else {
		y = Z(nx,1);
	}
	z = zi;
	zab = Z(abs(na-nb),1);

	if (nb > na) {
		if (nb == 1) {
			for (i = 1; i <= nx; i++) {
				y(i) = b(1)*x(i) + z(1);
				z = b2t*x(i);
			}
		} else {
			if (na == 0) {
				for (i = 1; i <= nx; i++) {
					y(i) = b(1)*x(i) + z(1);
					z = [[z(2:)]
						 [  0  ]] + b2t*x(i);
				}
			} else {
				for (i = 1; i <= nx; i++) {
					y(i) = b(1)*x(i) + z(1);

					z = [[z(2:)]
						 [  0  ]] + b2t*x(i) - [[a2t*y(i)]
												 [   zab  ]];
				}
			}
		}
	} else if (nb == na) {
		if (na == 0) {
			for (i = 1; i <= nx; i++) {
				y(i) = b(1)*x(i);
			}
		} else if (na == 1) {
			for (i = 1; i <= nx; i++) {
				y(i) = b(1)*x(i) + z(1);
				z = b2t*x(i) - a2t*y(i);
			}
		} else {
			for (i = 1; i <= nx; i++) {
				y(i) = b(1)*x(i) + z(1);

				z = [[z(2:)]
					 [  0  ]] + b2t*x(i) - a2t*y(i);
			}
		}
	} else {
		if (na == 1) {
			for (i = 1; i <= nx; i++) {
				y(i) = b(1)*x(i) + z(1);
				z = - a2t*y(i);
			}
		} else {
			if (nb == 0) {
				for (i = 1; i <= nx; i++) {
					y(i) = b(1)*x(i) + z(1);

					z = [[z(2:)]
						 [  0  ]] - a2t*y(i);
				}
			} else {
				for (i = 1; i <= nx; i++) {
					y(i) = b(1)*x(i) + z(1);

					z = [[z(2:)]
						 [  0  ]] + [[b2t*x(i)]
									 [   zab  ]] - a2t*y(i);
				}
			}
		}
	}
	
	if (Cols(x_) > Rows(x_)) {
		y = trans(y);
		z = trans(z);
	}

	return {y, z};
}