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

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

%% Generate the measurement
kvec = (0:499)';
vk = sin(2*pi*0.1*kvec) + sin(2*pi*0.25*kvec) + randn(500,2.25);

%% Compute periodogram
N = length(vk);
[Pxx,wvec] = periodogram(vk-mean(vk),[],N,1);
figure; plot(wvec,Pxx/sum(Pxx));

% Take Fourier Transform
vkf = fft(vk);
magvkf = abs(vkf(1:end/2)); phasevk = phase(vkf(1:end/2));

%% Estimate the signal
% Zero out the contributions (assumed to be) due to noise
magvkf2 = zeros(length(magvkf),1);
magvkf2(50:52) = magvkf(50:52); 
magvkf2(125:127) = magvkf(125:127);
vkfmod = magvkf2.*exp(i*phvkf);
vkfmod2 = [vkfmod ; 0 ; flipud(conj(vkfmod(2:end)))];
% Estimate the signal
vkhat = ifft(vkfmod2);

%% Plot against the true deterministic signal
xk = sin(2*pi*0.1*kvec) + sin(2*pi*0.25*kvec);
figure; plot((0:99),xk(1:100),'r-',(0:99),vkhat(1:100),'b--');