% COMPUTEPERIODOGRAM   Periodogram spectral estimation.
%   This function is used to calculate the Power Spectrum Sxx, and the
%   Cross Power Spectrum Sxy.
%
%   Sxx = COMPUTEPERIODOGRAM(X,WIN,NFFT) where x is a vector returns the
%   Power Spectrum over the whole Nyquist interval, [0, 2pi).
%
%   Sxy = COMPUTEPERIODOGRAM({X,Y},WIN,NFFT) returns the Cross Power
%   Spectrum over the whole Nyquist interval, [0, 2pi).
%
%   Inputs:
%    X           - Signal vector or a cell array of two elements containing
%                  two signal vectors.
%    WIN         - Window
%    NFFT        - Number of frequency points (FFT) or vector of
%    frequencies at whicch periodogram is desired (Goertzel)
%    WINCOMPFLAG - A string indicating the type of window compensation to
%                  be done. The choices are:
%                  'ms'  - compensate for Mean-square (Power) Spectrum;
%                          maintain the correct power peak heights.
%                  'psd' - compensate for Power Spectral Denstiy (PSD);
%                          maintain correct area under the PSD curve.
%
%   Output:
%    Sxx         - Power spectrum [Power] over the whole Nyquist interval.
%      or
%    Sxy         - Cross power spectrum [Power] over the whole Nyquist
%                  interval.
%
% A direct copy of MATLAB's function for LTPDA
%
% M Hewitson 08-05-08
%
% $Id: computeperiodogram.m,v 1.2 2008/08/01 13:19:42 ingo Exp $
%

%   Author(s): P. Pacheco
%   Copyright 1988-2006 The MathWorks, Inc.
%   $Revision: 1.2 $  $Date: 2008/08/01 13:19:42 $

function [P,f] = computeperiodogram(x,win,nfft,esttype,varargin)

  error(nargchk(3,5,nargin,'struct'));
  if nargin < 4,
    esttype = 'psd'; % Default, compenstate for window's power.
  end

  if nargin < 5 || isempty(varargin{1}),
    Fs = 2*pi;
  else
    Fs = varargin{1};
  end

  % Validate inputs and convert row vectors to column vectors.
  [x,Lx,y,is2sig,win,errid,errmsg] = validateinputs(x,win,nfft);
  if ~isempty(errmsg), error(errid, errmsg); end

  % Determine if FFT or Goertzel
  freqVectorSpecified = false;
  if (length(nfft)> 1), freqVectorSpecified = true; end

  % Window the data
  xw = x.*win;
  if is2sig, yw = y.*win; end

  % Evaluate the window normalization constant.  A 1/N factor has been
  % omitted since it will cancel below.
  if strcmpi(esttype,'ms'),
    % The window is convolved with every power spectrum peak, therefore
    % compensate for the DC value squared to obtain correct peak heights.
    U = sum(win)^2;
  else
    U = win'*win;  % compensates for the power of the window.
  end

  % Compute the periodogram power spectrum [Power] estimate
  % A 1/N factor has been omitted since it cancels

  [Xx,f] = ao.computeDFT(xw,nfft,Fs);
  if is2sig, [Yy,f] = ao.computeDFT(yw,nfft,Fs); end

  P = Xx.*conj(Xx)/U;      % Auto spectrum.
  if is2sig,
    P = Xx.*conj(Yy)/U;  % Cross spectrum.
  end

end

%--------------------------------------------------------------------------
function [x,Lx,y,is2sig,win,errid,errmsg] = validateinputs(x,win,nfft)
  % Validate the inputs to computexperiodogram and convert the input signal
  % vectors into column vectors.

  errid = '';
  errmsg = ''; % in case of early return.

  % Set defaults and convert to row vectors to columns.
  y     = [];
  Ly    = 0;
  is2sig= false;
  win   = win(:);
  Lw    = length(win);

  % Determine if one or two signal vectors was specified.
  Lx = length(x);
  if iscell(x),
    if length(x) > 1,
      y = x{2};
      y = y(:);
      is2sig = true;
    end
    x = x{1};
    Lx = length(x);
  end
  x = x(:);

  if is2sig,
    Ly  = length(y);
    if Lx ~= Ly,
      errid = generatemsgid('invalidInputSignalLength');
      errmsg = 'The length of the two input vectors must be equal to calculate the cross spectral density.';
      return;
    end
  end

  if Lx ~= Lw,
    errid = generatemsgid('invalidWindow');
    errmsg = 'The WINDOW must be a vector of the same length as the signal.';
    return;
  end

end