w2r.m

来自「A Matlab toolbox for exact linear time-i」· M 代码 · 共 56 行

% W2R - From time series to kernel representation of the MPUM.
%
% R = w2r(w,lmax,m,n,tol)
%
% W - time series
% LMAX  - upper bound for the system lag
% M - optional number of inputs 
% N - optional system order
% TOL - tolerance for order selection, default 1e-7
% R - 3d array of parameter of a kernel representation
%     R0 = R(:,:,1), ... ,Rlmax = R(:,:,lmax+1)

function R = w2r(w,lmax,m,n,tol)

dw  = size(w,2); % # of variables
lmax_1 = lmax+1;
T = size(w,1);
j = T - lmax;

if nargin < 5 | isempty(tol)
  tol = 1e-7; % default
end

% Left kernel of the Hankel matrix
r = triu(qr(blkhank(w,lmax_1,j)'))';
[U,S,V] = svd(r);
s = diag(S);
% Determine the order
if nargin < 4 | isempty(m) | isempty(n)
  rk = sum( s > tol);
  if exist('m') == 1 & exist('n') ~= 1 
    n = rk - lmax_1*m;
  elseif exist('m') ~= 1 & exist('n') == 1 
    m = (rk - n) / lmax_1;
    if ~isinteger(m)
      warning('The specified order is ignored.')
      m = floor(rk / lmax_1);
      n = rk - lmax_1 * m; 
    end
  else
    m = floor(rk / lmax_1);
    n = rk - lmax_1 * m; 
  end  
else
  rk = lmax_1*m + n;
end
R = U(:,rk+1:end)';
g = size(R,1);

% Convert R to a 3d array  
R_ = R;
R  = zeros(g,dw,lmax_1);
for i = 1:lmax_1
  R(:,:,i) = R_(:,(i-1)*dw+1:i*dw);
end