/* -*- MaTX -*-
 *
 *  NAME
 *      freqz() - Digital filter frequency response
 *
 *  SYNOPSIS
 *     {Gjw,w} = freqz(b,a,n)
 *        CoArray Gjw;
 *        Array w;
 *        Matrix b, a;
 *        Integer n;
 *
 *  DESCRIPTION
 *      freqz() returns the n-point frequency vector w and frequency
 *      response Gjw of the digital filter b(z)/a(z):
 *
 *                                     -1                -nb 
 *                   b(z)   b(1) + b(2)z + .... + b(nb+1)z
 *       G(e^(jw)) = ---- = ----------------------------
 *                                     -1                -na
 *                   a(z)    1   + a(2)z + .... + a(na+1)z
 *
 *      The response is evaluated at n points equally spaced 
 *      around the upper half of the unit circle.
 *
 *  EXAMPLE
 *      To plot magnitude and phase of a filter:
 *
 *        {Gjw,w} = freqz(b,a,n);
 *        mag = abs(Gjw);
 *        phase = angle(Gjw);
 *        mgplot_semilogx(1,w,mag); 
 *        mgplot_semilogx(2,w,phase)
 *
 *  SEE ALSO
 *      filter, fft, and freqs.
 *
 */

Func List freqz(b_, a_, n)
	Matrix b_, a_;
	Integer n;
{
	Integer na,nb,nn,N;
	Matrix a, b;
	CoArray Gjw;
	Array w;

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

	w = 2*PI/[0:n-1]/n;

	Gjw = fft([b, Z(1,n-nb)]) / fft([a, Z(1,n-na)]);
	Gjw = Gjw(1:n);

	return {Gjw,w};
}

/* -*- MaTX -*-
 *
 *  NAME
 *      freqzp() - Plot Digital filter frequency response
 *
 *  SYNOPSIS
 *      freqzp(b, a, n)
 *         Matrix b, a;
 *         Integer n;
 *
 *  SEE ALSO
 *      freqz and freqsp
 *
 */

Func void freqzp(b, a, n)
	Matrix b, a;
	Integer n;
{
	CoArray Gjw;
	Array w, mag, phase;

	{Gjw,w} = freqz(b, a, n);

    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(1,w,mag,{"magnitude"});
	mgplot(2,w,phase,{"phase"});
}