% Listing 11.4 of the book for Example 11.7
% "Principles of System Identification: Theory and Practice"

% Arun K. Tangirala
% April 02, 2015
% Freely distributed under no warranty

%% Generate data
N = 1000; uk = randn(N,1);

% Simulate process
Gd = tf([2 0],[1 -0.5],1);
yk = lsim(Gd,uk,(0:N-1)) + 0.5*randn(N,1);

%% Estimate FRF
[gamma_yu,w] = cpsd(yk,uk);
[gamma_uu,w] = cpsd(uk,uk);
Gwhat = gamma_yu./gamma_uu;

%% Compare estimated and true FRFs

Gw = squeeze(freqresp(Gd,w));
figure; plot(w,[abs(Gw) abs(Gwhat)]);