% COMPUTEDFTPERIODOGRAM compute periodogram with dft
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
% DESCRIPTION
% 
% Define one sided periodogram as:
% Xp(f) = 2*T*|DFT(w(t)*x(t))/T|^2
% Where w(t) are window sameples, which are supposed to be square
% integrable: sum(w(t)^2)=1
% 
% CALL 
% Sf = computeDftPeriodogram(x,fs,f,order,win,psll)
% Sf = computeDftPeriodogram(x,fs,f,order,win,[])
% 
% 
% INPUT
% 
% - x, data series, Nx1 double
% - fs, sampling frequency Hz, 1x1 double
% - f, frequenci vector Hz, 1x1 double
% - order, detrend order, -1,0,1,2,3,4,...
% - win, window name. e.g 'BH92'
% - psll, for Kaiser window only
% 
% REFERENCES
% 
% D. B. Percival and A. T. Walden, Spectral Analysis for Physical
% Applications (Cambridge University Press, Cambridge, 1993) p 291.
% 
% L Ferraioli 28-03-2011
%
% $Id: computeDftPeriodogram.m,v 1.2 2011/03/29 15:36:43 luigi Exp $
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Sf = computeDftPeriodogram(x,fs,f,order,win,psll)

  
  % detrend segment
  if order < 0
    xd = x; % no detrend
  else
    [xd,coeffs] = ltpda_polyreg(x,order);
  end
  
  N = numel(x);
  
  % get window
  switch lower(win)
    case 'kaiser'
      wnd = specwin('Kaiser', N, psll);
    otherwise
      wnd = specwin(win, N);
  end
  w = wnd.win;
  w = w./sqrt(wnd.ws2);  % compensates for the power of the window.
  
  if size(w)~=size(xd)
    w = w.';
  end
  
  % apply window
  xdw = xd.*w;
  
  T = 1/fs;
  
  Nf = numel(f);
  Sf = zeros(Nf,1);
  for ii=1:Nf
    xx = utils.math.dft(xdw,f(ii),T);
    xx = xx/T;
    Sf(ii) = 2*T*xx*conj(xx);
  end

  

end