% Listing 4.1 of the book
% For Figures 4.9 and 4.10

% 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(1275,'prbs',[0 1/5],[-1 1]);
xk = sim(mod_p,uk); mod_p.Noisevariance = var(xk)/10;
yk = sim(mod_p,uk,simOptions('AddNoise',true));
%% Collect the input-output data
Z = iddata(yk,uk,1); Ztrain = detrend(Z,0);
%% Estimate the impulse response 
mod_fir = impulseest(detrend(Ztrain,0));
[irest,kvec,~,sd_ir] = impulse(mod_fir,10);
figure; stem(kvec,irest); hold on
plot(kvec,sd_ir,'k--',kvec,-sd_ir,'g--');

% Compute the true values and plot them
ir_act = filter([0 0.6 -0.2],[1 -0.5],[1 zeros(1,9)]);
stem((0:9),ir_act,'ro','filled');

%% Compute step response estimates
[stepres,kvec] = step(mod_fir,15);
figure; plot(kvec,stepres,'b-'); hold on

% Compute the true values and plot them
sr_act = filter([0 0.6 -0.2],[1 -0.5],ones(15,1));
plot((0:14),sr_act,'r--');
plot((0:14),sr_act,'ro','MarkerFaceColor','red');