function [A,B,C,D,fs_out,emac,msv,H0,H1] = era(G,nlines,fs,nx,q,d)
%
%  ERA
%
%  [A,B,C,D,fs_out,emac,msv] = era(G,nlines,fs,nx,q,d)
%
%  Applies the Eigensystem Realization Algorithm (ERA) to the 
%  impulse response from the FRF matrix 'G'
%
%  nx is the number of states to use (2 states per mode)
%  q and d are the ERA Hankel parameters
%
%  Returns the state-space model matrics 'A,B,C,D' and the
%   'revised' sample rate fs_out (due to zero-padding)
%
%  Also returns the observability matrix 'V' and the
%   controllability matrix 'W'
%
%  If nx < 0, will plot Hankel singular values and
%   prompt for a cut-off

%  Reference for ERA: "An Eigensystem Realization Algorithm for
%    Modal Parameter Identification and Model Reduction," 
%    Juang, J-N, and Pappa, R.S., Journal of Guidance, Control, and
%    Dynamics, Vol. 8, No. 5, pp. 620-627, 1985.

%  Version SWD961118
%
%  Copyright (c) 1996 by Scott W. Doebling
%	Los Alamos National Laboratory
%	ESA-EA, doebling@lanl.gov
%

%  Define variables and allocate space for matrices

nresp = size(G,1) / nlines;
nref = size(G,2);

GBlk = zeros(size(G));

Y = zeros(nresp*(nlines*2),nref);

A = zeros(nx,nx);
B = zeros(nx,nref);
C = zeros(nresp,nx);
D = zeros(nresp,nref);

H = zeros(q*nresp,d*nref);

padex =  ceil(log(nlines)/log(2));

npad = 2 * (2^padex - nlines ) + 1;

nlines_out = 2^(padex+1);

fs_out = fs/((nlines-1) * 2) * nlines_out;


%  Convert G to Block form with Full Spectrum

disp('Putting FRF in Block Form...')

for j = 1:nref,

   for i = 1:nresp,

       Gij =  G(((i-1)*nlines + 1):(i*nlines),j);

       GBlk(  i:nresp:(nlines_out-1)*nresp+i , j ) = [Gij;zeros(npad,1);conj(Gij((nlines-1):-1:1))];

   end

end


%  Convert GBlk to Block-Form Impulse Response 'Y'

disp('Converting FRF to Impulse Response...')

for j = 1:nref,

   for i = 1:nresp,

       Y(  i:nresp:(nlines_out-1)*nresp+i , j ) = real(ifft(GBlk( i:nresp:(nlines_out-1)*nresp+i, j)));

   end

end

%  Form Hankel Matrix at k=0

disp('Forming Hankel Matrix at k=0...')

H0 = hankera(Y,nlines_out,0,q,d);

disp('Decomposing Hankel Matrix at k=0...')

[U0,S0,V0] = svd(H0);

if nx<0,
   semilogy(diag(S0))
   nx = input('Please enter number of states to save:');
end

U0N = U0(:,1:nx);
S0N = S0(1:nx,1:nx);
V0N = V0(:,1:nx);

%  Form Hankel Matrix at k=1

disp('Forming Hankel Matrix at k=1...')

H1 = hankera(Y,nlines_out,1,q,d);


%  Form A, B, C, D, V, W

disp('Forming A,B,C,D...')
SNsqrt = sqrt(S0N);
SNinv = sqrt(inv(S0N));

A = SNinv * U0N' * H1 * V0N * SNinv;
B = SNsqrt * V0N(1:nref,:)';
C = U0N(1:nresp,:) * SNsqrt;
D = Y(1:nresp,:);

% V = U0N(1:nresp,:) * SNsqrt;
% W = SNsqrt *  V0N(1:nref,:)';

V = U0N * SNsqrt;
W = SNsqrt *  V0N';

if nargout > 5, 

  disp('Computing EMAC...')

  [emac,msv] = emaccomp(A,V,W,q,d,fs_out);

end

return
