A = .5;
L = 2^8
k=12;
w0 = 2*pi*(k/L);

ndo = 100;
peak_was = zeros(1,ndo);
fprintf('x(t) = A*cos(w0*n+phi) + w\n')
fprintf('   A = .5\n')
fprintf('  w0 = 2*pi*(12/L) = 2*pi*0.0469\n')
fprintf(' phi = 2*pi*rand(1,1)\n')
fprintf('----------------------------------\n')
fprintf('In figure(1), we see many periodograms - one after another.\n')
for ntrial = 1:ndo
    % Below, compare sample mean with phi = 0 (deterministic), too...
    phi = 2*pi*rand(1,1);  % uniform random variable between 0 and 2*pi
    n = 1:L;
    s = 1; % sqrt of the VARIANCE, s = (var)^.5
    en = normrnd(0,s,1,L);  % white Gaussian noise...
    xn = A*cos(w0*n+phi) + en;
    
    % Calculate power density of the PROCESS (many, many, many trials):
    w = (2*pi/L)*[-(L/2):(L/2)-1];
    % Below, note that (w==w0) yields a DELTA function at w0...
    % ... also, s^2 is added on at ALL frequencies:
    Pxx = A^2*L/4*((w==w0)+(w==-w0)) + s^2;
    
    figure(1); clf;
    subplot(211)
    Fxn = fft(xn);  % plot such that w=0 is CENTERED on axes
    Fxn_shift = fftshift(Fxn);  % plot such that w=0 is CENTERED on axes
    plot(w,Pxx,'r-','LineWidth',1); hold on
    plot(w,(1/L)*(abs(Fxn_shift)).^2,'k.-'); hold on
    
    title('Power Spectral Density - For one trial')
    peak_was(ntrial) = (1/L)*abs(Fxn(k+1))^2;
    axis([-4 4 0 .125*L])
    plot([-10 10],A^2*L/4*[1 1],'g--')
    plot([-10 10],A^2*L/4*[1 1]+s^2,'b--')
    legend('Expected spectral density','AVERAGED periodogram','Power at carrier frequency','Peak spectral power expected')
    xlabel('Frequency, w')
    ylabel('(Power) Spectral Density, S(w)')
    
    subplot(212)
    plot(n,xn,'b.'); hold on
    title('Time samples')
    
    
    
    if ntrial==1
        FXNtotal = abs(Fxn_shift).^2;
        Xtotal = xn;
    else
        FXNtotal = FXNtotal + abs(Fxn_shift).^2;
        Xtotal = Xtotal + xn;
    end
    if 1
        pause(.05)
    else
        input('enter...')
    end
end

input('Hit enter to see ENSEMBLE average, for 100 trials of the process...')

figure(2); clf
subplot(211)
%Fxn = fft(xn);  % plot such that w=0 is CENTERED on axes
%Fxn_shift = fftshift(Fxn);  % plot such that w=0 is CENTERED on axes

plot(w,Pxx,'r-','LineWidth',2); hold on
plot(w,(1/L)*(1/ntrial)*FXNtotal,'k.-'); hold on
title('Power Spectral Density - Averaged for all trials')
axis([-4 4 0 .125*L])
plot([-10 10],A^2*L/4*[1 1],'g--')
plot([-10 10],A^2*L/4*[1 1]+s^2,'b--')
legend('Expected spectral density','AVERAGED periodogram','Power at carrier frequency','Peak spectral power expected')
xlabel('Frequency, w')
ylabel('(Power) Spectral Density, S(w)')

subplot(212)
plot(n,(1/ntrial)*Xtotal,'b.'); hold on
title('Averaged sample values - for all trials')



