tt = input('white noise input [1], step [0] ');
if tt==1, u = randn(1000,1); end
if tt==0, u = ones(1000,1); end
v = input('noise variance [.1] ');
y = dlsim([0 1 -1],[1 -1.5 0.7],u);
e = randn(size(y));
y = y + v*e;
figure(1),
subplot(211),plot(u),title('input');
subplot(212),plot(y),title('output');
[N,m]=size(u);
[N,l]=size(y);
s=input('number of block rows ');
% Construct Hankel matrices
U=u(1:N-s+1,:)';
Y=y(1:N-s+1,:)';
for i=2:s,
    U = [U;u(i:N-s+i)'];
    Y = [Y;y(i:N-s+i)'];
end
% compute RQ factorization 
r=triu(qr([U' Y']))';
R22 = r(s*m+1:s*(m+l),s*m+1:s*(m+l));
% compute svd
[uu,ss,vv]=svd(R22); figure(2),semilogy(diag(ss),'xr','LineWidth',2,'MarkerSize',12),
n = input('what is your choice of order: ');
Ct = uu(1:l,1:n);
% compute least squares 
At = uu(1:l*(s-1),1:n)\uu(l+1:s*l,1:n);
eig(At),roots([1 -1.5 0.7])