function [totPwr, initIdx, finalIdx, initPwr, finalPwr] = computeChannelPower(Pxx, f, f1, f2)
% CHANPOWER Compute a channel power measurement
%   Perform a rectangular approximation over the frequency range (F1, F2)
%   where each power density value (Pxx) is centered over the corresponding
%   frequency F:
%
%                          +-Pxx(2)-+
%                          |        |
%                 +-Pxx(1)-+        |
%                 |        |        +-Pxx(3)-+ ...
%                 |        |        |        |
%                 |<--df-->|<--df-->|<--df-->|
%                 +---()---+---()---+---()---+ ...
%                    f(1)     f(2)     f(3)      
%
%   Outputs:
%
%      totPwr  - the total integrated power between f1 and f2
%
%      initIdx  - the index into Pxx that contains F1
%      finalIdx - the index into Pxx that contains F2 
%
%      initPwr  - the power in Pxx(initIdx) not contained in totPwr
%      finalPwr - the power in Pxx(finalIdx) not contained in totPwr
%
%                     |          ... +---|-Pxx(finalIdx)----+
%                     |              |   |//////////////////|
%  +-----Pxx(initIdx)-|---+  ...     |   |/////finalPwr/////|
%  |//////////////////|   |          |   |//////////////////|
%  |////initPwr///////|   |          |   |//////////////////|
%  |<----------df-----|-->|          |<--|-------df-------->|
%  +-----------()-----|---+ ...  ... +---|------------------
%         f(initIdx)  |                  |  f(finalIdx)
%                    f1                  f2
%
%  This file is for internal use only.

%   Copyright 2013 The MathWorks, Inc.
%#codegen

df = mean(diff(f));
hdf = df/2;

% prevent errors when either f1 or f2 is NaN or out of range
if ~(f1 <= f2) || isempty(f)|| isnan(df) || f1 < f(1)-hdf || f(end)+hdf < f2
  totPwr = NaN;
  initIdx = [];
  finalIdx = [];
  initPwr = NaN;
  finalPwr = NaN;
  return
end

% extract initial power to discard
initIdx = find(f-hdf<=f1, 1, 'last');
initPwr = Pxx(initIdx) * (f1-(f(initIdx)-hdf));

% extract final power to discard
finalIdx = find(f2<=f+hdf, 1, 'first');
finalPwr = Pxx(finalIdx) * (f(finalIdx)+hdf-f2);

% compute total sum over region
allPwr = df*sum(Pxx(initIdx(1):finalIdx(1),1));

% discard the power before and after the selected region
totPwr = allPwr - initPwr - finalPwr;  
