function [z,x,varargout]=simlinsys(F,G,H,D,L,Q,R,M,u,xinit,NT,varargin)
% [z,x] = simlinsys(F,G,H,D,L,Q,R,M,u,xinit,NT,...)
% [z,x,statenoise,obsnoise] = simlinsys(F,G,H,D,L,Q,R,M,u,xinit,NT,...)
%
% Sample a linear discrete system
%
%  x(k+1) = F(k) x(k) + G u(k) + L(k) w(k)
%  z(k) = H(k) x(k) + D u(k) + n(k)
%
% with time-dependent state-space matrices F(k), G(k), H(k), D(k),
% L(k). The noise statistics are given by
%
%  w ~ N(0,Q)
%  n ~ N(0,R)
%  E[w*n'] = M
%
% Note that time-dependent cross correlations in M are not sampled
% correctly!
%
% The matrix u must have as many rows as time samples and as many
% columns as system inputs. If either u or x0 are empty, zero
% input/initial conditions are assumed. Note that the output vectors x
% and z start from time t=0, i.e.,
%  x(0) = x0;
%  z(0) = H*0;
%
% Output arguments statenoise and obsnoise give the realization of the
% noise processes w(k) and n(k) used to create the sample.
%
% All optional arguments are passed to simgauss.
%
% See also SIMGAUSS, SIMLTISYS

% Notes:
% - For practical reasons all the vectors are denoted as column
%   vectors during calculation (i.e. time increases in 2nd dimension)
%   and are transposed before output.
% - xinit is recorded as the first element in x, i.e.,
%    x(:,1) = x0;
%    z(:,1) = H*x0;
% - The observation input is instantaneous, i.e.,
%    z(:,n) = H*x(:,n) + D*u(:,n)
% - The state input is delayed by one timestep, i.e.,
%    x(:,n+1) = F*x(:,n) + G*u(:,n)
% All of the above is (tested to be!) consistent with the behaviour of
% lsim and initial.

% Built-in sampling functions use ltitr:
%
% [...]
%
% LTITR(A,B,U,X0) can be used if initial conditions exist.
% Here is what it implements, in high speed:
%
% for i=1:n
%   x(:,i) = x0;
%   x0 = a * x0 + b * u(i,:).';
% end
% x = x.';


nArg = 2;                               % number of fixed output arguments

%% parse and test input
vargs = varargs2cell(varargin);  % this is directly passed to simgauss

if size(M,3)~=1 || size(Q,3)~=1 || size(R,3)~=1
    error('time dependent noise correlations are not supported')
end

% check model
% h = checkdim(F,G,H,D,L,Q,R,M,z,u,xinit,pinit)
if ~checkdim(F,G,H,D,L,Q,R,M,[],u,xinit,[])
    error('Model dimensions are not consistent');
end

% check if inputs are time-dependent
% x_TD=true <> x is time dependent!
if size(F,3)==1, F_TD=false; f=F; else F_TD=true; end
if size(G,3)==1, D_TD=false; g=G; else G_TD=true; end
if size(H,3)==1, H_TD=false; h=H; else H_TD=true; end
if size(D,3)==1, G_TD=false; d=D; else D_TD=true; end
if size(L,3)==1, L_TD=false; l=L; else L_TD=true; end

%% preallocate and initialize vectors
[os ss ~]=size(H);                      % state dimensions
si = size(G,2);                         % input dimension
sn = size(Q,1);                         % statenoise dimension
on = size(R,1);                         % obsnoise dimension

% set defaults for empty vectors
if isempty(xinit)
    xinit = zeros(1,ss);
end
if isempty(u)
    u = sparse(NT,si);
end
% make xinit into column vector
if size(xinit,1) < size(xinit,2)
    xinit = xinit';
end

%% sample timeseries

% create noise input if necessary
if ~(isempty(Q) || isempty(R) || isempty(L))
    % for noisy systems
    if isempty(M)
        M = zeros(size(Q,1),size(R,2));
    end
    % full covariance
    fullcov = [[Q,M];[M',R]];
    fullnoise = simgauss(0,fullcov,NT,'testpos',true,vargs{:});
    % separate noise processes
    statenoise = fullnoise(1:sn,:);
    obsnoise = fullnoise(sn+1:on+sn,:);
else
    % no noise
    statenoise = sparse(max(1,sn),NT);
    obsnoise = sparse(max(1,on),NT);
end

% integrate process (analog to ltitr)
x = zeros(ss,NT);                       % state output
z = zeros(os,NT);                       % observation output

x0 = xinit;
for tt = 1:NT
    if F_TD, f=F(:,:,tt); end;
    if G_TD, g=G(:,:,tt); end;
    if H_TD, h=H(:,:,tt); end;
    if D_TD, d=D(:,:,tt); end;
    if L_TD, l=L(:,:,tt); end;

    x(:,tt) = x0;
    z(:,tt) = h*x0 + d*u(tt,:).' + obsnoise(:,tt);
    x0 = f*x0 + l*statenoise(:,tt) + g*u(tt,:).';
end
x = x.'; z = z.';

%% setting output arguments
nout = nargout-nArg;
if nout == 1
    varargout(1) = {statenoise.'};
elseif nout == 2
    varargout(1:2) = {statenoise.',obsnoise.'};
end