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

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

%% Generate signal

f0 = 0.1; f1 = 0.16; N = 240; kvec = (0:N-1)';
xk = sin(2*pi*f0*kvec) + 0.27*sin(2*pi*f1*kvec);

% Add measurement noise s.t. SNR = 2
yk = xk + randn(N,1)*std(xk)/sqrt(2); 
yk = yk - mean(yk);

%% Compute power spectral density

% FFT
ykf = fft(yk); yk_psd = abs(ykf(1:end/2)).^2;

% Plot the PSD of measurement
fvec = (0:1/N:0.5-1/N);
figure; stem(fvec,yk_psd/N,'Markersize',0);

% Window the measurement and plot the PSD
% Hann
win_hann = hanning(N); ykw1 = win_hann.*yk;
ykw1f = fft(ykw1); ykw1_psd = abs(ykw1f(1:end/2)).^2;
figure; stem(fvec,ykw1_psd/sum(win_hann.^2),'Markersize',0);
% Blackman
win_black = blackman(N); ykw3 = win_black.*yk;
ykw3f = fft(ykw3); ykw3_psd = abs(ykw3f(1:end/2)).^2;
figure; stem(fvec,ykw3_psd/sum(win_black.^2),'Markersize',0);