EE 451

Converting a transfer function to cascade form.

function [bcas,acas,g]=tocas(b,a)
%TOCAS Convert direct-form transfer function to cascade of biquad sections
%      
%      [BCAS,ACAS,G] = TOCAS(B,A) converts the direct-form transfer function
%                                    -1                -nb 
%                 B(z)   b(1) + b(2)z + .... + b(nb+1)z
%          H(z) = ---- = ----------------------------
%                                    -1                -na
%                 A(z)    1   + a(2)z + .... + a(na+1)z
%
%      to the cascade form
%                              -1       -2              -1       -2
%                   (1 + b(12)z + b(12)z  )  (1 + b(22)z + b(22)z  )
%          H(z) = G -----------------------  --------------------------- ...
%                              -1       -2              -1       -2
%                   (1 + a(12)z + a(12)z  )  (1 + a(22)z + a(22)z  )
%
%      BCAS will be a matrix of the numerator coefficients, 
%      ACAS will be a matrix of the denominator coefficients, 
%      G will be the gain.

p=cplxpair(roots(a));
z=cplxpair(roots(b));

M=length(z);
k = find(b ~= 0);
ninf = k(1)-1;  % Number of zeros at infinity (leading zeros in b)
N=length(p);

g=b(ninf+1)/a(1);

nsec = floor((max(M+ninf,N)+1)/2);
bcas=[ones(nsec,1) zeros(nsec,2)];
acas=[ones(nsec,1) zeros(nsec,2)];

for k=2:2:N
  acas(k/2,:)=real(poly([p(k-1),p(k)]));
end

if (rem(N,2)==1) 
    acas((N+1)/2,1:2) = poly(p(N));
end

for k=2:2:M
  bcas(k/2,:)=real(poly([z(k-1),z(k)]));
end

if (rem(M,2)==1) 
    bcas((M+1)/2,1:2) = poly(z(M));
end

k = floor((M+3)/2);

while (ninf >= 2)
  bcas(k,:) = [0 0 1];
  ninf = ninf - 2;
  k = k + 1;
end;

if (ninf == 1)
  if (rem(M,2) == 0)
    bcas(k,:) = [0 1 0];
  else
    bcas((M+1)/2,:) = [0 poly(z(M))];
  end
end