%% Periodogram
% In general terms, one way of estimating the PSD of a process is to simply find the discrete-time Fourier transform of the samples of the process (usually done on a grid with an FFT) and appropriately scale the magnitude squared of the result. This estimate is called the periodogram.
% The periodogram estimate of the PSD of a signal x_L(n) of length L is
%  P_{xx}(f) =
% {1\over{LF_s}}\left|\sum_{n=0}^{L-1}x_L(n)e^{-j2\pi fn/F_s}\right|^2,
% where Fs is the sampling frequency.
% In practice, the actual computation of P_{xx}(f) can be performed only at a finite number of frequency points, and usually employs an FFT. Most implementations of the periodogram method compute the N-point PSD estimate at the frequencies
%  f_k={{kF_s}\over N},\;\;\;k=0,1,\dots,N-1.
% In some cases, the computation of the periodogram via an FFT algorithm is more efficient if the number of frequencies is a power of two. Therefore it is not uncommon to pad the input signal with zeros to extend its length to a power of two.
% As an example of the periodogram, consider the following 1001-element signal xn, which consists of two sinusoids plus noise:

fs = 1000;                % Sampling frequency
t = (0:fs)/fs;            % One second worth of samples
A = [1 2];                % Sinusoid amplitudes (row vector)
f = [150;140];            % Sinusoid frequencies (column vector)
xn = A*sin(2*pi*f*t) + 0.1*randn(size(t));
% The three last lines are equivalent to
% xn = sin(2*pi*150*t) + 2*sin(2*pi*140*t) + 0.1*randn(size(t));

%
% The periodogram estimate of the PSD can be computed using periodogram. In this case, the data vector is multiplied by a Hamming window to produce a modified periodogram.

[Pxx,F] = periodogram(xn,hamming(length(xn)),length(xn),fs);
plot(F,10*log10(Pxx))
xlabel('Hz')
ylabel('dB')
title('Modified Periodogram Power Spectral Density Estimate')

% Algorithm
% Periodogram computes and scales the output of the FFT to produce the power vs. frequency plot as follows.
% If the input signal is real-valued, the magnitude of the resulting FFT is symmetric with respect to zero frequency (DC). For an even-length FFT, only the first (1 + nfft/2) points are unique. Determine the number of unique values and keep only those unique points.
% Take the squared magnitudes of the unique FFT values. Scale the squared magnitudes (except for DC) by 2/(F_sN), where N is the length of signal prior to any zero padding. Scale the DC value by 1/(F_sN).
% Create a frequency vector from the number of unique points, the nfft and the sampling frequency.
% Plot the resulting magnitude squared FFT against the frequency.
% Copyright 2012 The MathWorks, Inc.