%% 802.11ax Compressed Beamforming Packet Error Rate Simulation
%
% This example shows how to measure the packet error rate of a beamformed
% IEEE(R) 802.11ax(TM) high efficiency single user (HE SU) format link with
% different beamforming feedback quantization levels.

% Copyright 2018-2021 The MathWorks, Inc.

%% Introduction
% Transmit beamforming focuses energy towards a receiver to improve the SNR
% of a link. In this scheme, the transmitter is called a beamformer and the
% receiver is called a beamformee. A steering matrix is used by the
% beamformer to direct the energy to the beamformee. The steering matrix is
% calculated using channel state information obtained through channel
% measurements. These measurements are obtained by sounding the channel
% between beamformer and beamformee. To sound the channel, the beamformer
% sends a null data packet (NDP) to the beamformee. The beamformee measures
% the channel information during sounding to calculate a feedback matrix.
% This matrix is compressed in the form of quantized angles (phi and psi)
% and fed back to the beamformer. The beamformer can then calculate the
% feedback matrix from the quantized angles to create a steering matrix and
% beamform transmissions to the beamformee. The process of forming steering
% matrix is shown in <docid:wlan_ug#example-VHTBeamformingExample
% 802.11ac Transmit Beamforming>.
%
% In this example, a 4x2 MIMO configuration is considered between a
% transmitter and receiver, with two space-time streams used for each data
% packet transmission. An end-to-end simulation is used to determine the
% packet error rate (PER) for an 802.11ax [ <#12 1> ] single user (SU)
% format link with compressed beamforming feedback quantization for
% different quantization levels and a selection of SNR points. A plot is
% generated showing the PER vs SNR curve for each quantization resolution.
% This example does not consider grouping of sub-carriers (see Section
% 9.4.1.65 in [ <#12 1> ]).
%
% <<../HECompressedBeamformingExampleDiagram.png>>
%
%% Waveform Configuration
% An HE-SU packet is a full-band transmission to a single user. The
% transmit parameters for the HE-SU format are configured using a
% <docid:wlan_ref#mw_12ce61d2-4b46-4827-9281-e65b2413e9ba wlanHESUConfig>
% object. The properties of the object contain the configuration. In this
% example, the object is configured for a 20 MHz channel bandwidth, 4
% transmit antennas, 2 space-time streams and 16-QAM rate-1/2 (MCS 3).

NumTxAnts = 4; % Number of transmit antennas
NumSTS = 2;    % Number of space-time streams
NumRxAnts = 2; % Number of receive antennas
cfgHEBase = wlanHESUConfig;
cfgHEBase.ChannelBandwidth = 'CBW20';      % Channel bandwidth
cfgHEBase.NumSpaceTimeStreams = NumSTS;    % Number of space-time streams
cfgHEBase.NumTransmitAntennas = NumTxAnts; % Number of transmit antennas
cfgHEBase.APEPLength = 1e3;                % Payload length in bytes
cfgHEBase.ExtendedRange = false;           % Do not use extended range format
cfgHEBase.Upper106ToneRU = false;          % Do not use upper 106 tone RU
cfgHEBase.PreHESpatialMapping = false;     % Spatial mapping of pre-HE fields
cfgHEBase.GuardInterval = 0.8;             % Guard interval duration
cfgHEBase.HELTFType = 4;                   % HE-LTF compression mode
cfgHEBase.ChannelCoding = 'LDPC';          % Channel coding
cfgHEBase.MCS = 3;                         % Modulation and coding scheme
cfgHEBase.SpatialMapping = 'Custom';       % Custom for beamforming

%% Null Data Packet (NDP) Configuration
% The NDP transmission is configured to have data length of zero. Since
% the NDP is used to obtain the channel state information, the number of
% space-time streams is equal to number of transmit antennas. This results
% in a direct mapping of each space-time stream to a transmit antenna.

cfgNDP = cfgHEBase;
cfgNDP.APEPLength = 0;                  % NDP has no data
cfgNDP.NumSpaceTimeStreams = NumTxAnts; % For feedback matrix calculation
cfgNDP.SpatialMapping = 'Direct';       % Each TxAnt carries a STS

%% Channel Configuration
% In this example, a TGax NLOS indoor channel model is used with delay
% profile Model-B. The Model-B profile is considered NLOS when the distance
% between transmitter and receiver is greater than or equal to 5 meters.
% This is described further in <docid:wlan_ref#mw_43b5900e-69e1-4636-b084-1e72dbd46293 wlanTGaxChannel>.
% A 4x2 MIMO channel is simulated in this example.

% Create and configure the TGax channel
chanBW = cfgHEBase.ChannelBandwidth;
tgaxChannel = wlanTGaxChannel;
tgaxChannel.DelayProfile = 'Model-B';
tgaxChannel.NumTransmitAntennas = NumTxAnts;
tgaxChannel.NumReceiveAntennas = NumRxAnts;
tgaxChannel.TransmitReceiveDistance = 5; % Distance in meters for NLOS
tgaxChannel.ChannelBandwidth = chanBW;
tgaxChannel.LargeScaleFadingEffect = 'None';
tgaxChannel.NormalizeChannelOutputs = false;
fs = wlanSampleRate(cfgHEBase);
tgaxChannel.SampleRate = fs;

%% Simulation Parameters
% This example compares the performance of beamforming with two different
% resolutions of compression quantization, and without compression. For
% each quantization resolution, an end to end simulation with various SNR
% values is run to determine the packet error rate. 802.11ax Draft 4.1
% specifies only two sets of quantization resolution for single user
% beamforming (Table 9-31a in [ <#12 1> ]). The value of |codeBookSize|
% determines the number of bits used to quantize the beamforming feedback
% angles (phi and psi) in this simulation. When |codeBookSize| is Inf, no
% compression is performed. The quantization levels selected by
% |codeBookSize| are shown in the table below:
%
%  codeBookSize            Compression Configuration
% --------------------------------------------------------
%      0                   NumBitsphi = 4; NumBitspsi = 2
%      1                   NumBitsphi = 6; NumBitspsi = 4
%     Inf                  No compression
% --------------------------------------------------------

codeBookSize = [0 1 Inf];

%%
% A number of packets are generated, passed through a channel and
% demodulated to determine the packet error rate for each compression
% configuration at each SNR (dB) value in the |snr| vector.

snr = 10:2:18;

%%
% The number of packets tested at each SNR point is limited to a
% |maxNumErrors| or |maxNumPackets|:
%
% # |maxNumErrors| is the maximum number of packet errors simulated at each
% SNR point. When the number of packet errors reaches this limit, the
% simulation at this SNR point is complete.
% # |maxNumPackets| is the maximum number of packets simulated at each SNR
% point and limits the length of the simulation if the packet error limit
% is not reached.
%
% The numbers chosen in this example will lead to a very short simulation.
% For statistically meaningful results we recommend increasing these
% numbers.

maxNumErrors = 10;   % The maximum number of packet errors at an SNR point
maxNumPackets = 100; % The maximum number of packets at an SNR point

%% Processing SNR Points
% For each SNR point, a number of packets are tested and the packet error
% rate calculated. The pre-HE preamble of 802.11ax is backwards compatible
% with 802.11ac(TM), therefore in this example the front-end
% synchronization components for a VHT waveform are used to synchronize the
% HE waveform at the receiver. For each packet, the following processing
% steps occur.
%
% The beamformer obtains the steering matrix by transmitting an NDP which
% is processed by the beamformee to create a feedback matrix:
%
% # An NDP waveform is transmitted through an indoor TGax channel model.
% Different channel realizations are modeled for different packets.
% # AWGN is added to the received waveform to create the desired average
% SNR per active subcarrier after OFDM demodulation.
% # The packet is detected at the beamformee.
% # Coarse carrier frequency offset is estimated and corrected.
% # 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.
% # Fine carrier frequency offset is estimated and corrected.
% # The HE-LTF is extracted from the synchronized received waveform. The
% HE-LTF is OFDM demodulated and channel estimation is performed.
% # Singular value decomposition is performed on the estimated channel and
% the beamforming feedback matrix, V is calculated.
% # If there is no compression, this feedback matrix, V will be used as the
% steering matrix by the beamformer.
% # If compression is used, the feedback matrix, V will be compressed and
% quantized to create a set of angles as specified in the standard.
%
% The beamformer transmits a data packet using the recovered steering
% matrix and the beamformee decodes the beamformed data transmission to
% recover the PSDU:
%
% # Since the current example assumes zero delay in getting the beamforming
% feedback from the beamformee, the quantized angles are converted back to
% the beamforming feedback matrix, V.
% # A PSDU is created and encoded to create a single packet waveform with
% the steering matrix set to the beamforming feedback matrix, V.
% # The waveform is passed through the same indoor TGax channel realization
% as the NDP transmission.
% # AWGN is added to the received waveform.
% # As with NDP, synchronization and HE channel estimation are performed.
% # The data field is extracted from the synchronized received waveform
% and OFDM demodulated.
% # Common phase error pilot tracking is performed to track any residual
% carrier frequency offset.
% # Noise estimation is performed using the demodulated data field pilots
% and single-stream channel estimate at pilot subcarriers.
% # The phase corrected OFDM symbols are equalized with the channel
% estimate.
% # The equalized symbols are demodulated and decoded to recover the PSDU.
%
% A <docid:matlab_ref#f71-813245 parfor> loop can be used to parallelize
% processing of the SNR points. To enable the use of parallel computing
% for increased speed comment out the 'for' statement and uncomment the
% 'parfor' statement below.

numQuant = numel(codeBookSize);
numSNR = numel(snr); % Number of SNR points
packetErrorRate = zeros(numQuant,numSNR);

% Get occupied subcarrier indices and OFDM parameters
ofdmInfo = wlanHEOFDMInfo('HE-Data',cfgHEBase);

% Indices to extract fields from the PPDU
ind = wlanFieldIndices(cfgHEBase);
indSound = wlanFieldIndices(cfgNDP);

for ibf = 1:numQuant
    switch codeBookSize(ibf) % See P802.11ax/D4.1 Section 9.4.1.64
        case 0
            NumBitsPsi = 2; % Number of bits for psi
            NumBitsPhi = 4; % Number of bits for phi
            disp('End-to-End simulation with compressed beamforming quantization with');
            disp(['Number of Bits for phi = ' num2str(NumBitsPhi) ...
                ' and Number of Bits for psi = ' num2str(NumBitsPsi)]);
        case 1
            NumBitsPsi = 4; % Number of bits for psi
            NumBitsPhi = 6; % Number of bits for phi
            disp('End-to-End simulation with compressed beamforming quantization with');
            disp(['Number of Bits for phi = ' num2str(NumBitsPhi) ...
                ' and Number of Bits for psi = ' num2str(NumBitsPsi)]);
        otherwise
            disp('End-to-End simulation with non-compressed beamforming');
    end
    
    %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',100);
        stream.Substream = isnr;
        RandStream.setGlobalStream(stream);

        % Account for noise energy in nulls so the SNR is defined per
        % active subcarrier
        packetSNR = snr(isnr)-10*log10(ofdmInfo.FFTLength/ofdmInfo.NumTones);
        
        % Create an instance of the HE configuration object per SNR point
        % simulated. This will enable to use parfor
        cfgHE = cfgHEBase;
        
        % Loop to simulate multiple packets
        numPacketErrors = 0;
        numPkt = 1; % Index of packet transmitted
        while numPacketErrors<=maxNumErrors && numPkt<=maxNumPackets
            % Null data packet transmission
            tx = wlanWaveformGenerator([],cfgNDP);
            
            % Add trailing zeros to allow for channel delay
            txPad = [tx; zeros(50,cfgNDP.NumTransmitAntennas)];
            
            % Pass through a fading indoor TGax channel
            reset(tgaxChannel); % Reset channel for different realization
            rx = tgaxChannel(txPad);

            % Pass the waveform through AWGN channel
            rx = awgn(rx,packetSNR);
            
            % Calculate the steering matrix at the beamformee
            V = heUserBeamformingFeedback(rx,cfgNDP,true);
            
            if isempty(V)
                % User feedback failed, packet error
                numPacketErrors = numPacketErrors+1;
                numPkt = numPkt+1;
                continue; % Go to next loop iteration 
            end
            
            if ~isinf(codeBookSize(ibf))
                % Find quantized angles of the beamforming feedback matrix
                angidx = bfCompressQuantize(V(:,1:NumSTS,:),NumBitsPhi,NumBitsPsi);
                
                % Calculate steering matrix from the quantized angles at
                % beamformer:
                % Assuming zero delay in transmitting the quantized angles
                % from beamformee to beamformer, the steering matrix is
                % calculated from the quantized angles and is used in the
                % data transmission of beamformer.
                
                [~,Nc,Nr] = size(V(1,1:NumSTS,:));
                V = bfDecompress(angidx,Nr,Nc,NumBitsPhi,NumBitsPsi);
            end

            steeringMat = V(:,1:NumSTS,:);
            
            % Beamformed data transmission
            psduLength = getPSDULength(cfgHE); % PSDU length in bytes
            txPSDU = randi([0 1],psduLength*8,1); % Generate random PSDU
            cfgHE.SpatialMappingMatrix = steeringMat;
            tx = wlanWaveformGenerator(txPSDU,cfgHE);
            
            % Add trailing zeros to allow for channel delay
            txPad = [tx; zeros(50,cfgHE.NumTransmitAntennas)];
            
            % Pass through a fading indoor TGax channel
            rx = tgaxChannel(txPad);
            
            % Pass the waveform through AWGN channel
            rx = awgn(rx,packetSNR);
            
            % Packet detect and determine coarse packet offset
            coarsePktOffset = wlanPacketDetect(rx,chanBW);
            if isempty(coarsePktOffset) % If empty no L-STF detected; packet error
                numPacketErrors = numPacketErrors+1;
                numPkt = numPkt+1;
                continue; % Go to next loop iteration
            end
            
            % Extract L-STF and perform coarse frequency offset correction
            lstf = rx(coarsePktOffset+(ind.LSTF(1):ind.LSTF(2)),:);
            coarseFreqOff = wlanCoarseCFOEstimate(lstf,chanBW);
            rx = helperFrequencyOffset(rx,fs,-coarseFreqOff);
            
            % Extract the non-HT fields and determine fine packet offset
            nonhtfields = rx(coarsePktOffset+(ind.LSTF(1):ind.LSIG(2)),:);
            finePktOffset = wlanSymbolTimingEstimate(nonhtfields,chanBW);
            
            % Determine final packet offset
            pktOffset = coarsePktOffset+finePktOffset;
            
            % If packet detected outwith the range of expected delays from
            % the channel modeling; packet error
            if pktOffset>50
                numPacketErrors = numPacketErrors+1;
                numPkt = numPkt+1;
                continue; % Go to next loop iteration
            end
            
            % Extract L-LTF and perform fine frequency offset correction
            rxLLTF = rx(pktOffset+(ind.LLTF(1):ind.LLTF(2)),:);
            fineFreqOff = wlanFineCFOEstimate(rxLLTF,chanBW);
            rx = helperFrequencyOffset(rx,fs,-fineFreqOff);
            
            % HE-LTF demodulation and channel estimation
            rxHELTF = rx(pktOffset+(ind.HELTF(1):ind.HELTF(2)),:);
            heltfDemod = wlanHEDemodulate(rxHELTF,'HE-LTF',cfgHE);
            [chanEst,pilotEst] = heLTFChannelEstimate(heltfDemod,cfgHE);
            
            % Data demodulate
            rxData = rx(pktOffset+(ind.HEData(1):ind.HEData(2)),:);
            demodSym = wlanHEDemodulate(rxData,'HE-Data',cfgHE);
            
            % Pilot phase tracking
            % Average single-stream pilot estimates over symbols (2nd dimension)
            pilotEstTrack = mean(pilotEst,2);
            demodSym = heCommonPhaseErrorTracking(demodSym,pilotEstTrack,cfgHE);
            
            % Estimate noise power in HE fields
            nVarEst = heNoiseEstimate(demodSym(ofdmInfo.PilotIndices,:,:),pilotEstTrack,cfgHE);
            
            % Extract data subcarriers from demodulated symbols and channel
            % estimate
            demodDataSym = demodSym(ofdmInfo.DataIndices,:,:);
            chanEstData = chanEst(ofdmInfo.DataIndices,:,:);
            
            % Equalization and STBC combining
            [eqDataSym,csi] = heEqualizeCombine(demodDataSym,chanEstData,nVarEst,cfgHE);
            
            % Recover data
            rxPSDU = wlanHEDataBitRecover(eqDataSym,nVarEst,csi,cfgHE,'LDPCDecodingMethod','layered-bp');
            
            % Determine if any bits are in error, i.e. a packet error
            packetError = ~isequal(txPSDU,rxPSDU);
            numPacketErrors = numPacketErrors+packetError;
            numPkt = numPkt+1;
        end
        
        % Calculate packet error rate (PER) at SNR point
        packetErrorRate(ibf,isnr) = numPacketErrors/(numPkt-1);
        disp(['MCS ' num2str(cfgHE.MCS) ','...
            ' SNR ' num2str(snr(isnr)) ...
            ' completed after ' num2str(numPkt-1) ' packets,'...
            ' PER:' num2str(packetErrorRate(ibf,isnr))]);
    end
    disp(newline);
end

%% Plot Packet Error Rate vs Signal to Noise Ratio

figure;
lineTypes = ["k-o" "b-s" "r-*"];
semilogy(snr,packetErrorRate(1,:),lineTypes(1));
hold on;
grid on;
xlabel('SNR (dB)');
ylabel('PER');
for ibf = 2:numQuant
    semilogy(snr,packetErrorRate(ibf,:),lineTypes(ibf));
end
dataStr = [string(['Compressed Beamforming, ' newline ...
                   'NumBitsPhi = 4, NumBitsPsi = 2' newline])...
           string(['Compressed Beamforming, ' newline ...
                   'NumBitsPhi = 6, NumBitsPsi = 4' newline]) ...
                   "Non-Compressed Beamforming"];
legend(dataStr);
title(sprintf('802.11ax Beamforming PER for Channel %s, %s, %s',tgaxChannel.DelayProfile,cfgHEBase.ChannelBandwidth,cfgHEBase.ChannelCoding));

%%
% The number of packets tested at each SNR point is controlled by two
% parameters: |maxNumErrors| and |maxNumPackets|. For meaningful results,
% these values should be larger than those presented in this example. As an
% example, the figure below was created by running a longer simulation with
% |maxNumErrors|:1e3 and |maxNumPackets|:1e4.
%
% <<../HECompressedBeamformingExamplePER.png>>

%% Appendix
% This example uses the following helper functions:
%
% * <matlab:edit('bfCompressQuantize.m') bfCompressQuantize.m>
% * <matlab:edit('bfDecompress.m') bfDecompress.m>
% * <matlab:edit('heCommonPhaseErrorTracking.m') heCommonPhaseErrorTracking.m>
% * <matlab:edit('heEqualizeCombine.m') heEqualizeCombine.m>
% * <matlab:edit('helperFrequencyOffset.m') helperFrequencyOffset.m>
% * <matlab:edit('heLTFChannelEstimate.m') heLTFChannelEstimate.m>
% * <matlab:edit('heNoiseEstimate.m') heNoiseEstimate.m>
% * <matlab:edit('heUserBeamformingFeedback.m') heUserBeamformingFeedback.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.
