% Listing for Example 22.5 of the book on
% "Principles of System Identification: Theory and Practice"

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

% SCRIPT TO DEMONSTRATE THE IMPORTANCE OF PRE-WHITENING IN RESIDUAL
% CORRELATION ANALYSIS
%
% Chapter 22 of the text
%

% Arun K. Tangirala
% May 16, 2014

%% Data generating process
% DGP
mod_dgp = idpoly(1,[0 0 0.3 1 0.5 0.2],1,1,1,'Noisevariance',1,'Ts',1);
% Generate input
%uk = idinput(2044,'prbs',[0 1/4],[-1 1]); N = length(uk);
uk = idinput(2555,'prbs',[0 1/5],[-1 1]); N = length(uk);
% Simulate the process
xk = sim(mod_dgp,uk);
mod_dgp.Noisevariance = var(xk)/10;
yk = sim(mod_dgp,uk,simOptions('Addnoise',true));

%% Data analysis and identification

% Plot the input-output data
Zk = iddata(yk,uk,1);
Ztrain = detrend(Zk,0);

mod_fir2 = arx(Ztrain,[0 2 2]);
mod_fir4 = arx(Ztrain,[0 4 2]);

%% Residual analysis

% Without pre-whitening the input
resid(Ztrain,mod_fir2);

% ACF of the input
acf(uk,20);

% After pre-whitening the input
Udata = iddata(uk,[],1);
mod_tsuk = armax(Udata,[0 4]);
figure
resid(mod_tsuk,Udata);
uk_resid = pe(mod_tsuk,uk);
perr = pe(Ztrain,mod_fir2);
Lmax = 20;
[ccf_epwu,lags] = xcov([perr.y uk_resid],Lmax,'coeff');

% Plot the CCF and the significance levels
figure
stem(lags,ccf_epwu(:,2),'Markerfacecolor','blue');
hold on
plot([-Lmax Lmax],ones(1,2)*2.58/sqrt(N),'r--');
plot([-Lmax Lmax],-ones(1,2)*2.58/sqrt(N),'r--');
ylabel('CCF','fontsize',16,'fontweight','bold'); 
xlabel('Lags','fontsize',16,'fontweight','bold');
%title('With pre-whitened series','fontsize',16,'fontweight','bold');
prepare_figure

% Repeat for correct model
perr4 = pe(Ztrain,mod_fir4);
[ccf_epwu4,lags] = xcov([perr4.y uk_resid],Lmax,'coeff');
% Plot the CCF and the significance levels
figure
stem(lags,ccf_epwu4(:,2),'Markerfacecolor','blue');
hold on
plot([-Lmax Lmax],ones(1,2)*2.58/sqrt(N),'r--');
plot([-Lmax Lmax],-ones(1,2)*2.58/sqrt(N),'r--');
ylabel('CCF between residuals and pre-whitened input','fontsize',16,'fontweight','bold'); 
xlabel('Lags','fontsize',16,'fontweight','bold');
prepare_figure


