%% Recovery Procedure for an 802.11ax Packet
%
% This example shows how to detect a packet and decode payload bits in a
% received IEEE(R) 802.11ax(TM) waveform. The receiver recovers the packet
% format parameters from the preamble fields to decode the data field and
% the MAC frame. 
  
% Copyright 2018-2021 The MathWorks, Inc.

%% Introduction
% In an 802.11ax packet the transmission parameters are signaled
% to the receiver using the L-SIG, HE-SIG-A, and HE-SIG-B preamble fields [ <#16 1> ]:
% 
% * The L-SIG field contains information to allow the receiver to determine
% the transmission time of the packet.
% * The HE-SIG-A field contains common transmission parameters for HE-MU
% users and all transmission parameters for HE-SU and HE-EXT-SU packets.
% * The combination of length information in the L-SIG field, the
% modulation type, and the number of OFDM symbols in the HE-SIG-A field
% determines the HE packet format.
% * The HE-SIG-B field contains Resource unit (RU) allocation information
% and the transmission parameters for the users in an HE-MU packet.
%
% In this example we detect and decode an HE-MU packet within a generated
% waveform. This example can also recover HE-SU and HE-EXT-SU packets.
% All transmission parameters apart from the channel bandwidth are assumed
% to be unknown and are therefore retrieved from the decoded L-SIG,
% HE-SIG-A, and HE-SIG-B preamble fields. The recovered transmission
% parameters are used to decode the HE-Data field. Additionally, the
% following analysis is performed:
%
% * The waveform of the detected packet is recovered and displayed.
% * The spectrum of the detected packet is recovered and displayed.
% * The constellation of the equalized data symbols for all spatial streams
% is displayed.
% * The Error Vector Magnitude (EVM) of each field is measured.
% * An A-MPDU is detected and Frame Check Sequence (FCS) is determined for
% the recovered MAC frame.
% * The EVM per data symbol and spatial stream averaged over subcarriers is
% displayed.
% * The EVM per data subcarrier and spatial stream averaged over symbols is
% displayed.

%% Waveform Transmission
% In this example an 802.11ax HE-MU waveform is synthesized but you can use
% a captured waveform. You can use MATLAB(R) to acquire I/Q data from a
% wide range of instruments with the Instrument Control Toolbox(TM) and
% software defined radio platforms.
%
% The synthesized waveform is impaired by a 2x2 TGax indoor fading
% channel, additive white Gaussian noise, and carrier frequency offset. To
% generate an HE-MU waveform we configure an HE-MU format configuration 
% object <matlab:edit('wlanHEMUConfig.m') wlanHEMUConfig>. Note that the <matlab:edit('wlanHEMUConfig.m') wlanHEMUConfig> 
% configuration object is used at the transmitter side only. The receiver
% will formulate an HE recovery configuration object <matlab:edit('wlanHERecoveryConfig.m') wlanHERecoveryConfig>. 
% The unknown properties of the HE recovery configuration object are set
% after decoding the information bits in the L-SIG, HE-SIG-A, and HE-SIG-B
% fields. The helper function <matlab:edit('heSigRecGenerateWaveform.m') heSigRecGenerateWaveform> generates the impaired 
% waveform. The following processing steps are performed:
%  
% * A random payload of MSDUs is created for the MAC frame, which is
% encoded into an HE-MU packet.
% * The waveform is passed through a TGax indoor fading channel model.
% * Carrier frequency offset (CFO) and Additive white Gaussian noise (AWGN)
% are added to the waveform.

% A mixed OFDMA and MU-MIMO configuration is defined for an HE-MU packet.
% The allocation index 17 defines two 52-tone RUs, with one user in each
% RU, and one 106-tone RU. The 106-tone RU has two users in a MU-MIMO
% configuration.
cfgMU = wlanHEMUConfig(17);
cfgMU.NumTransmitAntennas = 2;

% Configure RU 1 and user 1
cfgMU.RU{1}.SpatialMapping = 'Direct';
cfgMU.User{1}.STAID = 1;
cfgMU.User{1}.APEPLength = 1e3;
cfgMU.User{1}.MCS = 5;
cfgMU.User{1}.NumSpaceTimeStreams = 2;
cfgMU.User{1}.ChannelCoding = 'LDPC';

% Configure RU 2 and user 2
cfgMU.RU{2}.SpatialMapping = 'Fourier';
cfgMU.User{2}.STAID = 2;
cfgMU.User{2}.APEPLength = 500;
cfgMU.User{2}.MCS = 4;
cfgMU.User{2}.NumSpaceTimeStreams = 1;
cfgMU.User{2}.ChannelCoding = 'BCC';

% Configure RU 3 and user 1
cfgMU.RU{3}.SpatialMapping = 'Fourier';
cfgMU.User{3}.STAID = 3;
cfgMU.User{3}.APEPLength = 100;
cfgMU.User{3}.MCS = 2;
cfgMU.User{3}.NumSpaceTimeStreams = 1;
cfgMU.User{3}.ChannelCoding = 'BCC';

% Configure RU 3 and user 2
cfgMU.User{4}.STAID = 4;
cfgMU.User{4}.APEPLength = 500;
cfgMU.User{4}.MCS = 3;
cfgMU.User{4}.NumSpaceTimeStreams = 1;
cfgMU.User{4}.ChannelCoding = 'LDPC';

% Specify propagation channel
numRx = 2; % Number of receive antennas
delayProfile = 'Model-D'; % TGax channel delay profile

% Specify impairments
noisePower = -40; % Noise power to apply in dBW
cfo = 62e3; % Carrier frequency offset Hz

% Generate waveform
rx = heSigRecGenerateWaveform(cfgMU,numRx,delayProfile,noisePower,cfo);

%% Packet Recovery
% The signal to process is stored in the variable |rx|. The processing
% steps to recover a packet are:
%
% * An HE recovery configuration object, <matlab:edit('wlanHERecoveryConfig.m') wlanHERecoveryConfig> is created. 
% The object properties are updated as preamble fields are decoded.
% * The packet is detected and synchronized.
% * The L-LTF is extracted and demodulated. The demodulated L-LTF symbols
% do not include tone rotation for each 20 MHz segment as described in
% [ <#16 2> ], section 21.3.7.5.
% * The demodulated L-LTF symbols are used for channel and noise estimates.
% * The time-domain signal containing samples equivalent to four OFDM
% symbols immediately following the L-LTF are used to determine the HE
% packet format. The packet format is updated in the <matlab:edit('wlanHERecoveryConfig.m') wlanHERecoveryConfig> object.
% * The L-LTF is demodulated. The demodulated L-LTF symbols include tone
% rotation for each 20 MHz segment as described in [ <#16 2> ], section
% 21.3.7.5. The L-LTF channel estimates (with tone rotation) are used to
% decode the pre-HE-LTF.
% * The L-SIG and RL-SIG fields are extracted. The channel is estimated on
% an extra four subcarriers per subchannel in the L-SIG and RL-SIG fields.
% The L-LTF channel estimates are updated to include the channel estimates
% on the extra subcarriers.
% * The information bits in the L-SIG field are recovered to determine the
% length of the packet in microseconds.
% * After HE-SIG-A decoding, the recovery configuration object is updated
% with common transmission parameters for an HE-MU packet and all
% transmission parameters for HE-SU and HE-EXT-SU packets.
% * For an HE-MU packet format the HE-SIG-B field is decoded. For a
% non-compressed SIGB waveform the HE-SIG-B common field is decoded first
% followed by HE-SIG-B user field. For a compressed SIGB waveform only the
% HE-SIG-B user field is decoded.
% * For an HE-MU format without SIGB compression the RU allocation and user
% transmission parameters are recovered from the HE-SIG-B field. For a
% compressed SIGB waveform the RU allocation information is inferred from
% HE-SIG-A field and the user transmission parameters are determined from
% the HE-SIG-B user field bits.
% * The <matlab:edit('wlanHERecoveryConfig.m') wlanHERecoveryConfig> object
% is created using the recovered transmission parameters for each user
% after HE-SIG-B decoding.
% * The HE-LTF field is extracted and demodulated. The demodulated symbols
% are used for channel estimation of the subcarriers allocated to the user
% of interest. The MIMO channel estimates are used to decode the HE-Data
% field.
% * The HE-Data field is extracted and the PSDU bits are recovered using
% the <matlab:edit('wlanHERecoveryConfig.m') wlanHERecoveryConfig> object for each user.
% * Detect A-MPDU within the recovered PSDU and check the FCS for the
% recovered MAC frame.

%% Setup Waveform Recovery Parameters
% In this example all transmission parameters apart from the channel
% bandwidth are assumed to be unknown and will be recovered. A recovery 
% configuration object, <matlab:edit('wlanHERecoveryConfig.m') wlanHERecoveryConfig>, is created to 
% store the recovered information in the L-SIG, HE-SIG-A, and HE-SIG-B
% preamble fields. The transmission properties in <matlab:edit('wlanHERecoveryConfig.m') wlanHERecoveryConfig>
% are updated after subsequent decoding of the preamble fields. The
% following code configures objects and variables for processing.

chanBW = cfgMU.ChannelBandwidth; % Assume channel bandwidth is known
sr = wlanSampleRate(cfgMU); % Sample rate

% Specify pilot tracking method for recovering the data field. This can be:
%  'Joint' - use joint common phase error and sample rate offset tracking
%  'CPE' - use only common phase error tracking
% When recovering 26-tone RUs only CPE tracking is used as the joint
% tracking algorithm is susceptible to noise.
pilotTracking = 'Joint';

% Create an HE recovery configuration object and set the channel bandwidth
cfgRx = wlanHERecoveryConfig;
cfgRx.ChannelBandwidth = chanBW;

% The recovery configuration object is used to get the start and end
% indices of the pre-HE-SIG-B field.
ind = wlanFieldIndices(cfgRx);

% Setup plots for the example
[spectrumAnalyzer,timeScope,ConstellationDiagram,EVMPerSubcarrier,EVMPerSymbol] = heSigRecSetupPlots(sr);

% Minimum packet length is 10 OFDM symbols
lstfLength = double(ind.LSTF(2));
minPktLen = lstfLength*5; % Number of samples in L-STF

rxWaveLen = size(rx,1); 

%% Front-End Processing
% The front-end processing consists of packet detection, coarse carrier
% frequency offset correction, timing synchronization, and fine carrier
% frequency offset correction. A <docid:matlab_ref#bub68r7-1 while> loop is
% used to detect and synchronize a packet within the received waveform. The
% sample offset |searchOffset| is used to index into |rx| to detect a
% packet. The first packet within |rx| is detected and processed. If the
% synchronization fails for the detected packet, the sample index offset
% |searchOffset| is incremented to move beyond the processed packet in
% |rx|. This is repeated until a packet has been successfully detected and
% synchronized.

searchOffset = 0; % Offset from start of waveform in samples
while (searchOffset + minPktLen) <= rxWaveLen
    % Packet detection
    pktOffset = wlanPacketDetect(rx,chanBW,searchOffset);

    % Adjust packet offset
    pktOffset = searchOffset + pktOffset;
    if isempty(pktOffset) || (pktOffset + ind.LSIG(2) > rxWaveLen)
        error('** No packet detected **');
    end

    % Coarse frequency offset estimation and correction using L-STF
    rxLSTF = rx(pktOffset+(ind.LSTF(1):ind.LSTF(2)), :);
    coarseFreqOffset = wlanCoarseCFOEstimate(rxLSTF,chanBW);
    rx = helperFrequencyOffset(rx,sr,-coarseFreqOffset);

    % Symbol timing synchronization
    searchBufferLLTF = rx(pktOffset+(ind.LSTF(1):ind.LSIG(2)),:);
    pktOffset = pktOffset+wlanSymbolTimingEstimate(searchBufferLLTF,chanBW);

    % Fine frequency offset estimation and correction using L-STF
    rxLLTF = rx(pktOffset+(ind.LLTF(1):ind.LLTF(2)),:);
    fineFreqOffset = wlanFineCFOEstimate(rxLLTF,chanBW);
    rx = helperFrequencyOffset(rx,sr,-fineFreqOffset);

    % Timing synchronization complete: packet detected
    fprintf('Packet detected at index %d\n',pktOffset + 1);

    % Display estimated carrier frequency offset
    cfoCorrection = coarseFreqOffset + fineFreqOffset; % Total CFO
    fprintf('Estimated CFO: %5.1f Hz\n\n',cfoCorrection);

    break; % Front-end processing complete, stop searching for a packet
end

% Scale the waveform based on L-STF power (AGC)
gain = 1./(sqrt(mean(rxLSTF.*conj(rxLSTF)))); 
rx = rx.*gain;

%% Packet Format Detection
% The time-domain samples equivalent to four OFDM symbols immediately
% following the L-LTF are used to determine the HE packet format [ <#16 1>
% Figure. 27-63]. The L-LTF is extracted and demodulated. For format
% detection, the demodulated L-LTF symbols must not include tone rotation
% for each 20 MHz segment as described in [ <#16 2> ], section 21.3.7.5.
% The demodulated L-LTF is used for channel and noise estimation. The L-LTF
% channel (without tone rotation) and noise power estimates are used to
% detect the packet format.

rxLLTF = rx(pktOffset+(ind.LLTF(1):ind.LLTF(2)),:);
lltfDemod = wlanLLTFDemodulate(rxLLTF,chanBW);
lltfChanEst = wlanLLTFChannelEstimate(lltfDemod,chanBW);
noiseVar = helperNoiseEstimate(lltfDemod);

disp('Detect packet format...');
rxSIGA = rx(pktOffset+(ind.LSIG(1):ind.HESIGA(2)),:);
pktFormat = wlanFormatDetect(rxSIGA,lltfChanEst,noiseVar,chanBW);
fprintf('  %s packet detected\n\n',pktFormat);

% Set the packet format in the recovery object and update the field indices
cfgRx.PacketFormat = pktFormat;
ind = wlanFieldIndices(cfgRx);

%% L-LTF Channel Estimate
% Demodulate the L-LTF and perform channel estimation. The demodulated
% L-LTF symbols include tone rotation for each 20 MHz segment as described
% in [ <#16 2> ], section 21.3.7.5. The L-LTF channel estimates (with tone
% rotation) are used to equalize and decode the pre-HE-LTF fields.

lltfDemod = wlanHEDemodulate(rxLLTF,'L-LTF',chanBW);
lltfChanEst = wlanLLTFChannelEstimate(lltfDemod,chanBW);

%% L-SIG and RL-SIG Decoding
% The L-SIG field is used to determine the receive time, or RXTIME, of the
% packet. The RXTIME is calculated using the length bits of the L-SIG
% payload. The L-SIG and RL-SIG fields are recovered to perform the channel
% estimate on the extra subcarriers in the L-SIG and RL-SIG fields. The
% |lltfChanEst| channel estimates are updated to include the channel
% estimates on extra subcarriers in the L-SIG and RL-SIG fields. The L-SIG
% payload is decoded using an estimate of the channel and noise power
% obtained from the L-LTF field. The L-SIG length property in
% <matlab:edit('wlanHERecoveryConfig.m') wlanHERecoveryConfig> is updated after L-SIG decoding.

disp('Decoding L-SIG... ');
% Extract L-SIG and RL-SIG fields
rxLSIG = rx(pktOffset+(ind.LSIG(1):ind.RLSIG(2)),:);

% OFDM demodulate
helsigDemod = wlanHEDemodulate(rxLSIG,'L-SIG',chanBW);

% Estimate CPE and phase correct symbols
helsigDemod = preHECommonPhaseErrorTracking(helsigDemod,lltfChanEst,'L-SIG',chanBW);

% Estimate channel on extra 4 subcarriers per subchannel and create full
% channel estimate
preheInfo = wlanHEOFDMInfo('L-SIG',chanBW);
preHEChanEst = preHEChannelEstimate(helsigDemod,lltfChanEst,preheInfo.NumSubchannels);

% Average L-SIG and RL-SIG before equalization
helsigDemod = mean(helsigDemod,2);

% Equalize data carrying subcarriers, merging 20 MHz subchannels
[eqLSIGSym,csi] = preHESymbolEqualize(helsigDemod(preheInfo.DataIndices,:,:), ...
    preHEChanEst(preheInfo.DataIndices,:,:),noiseVar,preheInfo.NumSubchannels);

% Decode L-SIG field
[~,failCheck,lsigInfo] = wlanLSIGBitRecover(eqLSIGSym,noiseVar,csi);

if failCheck
    disp(' ** L-SIG check fail **');
else
    disp(' L-SIG check pass');
end
% Get the length information from the recovered L-SIG bits and update the
% L-SIG length property of the recovery configuration object
lsigLength = lsigInfo.Length;
cfgRx.LSIGLength = lsigLength;

% Measure EVM of L-SIG symbols
EVM = comm.EVM;
EVM.ReferenceSignalSource = 'Estimated from reference constellation';
EVM.Normalization = 'Average constellation power';
EVM.ReferenceConstellation = wlanReferenceSymbols('BPSK');
rmsEVM = EVM(eqLSIGSym);
fprintf(' L-SIG EVM: %2.2fdB\n\n',20*log10(rmsEVM/100));

% Calculate the receive time and corresponding number of samples in the
% packet
RXTime = ceil((lsigLength + 3)/3) * 4 + 20; % In microseconds
numRxSamples = round(RXTime * 1e-6 * sr);   % Number of samples in time

fprintf(' RXTIME: %dus\n',RXTime);
fprintf(' Number of samples in the packet: %d\n\n',numRxSamples);

%%
% The waveform and spectrum of the detected packet within |rx| is
% displayed given the calculated RXTIME and corresponding number of
% samples.

sampleOffset = max((-lstfLength + pktOffset),1); % First index to plot
sampleSpan = numRxSamples + 2*lstfLength; % Number samples to plot
% Plot as much of the packet (and extra samples) as we can
plotIdx = sampleOffset:min(sampleOffset + sampleSpan,rxWaveLen);

% Configure timeScope to display the packet
timeScope.TimeSpan = sampleSpan/sr;
timeScope.TimeDisplayOffset = sampleOffset/sr;
timeScope.YLimits = [0 max(abs(rx(:)))];
timeScope(abs(rx(plotIdx,:)));
release(timeScope);

% Display the spectrum of the detected packet
spectrumAnalyzer(rx(pktOffset + (1:numRxSamples),:));
release(spectrumAnalyzer);

%% HE-SIG-A Decoding
% The HE-SIG-A field contains the transmission configuration of an HE
% packet. An estimate of the channel and noise power obtained from the
% L-LTF is required to decode the HE-SIG-A field.

disp('Decoding HE-SIG-A...')
rxSIGA = rx(pktOffset+(ind.HESIGA(1):ind.HESIGA(2)),:);
sigaDemod = wlanHEDemodulate(rxSIGA,'HE-SIG-A',chanBW);
hesigaDemod = preHECommonPhaseErrorTracking(sigaDemod,preHEChanEst,'HE-SIG-A',chanBW);

% Equalize data carrying subcarriers, merging 20 MHz subchannels
preheInfo = wlanHEOFDMInfo('HE-SIG-A',chanBW);
[eqSIGASym,csi] = preHESymbolEqualize(hesigaDemod(preheInfo.DataIndices,:,:), ...
                                      preHEChanEst(preheInfo.DataIndices,:,:), ...
                                      noiseVar,preheInfo.NumSubchannels);
% Recover HE-SIG-A bits
[sigaBits,failCRC] = wlanHESIGABitRecover(eqSIGASym,noiseVar,csi);

% Perform the CRC on HE-SIG-A bits
if failCRC
    disp(' ** HE-SIG-A CRC fail **');
else
    disp(' HE-SIG-A CRC pass');
end

% Measure EVM of HE-SIG-A symbols
release(EVM);
if strcmp(pktFormat,'HE-EXT-SU')
    % The second symbol of an HE-SIG-A field for an HE-EXT-SU packet is
    % QBPSK.
    EVM.ReferenceConstellation = wlanReferenceSymbols('BPSK',[0 pi/2 0 0]);
    % Account for scaling of L-LTF for an HE-EXT-SU packet
    rmsEVM = EVM(eqSIGASym*sqrt(2));
else
    EVM.ReferenceConstellation = wlanReferenceSymbols('BPSK');
    rmsEVM = EVM(eqSIGASym);
end
fprintf(' HE-SIG-A EVM: %2.2fdB\n\n',20*log10(mean(rmsEVM)/100));

%% Interpret Recovered HE-SIG-A bits
% The <matlab:edit('wlanHERecoveryConfig.m') wlanHERecoveryConfig> object
% is updated after interpreting the recovered HE-SIG-A bits. 

cfgRx = interpretHESIGABits(cfgRx,sigaBits);          
ind = wlanFieldIndices(cfgRx); % Update field indices
%%
% Display the common transmission configuration obtained from HE-SIG-A
% field for an HE-MU packet. The properties indicated by -1 are unknown or
% undefined. The unknown user-related properties are updated after
% successful decoding of the HE-SIG-B field.
disp(cfgRx)

%% HE-SIG-B Decoding
% For an HE-MU packet the HE-SIG-B field contains:
%
% * The RU allocation information for a non-compressed SIGB waveform is
% inferred from HE-SIG-B Common field [ <#16 1> Table. 27-23]. For a
% compressed SIGB waveform the RU allocation information is inferred from
% the recovered HE-SIG-A bits.
% * For a non-compressed SIGB waveform the number of HE-SIG-B symbols are
% updated in the <matlab:edit('wlanHERecoveryConfig.m') wlanHERecoveryConfig> object. 
% The symbols are only updated if the number of HE-SIG-B symbols indicated
% in the HE-SIG-A field is set to 15 and all content channels pass the CRC.
% The number of HE-SIG-B symbols indicated in the HE-SIG-A field are not
% updated if any HE-SIG-B content channel fails the CRC.
% * The user transmission parameters for both SIGB compressed and
% non-compressed waveforms are inferred from the HE-SIG-B user field [ <#16 1>
% Table. 27-25, 27-26].
%
% An estimate of the channel and noise power obtained from the L-LTF is
% required to decode the HE-SIG-B field.

if strcmp(pktFormat,'HE-MU')
    fprintf('Decoding HE-SIG-B...\n');
    if ~cfgRx.SIGBCompression
        fprintf(' Decoding HE-SIG-B common field...\n');
        s = getSIGBLength(cfgRx);
        % Get common field symbols. The start of HE-SIG-B field is known
        rxSym = rx(pktOffset+(double(ind.HESIGA(2))+(1:s.NumSIGBCommonFieldSamples)),:);
        % Decode HE-SIG-B common field
        [status,cfgRx] = heSIGBCommonFieldDecode(rxSym,preHEChanEst,noiseVar,cfgRx);

        % CRC on HE-SIG-B content channels
        if strcmp(status,'Success')
            fprintf('  HE-SIG-B (common field) CRC pass\n');
        elseif strcmp(status,'ContentChannel1CRCFail')
            fprintf('  ** HE-SIG-B CRC fail for content channel-1\n **');
        elseif strcmp(status,'ContentChannel2CRCFail')
            fprintf('  ** HE-SIG-B CRC fail for content channel-2\n **');
        elseif any(strcmp(status,{'UnknownNumUsersContentChannel1','UnknownNumUsersContentChannel2'}))
            error('  ** Unknown packet length, discard packet\n **');
        else
            % Discard the packet if all HE-SIG-B content channels fail
            error('  ** HE-SIG-B CRC fail **');
        end
        % Update field indices as the number of HE-SIG-B symbols are
        % updated
        ind = wlanFieldIndices(cfgRx);
    end

    % Get complete HE-SIG-B field samples
    rxSIGB = rx(pktOffset+(ind.HESIGB(1):ind.HESIGB(2)),:);
    fprintf(' Decoding HE-SIG-B user field... \n');
    % Decode HE-SIG-B user field
    [failCRC,cfgUsers] = heSIGBUserFieldDecode(rxSIGB,preHEChanEst,noiseVar,cfgRx);

    % CRC on HE-SIG-B users
    if ~all(failCRC)
        fprintf('  HE-SIG-B (user field) CRC pass\n\n');
        numUsers = numel(cfgUsers);
    elseif all(failCRC)
        % Discard the packet if all users fail the CRC
        error('  ** HE-SIG-B CRC fail for all users **');
    else
        fprintf('  ** HE-SIG-B CRC fail for at least one user\n **');
        % Only process users with valid CRC
        numUsers = numel(cfgUsers);
    end

else % HE-SU, HE-EXT-SU
    cfgUsers = {cfgRx};
    numUsers = 1;
end

%% HE-Data Decoding
% The updated <matlab:edit('wlanHERecoveryConfig.m') wlanHERecoveryConfig>
% object for each user can then be used to recover the PSDU bits for each
% user in the HE-Data field.

cfgDataRec = trackingRecoveryConfig;
cfgDataRec.PilotTracking = pilotTracking;

fprintf('Decoding HE-Data...\n');
for iu = 1:numUsers
    % Get recovery configuration object for each user
    user = cfgUsers{iu};
    if strcmp(pktFormat,'HE-MU')
        fprintf(' Decoding User:%d, STAID:%d, RUSize:%d\n',iu,user.STAID,user.RUSize);
    else
        fprintf(' Decoding RUSize:%d\n',user.RUSize);
    end

    heInfo = wlanHEOFDMInfo('HE-Data',chanBW,user.GuardInterval,[user.RUSize user.RUIndex]);

    % HE-LTF demodulation and channel estimation  
    rxHELTF = rx(pktOffset+(ind.HELTF(1):ind.HELTF(2)),:);
    heltfDemod = wlanHEDemodulate(rxHELTF,'HE-LTF',chanBW,user.GuardInterval, ...
        user.HELTFType,[user.RUSize user.RUIndex]);
    [chanEst,pilotEst] = heLTFChannelEstimate(heltfDemod,user);

    % Number of expected data OFDM symbols
    symLen = heInfo.FFTLength+heInfo.CPLength;
    numOFDMSym = double((ind.HEData(2)-ind.HEData(1)+1))/symLen;

    % HE-Data demodulation with pilot phase and timing tracking
    % Account for extra samples when extracting data field from the packet
    % for sample rate offset tracking. Extra samples may be required if the
    % receiver clock is significantly faster than the transmitter.
    maxSRO = 120; % Parts per million
    Ne = ceil(numRxSamples*maxSRO*1e-6); % Number of extra samples
    Ne = min(Ne,rxWaveLen-numRxSamples); % Limited to length of waveform
    numRxSamplesProcess = numRxSamples+Ne;
    rxData = rx(pktOffset+(ind.HEData(1):numRxSamplesProcess),:);
    if user.RUSize==26
        % Force CPE only tracking for 26-tone RU as algorithm susceptible
        % to noise
        cfgDataRec.PilotTracking = 'CPE';
    else
        cfgDataRec.PilotTracking = pilotTracking;
    end
        [demodSym,cpe,peg] = heTrackingOFDMDemodulate(rxData,chanEst,numOFDMSym,user,cfgDataRec);

    % Estimate noise power in HE fields
    demodPilotSym = demodSym(heInfo.PilotIndices,:,:);
    nVarEst = heNoiseEstimate(demodPilotSym,pilotEst,user);

    % Equalize
    [eqSym,csi] = heEqualizeCombine(demodSym,chanEst,nVarEst,user);

    % Discard pilot subcarriers
    eqSymUser = eqSym(heInfo.DataIndices,:,:);
    csiData = csi(heInfo.DataIndices,:);

    % Demap and decode bits
    rxPSDU = wlanHEDataBitRecover(eqSymUser,nVarEst,csiData,user,'LDPCDecodingMethod','layered-bp');

    % Deaggregate the A-MPDU
    [mpduList,~,status] = wlanAMPDUDeaggregate(rxPSDU,wlanHESUConfig);
    if strcmp(status,'Success')
        fprintf('  A-MPDU deaggregation successful \n');
    else
        fprintf('  A-MPDU deaggregation unsuccessful \n');
    end

    % Decode the list of MPDUs and check the FCS for each MPDU
    for i = 1:numel(mpduList)
        [~,~,status] = wlanMPDUDecode(mpduList{i},wlanHESUConfig,'DataFormat','octets');
        if strcmp(status,'Success')
            fprintf('  FCS pass for MPDU:%d\n',i);
        else
            fprintf('  FCS fail for MPDU:%d\n',i);
        end
    end
 
    % Plot equalized constellation of the recovered HE data symbols for all
    % spatial streams per user
    hePlotEQConstellation(eqSymUser,user,ConstellationDiagram,iu,numUsers);

    % Measure EVM of HE-Data symbols
    release(EVM);
    EVM.ReferenceConstellation = wlanReferenceSymbols(user);
    rmsEVM = EVM(eqSymUser(:));
    fprintf('  HE-Data EVM:%2.2fdB\n\n',20*log10(rmsEVM/100));
    
    % Plot EVM per symbol of the recovered HE data symbols
    hePlotEVMPerSymbol(eqSymUser,user,EVMPerSymbol,iu,numUsers);
    
    % Plot EVM per subcarrier of the recovered HE data symbols
    hePlotEVMPerSubcarrier(eqSymUser,user,EVMPerSubcarrier,iu,numUsers);
end


%% Appendix
% This example uses the following helper functions:
%
% * <matlab:edit('heCommonPhaseErrorTracking.m') heCommonPhaseErrorTracking.m>
% * <matlab:edit('heCPECorrection.m') heCPECorrection.m>
% * <matlab:edit('heEqualizeCombine.m') heEqualizeCombine.m> 
% * <matlab:edit('helperFrequencyOffset.m') helperFrequencyOffset.m>
% * <matlab:edit('helperNoiseEstimate.m') helperNoiseEstimate.m>
% * <matlab:edit('helperPlotPositions.m') helperPlotPositions.m>
% * <matlab:edit('helperSymbolEqualize.m') helperSymbolEqualize.m>
% * <matlab:edit('heLTFChannelEstimate.m') heLTFChannelEstimate.m>
% * <matlab:edit('heNoiseEstimate.m') heNoiseEstimate.m>
% * <matlab:edit('hePlotEQConstellation.m') hePlotEQConstellation.m>
% * <matlab:edit('hePlotEVMPerSubcarrier.m') hePlotEVMPerSubcarrier.m>
% * <matlab:edit('hePlotEVMPerSymbol.m') hePlotEVMPerSymbol.m>
% * <matlab:edit('heSigRecSetupPlots.m') heSigRecSetupPlots.m>
% * <matlab:edit('heSigRecGenerateWaveform.m') heSigRecGenerateWaveform.m>
% * <matlab:edit('heSIGBCommonFieldDecode.m') heSIGBCommonFieldDecode.m>
% * <matlab:edit('heSIGBMergeSubchannels.m') heSIGBMergeSubchannels.m>
% * <matlab:edit('heSIGBUserFieldDecode.m') heSIGBUserFieldDecode.m>
% * <matlab:edit('heTrackingOFDMDemodulate.m') heTrackingOFDMDemodulate.m>
% * <matlab:edit('preHEChannelEstimate.m') preHEChannelEstimate.m>
% * <matlab:edit('preHESymbolEqualize.m') preHESymbolEqualize.m>
% * <matlab:edit('trackingRecoveryConfig.m') trackingRecoveryConfig.m>

%% Selected Bibliography
% # IEEE P802.11ax(TM)/D4.1 Draft Standard for Information technology
% - Telecommunications and information exchange between systems - Local and
% metropolitan area networks - Specific requirements - Part 11: Wireless
% LAN Medium Access Control (MAC) and Physical Layer (PHY) Specifications -
% Amendment 6: Enhancements for High Efficiency WLAN.
% # IEEE Std 802.11(TM)-2016 IEEE Standard for Information technology -
% Telecommunications and information exchange between systems - Local and
% metropolitan area networks - Specific requirements - Part 11: Wireless
% LAN Medium Access Control (MAC) and Physical Layer (PHY) Specifications.