/* -*- MaTX -*-
 *
 *  NAME
 *      bilinear() - Bilinear transformation
 *
 *  SYNOPSIS
 *      {Ad,Bd,Cd,Dd} = bilinear(A,B,C,D,fs)
 *         Matrix Ad,Bd,Cd,Dd;
 *         Matrix A,B,C,D;
 *         Real fs;
 *
 *      {Ad,Bd,Cd,Dd} = bilinear(A,B,C,D,fs,fp)
 *         Matrix Ad,Bd,Cd,Dd;
 *         Matrix A,B,C,D;
 *         Real fs;
 *         Real fp;
 *
 *  DESCRIPTION
 *      bilinear(A,B,C,D,fs) transforms the continuous-time transfer function
 *      to a discrete-time transfer function obtained by the bilinear 
 *      transformation:
 *
 *        G(z) = G(s) |
 *                    | s = 2*fs*(z-1)/(z+1)
 *
 *      where fs is the sample frequency in Hz.
 *
 *      bilinear(...,fp) uses prewarping before the bilinear 
 *      transformation so that the frequency responses before and 
 *      after mapping match exactly at the frequency fp (Hz).
 *
 *  SEE ALSO
 *      c2d and d2c
 */

Func List bilinear(A, B, C, D, fs, fp, ...)
	Matrix A, B, C, D;
	Real fs, fp;
{
	Integer n;
	Real r,t;
	Matrix Ad,Bd,Cd,Dd,T1,T2;

	error(nargchk(5, 6, nargs, "bilinear"));
	
	if (nargs == 6) {
		fp = 2*PI*fp;
		fs = fp/tan(fp/fs/2)/2;
	}

	// Bilinear transformation
	n = Rows(A);
	t = 1/fs;
	r = sqrt(t);

	T1 = I(n) + A*t/2;
	T2 = I(n) - A*t/2;
	Ad = T2 \ T1;
	Bd = t/r*(T2 \ B);
	Cd = r*C/T2;
	Dd = C/T2*B*t/2 + D;

	return {Ad,Bd,Cd,Dd};
}