% Listing 21.3 for estimating an OE model
% "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_oe = idpoly(1,[0 0 0.6 -0.2],1,1,[1 -0.5],'Noisevariance',0.1);

% Create input sequence 
uk = idinput(2555,'prbs',[0 0.2],[-1 1]); 

% Simulate the process 
yk = sim(p_oe,uk,simOptions('Addnoise',true));

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

%% Estimate the OE model

% Compute IR for time-delay estimation 
mod_fir = impulseest(Ztrain);
figure; impulseplot(mod_fir,'sd',3);
% Time-delay = 2 samples
% Estimate OE model 
nb = 2; nf = 1; nk = 2; 
mod_oe = oe(Ztrain,[nb nf nk]);

%% Residual analysis

% Check the residual plot	 
figure ; resid(Ztrain, mod_oe);
% Analyze the model
present(mod_oe) 