function [bits, eqDataSym, varargout] = trackingNonHTDataRecover( ...
    rxNonHTData, chanEst, noiseVarEst, cfgNonHT, varargin)
%trackingNonHTDataRecover Recover information bits from non-HT Data field signal with pilot tracking
%
%   BITS = trackingNonHTDataRecover(RXNONHTDATA, CHANEST, NOISEVAREST,
%   CFGNONHT) recovers the information bits in the non-HT Data field for a
%   non-HT OFDM format transmission with joint sample rate offset and
%   residual carrier frequency offset tracking.
%
%   BITS is an int8 column vector of length 8*CFGNONHT.PSDULength
%   containing the recovered information bits.
%
%   RXNONHTDATA is the received time-domain non-HT Data field signal. It is
%   a Ns-by-Nr matrix of real or complex values, where Ns represents the
%   number of time-domain samples in the non-HT Data field and Nr
%   represents the number of receive antennas. Ns can be greater than the
%   non-HT Data field length; in this case additional samples at the end of
%   RXNONHTDATA, if not required, are not used. When sample rate offset
%   tracking is enabled using the optional CFGREC argument, additional
%   samples may be required in RXNONHTDATA. This is to allow for the
%   receiver running at a higher sample rate than the transmitter and
%   therefore more samples being required.
%
%   CHANEST is the estimated channel at data and pilot subcarriers based on
%   the L-LTF. It is a real or complex array of size Nst-by-1-by-Nr, where
%   Nst represents the total number of occupied subcarriers. The singleton
%   dimension corresponds to the single transmitted stream in the L-LTF
%   which includes the combined cyclic shifts if multiple transmit antennas
%   are used.
%
%   NOISEVAREST is the noise variance estimate. It is a real, nonnegative
%   scalar.
%
%   CFGNONHT is the format configuration object of type <a href="matlab:help('wlanNonHTConfig')">wlanNonHTConfig</a> 
%   that specifies the non-HT format parameters. Only OFDM modulation is
%   supported.
% 
%   BITS = trackingNonHTDataRecover(..., CFGREC) allows different algorithm
%   options for data recovery via the input CFGREC, which is a
%   <a href="matlab:help('trackingRecoveryConfig')">trackingRecoveryConfig</a> configuration object. When the CFGREC input is not 
%   specified, the default property values of the <a href="matlab:help('trackingRecoveryConfig')">trackingRecoveryConfig</a> object
%   are adopted in the recovery. Joint sample rate offset and residual
%   carrier frequency offset tracking is enabled by default.
%
%   [..., EQDATASYM, CPE, PEG, PILOTGAIN, EQPILOTSYM] =
%   trackingNonHTDataRecover(...) also returns the equalized data
%   subcarriers, common phase error, phase error gradient, pilot gain and
%   equalized pilot subcarriers.
%
%   EQDATASYM is a complex 48-by-Nsym matrix containing the equalized
%   symbols at data carrying subcarriers. There are 48 data carrying
%   subcarriers in the non-HT Data field. Nsym represents the number of
%   OFDM symbols in the non-HT Data field.
%
%   CPE is a column vector of length Nsym containing the common phase error
%   between each received and expected OFDM symbol.
%   
%   PEG is a column vector of length Nsym containing the phase error
%   gradient per OFDM symbol in degrees per subcarrier. This error is
%   caused by a sample rate offset between transmitter and receiver.
%
%   PILOTGAIN is a Nsym-by-4 array containing the pilot gains.
%
%   EQPILOTSYM is a complex 4-by-Nsym matrix containing the equalized
%   symbols at pilot carrying subcarriers. There are 4 pilot carrying
%   subcarriers in the non-HT Data field.
%
%   Example: 
%   %  Recover a non-HT Data field signal through a SISO AWGN channel
%   %  using ZF equalization.
%  
%     cfgNonHT = wlanNonHTConfig('PSDULength', 1024);  % non-HT OFDM 
%     txBits = randi([0 1], 8*cfgNonHT.PSDULength, 1); % PSDU bits
%     tx = wlanNonHTData(txBits, cfgNonHT);      % non-HT Data field signal
% 
%     % Add AWGN, with noise variance of 1
%     rx = awgn(tx, 1, 1);
%
%     % Configure recovery object
%     cfgRec = trackingRecoveryConfig('EqualizationMethod', 'ZF'); 
%     % Recover PSDU bits 
%     rxBits = trackingNonHTDataRecover(rx, ones(52,1), 1, cfgNonHT, cfgRec);
%   
%     [numerr, ber] = biterr(rxBits, txBits); % Compare bits
%     disp(ber)
%     
%   See also trackingRecoveryConfig. 

%   Copyright 2016-2019 The MathWorks, Inc.

%#codegen

narginchk(4,5);
nargoutchk(0,6);

% Calculate CPE and/or PEG if requested
calculateCPE = false;
calculatePEG = false;
if nargout>2
    calculateCPE = true;
    if nargout>3
        calculatePEG = true;
    end
end

% Non-HT configuration input self-validation
validateattributes(cfgNonHT, {'wlanNonHTConfig'}, {'scalar'}, mfilename, 'format configuration object');
% Only applicable for OFDM and DUP-OFDM modulations
coder.internal.errorIf( ~strcmp(cfgNonHT.Modulation, 'OFDM'), 'wlan:wlanNonHTDataRecover:InvalidModulation');
s = validateConfig(cfgNonHT);

% Validate rxNonHTData
validateattributes(rxNonHTData, {'double'}, {'2d','finite'}, 'rxHTData', 'Non-HT OFDM Data field signal'); 
% Validate chanEst
validateattributes(chanEst, {'double'}, {'3d','finite'}, 'chanEst', 'channel estimates'); 
% Validate noiseVarEst
validateattributes(noiseVarEst, {'double'}, {'real','scalar','nonnegative','finite'}, 'noiseVarEst', 'noise variance estimate'); 

% Optional recovery configuration input validation
if nargin == 5
    validateattributes(varargin{1}, {'trackingRecoveryConfig'}, {'scalar'}, mfilename, 'recovery configuration object');

    symOffset = varargin{1}.OFDMSymbolOffset;
    eqMethod  = varargin{1}.EqualizationMethod; 
    pilotTracking = varargin{1}.PilotTracking;
    pilotTrackingWindow = varargin{1}.PilotTrackingWindow;
    pilotGainTracking = varargin{1}.PilotGainTracking;
else % use defaults
    symOffset = 0.75;
    eqMethod  = 'MMSE'; 
    pilotTracking = 'Joint';
    pilotTrackingWindow = 9;
    pilotGainTracking = false;
end
cfgTrack = struct('calculateCPE', calculateCPE, 'pilotTracking', pilotTracking, 'pilotTrackingWindow', pilotTrackingWindow, 'pilotGainTracking', pilotGainTracking);

numRx = size(rxNonHTData, 2);

mcsTable = wlan.internal.getRateTable(cfgNonHT);

numOFDMSym = s.NumDataSymbols;

% Get OFDM configuration
[cfgOFDM, dataInd, pilotInd] = wlan.internal.wlanGetOFDMConfig(cfgNonHT.ChannelBandwidth, 'Long', 'Legacy');

% Cross validate inputs
numST = numel([dataInd; pilotInd]); % Total number of occupied subcarriers
coder.internal.errorIf(size(chanEst, 1) ~= numST, 'wlan:wlanNonHTDataRecover:InvalidNHTChanEst1D', numST);
coder.internal.errorIf(size(chanEst, 2) ~= 1, 'wlan:wlanNonHTDataRecover:InvalidNHTChanEst2D');
coder.internal.errorIf(size(chanEst, 3) ~= numRx, 'wlan:wlanNonHTDataRecover:InvalidNHTChanEst3D');

% Extract pilot subcarriers from channel estimate
chanEstPilots = chanEst(pilotInd,:,:);

% Cross-validation between inputs
minInputLen = numOFDMSym*(cfgOFDM.FFTLength+cfgOFDM.CyclicPrefixLength);
coder.internal.errorIf(size(rxNonHTData, 1) < minInputLen, 'wlan:wlanNonHTDataRecover:ShortNHTDataInput', minInputLen);

% Get reference pilots, from IEEE Std 802.11-2012, Eqn 18-22
z = 1; % Offset by 1 to account for L-SIG pilot symbol
fnRefPilots = @()wlan.internal.nonHTPilots(numOFDMSym, z);

% OFDM demodulation and optional pilot tracking
[ofdmDemod,cpe,peg,pilotgain] = trackingOFDMDemodulate(rxNonHTData, chanEstPilots, fnRefPilots, numOFDMSym, symOffset, cfgOFDM, cfgTrack);
if calculateCPE
    varargout{1} = cpe; 
end
if calculatePEG
    varargout{2} = peg;
end
varargout{3} = pilotgain;

% Equalization
[eqSym, csi] = wlan.internal.wlanEqualize(ofdmDemod, chanEst, eqMethod, noiseVarEst);

eqDataSym = eqSym(dataInd,:,:);
csiData = csi(dataInd,:);
eqPilotSym = eqSym(pilotInd,:,:);

if nargout>5
   varargout{4} = eqPilotSym;
end

% Constellation demapping
qamDemodOut = wlanConstellationDemap(eqDataSym, noiseVarEst, mcsTable.NBPSCS);

% Apply bit-wise CSI and concatenate OFDM symbols in the first dimension
qamDemodOut = bsxfun(@times, ...
              reshape(qamDemodOut, mcsTable.NBPSCS, [], numOFDMSym), ...
              reshape(csiData, 1, [])); % [Nbpscs Nsd Nsym]
qamDemodOut = reshape(qamDemodOut, [], 1);

% Deinterleave 
deintlvrOut = wlanBCCDeinterleave(qamDemodOut, 'Non-HT', mcsTable.NCBPS);

% Channel decoding
decBits = wlanBCCDecode(deintlvrOut, mcsTable.Rate);

% Derive initial state of the scrambler
scramInit = wlan.internal.scramblerInitialState(decBits(1:7));

% Remove pad and tail bits, and descramble
if all(scramInit==0)
    % Scrambler initialization invalid (0), therefore do not descramble
    descramDataOut = decBits(1:(16+8*cfgNonHT.PSDULength));
else
    descramDataOut = wlanScramble(decBits(1:(16+8*cfgNonHT.PSDULength)), scramInit);
end

% Remove the 16 service bits
bits = descramDataOut(17:end);   

end