%% 802.11ax Packet Error Rate Simulation for Uplink Trigger-Based Format
%
% This example shows how to measure the packet error rate of an IEEE(R)
% 802.11ax(TM) high efficiency (HE) uplink, trigger-based (TB) format.

% Copyright 2018-2019 The MathWorks, Inc.

%% Introduction
% The 802.11ax [ <#14 1> ] HE trigger-based (HE TB) format allows for OFDMA
% or MU-MIMO transmission in the uplink. An HE TB transmission is
% controlled entirely by an access point (AP). All the parameters required
% for the transmission are provided in a trigger frame to all STAs
% participating in the HE TB transmission. Each station (STA) transmits an
% HE TB packet simultaneously, when triggered by the AP as shown in the
% following diagram.
%
% <<../heMUTBScenario.png>>
%
% In this example an end-to-end simulation is used to determine the packet
% error rate of an HE TB link for four STAs in a MU-MIMO configuration. At
% each SNR point multiple packets are transmitted with no impairments apart
% from channel and noise. The received packets are demodulated and the
% PSDUs recovered for each STA. The PSDUs are compared to those transmitted
% to determine the number of packet errors and hence the packet error rate
% for all users. Packet detection, timing synchronization and symbol
% equalization is performed by the receiver. No frequency offset correction
% is performed in this example. The processing of HE TB processing chain is
% shown in the following diagram.
%
% <<../heMUTBProcessing.png>>
%
% The receiver performs a minimum-mean-square-error-based ordered
% successive interference cancellation (MMSE-SIC) process for data
% equalization [ <#14 2> ]. To avoid error propagation in the cancellation
% stage, data streams for all STAs are sorted in descending order based on
% the channel state information and equalized sequentially. This diagram
% shows the procedure of MMSE-SIC equalization.
%
% <<../MMSESICEqualization.png>>

%% Equalization Method
% In this example, you can specify the equalization method as 'mmse' or
% 'mmse-sic'. The default equalizer is 'mmse-sic'.

equalizationMethod = 'mmse-sic';

%% User Configuration
% In this example the allocation information and transmit parameters for
% multiple uplink STAs are configured using an |heTBSystemConfig|
% object.

allocationIndex = 195; % Four uplink users in a MU-MIMO configuration
cfgSys = heTBSystemConfig(allocationIndex);

%%
% In a trigger-based transmission some parameters are the same for all
% uplink users, while some can differ. The |User| property of |cfgSys|
% contains a cell array of user configurations. Each element of the cell
% array is an object which can be configured to set the parameters of
% individual users. In this example, all users have the same transmission
% parameters.

% These parameters are same for all users in the MU-MIMO system
cfgSys.HELTFType = 4;          % HE-LTF compression mode
cfgSys.GuardInterval = 3.2;    % Guard interval type
cfgSys.SingleStreamPilots = 1; % Single stream pilot transmission of HE-LTF
numRx = 8;                     % Number of receive(AP) antennas

% The individual parameters for each user are specified below
allocInfo = ruInfo(cfgSys);
numUsers = allocInfo.NumUsers; % Number of uplink users

for userIdx = 1:numUsers
    cfgSys.User{userIdx}.NumTransmitAntennas = 1;
    cfgSys.User{userIdx}.NumSpaceTimeStreams = 1;
    cfgSys.User{userIdx}.SpatialMapping = 'Direct';
    cfgSys.User{userIdx}.MCS = 7;
    cfgSys.User{userIdx}.APEPLength = 1e3;
    cfgSys.User{userIdx}.ChannelCoding = 'LDPC';
end

%%
% A trigger-based transmission for a single user within the system is
% configured with an |wlanHETBConfig| object. The transmission configurations
% for all users are generated using the method |getUserConfig|. A cell
% array of four HE TB objects is created to describe the transmission of
% four users.
cfgTB = getUserConfig(cfgSys);

%% Simulation Parameters
% For each SNR point (dB) in the |snr| vector a number of packets are
% generated, passed through a channel and demodulated to determine the
% packet error rate.

snr = 20:2:24;

% The sample rate and field indices for the HE TB packet is same for all
% users. Here the trigger configuration of the first user is used to get
% the sample rate and field indices of the HE TB PPDU.
fs = wlanSampleRate(cfgTB{1});    % Same for all users
ind = wlanFieldIndices(cfgTB{1}); % Same for all users

%% Channel Configuration
% In this example, a TGax NLOS indoor channel model is used with delay
% profile Model-B. Model-B is considered NLOS when the distance between the
% transmitter and receiver is greater than or equal to 5 meters. This is
% described further in <matlab:doc('wlanTGaxChannel') wlanTGaxChannel>. In
% this example all the STAs are assumed to be at the same distance from the
% AP.

tgaxBase = wlanTGaxChannel;
tgaxBase.SampleRate = fs;
tgaxBase.TransmissionDirection = 'Uplink';
tgaxBase.TransmitReceiveDistance = 10;
chanBW = cfgSys.ChannelBandwidth;
tgaxBase.ChannelBandwidth = chanBW;
tgaxBase.NumReceiveAntennas = numRx;

%%
% An individual channel is created for each of the four users. Each channel
% is a clone of |tgaxBase|, but with a different |UserIndex| property, and
% is stored in a cell array |tgax|. The |UserIndex| property of each
% individual channel is set to provide a unique channel for each user. In
% this example a random channel realization is used for each packet by
% randomly varying the |UserIndex| property for each transmitted packet.

% A cell array stores the channel objects, one per user
tgax = cell(1,numUsers);
for userIdx = 1:numUsers
    tgax{userIdx} = clone(tgaxBase);
    tgax{userIdx}.NumTransmitAntennas = cfgSys.User{userIdx}.NumTransmitAntennas;
    tgax{userIdx}.UserIndex = userIdx;
end

%% Processing SNR Points
% For each SNR point a number of packets are tested and the packet error
% rate is calculated. The pre-HE preamble of 802.11ax is backwards
% compatible with 802.11ac(TM), therefore in this example the timing
% synchronization components for a VHT waveform are used to synchronize the
% HE waveform at the receiver. For each user, the following processing
% steps occur to create a waveform at the receiver containing all four
% users:
%
% # To create an HE TB waveform, a PSDU is created and encoded for each
% user based on predefined user parameters.
% # The waveform for each user is passed through an indoor TGax channel
% model. Different channel realizations are modeled for different users
% and packets, by randomly varying the |UserIndex| property of the channel.
% This results in same spatial correlation properties for all users.
% # The waveforms for all HE TB users are scaled and combined to ensure
% same SNR for each user after the addition of noise.
% # AWGN is added to the received waveform to create the desired average
% SNR per subcarrier after OFDM demodulation for each user.
% <matlab:doc('comm.AWGNChannel') comm.AWGNChannel> is configured to
% provide the correct SNR. The configuration accounts for normalization
% within the channel by the number of receive antennas, and the noise
% energy in unused subcarriers which is removed during OFDM demodulation.

%%
% At the receiver (AP) the following processing steps occur:
%
% # The packet is detected.
% # Fine timing synchronization is established. The L-STF, L-LTF and L-SIG
% samples are provided for fine timing to allow for packet detection at the
% start or end of the L-STF.
% # The HE-LTF and HE-Data fields for all users are extracted from the
% synchronized received waveform. The HE-LTF and HE-Data fields are OFDM
% demodulated.
% # The demodulated HE-LTF is extracted for each RU and channel estimation
% is performed.
% # Noise estimation is performed using the demodulated data field pilots
% for each RU.
% # The data field is extracted and equalized for all users within an RU,
% from the demodulated data field.
% # For each RU, and user within the RU, the spatial streams for a user are
% demodulated and decoded to recover the transmitted PSDU.
%
% A parfor loop can be used to parallelize processing of the SNR points,
% therefore for each SNR point an AWGN channel is created and configured
% with comm.AWGNChannel. To enable the use of parallel computing for
% increased speed comment out the 'for' statement and uncomment the
% 'parfor' statement below.

ofdmInfo = wlanHEOFDMInfo('HE-Data',cfgSys.ChannelBandwidth,cfgSys.GuardInterval);
numSNR = numel(snr); % Number of SNR points
numPackets = 50; % Number of packets to simulate
packetErrorRate = zeros(numUsers,numSNR);
txPSDU = cell(numUsers);

% parfor isnr = 1:numSNR % Use 'parfor' to speed up the simulation
for isnr = 1:numSNR
    % Set random substream index per iteration to ensure that each
    % iteration uses a repeatable set of random numbers
    stream = RandStream('combRecursive','Seed',0);
    stream.Substream = isnr;
    RandStream.setGlobalStream(stream);

    % Create an instance of the AWGN channel per SNR point simulated
    awgn = comm.AWGNChannel;
    awgn.NoiseMethod = 'Signal to noise ratio (SNR)';
    awgn.SignalPower = 1/numRx;
    sysInfo = ruInfo(cfgSys);

    % Simulate multiple packets
    numPacketErrors = zeros(numUsers,1);
    for pktIdx = 1:numPackets

        % Transmit processing
        rxWaveform = 0;
        packetError = zeros(numUsers,1);
        txPSDU = cell(1,numUsers);

        % Generate random channel realization for each packet by varying
        % the UserIndex property of the channel. This assumes all users
        % have the same number of transmit antennas.
        chPermutations = randperm(numUsers);
        for userIdx = 1:numUsers
            % HE TB config object for each user
            cfgUser = cfgTB{userIdx};

            % Generate a packet with random PSDU
            txPSDU{userIdx} = randi([0 1],getPSDULength(cfgUser)*8,1,'int8');

            % Generate HE TB waveform, containing payload for single user
            txTrig = wlanWaveformGenerator(txPSDU{userIdx},cfgUser);

            % Pass waveform through a random TGax Channel
            channelIdx = chPermutations(userIdx);
            reset(tgax{channelIdx}); % New channel realization
            rxTrig = tgax{channelIdx}([txTrig; zeros(15,size(txTrig,2))]);

            % Scale the transmit power of the user within an RU. This is to
            % ensure same SNR for each user after the addition of noise.
            ruNum = cfgSys.User{userIdx}.RUNumber;
            SF = sqrt(1/sysInfo.NumUsersPerRU(ruNum))*sqrt(cfgUser.RUSize/(sum(sysInfo.RUSizes)));

            % Combine uplink users into one waveform
            rxWaveform = rxWaveform+SF*rxTrig;
        end

        % Pass the waveform through AWGN channel. Account for noise
        % energy in nulls so the SNR is defined per active subcarriers.
        awgn.SNR = snr(isnr)-10*log10(ofdmInfo.FFTLength/sum(sysInfo.RUSizes));
        rxWaveform = awgn(rxWaveform);

        % Receive processing
        % Packet detect and determine coarse packet offset
        coarsePktOffset = wlanPacketDetect(rxWaveform,chanBW,0,0.05);
        if isempty(coarsePktOffset) % If empty no L-STF detected; packet error
            numPacketErrors = numPacketErrors+1;
            continue; % Go to next loop iteration
        end

        % Extract the non-HT fields and determine fine packet offset
        nonhtfields = rxWaveform(coarsePktOffset+(ind.LSTF(1):ind.LSIG(2)),:);
        finePktOffset = wlanSymbolTimingEstimate(nonhtfields,chanBW);

        % Determine final packet offset
        pktOffset = coarsePktOffset+finePktOffset;

        % If packet detected out with the range of expected delays from
        % the channel modeling; packet error
        if pktOffset>50
            numPacketErrors = numPacketErrors+1;
            continue; % Go to next loop iteration
        end

        % Extract HE-LTF and HE-Data fields for all RUs
        rxLTF = rxWaveform(pktOffset+(ind.HELTF(1):ind.HELTF(2)),:);
        rxData = rxWaveform(pktOffset+(ind.HEData(1):ind.HEData(2)),:);

        for ruIdx = 1:allocInfo.NumRUs
            % Demodulate HE-LTF and HE-Data field for the RU of interest
            ru = [allocInfo.RUSizes(ruIdx) allocInfo.RUIndices(ruIdx)];
            demodHELTFRU = wlanHEDemodulate(rxLTF,'HE-LTF',chanBW,cfgSys.GuardInterval,cfgSys.HELTFType,ru);
            demodHEDataRU = wlanHEDemodulate(rxData,'HE-Data',chanBW,cfgSys.GuardInterval,ru);

            % Channel estimate
            [chanEst,ssPilotEst] = heLTFChannelEstimate(demodHELTFRU,cfgSys,ruIdx);

            % Get indices of data and pilots within RU (without nulls)
            ruOFDMInfo = wlanHEOFDMInfo('HE-Data',cfgSys.ChannelBandwidth,cfgSys.GuardInterval, ...
                [allocInfo.RUSizes(ruIdx) allocInfo.RUIndices(ruIdx)]);

            % Estimate noise power in HE fields of each user
            nVarEst = heNoiseEstimate(demodHEDataRU(ruOFDMInfo.PilotIndices,:,:),ssPilotEst,cfgSys,ruIdx);

            % Discard pilot subcarriers
            demodDataSym = demodHEDataRU(ruOFDMInfo.DataIndices,:,:);
            chanEstData = chanEst(ruOFDMInfo.DataIndices,:,:);

            % Equalize
            if strcmpi(equalizationMethod,'mmse-sic')
                [eqSym,csi] = heSuccessiveEqualize(demodDataSym,chanEstData,nVarEst,cfgSys,ruIdx);
            else
                [eqSym,csi] = heEqualizeCombine(demodDataSym,chanEstData,nVarEst,cfgSys);
            end

            for userIdx = 1:allocInfo.NumUsersPerRU(ruIdx)
                % Get TB config object for each user
                userNum = cfgSys.RU{ruIdx}.UserNumbers(userIdx);
                cfgUser = cfgTB{userNum};

                % Get space-time stream indices for the current user
                stsIdx = cfgUser.StartingSpaceTimeStream-1+(1:cfgUser.NumSpaceTimeStreams);

                % Demap and decode bits
                rxPSDU = wlanHEDataBitRecover(eqSym(:,:,stsIdx),nVarEst,csi(:,stsIdx),cfgUser,'LDPCDecodingMethod','layered-bp');

                % PER calculation
                packetError(userNum) = any(biterr(txPSDU{userNum},rxPSDU));
            end
        end
        numPacketErrors = numPacketErrors+packetError;
    end

    % Calculate packet error rate (PER) at SNR point
    packetErrorRate(:,isnr)= numPacketErrors/numPackets;
    disp(['SNR ' num2str(snr(isnr)) ...
          ' completed for ' num2str(numUsers) ' users']);

end

%% Plot Packet Error Rate vs SNR
markers = 'ox*sd^v><ph+ox*sd^v';
color = 'bmcrgbrkymcrgbrkymc';
figure;

for nSTA = 1:numUsers
    semilogy(snr,packetErrorRate(nSTA,:).',['-' markers(nSTA) color(nSTA)]);
    hold on;
end

grid on;
xlabel('SNR (dB)');
ylabel('PER');
dataStr = arrayfun(@(x)sprintf('STA- %d',x),1:numUsers,'UniformOutput',false);
legend(dataStr);
title('PER for uplink 802.11ax link');

%%
% The number of packets tested at each SNR point is controlled by
% |numPackets|. For meaningful results, this value should be larger than
% those presented in this example. The figure below was created by running
% a longer simulation with |numPackets|:1e4 and |snr|:20:2:28, which shows
% the packet error rate of both MMSE equalizer and MMSE-SIC equalizer.
%
% <<../HETBPER.png>>

%% Appendix
% This example uses the following helper functions and objects:
%
% * <matlab:edit('heEqualizeCombine.m') heEqualizeCombine.m>
% * <matlab:edit('heLTFChannelEstimate.m') heLTFChannelEstimate.m>
% * <matlab:edit('heNoiseEstimate.m') heNoiseEstimate.m>
% * <matlab:edit('heSuccessiveEqualize.m') heSuccessiveEqualize.m>
% * <matlab:edit('heTBRU.m') heTBRU.m>
% * <matlab:edit('heTBSystemConfig.m') heTBSystemConfig.m>
% * <matlab:edit('heTBUser.m') heTBUser.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.
% # M. Debbah, B. Muquet, M. de Courville, M. Muck, S. Simoens, and P.
% Loubaton. A MMSE successive interference cancellation scheme for a new
% adjustable hybrid spread OFDM system. IEEE 51st Vehicular Technology
% Conference Proceedings, pp. 745-749, vol. 2, 2000.
