function [totPwr, startFreq, stopFreq] = OBW(Pxx, f, f1, f2, pctObw)
% OBW Compute power and border of the occupied bandwidth over a specified band
%#codegen

% obtain the total power between the channel, the indices which 
% contain the frequency boundaries and the excess (unwanted) power
% below the initial frequency
[totPwr, initIdx, finalIdx, initPwr] = computeChannelPower(Pxx, f, f1, f2);

% prevent errors when either f1 or f2 is NaN or out of range
if isnan(totPwr(1))
  startFreq = NaN;
  stopFreq = NaN;
  return
end

% obtain average bin width
df = mean(diff(f));
hdf = df/2;

% ensure we always have at least two points to interpolate
cumPwr = [0; df*cumsum(Pxx(initIdx(1):finalIdx(1),1))]-initPwr(1);  
fOffset = [f(initIdx)-hdf; f(initIdx(1):finalIdx(1))+hdf];

% distribute the unoccupied bandwidth evenly between the borders
startThresh = (100-pctObw)/200 * totPwr;
stopThresh = (100+pctObw)/200 * totPwr;

% now compute the borders of the occupied bandwidth

startFreq = interpFreq(cumPwr, fOffset, startThresh);
stopFreq  = interpFreq(cumPwr, fOffset, stopThresh);


function f = interpFreq(cumPwr, fOffset, pwrThresh)
idx = find(cumPwr > pwrThresh(1), 1, 'first');
if ~isempty(idx) && idx(1)>1
  i = idx(1);
  f = signal.internal.linterp(fOffset(i-1), fOffset(i), ...
                              cumPwr(i-1), cumPwr(i), pwrThresh(1));
else
  f = NaN;
end
