function iqSamples = bleCTEIQSample(cteSamples,varargin)
%bleCTEIQSample Perform IQ sampling on CTE field of Bluetooth LE packet
%
%   IQSAMPLES = bleCTEIQSample(CTESAMPLES) returns in-phase and quadrature
%   (IQ) samples, IQSAMPLES, by performing IQ sampling at each microsecond
%   of the reference period and each sample slot of the constant tone
%   extension (CTE), CTESAMPLES.
%
%   IQSAMPLES is a complex-valued column vector of type double. This value
%   contains IQ samples corresponding to 8 microseconds reference period
%   and sample slots of CTE.
%
%   CTESAMPLES is a complex-valued time-domain signal with size Ns-by-1,
%   where Ns represents the number of CTE samples received. This table
%   shows how the 'Mode' and 'SamplesPerSymbol' values impacts Ns.
%   -------------------------------------------------------
%   |Value of Mode|      Value of Ns       |  Multiple of  |
%   -------------------------------------------------------
%   |     LE1M    |   [16*sps, 160*sps]    |     8*sps     |
%   |     LE2M    |   [32*sps, 320*sps]    |     16*sps    |
%   -------------------------------------------------------
%
%   IQSAMPLES = bleCTEIQSample(...,Name=Value) specifies additional
%   name-value arguments described below. When a name-value argument is not
%   specified, the function uses the default value.
%
%   Mode                  Specify the PHY reception mode as 'LE1M' or
%                         'LE2M'. The default value is 'LE1M'.
%
%   SamplesPerSymbol      Specify samples per symbol as a positive
%                         integer. The default value is 8.
%
%   SlotDuration          Specify the switch and sample slot duration (in
%                         microseconds) as 1 or 2. The default value is 2.
%
%   SampleOffset          Specify the sample offset as a positive integer
%                         in the range [sps/8,7*sps/8] for LE1M and
%                         [sps/4,7*sps/4] for LE2M. The default value is 4
%                         (sps/2).
%
%   % Example:
%
%   % Generate Bluetooth LE waveform by enabling CTE for connection-
%   % oriented scenario with 2 microseconds slots and the CTE time of 16
%   % microseconds. Add white Gaussian noise to the waveform. Perform IQ
%   % sampling on the CTE portion of the received Bluetooth LE waveform.
%
%   % Generate a connection data channel protocol data unit (PDU) for
%   % angle of departure (AoD) CTE with 2 microseconds slots and CTE time
%   % of 16 microseconds
%   pduHex = '1B820127'; % Valid PDU in hexadecimal
%   pdu = int2bit(hex2dec(pduHex),32,false);
%   crcGen = comm.CRCGenerator('z^24+z^10+z^9+z^6+z^4+z^3+z+1',...
%   InitialConditions=int2bit(hex2dec('555551'),24)',...
%   DirectMethod=true); % Create a cyclic redundancy check (CRC) generator
%   pduCRC = crcGen(pdu); % Append CRC to the PDU
%   sps = 16; % Samples per symbol
%   txWaveform = bleWaveformGenerator(pduCRC,ChannelIndex=36,...
%                DFPacketType='ConnectionCTE',SamplesPerSymbol=sps);
%   snr = 30; % In dB
%   rxWaveform = awgn(txWaveform,snr); % Pass the waveform through the AWGN channel
%
%   % Perform IQ sampling for a slot duration of 2 microseconds and take
%   % the 8th sample after an antenna switch has happened
%   slotDuration = 2; % Slot duration
%   sampleOffset = 8; % Sample offset
%   cteSamples = rxWaveform(end-16*sps+1:end); % Waveform corresponding to CTE
%   iqSamples = bleCTEIQSample(cteSamples,SamplesPerSymbol=sps,...
%               SlotDuration=slotDuration,SampleOffset=sampleOffset);
%
%   See also bleWaveformGenerator, bleIdealReceiver, bleAngleEstimate,
%   bleAngleEstimateConfig.

%   Copyright 2021 The MathWorks, Inc.

%#codegen

% Check the number of input arguments
narginchk(1,9);

% Validation of input arguments
[phyMode,sps,slotDuration,sampleOffset] = ....
                                    validateInputArgs(cteSamples,varargin{:});
inputClass = class(cteSamples);

% Perform Gaussian matched filtering
bt = 0.5; % Bandwidth-bit period product
filtLen = 2*sps+1;
t = linspace(-1,1,filtLen);
alpha = sqrt(log(2)/2)/(bt);
h = (sqrt(pi)/alpha)*exp(-(t*pi/alpha).^2); % Compute the filter coefficients
hNorm = h./sum(h); % Normalize the filter coefficients
cteSamplesConv = conv(cteSamples,hNorm(:),'same');

% Perform IQ sampling on CTE
guardLen = 4; % Guard period
referenceLen = 8; % Reference period
multFactor = 1+strcmp(phyMode,'LE2M');
refEndSample = (guardLen+referenceLen-1)*sps*multFactor+sampleOffset;
refSamplesTimeSteps = round(guardLen*sps*multFactor+...
                      sampleOffset:sps*multFactor:refEndSample);
remSamplesTimeSteps = round(refEndSample+2*slotDuration*sps*multFactor:...
                      2*slotDuration*sps*multFactor:length(cteSamplesConv));
iqSamplesRef = cteSamplesConv(refSamplesTimeSteps); % Reference period IQ samples
iqSamplesRem = cteSamplesConv(remSamplesTimeSteps); % Sample slot IQ samples
% Append reference IQ samples and samples corresponding to sample slots
iqSamplesRefRem = [iqSamplesRef;iqSamplesRem];

% DC offset correction
iqSamplesDCFree = iqSamplesRefRem-mean(iqSamplesRef);

% Remove the frequency offset due to modulated ones
modeFreqOffset = 250e3*multFactor;
timeSteps = (0:referenceLen-1).';
freqVec = modeFreqOffset*timeSteps/1e6;
iqSamplesModeFreqFree = iqSamplesDCFree(1:referenceLen).*exp(1i*2*pi*(-freqVec));

% Estimate frequency offset
estFreqOffset = mean(diff(unwrap(angle(iqSamplesModeFreqFree))))*1e6/(2*pi);

% Compensate frequency offset
inputLength = length(iqSamplesDCFree);
timeSteps = [refSamplesTimeSteps remSamplesTimeSteps];
symbolRate = 1e6*multFactor; % Symbol rate
freqVec = (estFreqOffset*timeSteps /(sps*symbolRate)).';
iqSamples = iqSamplesDCFree.*exp(1i*2*pi*(-cast(freqVec(1:inputLength),inputClass)));

% Validate the input arguments
    function [phyMode,sps,slotDuration,sampleOffset] = ....
                                    validateInputArgs(cteSamples,varargin)

if(nargin==1)
    % For 'bleCTEIQSample(cteSamples)'
    phyMode = 'LE1M';
    spsNumeric = 8;
    slotDurationNumeric = 2;
    sampleOffsetNumeric = 4;
else
    % For 'bleCTEIQSample(cteSamples,N1=V1,...)'
    if isempty(coder.target) % Simulation path
        p = inputParser;
        p.PartialMatching = true;
        % Set default values for the optional arguments
        addParameter(p,'Mode','LE1M');
        addParameter(p,'SamplesPerSymbol',8);
        addParameter(p,'SlotDuration',2);
        addParameter(p,'SampleOffset',4);
        parse(p,varargin{:}); % Parse inputs
        res = p.Results;
        phyMode = res.Mode;
        spsNumeric = res.SamplesPerSymbol;
        slotDurationNumeric = res.SlotDuration;
        sampleOffsetNumeric = res.SampleOffset;
    else % Codegen path
        nvPairs = struct('Mode',uint32(0),'SamplesPerSymbol',uint32(0),...
            'SlotDuration',uint32(0),'SampleOffset',uint32(0));
        % Select parsing options
        popts = struct('PartialMatching', true,...
            'CaseSensitivity',false);
        % Parse inputs
        pStruct = coder.internal.parseParameterInputs(nvPairs,popts,varargin{:});

        % Get values for the N-V pair or set defaults for the optional arguments
        phyMode = coder.internal.getParameterValue(pStruct.Mode,...
            'LE1M',varargin{:});
        spsNumeric = coder.internal.getParameterValue(pStruct.SamplesPerSymbol,...
            8,varargin{:});
        slotDurationNumeric = coder.internal.getParameterValue(pStruct.SlotDuration,2,...
            varargin{:});
        sampleOffsetNumeric = coder.internal.getParameterValue(pStruct.SampleOffset,...
            4,varargin{:});
    end
phyMode = validatestring(phyMode,{'LE1M','LE2M'},mfilename,'Mode');
validateattributes(spsNumeric,{'numeric'},{'nonnan','finite','integer',...
                'positive','scalar'},mfilename,'SamplesPerSymbol');
validateattributes(slotDurationNumeric,{'numeric'},{'integer','positive',...
                    'scalar','<=',2}, mfilename,'SlotDuration');
validateattributes(sampleOffsetNumeric,{'numeric'},{'integer','positive',...
                    'scalar'}, mfilename,'SampleOffset');
end
% Validate the received CTE samples
validateattributes(cteSamples,{'double','single'},{'nonnan','finite','column'},...
                        mfilename,'CTESamples');
% Typecast sps, slotDuration and sampleOffset to double, if not a float
if ~isfloat(spsNumeric) || ~isfloat(slotDurationNumeric) || ~isfloat(sampleOffsetNumeric)
    sps = double(spsNumeric);
    slotDuration = double(slotDurationNumeric);
    sampleOffset = double(sampleOffsetNumeric);
else
    sps = spsNumeric;
    slotDuration = slotDurationNumeric;
    sampleOffset = sampleOffsetNumeric;
end

modeFactor = 1+strcmp(phyMode,'LE2M');
coder.internal.errorIf(length(cteSamples) < 16*sps*modeFactor || length(cteSamples) > 160*sps*modeFactor,...
    'bluetooth:bleCTEIQSample:InputLength',phyMode,16*modeFactor,160*modeFactor,length(cteSamples));
coder.internal.errorIf(rem(length(cteSamples),sps*8*modeFactor) ~= 0,...
    'bluetooth:bleCTEIQSample:InputLengthMultiple',phyMode,8*modeFactor,length(cteSamples));

% Validate the sample offset
coder.internal.errorIf(sampleOffset < round(sps*modeFactor/8) || sampleOffset > round(7*sps*modeFactor/8),...
    'bluetooth:bleCTEIQSample:SampleOffsetValue',phyMode,round(sps*modeFactor/8),round(7*sps*modeFactor/8),sampleOffset);
end

end