% Listing 21.4 for estimating a BJ model (Example 21.13)
% "Principles of System Identification: Theory and Practice"

% Arun K. Tangirala
% December 06, 2015
% Freely distributed under no warranty

%% Data generation

% Create the plant and noise model objects 
p_bj = idpoly(1,[0 0 0.6 -0.2],[1 0.2],[1 -0.4],[1 -0.5],'Noisevariance',0.1);

% Create input sequence 
uk = idinput(2046,'prbs',[0 0.4],[-1 1]); 
% Simulate the process 
yk = sim(p_bj,uk,simOptions('Addnoise',true));

% Build iddata objects and remove means
Zk = iddata(yk,uk,1); Ztrain = detrend(z,0); 

%% Estimate the BJ model

% Compute IR for time-delay estimation 
mod_fir = impulseest(Ztrain);
figure; impulseplot(mod_fir,'sd',3);
% Time-delay = 2 samples
% Estimate BJ model (assume known orders ) 
nb = 2; nc = 1; nd = 1; nf = 1; nk = 2; 
mod_bj = bj(zd,[nb nc nd nf nk]);
% Present the model
present(mod_bj) 

%% Residual analysis

% Check the residual plot	 
figure ; resid(Ztrain, mod_bj);