function [power, acpr] = ACPR(Pxx, f, chanCtr, offsets, nOffsets, adjWidth, filtShape, filtCoeff, chanWidth)
% Spectrum Analyzer callback routine for computing ACPR
%    Pxx       - input power spectrum
%    f         - frequency points corresponding to the input power spectrum
%    chanCtr   - channel center frequency (in Hz)
%    offsets   - adjacent channel offsets
%    nOffsets  - number of desired offsets
%    adjWidth  - width of adjacent channel
%    filtShape - 1 = None, 2 = Gaussian, 3 = RRC
%    filtCoeff - governs alpha parameter for Gaussian / RRC filter
%    chanWidth - width of main channel
%
%  This file is for internal use only.  It is subject to change.

%   Copyright 2011-2013 The MathWorks, Inc.
%#codegen


% initialize lower and upper sideband results
acpr = zeros(nOffsets, 2);

% compute the main channel power
mcp = fchanpower(Pxx, f, chanCtr, chanWidth, filtShape, filtCoeff, chanWidth);

% for each offset, compute the power ratio relative to the main channel
for i=1:nOffsets
  offset = offsets(i);
  % compute lower sideband
  acp = fchanpower(Pxx, f, chanCtr-offset, adjWidth, filtShape, filtCoeff, chanWidth);
  acpr(i, 1) = 10*log10(acp(1) / mcp(1));
  % compute upper sideband
  acp = fchanpower(Pxx, f, chanCtr+offset, adjWidth, filtShape, filtCoeff, chanWidth);
  acpr(i, 2) = 10*log10(acp(1) / mcp(1));
end

% always convert main channel power to dBm or dBW
power = 10*log10(mcp(1));


function totPwr = fchanpower(Pxx, f, fc, fw, filtFlag, coeff, symRate)
% FCHANPOWER Compute a filtered channel power measurement
% centered at fc, over width fw with bt/rolloff factor 
%   Pxx      - spectral power density
%   f        - corresponding frequencies 
%   fc       - center of (adjacent) channel
%   fw       - width of (adjacent) channel
%   filtFlag - 1 = None, 2 = Gaussian, 3 = RRC
%   coeff    - (BT product for gaussian, alpha for RRC)
%   symRate  - symbol rate (set to full main channel width)

if filtFlag == 2  %Gaussian
  %coeff is BT product
  Hxx = gaussFreqResp(f, fc, symRate, coeff);
  Pxx = Pxx .* Hxx.^2;
elseif filtFlag == 3 %RRC
  %coeff is roll off ratio
  Hxx = rrcFreqResp(f, fc, symRate, coeff);
  Pxx = Pxx .* Hxx.^2;
end
totPwr = computeChannelPower(Pxx, f, fc-fw/2, fc+fw/2);

function Hxx = gaussFreqResp(f, fc, symRate, BT)
%GAUSSFREQRESP construct gaussian frequency response (magnitude)
%   Reference:
%     Rappaport T.S., Wireless Communications Principles and Practice,  
%     2nd Ed., Prentice Hall, 2002.  

% see reference, eqns (5.53, 5.54)
alpha = sqrt(log(2)/2) / BT;
Hxx = exp(-(alpha*(f-fc)/symRate).^2);

function Hxx = rrcFreqResp(f, fc, symRate, alpha)
%RRCFREQRESP construct root raised cosine frequency response (magnitude)
%   Reference:
%     Feher K., Advanced Digital Communications Systems and Signal 
%     Processing Techniques, Prentice Hall, 1987.

% see reference, eqn (7.4):
%  replace omega with 2*pi*(f-fc)
%  replace 1/Ts with symRate
inner = abs(f-fc) <= (symRate/2)*(1-alpha);
outer = abs(f-fc) >= (symRate/2)*(1+alpha);
Hxx = ones(size(f));
%Hxx(inner) = 1;    
Hxx(outer) = 0;
Hxx(~inner & ~outer) = cos(pi/(4*alpha) * (2*abs(f(~inner & ~outer)-fc)/symRate-(1-alpha)));
