function [singular_values] = mkhank0(G,nlines,fs,q,d)
%
%  Name: mkhank0.m
%
%  Usage: function [singular_values] = mkhank0(G,nlines,fs,q,d)
%
%  Description:
%  Modified verision of era.m.
%  Allows for svd of H(0) for selecting model order for ERA.
%
%  Applies the Eigensystem Realization Algorithm (ERA) to the 
%  impulse response from the FRF matrix 'G'
%
%  q and d are the ERA Hankel parameters
%
%  Returns the Hankel Matrix H(0) singular values vector
%
%  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 EGS19960826
%
%  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(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
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,q,d);
fprintf(1,' Done. \n');

%  Decompose Hankel Matrix at k=0
fprintf(1,'Decomposing Hankel Matrix at k=0...');
[U0,S0,V0] = svd(H0);
singular_values = diag(S0);
fprintf(1,' Done. \n\n');

return
