/* -*- MaTX -*-
 *
 *  NAME
 *      freqs() - Analog filter frequency response
 *
 *  SYNOPSIS
 *      Gjw = freqs(b,a,w)
 *        CoArray Gjw;
 *        Matrix b,a;
 *        Array w;
 *
 *  DESCRIPTION
 *      freqs() returns the frequency response Gjw
 *      at frequencies in w of the analog filter b(s)/a(s):
 *
 *                         nb-1      nb-2
 *             b(s)   b(1)s  +  b(2)s   + ... +  b(nb)
 *      G(s) = ---- = --------------------------------
 *                         na-1      na-2
 *             a(s)   a(1)s  +  a(2)s   + ... +  a(na)
 *
 *
 *  EXAMPLE
 *      To plot magnitude and phase of a filter:
 *
 *       w = logspace(d1,d2);
 *       Gjw = freqs(b,a,w);
 *       mg = abs(Gjw);
 *       ph = angle(Gjw);
 *       mgplot_semilogx(1,w,mg);
 *       mgplot_semilogx(2,w,ph);
 *
 *  SEE ALSO
 *      logspace, freqsp, and freqz.
 */

Func CoArray freqs(b, a, w)
	Matrix b, a;
	Array w;
{
	Complex j;
	CoArray jw,Gjw;

	j = (0,1);
	jw = j*w;
	if (version() == 4) {
		Gjw = eval(Polynomial(fliplr(b)), jw)/eval(Polynomial(fliplr(a)), jw);
	} else {
		Gjw = eval(Polynomial(b), jw) / eval(Polynomial(a), jw);
	}
	return Gjw;
}

/* -*- MaTX -*-
 *
 *  NAME
 *      freqsp() - Plot Analog filter frequency response
 *
 *  SYNOPSIS
 *      freqsp(b, a)
 *        Matrix b, a;
 *
 *      freqsp(b, a, w)
 *        Matrix b, a;
 *        Array w;
 *
 *  SEE ALSO
 *      freqs, logspace, and freqz.
 */

Func void freqsp(b, a, w, ...)
	Matrix b, a;
	Array w;
{
	Array mag, phase;
	CoArray Gjw;
	
	if (nargs < 2 || 3 < nargs) {
		error("freqsp(): Incorrect number of arguments.\n");
	}
	
	if (nargs == 2) {
		w = logspace(-1.0, 1.0);
	}
    Gjw = freqs(b,a,w);
    mag = 20*log10(abs(Gjw));
	phase = unwrap(angle(Gjw)*180/PI);

	mgplot_grid(1,1);
	mgplot_grid(2,1);
	mgplot_xlabel(1, "Frequency (rad/sec)");
	mgplot_xlabel(2, "Frequency (rad/sec)");
	mgplot_ylabel(1, "Magnitude (dB)");
	mgplot_ylabel(2, "Phase (deg)");
    mgplot_semilogx(1,w,mag,{"magnitude"});
	mgplot_semilogx(2,w,phase,{"phase"});
}