function Gr=ps42_6242_2004(r,h)
% function Gr=ps42_6242_2004(r,h)
% 
% implements algorithm of problem 4.2 of 6.242/2004

if nargin<1, r=5; end
if nargin<2, h=10; end
if length(h)==1, h=[zeros(1,2*h+1) (1:h)/h ((h-1):-1:1)/h]; end

h=h(:);
N=length(h);
H=hankel(h);
W=zeros(N);
L=zeros(N);
W(1,1)=2;
L(1,1)=-1;
for k=2:N,
    W(k,k)=4;
    W(k,k-1)=1;
    W(k-1,k)=1;
    L(k,k)=-1;
    L(k-1,k)=1;
end
e=[1 zeros(1,N-1)];
S=chol(W)';
[V,D]=eig(S'*H*S);
[dd,ii]=sort(-abs(diag(D)));
g=S*V(:,ii(1:r));
q=V(:,ii(1:r))'*inv(S);
g=g*inv(q*H*g);
Gr=ss(q*L*H*g,q*h,e*H*g,0);


close(gcf)
tt=0:(N-1);
hr=squeeze(impulse(Gr,tt));
h0=squeeze(impulse(ss(L,h,e,0),tt));
subplot(2,1,1); plot(tt,h,tt,hr);grid
subplot(2,1,2); plot(tt,hr,tt,h0); grid
disp(['Model error:  ' num2str(max(abs(h-hr)))])
disp(['Match error:  ' num2str(max(abs(h0-hr)))])