% Create a format configuration object for a DMG OFDM transmission
cfgDMG = wlanDMGConfig;
cfgDMG.PSDULength = 8192; % bytes

% For DMG OFDM PHY, the valid range of MCS is 13-24(inclusive)
mcs = 21; % OFDM PHY, 16QAM, rate 13/16

% SNR ranges to use for AWGN
snrRanges = {...
    -1:0.5:1.5, ...  % MCS 13
    0:0.5:2.5, ...   % MCS 14
    1.5:0.5:4, ...   % MCS 15
    3:0.5:5.5, ...   % MCS 16
    4.5:0.5:7, ...   % MCS 17
    7.5:0.5:10, ...  % MCS 18
    9:0.5:11.5, ...  % MCS 19
    10.5:0.5:13, ... % MCS 20
    12:0.5:14.5, ... % MCS 21
    14.5:0.5:17, ... % MCS 22
    16.5:0.5:19, ... % MCS 23
    17.5:0.5:20, ... % MCS 24
    };

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

% OFDM information
ofdmInfo = wlanDMGOFDMInfo();

% Indices of data and pilot occupied subcarriers
cfgDMG.MCS = mcs(1); % Set OFDM MCS to get subcarrier indices
Nsd = numel(ofdmInfo.DataIndices); % Number of data carrying subcarriers

numSNR = numel(snrRanges{1}); % Number of SNR points
numMCS = numel(mcs);          % Number of MCS
packetErrorRate = zeros(numMCS,numSNR);

for imcs = 1:numMCS
    cfgDMG.MCS = mcs(imcs);
    if ~strcmp(phyType(cfgDMG),'OFDM')
        error('This example only supports DMG OFDM PHY simulation');
    end

    % Indices of fields within the packet
    fieldIndices = wlanFieldIndices(cfgDMG);

    % SNR points to simulate from MCS
    snr = snrRanges{cfgDMG.MCS-12};

    %parfor isnr = 1:numSNR % Use 'parfor' to speed up the simulation
    for isnr = 1:numSNR % 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 = isnr;
        RandStream.setGlobalStream(stream);

        % Create an instance of the AWGN channel per SNR point simulated
        awgnChannel = comm.AWGNChannel;
        awgnChannel.NoiseMethod = 'Signal to noise ratio (SNR)';
        awgnChannel.SignalPower = 1;
        % Account for noise energy in nulls so the SNR is defined per
        % active subcarrier
        awgnChannel.SNR = snr(isnr)-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 waveform
            txPSDU = randi([0 1],cfgDMG.PSDULength*8,1); % PSDULength in bytes
            tx = wlanWaveformGenerator(txPSDU,cfgDMG);

            % Pass the waveform through AWGN channel model
            rx = awgnChannel(tx);

            % Extract data field
            rxData = rx(fieldIndices.DMGData(1):fieldIndices.DMGData(2));

            % OFDM demodulate
            demodSym = wlanDMGOFDMDemodulate(rxData);
            dataSym = demodSym(ofdmInfo.DataIndices,:); % Discard pilots

            % Equalize
            chanSym = ones(Nsd,1); % Set channel gains to 1 as AWGN channel
            nVar = 10^(-snr(isnr)/10); % Noise variance
            [eqSym,csi] = helperSymbolEqualize(dataSym,chanSym,nVar);

            % Recover data
            rxPSDU = wlanDMGDataBitRecover(eqSym,nVar,csi,cfgDMG);

            % Determine if any bits are in error, i.e. a packet error
            packetError = any(biterr(txPSDU,rxPSDU));
            numPacketErrors = numPacketErrors+packetError;
            numPkt = numPkt+1;
        end

        % Calculate packet error rate (PER) at SNR point
        packetErrorRate(imcs,isnr) = numPacketErrors/(numPkt-1);
        disp(['MCS ' num2str(mcs(imcs)) ','...
              ' SNR ' num2str(snr(isnr)) ...
              ' completed after ' num2str(numPkt-1) ' packets,'...
              ' PER:' num2str(packetErrorRate(imcs,isnr))]);
    end
end

markers = 'ox*sd^v><ph+';
color = 'bmcrgbrkymcr';
figure;
for imcs = 1:numMCS
    semilogy(snrRanges{mcs(imcs)-12},packetErrorRate(imcs,:).',['-' markers(imcs) color(imcs)]);
    hold on;
end
grid on;
xlabel('SNR (dB)');
ylabel('PER');
dataStr = arrayfun(@(x)sprintf('MCS %d',x),mcs,'UniformOutput',false);
legend(dataStr);
title('PER for DMG OFDM-PHY with AWGN channel');