%% 802.11ah Packet Error Rate Simulation for 2x2 TGah Channel
%
% This example shows how to measure the packet error rate of an IEEE(R)
% 802.11ah(TM) S1G short preamble link using an end-to-end simulation with
% a fading TGah indoor channel model and additive white Gaussian noise.

% Copyright 2017-2018 The MathWorks, Inc.

%% Introduction
% In this example an end-to-end simulation is used to determine the packet
% error rate for an 802.11ah [ <#11 1> ] S1G short preamble link with a
% fading channel at a selection of SNR points. For each SNR point multiple
% packets are transmitted through a channel, demodulated and the PSDUs
% recovered. The PSDUs are compared to those transmitted to determine the
% number of packet errors and hence the packet error rate. The processing
% for each packet is summarized in the following diagram.
%
% <<../S1GPERExampleDiagram.png>>
%
% This example simulates the S1G-Short format with no impairment correction
% apart from channel estimation and equalization. The received signal is
% synchronized to the start of the packet by compensating for a known delay
% and the default OFDM demodulation symbol sampling offset. No frequency
% synchronization is performed. For information on how to automatically
% detect and synchronize to the received signal, see the following examples
% for 802.11n(TM) and 802.11ac(TM):
%
% * <docid:wlan_examples#example-HTMIMOPacketErrorRateExample 802.11n
% Packet Error Rate Simulation for 2x2 TGn Channel>
% * <docid:wlan_examples#example-VHTMIMOPacketErrorRateExample 802.11ac
% Packet Error Rate Simulation for 8x8 TGac Channel>
%
% This example also demonstrates how a <docid:matlab_ref#f71-813245 parfor>
% loop can be used instead of the <docid:matlab_ref#buhafgy-1 for> loop
% when simulating each SNR point to speed up a simulation.
% <docid:matlab_ref#f71-813245 parfor> as part of the Parallel Computing
% Toolbox(TM), executes processing for each SNR in parallel to reduce the
% total simulation time.

%% Waveform Configuration
% A single-user 802.11ah S1G short preamble transmission is simulated in
% this example. An S1G format configuration object contains the format
% specific configuration of the transmission. The object is created using
% the <docid:wlan_ref#bvb8630 wlanS1GConfig> function. The properties of
% the object contain the configuration. In this example the object is
% configured for a 2 MHz channel bandwidth, short preamble, 2 transmit
% antennas, 2 space-time streams, 256 bytes payload, and 64-QAM rate-5/6
% (MCS 7).

% Create S1G configuration object for single user S1G short preamble
% transmission with 2 transmit antennas and 2 space-time streams
cfgS1G = wlanS1GConfig;
cfgS1G.ChannelBandwidth = 'CBW2'; % 2 MHz channel bandwidth
cfgS1G.Preamble = 'Short';        % Short preamble
cfgS1G.NumTransmitAntennas = 2;   % 2 transmit antennas
cfgS1G.NumSpaceTimeStreams = 2;   % 2 space-time streams
cfgS1G.APEPLength = 256;          % APEP length in bytes
cfgS1G.MCS = 7;                   % 64-QAM rate-5/6

%% Channel Configuration
% In this example a TGah N-LOS indoor channel model is used with delay
% profile Model-B. For Model-B, when the distance between transmitter and
% receiver is greater than or equal to 5 meters, the model is NLOS. This is
% described further in <docid:wlan_ref#bviq0w4-1 wlanTGahChannel>. A 2x2
% MIMO channel is simulated in this example therefore 2 receive antennas
% are specified.

% Create and configure the TGah channel
tgahChannel = wlanTGahChannel;
tgahChannel.DelayProfile = 'Model-B';
tgahChannel.NumTransmitAntennas = cfgS1G.NumTransmitAntennas;
tgahChannel.NumReceiveAntennas = 2;
tgahChannel.TransmitReceiveDistance = 5; % Distance in meters for NLOS
tgahChannel.ChannelBandwidth = cfgS1G.ChannelBandwidth;
tgahChannel.LargeScaleFadingEffect = 'None';

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

snr = 25:10:45;

%%
% The number of packets tested at each SNR point is controlled by two
% parameters:
%
% # |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 meaningful results we recommend increasing the numbers.

maxNumErrors = 1e2;  % The maximum number of packet errors at an SNR point
maxNumPackets = 1e3; % Maximum number of packets at an SNR point

%%
% Set the remaining variables for the simulation.

% Indices for accessing each field within the time-domain packet
fieldInd = wlanFieldIndices(cfgS1G);

% OFDM information
ofdmInfo = wlanS1GOFDMInfo('S1G-Data',cfgS1G);

% Set the sampling rate of the channel
tgahChannel.SampleRate = wlanSampleRate(cfgS1G);

if ~strcmp(packetFormat(cfgS1G),'S1G-Short')
    error('This example only supports the S1G-Short packet format');
end

%% Processing SNR Points
% For each SNR point a number of packets are tested and the packet error
% rate calculated. At each SNR point:
%
% # Multiple data packets are transmitted through a 2x2 TGah channel with
% AWGN.
% # Each packet is time synchronized given a known delay.
% # The S1G-LTF1 and S1G-LTF2N fields are demodulated and channel
% estimation is performed.
% # The S1G-Data field is extracted from the synchronized received
% waveform and OFDM demodulated.
% # The data carrying subcarriers are equalized using the channel
% estimates.
% # The PSDU is recovered using the equalized data subcarriers, noise
% variance estimate, and channel state information (CSI).
% # The recovered PSDU of each packet is compared to those transmitted to
% determine the number of packet errors and hence the packet error rate.
%
% A <docid:matlab_ref#f71-813245 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 <docid:comm_ref#buiamu7-1
% comm.AWGNChannel>.

packetErrorRate = zeros(numel(snr),1);
%parfor i = 1:numel(snr) % Use 'parfor' to speed up the simulation
for i = 1:numel(snr) % Use 'for' to debug the simulation
    % 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 = i;
    RandStream.setGlobalStream(stream);
    
    % Create an instance of the AWGN channel per SNR point simulated
    awgnChannel = comm.AWGNChannel;
    awgnChannel.NoiseMethod = 'Signal to noise ratio (SNR)';
    % Normalization
    awgnChannel.SignalPower = 1/tgahChannel.NumReceiveAntennas;
    % Account for energy in nulls
    awgnChannel.SNR = snr(i)-10*log10(ofdmInfo.FFTLength/ofdmInfo.NumTones);
    
    % Loop to simulate multiple packets
    numPacketErrors = 0;
    numPkt = 1; % Index of packet transmitted
    while numPacketErrors<=maxNumErrors && numPkt<=maxNumPackets
        % Generate a packet for 802.11ah short preamble
        txPSDU = randi([0 1],cfgS1G.PSDULength*8,1);
        txWaveform = wlanWaveformGenerator(txPSDU,cfgS1G);
        
        % Add trailing zeros to allow for channel delay
        tx = [txWaveform; zeros(50,cfgS1G.NumTransmitAntennas)];
        
        % Pass through fading indoor TGah channel
        reset(tgahChannel); % Reset channel for different realization
        rx = tgahChannel(tx);
        
        % Add noise
        rx = awgnChannel(rx);

        % Synchronize 
        % The received signal is synchronized to the start of the packet by
        % compensating for a known delay and the default OFDM demodulation
        % symbol sampling offset.
        delay = 4;
        rxSync = rx(delay+1:end,:);
        
        % LTF demodulation and channel estimation
        % Demodulate S1G-LTF1
        rxLTF1 = rxSync(fieldInd.S1GLTF1(1):fieldInd.S1GLTF1(2),:);
        demodLTF1 = wlanS1GDemodulate(rxLTF1,'S1G-LTF1',cfgS1G);
        
        % If required, demodulate S1G-LTF2N, and perform channel estimation
        if cfgS1G.NumSpaceTimeStreams>1
            % Use S1G-LTF1 and S1G-LTF2N for channel estimation
            rxLTF2N = rxSync(fieldInd.S1GLTF2N(1):fieldInd.S1GLTF2N(2),:);
            demodLTF2N = wlanS1GDemodulate(rxLTF2N,'S1G-LTF2N',cfgS1G);
            chanEst = s1gLTFChannelEstimate([demodLTF1 demodLTF2N],cfgS1G);
        else
            % Use only S1G-LTF1 for channel estimation
            chanEst = s1gLTFChannelEstimate(demodLTF1,cfgS1G);
        end

        % Noise variance estimate from S1G-LTF1 demodulated symbols
        noiseVarEst = helperNoiseEstimate(demodLTF1);
        
        % Extract S1G-Data field
        rxData = rxSync(fieldInd.S1GData(1):fieldInd.S1GData(2),:);
        
        % OFDM demodulation
        demodSym = wlanS1GDemodulate(rxData,'S1G-Data',cfgS1G);
        
        % Extract data subcarriers from demodulated symbols and channel
        % estimate
        demodDataSym = demodSym(ofdmInfo.DataIndices,:,:);
        chanEstData = chanEst(ofdmInfo.DataIndices,:,:);
        
        % MMSE frequency domain equalization
        [eqDataSym,csi] = helperSymbolEqualize(demodDataSym,chanEstData,noiseVarEst); 

        % Recover PSDU bits
        rxPSDU = s1gDataBitRecover(eqDataSym,noiseVarEst,csi,cfgS1G);
      
        % Determine if any bits are in error, i.e. a packet error
        packetError = any(biterr(txPSDU,rxPSDU));
        numPacketErrors = numPacketErrors+packetError;
        numPkt = numPkt+1;
    end
    
    % Compute PER for this SNR point
    packetErrorRate(i) = numPacketErrors/(numPkt-1);
    disp(['SNR ' num2str(snr(i))...
          ' completed after ' num2str(numPkt-1) ' packets,'...
          ' PER: ' num2str(packetErrorRate(i))]);
end

%% Plot Packet Error Rate vs SNR Results

figure;
semilogy(snr,packetErrorRate,'-ob');
grid on;
xlabel('SNR (dB)');
ylabel('PER');
title(['802.11ah ' cfgS1G.Preamble ', ' num2str(cfgS1G.ChannelBandwidth(4:end)) ...
    'MHz, MCS' num2str(cfgS1G.MCS) ', ' ...
    num2str(tgahChannel.NumTransmitAntennas) 'x' num2str(tgahChannel.NumReceiveAntennas) ...
    ' TGah Channel ' num2str(tgahChannel.DelayProfile)]);

%% Further Exploration
% 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|, and the SNR range |snr =
% 20:8:60|.
%
% <<../S1GShortPER.png>>

%% Appendix
% This example uses the following helper functions:
%
% * <matlab:edit('helperNoiseEstimate.m') helperNoiseEstimate.m>
% * <matlab:edit('helperSymbolEqualize.m') helperSymbolEqualize.m>
% * <matlab:edit('s1gDataBitRecover.m') s1gDataBitRecover.m>
% * <matlab:edit('s1gLTFChannelEstimate.m') s1gLTFChannelEstimate.m>

%% Selected Bibliography
% # IEEE P802.11ah(TM) D5.0 IEEE 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 2: Sub 1 GHz License Exempt Operation.
