% Listing 20.2 for Example 20.3 of the book 
% "Principles of System Identification: Theory and Practice"

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

%% Data generation

% DGP
proc_armax = idpoly([1 -1.2 0.35],[0 0 2],[1 0.3],'Noisevariance',1);
% Simulate the process with a PRBS input
uk = idinput(2046,'prbs',[0 0.4],[-1 1]); % Full-length PRBS
yktrue = sim(proc_armax,uk); % Noise-free simulation
% Adjusting the noise variance s.t. SNR = 10
mult_fact = (1 + 0.3^2)*(1 + 0.35)/((1 + 0.35)*(1 - 0.35^2) - 1.2^2 + 0.35*1.2^2);
proc_armax.Noisevar = 0.1*var(yktrue)/mult_fact;
% Simulate with noise
yk = sim(proc_armax,uk,simOptions('Addnoise',true));

% Processing the data
dataset = iddata(yk,uk,1);
datatrain = detrend(dataset,0);

%% Estimating the FRF with ETFE
Gwhat = etfe(datatrain);
wvec = Gwhat.Frequency;
[mag_act,phase_act] = bode(proc_armax,wvec);

%% Plots

% Plot the estimated vs. true magnitude and phase response
figure
semilogx(wvec,squeeze(abs(Gwhat.Response)),'b-',wvec,squeeze(mag_act),'r--');
figure
semilogx(wvec,phase(squeeze(Gwhat.Response))*180/pi,'b-',wvec,squeeze(phase_act),'r--');
