mBilinear(b,a,Fs)

% use polynomial multiplication to computer bilinear transform
% the input is H(s) in decending order polynomial
% the output is H(z) in decending order polynomial
%
% Usage : [c,d] = mBilinear(a,b,Fs);
% with a:  vector in decending order of H(s) numerator
%        b:  vector in decending order of H(s) denumerator
%        c:  H(z) numerator in decending order 
%        d:  H(z) denumerator in decending order
%        Fs: sampling frequency
%
%  Note: vector a and b must have the same length. The highest order of
%           H(z)'s numerator and denumerator are scaled to be 1. This is
%           different to Matlab's bilinear function. For example:
%   [c,d] = mBilinear([0 1 1],[1 0 1],1) will return
%   c = [1.0000 0.6667 -0.3333] and d = [1.0000 -1.2000 1.0000]
%  while using Matlab's bilinear function, 
%   [p,q ] = bilinear([0 1 1],[1 0 1],1) will return
%   p = [0.6000 0.4000 -0.2000] and q = [1.0000 -1.2000 1.0000]
%   It is obvious that c = p/p(1);
%

function [c,d] = mBilinear(b,a,Fs)
n = length(b)-1;
btop = zeros(n,n+1);
atop = zeros(n,n+1);
line = 1;
for (k = n:-1:0) % kth term index
if(k == 0) ini = [1 1];
else ini = [1 -1];
end
for (i=k-1:-1:1)      % conv index
    ini=conv([1 -1],ini);
end
% see if need to conv with [1 1] j times
if (k == 0) len = n-k-1;
else len = n-k;
end
for (j=1:len)
    ini = conv([1 1],ini);
end
btop(line,:) = b(line)*2^k*ini*Fs;
atop(line,:) = a(line)*2^k*ini*Fs;
line = line + 1;
end
bf = sum(btop);
af = sum(atop);
c = bf / bf(1);
d = af/ af(1);