function H = orthonormal(num,den)
%
%H = ORTHONORMAL(NUM,DEN) returns the state-space description representing an orthonormal filter. 
%NUM and DEN are row vectors containing the coeffecients in the numerator and denominator, respectively of 
%the transfer function, starting with the highest order coefficient.
%
D=num(1);

%fill up with zeros in num
m = length(num);
n = length(den);
if (m~=n)
    num1 = num;
    for i=1:1:(n-m)
        num(i) = 0;
    end
    for i=(n-m)+1:1:n
        num(i) = num1(i-(n-m));
    end
end

%g = num
%p = den

%Check whether partial fraction expansion is needed
if num(1)~=0
    num=num-num(1)*den;
end;

N=length(den);
%separate odd and even denominator coefficients
if mod(N,2)==0
    for i=1:ceil(N/2)
        first(i)=den(2*i-1);
        second(i)=den(2*i);
    end;
else
    first(1)=den(1);
    for i=2:ceil(N/2)
        first(i)=den(2*i-1);
        second(i-1)=den(2*i-2);
    end;
end;


%loop for getting reactances
for i=1:(N-1)
    %i
    %first
    %second
    r(N-i)=first(1)/second(1);
    first=first(2:end)-r(N-i)*[second(2:end),zeros(1,(length(first(2:end))-length(second(2:end))))];
    temp=first;first=second;second=temp;
end;

%convert reactances in orthonormal coefficients
alpha(N-1)=1/r(N-1);
alpha(1:end-1)=1./sqrt((r(1:end-1).*r(2:end)));

%calculate F
beta1=sqrt(r(1)/pi);
F(1,N-1)=beta1*den(N);
F(2,N-2)=F(1,N-1)/alpha(1);
for i=3:N-1
    F(i,:)=([F(i-1,2:end),0]+alpha(i-2)*F(i-2,:))/alpha(i-1);
end;

%create state space system
A=diag(alpha(1:end-1),1)+diag(-alpha(1:end-1),-1);
A(N-1,N-1)=-alpha(N-1);
B(N-1)=sqrt(alpha(N-1)/pi);
B=B.';
C(N-1)=num(2)/F(N-1,1);
C(N-2)=num(3)/F(N-2,2);
for i=3:N-1     
    C(N-i)=(num(i+1)-C*F(:,i))/F(N-i,i);
end;

H=ss(A,B,C,D);