## Copyright (C) 1995-2013 Friedrich Leisch
## Copyright (C) 2010 Alois Schloegl
## Copyright (C) 2014 Drew Abbot
##
## This file is part of Octave.
##
## Octave is free software; you can redistribute it and/or modify it
## under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 3 of the License, or (at
## your option) any later version.
##
## Octave is distributed in the hope that it will be useful, but
## WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
## General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with Octave; see the file COPYING.  If not, see
## <http://www.gnu.org/licenses/>.

## -*- texinfo -*-
## @deftypefn {Function File} {[Pxx, @var{w}] =} periodogram (@var{x})
## For a data matrix @var{x} from a sample of size @var{n}, return the
## periodogram.  The angular frequency is returned in @var{w}.
##
## [Pxx,w] = periodogram (@var{x}).
##
## [Pxx,w] = periodogram (@var{x},win).
##
## [Pxx,w] = periodogram (@var{x},win,nfft).
##
## [Pxx,f] = periodogram (@var{x},win,nfft,Fs).
##
## [Pxx,f] = periodogram (@var{x},win,nfft,Fs,"range").
##
## @itemize
## @item x: data; if real-valued a one-sided spectrum is estimated,
## if complex-valued or range indicates @qcode{"@nospell{twosided}"}, the full
## spectrum is estimated.
##
## @item win: weight data with window, x.*win is used for further computation,
## if window is empty, a rectangular window is used.
##
## @item nfft: number of frequency bins, default max (256, 2.^ceil (log2 (length (x)))).
##
## @item Fs: sampling rate, default 1.
##
## @item range: @qcode{"@nospell{onesided}"} computes spectrum from [0..nfft/2+1].
## @qcode{"@nospell{twosided}"} computes spectrum from [0..nfft-1].  These
## strings can appear at any position in the list input arguments after
## window.
##
## @item @nospell{Pxx}: one-, or two-sided power spectrum.
##
## @item w: angular frequency [0..2*pi) (two-sided) or [0..pi] one-sided.
##
## @item f: frequency [0..Fs) (two-sided) or [0..Fs/2] one-sided.
## @end itemize
## @end deftypefn

## Author: FL <address@hidden>
## Description: Compute the periodogram

function [pxx, f] = periodogram (x, varargin)

  ## check input arguments

  if (nargin < 1 || nargin > 5)
    print_usage ();
  endif

  nfft = []; fs = []; range = []; window = [];
  j = 1;
  for k = 1:length (varargin)
    if (ischar (varargin{k}))
      range = varargin{k};
    else
      switch (j)
        case 1
          window = varargin{k};
        case 2
          nfft   = varargin{k};
        case 3
          fs     = varargin{k};
        case 4
          range  = varargin{k};
      endswitch
      j++;
    endif
  endfor

  if (ischar (window))
    range = window;
    window = [];
  endif;
  if (ischar (nfft))
    range = nfft;
    nfft = [];
  endif;
  if (ischar (fs))
    range = fs;
    fs = [];
  endif;

  if (isrow (x))
    x = x.';
  endif
  
  n = size (x,1);

  if (! isempty (window))
    if (isrow (window))
      window = window.';
    endif
    if ( size(window,1) != n || ~isvector(window) )
      error('periodogram: the window must be a vector of the same length as the signal.')
    endif
    x .*= window;
  endif

  if (numel (nfft)>1)
    error ("nfft must be scalar");
  endif
  if (isempty (nfft))
    nfft = max (256, 2.^ceil (log2 (n)));
  endif

  if (strcmp (range, "onesided"))
    range = 1;
  elseif (strcmp (range, "twosided"))
    range = 2;
  else
    range = 2-isreal (x);
  endif

  ## compute periodogram

  if (n>nfft)
    Pxx = 0;
    rr = rem (length (x), nfft);
    if (rr)
      x = [x(:); (zeros (nfft-rr, 1))];
    endif
    x = sum (reshape (x, nfft, []), 2);
  endif

  if (! isempty (window))
    n = sumsq (window);
  endif;
  Pxx = (abs (fft (x, nfft))) .^ 2 / n ;

  if (nargin<4)
    Pxx /= 2*pi;
  elseif (! isempty (fs))
    Pxx /= fs;
  endif

  ## generate output arguments

  if (range == 1)  # onesided
    if ( ~ rem(nfft,2) ) # nfft is even
      psd_len = nfft/2+1;
      Pxx = Pxx(1:psd_len) + [0; Pxx(nfft:-1:psd_len+1); 0];
    else # nfft is odd
      psd_len = (nfft+1)/2;
      Pxx = Pxx(1:psd_len) + [0; Pxx(nfft:-1:psd_len+1)];
    endif
  endif

  if (nargout != 1)
    if (range == 1)
      f = (0:nfft/2)'/nfft;
    elseif (range == 2)
      f = (0:nfft-1)'/nfft;
    endif
    if (nargin<4)
      f *= 2*pi; # generate w=2*pi*f
    elseif (! isempty (fs))
      f *= fs;
    endif
  endif

  if (nargout == 0)
    if (nargin<4)
      plot (f/(2*pi), 10*log10 (Pxx));
      xlabel ("normalized frequency [x pi rad]");
      ylabel ("Power density [dB/rad/sample]");
    else
      plot (f, 10*log10 (Pxx));
      xlabel ("frequency [Hz]");
      ylabel ("Power density [dB/Hz]");
    endif
    grid on;
    title ("Periodogram Power Spectral Density Estimate");
  else
    pxx = Pxx;
  endif

endfunction