% Listing 16.9 for Example 16.7 of the book 
% "Principles of System Identification: Theory and Practice"

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

%% Generate WN sequence and compute theoretical PSD
N = 1024;
ek = randn(N,1); ekd = ek - mean(ek);
th_psd = 1/(2*pi);

%% Periodogram and their modifications

% Periodogram
[psd_hat,wv] = periodogram(ekd,rectwin(N),N,'twosided');
psd_per = psd_hat(1:end/2);

% Welch's averaged periodogram
[psd_hatw,wv2] = pwelch(ekd,hanning(128),64,[],'twosided');
psd_welch = psd_hatw(1:end/2);

% B-T method
L = 127; 
[acf_vv,lags] = xcov(ekd,L,'biased');
win_fn = bartlett(length(lags));
wacf_vv = win_fn.*acf_vv;
psd_bt = fft(wacf_vv(128:end),256) + 256*ifft(wacf_vv(128:-1:1),256) - wacf_vv(128);

% Neglect imaginary parts due to rounding errors
psd_bt = real(psd_bt(1:end/2))/(2*pi); 

% MTM 
[psd_hatm,wv3] = pmtm(ekd,[],[],'twosided');
psd_mtm = psd_hatm(1:end/2);

figure
subplot(2,2,1)
plot(wv(1:end/2),psd_per,wv([1 end/2]),th_psd*[1 1],'r--'); 
subplot(2,2,2)
plot(wv2(1:end/2),psd_welch,wv2([1 end/2]),th_psd*[1 1],'r--')
subplot(2,2,3)
plot(wv2(1:end/2),psd_bt,wv2([1 end/2]),th_psd*[1 1],'r--')
subplot(2,2,4)
plot(wv3(1:end/2),psd_mtm,wv3([1 end/2]),th_psd*[1 1],'r--')