function [H0,Y,U0,S0,V0,hERAconva,fs_out,nresp,nref] = mkCEhank(G,nlines,fs,samples)
%
%  Name: mkCEhank.m
%
%  Usage: function [H0,Y,U0,S0,V0,heraconva,fs_out,nresp,nref] = mkCEhank(G,nlines,fs,samples)
%
%  Description:
%   Modified verision of era.m.
%   Allows for svd of H(0) for selecting model order for CE.
%
%   Applies the Eigensystem Realization Algorithm (ERA) to the 
%   impulse response from the FRF matrix 'G'
%
%   samples is the number of time samples to use from the impulse response
%
%   Returns the Hankel Matrices H0 and the
%   'revised' sample rate fs_out (due to zero-padding)
%
%  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 EGS19960815
%
%  Author(s): Scott Doebling, doebling@lanl.gov
%             Phillip Cornwell, cornwell@lanl.gov
%	      Erik Straser, estraser@lanl.gov
%
%  Institution: Los Alamos National Laboratory, ESA-EA
%
%  Copyright (c) 1996, The Regents of the University
%       of California. (See top level of code for
%       full copyright statement.)


%  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);

H0 = zeros(samples*nresp,samples*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

fprintf(1,'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

fprintf(1,' Done. \n');


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

fprintf(1,'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

fprintf(1,' Done. \n');


%  Form Hankel Matrix at k=0

fprintf(1,'Forming Hankel Matrix at k=0...');

H0 = hankera(Y,nlines_out,0,samples,samples);

fprintf(1,' Done. \n');


%  Decompose Hankel Matrix at k=0

fprintf(1,'Decomposing Hankel Matrix at k=0...');

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

fprintf(1,' Done. \n');


%  Plot singluar values of H(0) versus pole number

hERAconva=axes('position',[.1 .2 .7 .7]);
   semilogy(diag(S0))

fprintf(1,'\n');

return
