function [S,Sconf,per,tau,exctn]=armode(A,C,th)
%ARMODE	Spectral decomposition of AR process.
%
% 	[S,Sconf,per,tau,exctn]=ARMODE(A,C,th) produces spectral
%	information for an AR(p) model that was fitted using AR. 
%	The required input arguments are all output of AR.
%
%	On output, the columns of S contain the model's normal modes,
%	and Sconf contains 68.27% confidence intervals for the
%	components of S.
%
%	per and tau are two-row matrices, in the first row containing
%	the period per(1,k) and damping rate tau(1,k) associated with 
%	the mode S(:,k). 68.27% confidence intervals for the estimated 
%	periods and damping rates are returned as the second row of
%	per and tau, respectively. 
%
%	exctn gives fractions indicating the relative importance 
%	(`excitation') of each mode. The sum of the excitations of all
%	modes equals one.
%
%	See also: AR, ARCONF

%       Modified 2-Sep-97
%       Author: Tapio Schneider
%		AOS Program, Princeton University
%		Princeton, NJ 08544-0710
%		tapio@splash.Princeton.EDU

% Dimension of state space
m = size(C,1);
p = size(A,2)/m; 			% order of process
if p <= 0 
   error('Order must be greater 0.'); 
end

% Only consider real parts of input arguments. These correspond to 
% parameter estimates, while imaginary parts may be confidence intervals
% if ARCONF has been called preciously
C = real(C);
A = real(A);

gamma = 1; 				% compute 68.27% confidence intervals
sqrt_gamma = sqrt(gamma);

% Assemble coefficient matrix of equivalent AR(1) process
A1 = [A; eye((p-1)*m) zeros((p-1)*m,m)];

% Eigenvalues and eigenvectors of AR(1)-matrix
[BigS,d] = eig(A1);  			% columns of BigS are eigenvectors
lambda   = diag(d);    			% vector containing eigenvalues

% Fix phase of eigenvectors such that each vector's
% real part and imaginary part are orthogonal
BigS = adjph(BigS);
% Return only last m components of each eigenvector
S = BigS((p-1)*m+1:p*m,:);

% Compute inverse of BigS for later use
BigS_inv = inv(BigS);

% Recover matrix Uinv which appears in the estimates' asymptotic 
% covariance matrix (Uinv is output of AR)
if (size(th,2) == m*p+1)
%  The intercept vector has been fitted by AR; to compute confidence
%  intervals for the modes, this vector is irrelevant. So, we don't need 
%  the first row and first column in Uinv, which correspond to elements 
%  in the intercept vector.   
   Uinv = th(3:size(th,1),2:size(th,2));

elseif (size(th,2) == m*p)
%  no intercept vector has been fitted
   Uinv = th(2:size(th,1),:);
else
   error('Input arguments of ARMODE must be output of AR.')
end

% Assemble asymptotic covariance matrix of estimator for coefficient matrix A
Sigma_A = kron(Uinv, C);

% Assemble noise covariance matrix of equivalent AR(1) process
BigCov = [C zeros(m,(p-1)*m); zeros((p-1)*m, p*m)];

% Get noise covariance matrix of decoupled system of relaxators and 
% oscillators
CovDcpld = BigS_inv*BigCov*BigS_inv';

% Compute period per, damping rate tau, and `excitation' exctn for each mode; 
% also get confidence intervals for per and tau
for j=1:m*p				% j counts the mode number
   a = real(lambda(j)); 		% real part of eigenvalue j
   b = imag(lambda(j)); 		% imaginary part of eigenvalue j
   abs_lambda_sq = abs(lambda(j))^2;  	% absolute value squared of eigenvalue j
   tau(1,j) = -2/log(abs_lambda_sq);    % damping rate of mode j

   % Get period of mode j (the following if-statements avoid warning messages 
   % due to possible division by zero)  
   if (b == 0 & a >= 0)
      per(1,j) = Inf;   		% purely real, nonnegative eigenvalue
   elseif (b == 0 & a < 0)
      per(1,j) = 2;     		% purely real, negative eigenvalue
   else      
      per(1,j) = 2*pi/abs(atan2(b,a)); 	% complex eigenvalue
   end

   % Compute excitation for mode j 
   exctn(j) = abs_lambda_sq^(p-1) * real(CovDcpld(j,j)/(1-abs_lambda_sq)); 

   % Assemble derivative of eigenvalue with respect to parameters in 
   % the coefficient matrix A 
   for k=1:m
       dot_lam(k:m:k+(m*p-1)*m) = BigS_inv(j,k).*BigS(1:m*p,j);
   end
   dot_a = real(dot_lam)'; 		% derivative real part of lambda(j)
   dot_b = imag(dot_lam)'; 		% derivative imaginary part of lambda(j)

   % Get confidence intervals for per and tau and store them in the 
   % respective second rows.
   % --------------------------------------------------------------
   % For tau: let phi be the derivative of the damping rate w.r. to 
   % parameters in A
   phi = tau(1,j)^2/(abs_lambda_sq)*(a*dot_a+b*dot_b);
   tau(2,j) = sqrt_gamma * sqrt(phi'*Sigma_A*phi);

   % Similarly for per: let phi be the derivative of the period w.r. to 
   % parameters in A
   phi = per(1,j)^2/(2*pi*abs_lambda_sq)*(b*dot_a-a*dot_b);
   per(2,j) = real( sqrt_gamma * sqrt(phi'*Sigma_A*phi));
end
% Give the excitation as `relative importance' that sums to one
exctn = exctn/sum(exctn);

% Compute confidence intervals for modes 
% --------------------------------------
% Shorthands for matrix products
XTX = real(BigS)'*real(BigS);
YTY = imag(BigS)'*imag(BigS);
XTY = real(BigS)'*imag(BigS);

% Need confidence intervals only for last m rows of BigS
row1 = (p-1)*m+1;			% first row for which confidence 
					% interval is needed

mp = m*p;                           	% dimension of equivalent AR(1) model
for k=1:mp        			% loop over columns of S
   for l=row1:mp  			% loop over rows of S
      % evaluate gradient of S_{lk}
      for ii=1:m
         for jj=1:mp
            % compute derivative w.r. to A(ii,jj)
            zsum = 0;
            zkkr = 0; 			% real part of Z_{kk}
            zkki = 0; 			% imaginary part of Z_{kk}
            for j=1:mp
               if (j ~= k) 		% sum up elements that appear in Z_{kk}
                  zjk  = BigS_inv(j,ii)*BigS(jj,k)/(lambda(k)-lambda(j));
                  zjkr = real(zjk);
                  zjki = imag(zjk);
                  zkkr = zkkr+zjki*(XTY(k,j)-XTY(j,k))-zjkr*(XTX(k,j)+YTY(k,j));
                  zkki = zkki+zjki*(YTY(k,j)-XTX(k,j))-zjkr*(XTY(k,j)+XTY(j,k));
                  zsum = zsum+BigS(l,j)*zjk;
               end
            end
            % now add Z_{kk}
            zkki = zkki/(XTX(k,k)-YTY(k,k));
            grad_S((jj-1)*m+ii) = zsum+BigS(l,k)*(zkkr+i*zkki);
         end     
      end 
      Sconf(l,k) = sqrt_gamma * ( sqrt( real(grad_S)*Sigma_A*real(grad_S)') ...
                              + i*sqrt( imag(grad_S)*Sigma_A*imag(grad_S)') );
   end
end

% Only return last m*p rows of Sconf
Sconf = Sconf(row1:m*p,:);


