% Listing 16.8 for Example 16.6 of the book 
% "Principles of System Identification: Theory and Practice"

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

%% Generate realizations and their psds
Nreal = 1000; Nobs2 = 2000; Nobs1 = Nobs2/8;

vk_psd1 = zeros(Nobs1/2,Nreal);
vk_psd2 = zeros(Nobs2/2,Nreal);

for i = 1:Nreal
    vk = randn(Nobs2,1);
    [t_psd1,wvec1] = periodogram(vk(1:Nobs1),rectwin(Nobs1),Nobs1,'twosided');
    [t_psd2,wvec2] = periodogram(vk,rectwin(Nobs2),Nobs2,'twosided');
    vk_psd1(:,i) = t_psd1(1:end/2);
    vk_psd2(:,i) = t_psd2(1:end/2);
end

%% Plot the results
true_psd = 1/(2*pi);
figure
h1 = subplot(321);
plot(h1,wvec1(1:end/2),vk_psd1(:,34)); hold on
plot([wvec1(1) wvec1(end/2)],[1 1]/(2*pi),'r--','linewidth',2);
h2 = subplot(322);
plot(h2,wvec2(1:end/2),vk_psd2(:,34)); hold on
plot([wvec2(1) wvec2(end/2)],[1 1]/(2*pi),'r--','linewidth',2);
h3 = subplot(323);
plot(h3,wvec1(1:end/2),mean(vk_psd1,2)); hold on
plot([wvec1(1) wvec1(end/2)],[1 1]/(2*pi),'r--','linewidth',2);
axis tight
set(gca,'ylim',[0.1 0.2])
h4 = subplot(324);
plot(h4,wvec2(1:end/2),mean(vk_psd2,2)); hold on
plot([wvec2(1) wvec2(end/2)],[1 1]/(2*pi),'r--','linewidth',2);

% Plot distributions of PSD estimates
h5 = subplot(325);
psdhat0 = vk_psd2(1,:)/true_psd;
hist(h5,psdhat0)
h6 = subplot(326);
psdhat1 = 2*vk_psd2(401,:)/true_psd;
hist(h6,psdhat1)

%% Fit chi-square distributions to two select estimates      
psdfit0 = fitdist(psdhat0','gamma');
psdfit1 = fitdist(psdhat1','gamma');