function [Sx, ff] = estimatePSD(sig, win, fs)
% estimatePSD - periodogram estimate of the PSD from sampled signal
%
% Usage:
%   [Sx, ff] = estimatePSD(sig, win, fs)

%% Initialization
% figure out how many segments we can average
lenWin = length(win);
Nsegs = floor(length(sig)/lenWin);

% signal power
Psig = sum(abs(sig).^2)/(length(sig)/fs);

if Nsegs == 0
    error('Not enough samples: sig must be longer than win.')
end

% compute the PSD estimate as the average of periodogram
Sx = zeros(1,lenWin);

for nn=1:Nsegs
    sigSeg = sig((nn-1)*lenWin+1:nn*lenWin);
    % compute periodogram
    Sx = Sx + 1/Nsegs*1/lenWin*abs(fft(sigSeg.*win)).^2;
end

% make sure DC is in the center
Sx = fftshift(Sx)/(lenWin/fs);

% scale Sx such that df*sum(Sx) = Psig
Px = sum(Sx)/(lenWin/fs);
Sx = Sx*Psig/Px;


% construct a time axis from -fs/2 to fs/2
ff = (0:lenWin-1)/lenWin*fs - fs/2;