function [x,y] = xxcorr(num,den,t)
% XXCORR is used to find the theoretical autocorrelation 
%      function of the output of a linear system given by
%       NUM(s)/DEN(s).  The time range specified is t.
%
%   [x,y] = xxcorr(num,den,t)
%

% Perform spactral factorisation
[b,a] = tffact(num,den);
[a,b,c,d] = tf2ss(b,a);
% Find the analytical solution of the system
[S,Rts,K] = analyt(a,c,b);
V = c*S;
ii = find(t>=0); x = t(ii); L = length(x); y = zeros(size(x)); l = 0;
for i=1:length(K)
   if abs(imag(Rts(i))) <= 1e-5
      for j=1:K
         l = l+1; y = y+V(l)*x.^(j-1).*exp(Rts(i)*x); 
      end
   else
      y = y+V(l+1)*exp(real(Rts(i))*x).*cos(imag(Rts(i))*x);
      y = y+V(l+2)*exp(real(Rts(i))*x).*sin(imag(Rts(i))*x);
      l = l+2;
   end
end
%  Extending to negative axis
x = [-x(L:-1:2) x]; y = [y(L:-1:2) y]; 


function [B,A] = tffact(num,den)
%  TFFACT is used to perform power spectrum factorisation
%     for given transfer function NUM(s)/DEN(s), 
%  
%     [B,A] = TFFACT(NUM,DEN)
%   
%     NUM(s)/DEN(s) = B(s)/A(s) + B(-s)/A(-s) 
%
%   where A(s) has all its poles on the left hand side of
%   the s-plane. 
if length(num)==length(den), num = num(2:length(num)); end 
Rts = roots(den); m = length(Rts);
for i = 1:m
   if real(Rts(i)>0), Rts(i) = -Rts(i); end
end
A = poly(Rts); A = real(A), Z = []; Z(length(num)) = 1;
for i=length(num)-1:-1:1; Z(i) = -Z(i+1);  end
num1 = num.*Z;
X = conv(num,num1),
X = X(1:2:length(X))';  X = 0.5*[zeros(m-length(X),1); X];
k = 0; Xx = 1;
if length(A)>2, 
   Xx = [(-1)^(m-1)*A(2) (-1)^m*A(1) zeros(1,m-2)]; 
end 
V0 = Xx  
for i=2:m
   V0 = [0 0 V0(1:m-2)]; k = k+2;
   if k < m+1,
      V0(2) = (-1)^m*A(k+1);     
      if k < m, V0(1) = (-1)^(m-1)*A(k+2); end  
   end
   Xx = [Xx; V0];
end

B = inv(Xx)*X; B = B';