## Copyright (C) 2015 GNU Octave Community <vintercepter@gmail.com>
##
## This program 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.
##
## This program 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
## this program; if not, see <http://www.gnu.org/licenses/>.

## -*- texinfo -*-
## @deftypefn {Function File} {} bandpower (@var{X})
## @deftypefnx {Function File} {} bandpower (@var{X}, @var{f}, [@var{upper} @var{lower}])
##
## bandpower (@var{X}) returns the average power of the signal @var{X}.
##
## bandpower (@var{X}, @var{f}, [@var{upper} @var{lower}]) returns the average 
## power of the signal @var{X} in the band of the frequencies @var{f} given by 
## the limits [@var{upper} @var{lower}]
##
## @seealso{periodogram}
## @end deftypefn



function x = bandpower (varargin)

  if (length (varargin) == 1)
    x = avg_pwr ([varargin{1}]);
    
  elseif (length (varargin) == 2)
    error ("bandpower: not enough inputs");
    print_usage ();
    
  elseif (length (varargin) == 3)
    x = pwr_freq (varargin{:});
    
  else
     print_usage ();
     
  endif

endfunction

function pwr = avg_pwr (x)

  pwr = (x * x') / length (x);
  
endfunction

function pw = pwr_freq (x, fs, limits)

    p = periodogram (x);  
    N = length (x)
    lower_limit = limits (1);
    upper_limit = limits (2);
    lower_limits_n = floor (lower_limit * N / fs);
    upper_limits_n = floor (upper_limit * N / fs);
    pw = 2 * sum (p (max ([1, lower_limits_n]):upper_limits_n)) / N;
    
endfunction
