function [A,D,OL]=idsys(Y,U,L,g,J)
% [A,D,OL]=idsys(Y,U,L,g,J)
% Purpose. Using the shift invariance method to estimate n, A and OL

[Uf,Yf,Up,Yp,YL,UL]=mat4oppg13(Y,U,L,g,J);
 
 %
  W=[Up;Yp;Uf];  % Matrix to remove noise
 K=size(W,2);   % Number of columns
 
 Zd=Yf*W'*pinv(W*W')*W;
 
 
 % Calculate projection Ufp such that Uf*Ufp=0
 Ufp=eye(K) - Uf'*pinv(Uf*Uf')*Uf;
 
 % Calculate the main projection Z
 Z=Zd*Ufp;
 
 % Using SVD
 [U,S,V]=svd(Z);
 diag(S)'  % Display singular values
 n=input('System order = ?');
 
 % Partition the SVD
 U1=U(:,1:n); S1=S(1:n,1:n); V1=V(:,1:n);
 
 O=U1; % The extended observability O_L+1
 
 % Use shift invariance technique to calculate A
 
 OL=O(1:L,:);
 OLA=O(2:end,:);
 
 A=inv(OL'*OL)*OL'*OLA;
 D=OL(1);