% Example 8.2.2 of Hayes, Statistical DSP
% Periodogram of a Sinusoid in Noise

phi=(rand-1/2)*2*pi;  % generating the random phase
sigma=1;
f0=50;  % hidden sine freq
T=1/(2*(f0+2*f0));  % sampling frequency 
w0=2*pi*f0;
A=1;
N=16;   % no of samples

X=[];
res=4*1024;  % dft resolution (implemented with zero padding), ensure res >= N
for i=1:500
x=A*sin(T*[0:N-1]*w0+phi)+randn(1,N)*sigma;

X=[X;(1/N)*abs(fft(x,res)).^2];

end

varX=var(X,1);
Xmean=mean(X,1);

% now plot estimated PSD for input sequence
plot([0:res/2-1]*(1/T)*(1/res),10*log10(Xmean(:,1:res/2)')) % plot until pi/T only
title('Estimated PSD')
xlabel('Frequency in Hz')
ylabel('Magnitude in dB (10*log10(|PSD|))')
axis([0   150   -30    40]);
% now plot the sine wave and the measured version in noise
%figure
%plot(x,'r-')  % measured signal
%hold on
%plot(A*sin(T*[0:N-1]*w0+phi),'-') %true signal
%title('True signal blue, Measured signal red')

