% Listing 4.3 of the book 
% "Principles of System Identification: Theory and Practice"

% Arun K. Tangirala
% April 01, 2015
% Freely distributed under no warranty

%% Create an IDPOLY object for the process
mod_p = idpoly(1,[0 0.6 -0.2],1,1,[1 -0.5],'Noisevariance',1);

% Generate input and simulate the process
uk = idinput(2555,'prbs',[0 0.2],[-1 1]);
xk = sim(mod_p,uk); mod_p.Noisevariance = var(xk)/10;
yk = sim(mod_p,uk,'Noise');

% Collect the input-output data
dataset = iddata(yk,uk,1);
datatilde = detrend(dataset,0);

%% Estimate a freely parametrized state-space model
Gss_hat = pem(datatilde);
% Print the state-space model
present(Gss_hat)

% Compare predictions
figure; compare(datatilde,Gss_hat);

% Fit an observability canonical form
Gss_can = pem(datatilde,'ssparameterization','canonical');
present(Gss_can);

%% Fit a structured state-space model

% Initialization
A = [0 1; 0 0.4]; B = [0.4 ; 0.3]; 
C = [1 0]; D = 0; K = [0 ; 0]; 
M = idss(A,B,C,D,K);

% Providing structural information
M.As = [0 1; 0 NaN]; M.Bs = [NaN ; NaN];
M.Cs = [1 0]; M.Ds = 0; M.Ks = [0 ; 0];
M.x0s = [0 ; 0];

% Estimating the structured model
Mhat = pem(datatilde,M);
present(Mhat)