function [a, b, c, d] = subspace(vmat, m, p, n, q)

% subspace-based identification for MIMO systems from frequency responce data
%  see the following ref. for detail:
%    @article{McKelvey:AC:1996,
%      author = "Tomas McKelvey and Huseyin Akcay and Lennart Ljung",
%      title = "Subspace-Based Multivariable System Identification from Frequency Response Data",
%      journal = ieeetac,
%      volume = "41",
%      number = "7",
%      pages = "960--979",
%      year = "1996"
%    }
%
% input
%  vmat : varying matrix of frequency responce data
%  m : input dimension of identifying plant
%  p : output dimension of identifying plant
%  n : order of the plant
%  q : design parameter which must satisfy M >= q + n
% output
%  [a, b, c, d] : state-space realization of identified plant
%
% Note: continuous-time frequency response is firstly discretized because the subspace based method is originally for discrete-time plant identification

w = getiv(vmat);
M = size(w,1);

T = 0.001; % sampling period for discretization
wd = 2 * atan(w * T / 2); % discretized angular frequency

Wm = [];
for l=0:q-1
  tmp = [];
  for k=1:M
    tmp = [tmp, exp(i*l*wd(k)) * eye(m)];
  end
  Wm = [Wm; tmp];
end
Wm = Wm / sqrt(M);

G_tmp = [];
for k=1:M
  G_tmp = mdiag(G_tmp, xtracti(vmat,k,1));
end

G = Wm * G_tmp;

tmp1 = [real(Wm), imag(Wm); real(G), imag(G)];
[Q, R] = qr(tmp1');
tmp = R';
R22 = tmp((q*p)+1:q*(p+m),(q*p)+1:2*m*M);
[U, S, V] = svd(R22);
Us = U(:,1:n);
Ss = S(1:n,1:n);

J1 = [eye((q-1)*p), zeros((q-1)*p,p)];
J2 = [zeros((q-1)*p,p), eye((q-1)*p)];
J3 = [eye(p), zeros(p,(q-1)*p)];

Ad = pinv(J1*Us) * J2*Us;
Cd = J3*Us;

X = [];
G0 = [];
for k=1:M
  tmp = [Cd/(exp(i*wd(k)) * eye(n) - Ad), eye(p)];
  X = [X; tmp];
  G0 = [G0; xtracti(vmat,k,1)];
end

tmp = pinv([real(X); imag(X)]) * [real(G0); imag(G0)];
Bd = tmp(1:n,:);
Dd = tmp(n+1:n+p,:);

[a, b, c, d] = bilin(Ad, Bd, Cd, Dd, -1, 'Tustin', T);
