function [a,b,d,k]=harxmr(y,u,n_arx,n)
% [a,b,d,k]=harxmr(y,u,n_arx,n)
% Purpose
% Model identification by first identification of a higher order
% OLS/ARX model followed by model reduction using Hankel matrix
% realization theory.

[N,nu]=size(u);
[N,ny]=size(y);

% identify a higher order ARX/OLS model of order n_arx >> n
th=arx([y u],[n_arx,n_arx,1]);
[A,B,D,E,K]=th2ss(th);            
%K=K*0;

% Calculate impulse response matrices for the higher order model
L=n_arx;
hw=zeros(2*L,nu+ny);
hw(1,:)=D*[B K];
for i=2:2*L
    hw(i,:)=D*(A^(i-1))*[B K];
end

% Make Hankel matrix Hn=OL*CJ
Hn=hank_m(hw,ny,nu+ny,L,L);

% SVD to find OL and CJ, and model matrices D,B and K
[U,S,V]=svd(Hn);
U1=U(:,1:n); S1=S(1:n,1:n); V1=V(:,1:n);
O=U1;
C=S1*V1';
d=O(1:ny,:); b=C(:,1:nu);
k=C(:,nu+1:nu+ny);

% Make Hankel matrix Ha=OL*A*CJ and calculate A
Ha=hank_m(hw(ny+1:2*L*ny,:),ny,nu+ny,L,L);
a=pinv(O'*O)*O'*Ha*C'*pinv(C*C');

% END harxmr