function [A,B,C,D] = ob2ss(poles,theta)
%   Given a model parameterised in orthonormal form as
%
%   G(z) = theta_1*B_1(z) + theta_2*B_2(z) + ... + theta_n*B_n(z)
%   
%   where the B_k(z) are the general orthonormal bases, defined as
%
%                  sqrt(1 - |psi_0|^2)
%       B_0 (z) = ---------------------
%                       z - psi_0
%
%                  sqrt(1 - |psi_k|^2)                     1 - z psi_i^*
%       B_k (z) = --------------------- PROD_[i=0]^[k-1] ----------------
%                       z - psi_k                          z - psi_i
%       for k >= 1
%  with
%  
%   theta:  theta = [theta_1;theta_2;...;theta_n]  (column) vector of parameters 
%           defining the model with respect to the orthonormal bases.
%
%   poles:  poles = [psi_1,psi_2,...,psi_n]  poles used to define the general
%           orthonormal basis functions parameterising the model.
%
%   The function computes the matrices A, B, C and D of a minimal state space 
%   realization.
%   References: 
%   Gomez, J.C.. "Analysis of Dynamic System Identification using Rational Orthonormal 
%   Bases", PhD. Thesis, University of Newcastle, Australia, February 1998.
%------------------------------------------------------------------------------------------
%   Usage:
%   >> poles = [0.5 0.4 0.3 0.2 0.1];
%   >> theta = [1 ;0.5; 2; 0.1;0.6]; 
%   >> [A,B,C,D] = ob2ss(poles,theta);
%------------------------------------------------------------------------------------------    
%
%   written by Juan Carlos Gomez.
%
%   Last Update 27/10/03.
%--------------------------------------------------------------------------
p = length(poles);
A = diag(poles);
B = [1;zeros(p-1,1)];
Ctilde = eye(p,p);
pp = -poles;
wun = ones(size(poles));
Lambda = sqrt(diag(wun - abs(poles).^2));

%  Run down first sub-diagonals putting in the correct entries.

for m=2:p
  for n = 1:m-1
    A(m,n) = 1-poles(m)*poles(m-1);
    Ctilde(m,n) = prod(pp(n:m-1));
  end;
end;

%  For the case of the A matrix re-do 3rd and lower sub-diagonals putting in
%  some product terms.

for m=3:p
  for n = 1:m-1
    A(m,n) = A(m,n)*prod(pp(n:m-2));
  end;
end;
C=theta'*Lambda*Ctilde;
D=0;