/* -*- MaTX -*-
 *
 *  NAME
 *      freqzw() - Digital filter frequency response
 *
 *  SYNOPSIS
 *      {h,w} = freqzw(b,a,w)
 *         CoArray h;
 *         Array w;
 *         Matrix b, a;
 *         Array w;
 *
 *  DESCRIPTION
 *      freqzw() returns the complex frequency response h
 *      at frequencies designated in vector w of the filter b/a:
 *                                -1                -nb 
 *         jw   B(z)   b(1) + b(2)z + .... + b(nb+1)z
 *       H(e) = ---- = ----------------------------
 *   	                       -1                -na
 *              A(z)    1   + a(2)z + .... + a(na+1)z
 *
 *  EXAMPLE
 *      To plot magnitude and phase of a filter:
 *
 *       w = logspace(d1,d2);
 *       {h} = freqzw(b,a,w);
 *       mg = abs(h); 
 *       ph = angle(h);
 *       mgplot_semilogx(1,w,mg);
 *       mgplot_semilogx(2,w,ph);
 *
 *  SEE ALSO
 *      fft and freqs
 */

Func List freqzw(b_, a_, w)
	Matrix b_, a_;
	Array w;
{
	Integer na,nb;
	Complex j;
	Matrix a, b;
	CoArray h,ejw;

	j = (0,1);

	a = makerowv(a_);
	b = makerowv(b_);
	na = length(a);
	nb = length(b);

	a = [a Z(1,nb-na)];
	b = [b Z(1,na-nb)];
	ejw = exp(j*w);

	if (version() == 4) {
		h = eval(Polynomial(fliplr(b)),ejw) / eval(Polynomial(fliplr(a)),ejw);
	} else {
		h = eval(Polynomial(b),ejw) / eval(Polynomial(a),ejw);
	}
	return {h, w};
}